|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00017 MODULE xc_lyp_adiabatic 00018 USE bibliography, ONLY: Lee1988,& 00019 cite_reference 00020 USE f77_blas 00021 USE input_section_types, ONLY: section_vals_type,& 00022 section_vals_val_get 00023 USE kinds, ONLY: dp 00024 USE mathconstants, ONLY: pi 00025 USE timings, ONLY: timeset,& 00026 timestop 00027 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00028 xc_dset_get_derivative 00029 USE xc_derivative_types, ONLY: xc_derivative_get,& 00030 xc_derivative_type 00031 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00032 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 00033 xc_rho_set_type 00034 #include "cp_common_uses.h" 00035 00036 IMPLICIT NONE 00037 PRIVATE 00038 00039 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00040 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_lyp_adiabatic' 00041 REAL(kind=dp), PARAMETER, PRIVATE :: a=0.04918_dp, b=0.132_dp, 00042 c=0.2533_dp,d=0.349_dp 00043 00044 PUBLIC :: lyp_adiabatic_lda_info, lyp_adiabatic_lsd_info, lyp_adiabatic_lda_eval, lyp_adiabatic_lsd_eval 00045 00046 CONTAINS 00047 00048 ! ***************************************************************************** 00060 SUBROUTINE lyp_adiabatic_lda_info(reference,shortform, needs, max_deriv, error) 00061 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00062 TYPE(xc_rho_cflags_type), 00063 INTENT(inout), OPTIONAL :: needs 00064 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00065 TYPE(cp_error_type), INTENT(inout) :: error 00066 00067 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_adiabatic_lda_info', 00068 routineP = moduleN//':'//routineN 00069 00070 IF ( PRESENT ( reference ) ) THEN 00071 reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LDA version}" 00072 END IF 00073 IF ( PRESENT ( shortform ) ) THEN 00074 shortform = "Lee-Yang-Parr correlation energy functional (LDA)" 00075 END IF 00076 IF (PRESENT(needs)) THEN 00077 needs%rho=.TRUE. 00078 needs%rho_1_3=.TRUE. 00079 needs%norm_drho=.TRUE. 00080 END IF 00081 IF (PRESENT(max_deriv)) max_deriv=1 00082 00083 END SUBROUTINE lyp_adiabatic_lda_info 00084 00085 ! ***************************************************************************** 00097 SUBROUTINE lyp_adiabatic_lsd_info(reference,shortform, needs, max_deriv, error) 00098 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00099 TYPE(xc_rho_cflags_type), 00100 INTENT(inout), OPTIONAL :: needs 00101 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00102 TYPE(cp_error_type), INTENT(inout) :: error 00103 00104 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_adiabatic_lsd_info', 00105 routineP = moduleN//':'//routineN 00106 00107 IF ( PRESENT ( reference ) ) THEN 00108 reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LSD version}" 00109 END IF 00110 IF ( PRESENT ( shortform ) ) THEN 00111 shortform = "Lee-Yang-Parr correlation energy functional (LSD)" 00112 END IF 00113 IF (PRESENT(needs)) THEN 00114 needs%rho_spin=.TRUE. 00115 needs%norm_drho_spin=.TRUE. 00116 needs%norm_drho=.TRUE. 00117 END IF 00118 IF (PRESENT(max_deriv)) max_deriv=1 00119 END SUBROUTINE lyp_adiabatic_lsd_info 00120 00121 ! ***************************************************************************** 00126 SUBROUTINE lyp_adiabatic_lda_eval(rho_set,deriv_set,grad_deriv,lyp_adiabatic_params,error) 00127 TYPE(xc_rho_set_type), POINTER :: rho_set 00128 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00129 INTEGER, INTENT(in) :: grad_deriv 00130 TYPE(section_vals_type), POINTER :: lyp_adiabatic_params 00131 TYPE(cp_error_type), INTENT(inout) :: error 00132 00133 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_adiabatic_lda_eval', 00134 routineP = moduleN//':'//routineN 00135 00136 INTEGER :: handle, npoints, stat 00137 INTEGER, DIMENSION(:, :), POINTER :: bo 00138 LOGICAL :: failure 00139 REAL(kind=dp) :: epsilon_norm_drho, 00140 epsilon_rho, lambda 00141 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, 00142 e_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, 00143 e_rho_rho, e_rho_rho_rho, norm_drho, rho, rho_1_3 00144 TYPE(xc_derivative_type), POINTER :: deriv 00145 00146 CALL timeset(routineN,handle) 00147 00148 failure=.FALSE. 00149 NULLIFY(bo) 00150 00151 CALL section_vals_val_get(lyp_adiabatic_params,"LAMBDA",r_val=lambda,error=error) 00152 CALL cite_reference(Lee1988) 00153 00154 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00155 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00156 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00157 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00158 IF (.NOT. failure) THEN 00159 CALL xc_rho_set_get(rho_set,rho_1_3=rho_1_3,rho=rho,& 00160 norm_drho=norm_drho,local_bounds=bo,rho_cutoff=epsilon_rho,& 00161 drho_cutoff=epsilon_norm_drho,error=error) 00162 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00163 00164 ! meaningful default for the arrays we don't need: let us make compiler 00165 ! and debugger happy... 00166 IF (cp_debug) THEN 00167 ALLOCATE(dummy(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 00168 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00169 ELSE 00170 dummy=> rho 00171 END IF 00172 00173 e_0 => dummy 00174 e_rho => dummy 00175 e_ndrho => dummy 00176 e_rho_rho => dummy 00177 e_ndrho_rho => dummy 00178 e_ndrho_ndrho => dummy 00179 e_rho_rho_rho => dummy 00180 e_ndrho_rho_rho => dummy 00181 e_ndrho_ndrho_rho => dummy 00182 00183 IF (grad_deriv>=0) THEN 00184 deriv => xc_dset_get_derivative(deriv_set,"",& 00185 allocate_deriv=.TRUE., error=error) 00186 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00187 END IF 00188 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 00189 deriv => xc_dset_get_derivative(deriv_set,"(rho)",& 00190 allocate_deriv=.TRUE.,error=error) 00191 CALL xc_derivative_get(deriv,deriv_data=e_rho,error=error) 00192 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",& 00193 allocate_deriv=.TRUE.,error=error) 00194 CALL xc_derivative_get(deriv,deriv_data=e_ndrho,error=error) 00195 END IF 00196 IF (grad_deriv>1.OR.grad_deriv<-1) THEN 00197 CALL cp_unimplemented_error(fromWhere=routineP, & 00198 message="derivatives bigger than 1 not implemented", & 00199 error=error, error_level=cp_failure_level) 00200 END IF 00201 00202 !$omp parallel default(none) & 00203 !$omp shared(rho, rho_1_3, norm_drho, e_0, e_rho, e_ndrho) & 00204 !$omp shared(e_rho_rho, e_ndrho_rho, e_ndrho_ndrho) & 00205 !$omp shared(e_rho_rho_rho, e_ndrho_rho_rho) & 00206 !$omp shared(e_ndrho_ndrho_rho, grad_deriv, npoints) & 00207 !$omp shared(epsilon_rho, epsilon_norm_drho, lambda, error) 00208 00209 CALL lyp_adiabatic_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho,& 00210 e_0=e_0,e_rho=e_rho,e_ndrho=e_ndrho,e_rho_rho=e_rho_rho,& 00211 e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, & 00212 e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho,& 00213 e_ndrho_ndrho_rho=e_ndrho_ndrho_rho,& 00214 grad_deriv=grad_deriv,& 00215 npoints=npoints,epsilon_rho=epsilon_rho,epsilon_norm_drho=epsilon_norm_drho,lambda=lambda,& 00216 error=error) 00217 00218 !$omp end parallel 00219 00220 IF (cp_debug) THEN 00221 DEALLOCATE(dummy,stat=stat) 00222 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00223 ELSE 00224 NULLIFY(dummy) 00225 END IF 00226 END IF 00227 CALL timestop(handle) 00228 END SUBROUTINE lyp_adiabatic_lda_eval 00229 00230 ! ***************************************************************************** 00235 SUBROUTINE lyp_adiabatic_lda_calc(rho, rho_1_3, norm_drho,& 00236 e_0,e_rho,e_ndrho,e_rho_rho,e_ndrho_rho,& 00237 e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho,& 00238 grad_deriv,npoints,epsilon_rho, epsilon_norm_drho,lambda,error) 00239 INTEGER, INTENT(in) :: npoints, grad_deriv 00240 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_rho, 00241 e_ndrho_rho_rho, e_rho_rho_rho, e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, 00242 e_ndrho, e_rho, e_0 00243 REAL(kind=dp), DIMENSION(1:npoints), 00244 INTENT(in) :: norm_drho, rho_1_3, rho 00245 REAL(kind=dp), INTENT(in) :: epsilon_rho, 00246 epsilon_norm_drho, lambda 00247 TYPE(cp_error_type), INTENT(inout) :: error 00248 00249 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_adiabatic_lda_calc', 00250 routineP = moduleN//':'//routineN 00251 00252 INTEGER :: ii 00253 REAL(kind=dp) :: cf, my_ndrho, my_rho, t10, t107, t11, t117, t12, t122, 00254 t125, t13, t14, t15, t153, t16, t17, t180, t189, t19, t195, t2, t20, 00255 t25, t28, t29, t3, t34, t36, t37, t38, t4, t40, t41, t42, t43, t45, 00256 t46, t47, t50, t51, t52, t57, t58, t59, t6, t63, t65, t7, t71, t77, 00257 t78, t87, t9, t94 00258 00259 cf=0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp) 00260 00261 !$omp do 00262 00263 DO ii=1,npoints 00264 my_rho = rho(ii) 00265 IF (my_rho>epsilon_rho) THEN 00266 IF( grad_deriv >=0 )THEN 00267 my_ndrho = norm_drho(ii) 00268 t2 = d * lambda 00269 t3 = my_rho ** (0.1e1_dp / 0.3e1_dp) 00270 t4 = 0.1e1_dp / t3 00271 t6 = 0.10e1_dp + t2 * t4 00272 t7 = 0.1e1_dp / t6 00273 t9 = a * b 00274 t10 = t9 * my_rho 00275 t11 = c * lambda 00276 t12 = t11 * t4 00277 t13 = EXP(-t12) 00278 t14 = t13 * t7 00279 t15 = my_ndrho ** 2 00280 t16 = my_rho ** 2 00281 t17 = t3 ** 2 00282 t19 = 0.1e1_dp / t17 / t16 00283 t20 = t15 * t19 00284 t25 = 0.30e1_dp + 0.70e1_dp * t12 + 0.70e1_dp * t2 * t4 * t7 00285 t28 = Cf - 0.1388888889e-1_dp * t20 * t25 00286 t29 = t14 * t28 00287 t34 = lambda ** 2 00288 t36 = t6 ** 2 00289 t37 = 0.1e1_dp / t36 00290 t38 = t37 * d 00291 t40 = t9 * t17 00292 t41 = c * t13 00293 t42 = t7 * t28 00294 t43 = t41 * t42 00295 t45 = t13 * t37 00296 t46 = t28 * d 00297 t47 = t45 * t46 00298 t50 = 0.1e1_dp / t17 / my_rho 00299 t51 = t9 * t50 00300 t52 = c * t4 00301 t57 = d ** 2 00302 t58 = t57 * lambda 00303 t59 = 0.1e1_dp / t17 00304 t63 = 0.70e1_dp * t52 + 0.70e1_dp * d * t4 * t7 - 0.70e1_dp * t58 * t59 * t37 00305 t65 = t14 * t15 * t63 00306 00307 e_0(ii) = e_0(ii) + 0.20e1_dp * lambda * (-a * my_rho * t7 - t10 * t29) + t34 * (a * t17 & 00308 * t38 + t40 * t43 + t40 * t47 + 0.13888888888888888889e-1_dp * t51 * & 00309 t65) 00310 00311 END IF 00312 IF(grad_deriv >=1 ) THEN 00313 t71 = a * t4 00314 t77 = lambda * t13 00315 t78 = t77 * t42 00316 t87 = t16 * my_rho 00317 t94 = 0.1e1_dp / t3 / my_rho 00318 t107 = 0.37037037037037037037e-1_dp * t15 / t17 / t87 * t25 - 0.1388888889e-1_dp & 00319 * t20 * (-0.2333333333e1_dp * t11 * t94 - 0.2333333333e1_dp * t2 & 00320 * t94 * t7 + 0.23333333333333333333e1_dp * t57 * t34 * t50 * t37) 00321 t117 = 0.1e1_dp / t36 / t6 00322 t122 = t9 * t4 00323 t125 = c ** 2 00324 t153 = 0.1e1_dp / t87 00325 t180 = 0.1e1_dp / t16 00326 t189 = 0.2e1_dp / 0.3e1_dp * t71 * t38 + 0.2e1_dp / 0.3e1_dp * a * t59 * t117 * & 00327 t57 * lambda + 0.2e1_dp / 0.3e1_dp * t122 * t43 + t9 * t59 * t125 * t78 & 00328 / 0.3e1_dp + 0.2e1_dp / 0.3e1_dp * t9 * t59 * c * t45 * t46 * lambda + t40 & 00329 * t41 * t7 * t107 + 0.2e1_dp / 0.3e1_dp * t122 * t47 + 0.2e1_dp / 0.3e1_dp * & 00330 t9 * t59 * t13 * t117 * t28 * t58 + t40 * t45 * t107 * d - 0.2314814815e-1_dp & 00331 * t9 * t19 * t65 + 0.46296296296296296297e-2_dp * t9 * t153 & 00332 * c * t77 * t7 * t15 * t63 + 0.46296296296296296297e-2_dp * t9 * t153 & 00333 * t13 * t37 * t15 * t63 * d * lambda + 0.13888888888888888889e-1_dp & 00334 * t51 * t14 * t15 * (-0.2333333333e1_dp * c * t94 - 0.2333333333e1_dp * & 00335 d * t94 * t7 + 0.70000000000000000000e1_dp * t57 * t50 * t37 * lambda & 00336 - 0.4666666667e1_dp * t57 * d * t34 * t180 * t117) 00337 00338 e_rho(ii) = e_rho(ii) + 0.20e1_dp * lambda * (-a * t7 - t71 * t38 * lambda / 0.3e1_dp - t9 * & 00339 t29 - t9 * t52 * t78 / 0.3e1_dp - t9 * t4 * t13 * t37 * t28 * t2 / 0.3e1_dp & 00340 - t10 * t14 * t107) + t34 * t189 00341 t195 = t14 * my_ndrho * t25 00342 00343 e_ndrho(ii) = e_ndrho(ii) + 0.55555555555555555556e-1_dp * lambda * a * b * t50 * t195 + t34 & 00344 * (-0.2777777778e-1_dp * t9 * t180 * c * t195 - 0.2777777778e-1_dp * t9 & 00345 * t180 * t13 * t37 * my_ndrho * t25 * d + 0.27777777777777777778e-1_dp * & 00346 t51 * t14 * my_ndrho * t63) 00347 00348 END IF 00349 END IF 00350 END DO 00351 00352 !$omp end do 00353 00354 END SUBROUTINE lyp_adiabatic_lda_calc 00355 00356 ! ***************************************************************************** 00361 SUBROUTINE lyp_adiabatic_lsd_eval(rho_set,deriv_set,grad_deriv,lyp_adiabatic_params,error) 00362 TYPE(xc_rho_set_type), POINTER :: rho_set 00363 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00364 INTEGER, INTENT(in) :: grad_deriv 00365 TYPE(section_vals_type), POINTER :: lyp_adiabatic_params 00366 TYPE(cp_error_type), INTENT(inout) :: error 00367 00368 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_adiabatic_lsd_eval', 00369 routineP = moduleN//':'//routineN 00370 00371 INTEGER :: handle, npoints, stat 00372 INTEGER, DIMENSION(:, :), POINTER :: bo 00373 LOGICAL :: failure 00374 REAL(kind=dp) :: epsilon_drho, epsilon_rho, 00375 lambda 00376 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, 00377 e_ndr_ndr, e_ndr_ra, e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, 00378 e_ndra_rb, e_ndrb, e_ndrb_ndrb, e_ndrb_ra, e_ndrb_rb, e_ra, e_ra_ra, 00379 e_ra_rb, e_rb, e_rb_rb, norm_drho, norm_drhoa, norm_drhob, rhoa, rhob 00380 TYPE(xc_derivative_type), POINTER :: deriv 00381 00382 CALL timeset(routineN,handle) 00383 failure=.FALSE. 00384 NULLIFY(deriv, bo) 00385 00386 CALL section_vals_val_get(lyp_adiabatic_params,"LAMBDA",r_val=lambda,error=error) 00387 CALL cite_reference(Lee1988) 00388 00389 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00390 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00391 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00392 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00393 IF (.NOT. failure) THEN 00394 CALL xc_rho_set_get(rho_set,& 00395 rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, & 00396 norm_drhob=norm_drhob, norm_drho=norm_drho, & 00397 rho_cutoff=epsilon_rho,& 00398 drho_cutoff=epsilon_drho, local_bounds=bo, error=error) 00399 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00400 00401 ! meaningful default for the arrays we don't need: let us make compiler 00402 ! and debugger happy... 00403 IF (cp_debug) THEN 00404 ALLOCATE(dummy(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 00405 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00406 ELSE 00407 dummy=> rhoa 00408 END IF 00409 e_0 => dummy 00410 e_ra => dummy 00411 e_rb => dummy 00412 e_ndra_ra => dummy 00413 e_ndra_rb => dummy 00414 e_ndrb_ra => dummy 00415 e_ndrb_rb => dummy 00416 e_ndr_ndr => dummy 00417 e_ndra_ndra => dummy 00418 e_ndrb_ndrb => dummy 00419 e_ndr => dummy 00420 e_ndra => dummy 00421 e_ndrb => dummy 00422 e_ra_ra => dummy 00423 e_ra_rb => dummy 00424 e_rb_rb => dummy 00425 e_ndr_ra => dummy 00426 e_ndr_rb => dummy 00427 00428 IF (grad_deriv>=0) THEN 00429 deriv => xc_dset_get_derivative(deriv_set,"",& 00430 allocate_deriv=.TRUE., error=error) 00431 CALL xc_derivative_get(deriv, deriv_data=e_0,error=error) 00432 END IF 00433 IF (grad_deriv==1.OR.grad_deriv==-1) THEN 00434 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)",& 00435 allocate_deriv=.TRUE.,error=error) 00436 CALL xc_derivative_get(deriv,deriv_data=e_ra,error=error) 00437 deriv => xc_dset_get_derivative(deriv_set,"(rhob)",& 00438 allocate_deriv=.TRUE.,error=error) 00439 CALL xc_derivative_get(deriv,deriv_data=e_rb,error=error) 00440 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",& 00441 allocate_deriv=.TRUE.,error=error) 00442 CALL xc_derivative_get(deriv,deriv_data=e_ndr,error=error) 00443 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)",& 00444 allocate_deriv=.TRUE.,error=error) 00445 CALL xc_derivative_get(deriv,deriv_data=e_ndra,error=error) 00446 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)",& 00447 allocate_deriv=.TRUE.,error=error) 00448 CALL xc_derivative_get(deriv,deriv_data=e_ndrb,error=error) 00449 END IF 00450 IF (grad_deriv>1.OR.grad_deriv<-1) THEN 00451 CALL cp_unimplemented_error(fromWhere=routineP, & 00452 message="derivatives bigger than 1 not implemented", & 00453 error=error, error_level=cp_failure_level) 00454 END IF 00455 00456 !$omp parallel default(none) & 00457 !$omp shared(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob) & 00458 !$omp shared(e_0, e_ra, e_rb, e_ndr, e_ndra, e_ndrb) & 00459 !$omp shared(grad_deriv, npoints) & 00460 !$omp shared(epsilon_rho, epsilon_drho, lambda, error) 00461 00462 CALL lyp_adiabatic_lsd_calc(& 00463 rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa,& 00464 norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb,& 00465 e_ndr=e_ndr,& 00466 e_ndra=e_ndra, e_ndrb=e_ndrb,& 00467 grad_deriv=grad_deriv, npoints=npoints, & 00468 epsilon_rho=epsilon_rho,epsilon_drho=epsilon_drho,lambda=lambda,error=error) 00469 00470 !$omp end parallel 00471 00472 IF (cp_debug) THEN 00473 DEALLOCATE(dummy,stat=stat) 00474 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00475 ELSE 00476 NULLIFY(dummy) 00477 END IF 00478 END IF 00479 CALL timestop(handle) 00480 END SUBROUTINE lyp_adiabatic_lsd_eval 00481 00482 ! ***************************************************************************** 00487 SUBROUTINE lyp_adiabatic_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob,& 00488 e_0, e_ra, e_rb,& 00489 e_ndr,& 00490 e_ndra, e_ndrb, & 00491 grad_deriv,npoints,epsilon_rho,epsilon_drho,lambda,error) 00492 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob, norm_drho, 00493 norm_drhoa, norm_drhob 00494 REAL(kind=dp), DIMENSION(*), 00495 INTENT(inout) :: e_0, e_ra, e_rb, e_ndr, 00496 e_ndra, e_ndrb 00497 INTEGER, INTENT(in) :: grad_deriv, npoints 00498 REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_drho, 00499 lambda 00500 TYPE(cp_error_type), INTENT(inout) :: error 00501 00502 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_adiabatic_lsd_calc', 00503 routineP = moduleN//':'//routineN 00504 00505 INTEGER :: ii 00506 REAL(KIND=dp) :: cf, my_ndrho, my_ndrhoa, my_ndrhob, my_rhoa, my_rhob, 00507 t1, t10, t100, t102, t103, t106, t108, t113, t115, t118, t119, t124, 00508 t125, t128, t129, t132, t135, t138, t14, t140, t141, t143, t145, t146, 00509 t15, t151, t153, t157, t16, t162, t165, t169, t17, t171, t174, t179, 00510 t18, t183, t186, t187, t188, t19, t194, t196, t199, t2, t200, t202, 00511 t21, t212, t216, t220, t222, t223, t225, t23, t231, t237, t24, t246, 00512 t25, t250, t259, t26, t266, t27, t270, t273, t276, t28, t280, t285, 00513 t288, t294, t3, t30, t300, t307, t31, t316, t32, t325, t348, t351, 00514 t355, t362, t387, t39, t394, t4, t41, t42, t421, t46, t47, t48, t49 00515 REAL(KIND=dp) :: t5, t51, t55, t58, t6, t62, t63, t65, t67, t7, t73, t74, 00516 t76, t77, t78, t80, t83, t84, t85, t86, t87, t9, t90, t91, t94, t95, 00517 t96, t97 00518 00519 cf=0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp) 00520 00521 !$omp do 00522 00523 DO ii=1,npoints 00524 my_rhoa=MAX(rhoa(ii),0.0_dp) 00525 my_rhob=MAX(rhob(ii),0.0_dp) 00526 IF (my_rhoa+my_rhob>epsilon_rho) THEN 00527 my_ndrhoa=norm_drhoa(ii) 00528 my_ndrhob=norm_drhob(ii) 00529 my_ndrho = norm_drho(ii) 00530 IF( grad_deriv >= 0 ) THEN 00531 t1 = a * my_rhoa 00532 t2 = my_rhoa + my_rhob 00533 t3 = 0.1e1_dp / t2 00534 t4 = my_rhob * t3 00535 t5 = d * lambda 00536 t6 = t2 ** (0.1e1_dp / 0.3e1_dp) 00537 t7 = 0.1e1_dp / t6 00538 t9 = 0.10e1_dp + t5 * t7 00539 t10 = 0.1e1_dp / t9 00540 t14 = a * b 00541 t15 = c * lambda 00542 t16 = t15 * t7 00543 t17 = EXP(-t16) 00544 t18 = t14 * t17 00545 t19 = t2 ** 2 00546 t21 = t6 ** 2 00547 t23 = 0.1e1_dp / t21 / t19 / t2 00548 t24 = t10 * t23 00549 t25 = my_rhoa * my_rhob 00550 t26 = my_rhoa ** 2 00551 t27 = my_rhoa ** (0.1e1_dp / 0.3e1_dp) 00552 t28 = t27 ** 2 00553 t30 = my_rhob ** 2 00554 t31 = my_rhob ** (0.1e1_dp / 0.3e1_dp) 00555 t32 = t31 ** 2 00556 t39 = t5 * t7 * t10 00557 t41 = 0.26111111111111111111e1_dp - 0.3888888889e0_dp * t16 - 0.3888888889e0_dp & 00558 * t39 00559 t42 = my_ndrho ** 2 00560 t46 = 0.25000000000000000000e1_dp - 0.5555555556e-1_dp * t16 - 0.5555555556e-1_dp & 00561 * t39 00562 t47 = my_ndrhoa ** 2 00563 t48 = my_ndrhob ** 2 00564 t49 = t47 + t48 00565 t51 = t16 + t39 - 0.110e2_dp 00566 t55 = my_rhoa * t3 * t47 + t4 * t48 00567 t58 = 0.12699208415745595798e2_dp * Cf * (t28 * t26 + t32 * t30) + t41 & 00568 * t42 - t46 * t49 - 0.1111111111e0_dp * t51 * t55 00569 t62 = 0.66666666666666666667e0_dp * t19 00570 t63 = t62 - t26 00571 t65 = t62 - t30 00572 t67 = t25 * t58 - 0.6666666667e0_dp * t19 * t42 + t63 * t48 + t65 * t47 00573 t73 = lambda ** 2 00574 t74 = t1 * my_rhob 00575 t76 = 0.1e1_dp / t6 / t2 00576 t77 = t9 ** 2 00577 t78 = 0.1e1_dp / t77 00578 t80 = t76 * t78 * d 00579 t83 = t14 * c 00580 t84 = t19 ** 2 00581 t85 = 0.1e1_dp / t84 00582 t86 = t85 * t17 00583 t87 = t10 * t67 00584 t90 = t78 * t85 00585 t91 = t67 * d 00586 t94 = t17 * t10 00587 t95 = t14 * t94 00588 t96 = t23 * my_rhoa 00589 t97 = c * t7 00590 t100 = d * t7 * t10 00591 t102 = d ** 2 00592 t103 = t102 * lambda 00593 t106 = t103 / t21 * t78 00594 t108 = -0.3888888889e0_dp * t97 - 0.3888888889e0_dp * t100 + 0.38888888888888888889e0_dp & 00595 * t106 00596 t113 = -0.5555555556e-1_dp * t97 - 0.5555555556e-1_dp * t100 + 0.55555555555555555556e-1_dp & 00597 * t106 00598 t115 = t97 + t100 - t106 00599 t118 = t108 * t42 - t113 * t49 - 0.1111111111e0_dp * t115 * t55 00600 t119 = my_rhob * t118 00601 00602 e_0(ii) = e_0(ii) + 0.20e1_dp * lambda * (-0.40e1_dp * t1 * t4 * t10 - t18 * t24 * t67) & 00603 + t73 * (0.40e1_dp * t74 * t80 + t83 * t86 * t87 + t18 * t90 * t91 - & 00604 t95 * t96 * t119) 00605 00606 END IF 00607 IF( grad_deriv == 1 .OR. grad_deriv == -1 ) THEN 00608 t124 = a * my_rhob 00609 t125 = t3 * t10 00610 t128 = 0.1e1_dp / t19 00611 t129 = my_rhob * t128 00612 t132 = 0.40e1_dp * t1 * t129 * t10 00613 t135 = 0.1e1_dp / t6 / t19 * t78 00614 t138 = 0.1333333333e1_dp * t74 * t135 * t5 00615 t140 = t84 * t2 00616 t141 = 0.1e1_dp / t140 00617 t143 = t141 * t17 * t87 00618 t145 = t14 * t15 * t143 / 0.3e1_dp 00619 t146 = t17 * t78 00620 t151 = t14 * t146 * t141 * t67 * t5 / 0.3e1_dp 00621 t153 = 0.1e1_dp / t21 / t84 00622 t157 = 0.11e2_dp / 0.3e1_dp * t18 * t10 * t153 * t67 00623 t162 = t15 * t76 00624 t165 = t5 * t76 * t10 00625 t169 = 0.1e1_dp / t21 / t2 00626 t171 = t102 * t73 * t169 * t78 00627 t174 = (0.12962962962962962963e0_dp * t162 + 0.12962962962962962963e0_dp & 00628 * t165 - 0.1296296296e0_dp * t171) * t42 00629 t179 = (0.18518518518518518519e-1_dp * t162 + 0.18518518518518518519e-1_dp & 00630 * t165 - 0.1851851852e-1_dp * t171) * t49 00631 t183 = 0.1111111111e0_dp * (-t162 / 0.3e1_dp - t165 / 0.3e1_dp + t171 / 0.3e1_dp) & 00632 * t55 00633 t186 = my_rhoa * t128 * t47 00634 t187 = t129 * t48 00635 t188 = t3 * t47 - t186 - t187 00636 t194 = 0.1333333333e1_dp * t2 * t42 00637 t196 = 0.13333333333333333333e1_dp * my_rhob 00638 t199 = 0.13333333333333333333e1_dp * my_rhoa 00639 t200 = t199 + t196 00640 t202 = my_rhob * t58 + t25 * (0.33864555775321588795e2_dp * Cf * t28 * my_rhoa & 00641 + t174 - t179 - t183 - 0.1111111111e0_dp * t51 * t188) - t194 + (-0.6666666667e0_dp & 00642 * my_rhoa + t196) * t48 + t200 * t47 00643 t212 = 0.5333333333e1_dp * t74 * t135 * d 00644 t216 = 0.1e1_dp / t77 / t9 00645 t220 = 0.26666666666666666667e1_dp * t74 / t21 / t19 * t216 * t103 00646 t222 = 4 * t83 * t143 00647 t223 = c ** 2 00648 t225 = 0.1e1_dp / t6 / t140 00649 t231 = t14 * t223 * t225 * lambda * t17 * t87 / 0.3e1_dp 00650 t237 = 0.2e1_dp / 0.3e1_dp * t14 * c * t225 * t146 * t91 * lambda 00651 t246 = 0.2e1_dp / 0.3e1_dp * t14 * t17 * t216 * t225 * t67 * t103 00652 t250 = 4 * t18 * t78 * t141 * t91 00653 t259 = t14 * t15 * t141 * t94 * t25 * t118 / 0.3e1_dp 00654 t266 = t14 * t146 * t141 * t25 * t118 * d * lambda / 0.3e1_dp 00655 t270 = 0.11e2_dp / 0.3e1_dp * t95 * t153 * my_rhoa * t119 00656 t273 = c * t76 00657 t276 = d * t76 * t10 00658 t280 = t102 * t169 * t78 * lambda 00659 t285 = t102 * d * t73 * t128 * t216 00660 t288 = (0.12962962962962962963e0_dp * t273 + 0.12962962962962962963e0_dp & 00661 * t276 - 0.3888888889e0_dp * t280 + 0.25925925925925925926e0_dp * t285) & 00662 * t42 00663 t294 = (0.18518518518518518519e-1_dp * t273 + 0.18518518518518518519e-1_dp & 00664 * t276 - 0.5555555556e-1_dp * t280 + 0.37037037037037037037e-1_dp * t285) & 00665 * t49 00666 t300 = 0.1111111111e0_dp * (-t273 / 0.3e1_dp - t276 / 0.3e1_dp + t280 - 0.2e1_dp & 00667 / 0.3e1_dp * t285) * t55 00668 t307 = 0.40e1_dp * t124 * t80 - t212 + t220 - t222 + t231 + t237 + t83 & 00669 * t86 * t10 * t202 + t246 - t250 + t18 * t90 * t202 * d - t259 - & 00670 t266 + t270 - t18 * t24 * t119 - t95 * t96 * my_rhob * (t288 - t294 - & 00671 t300 - 0.1111111111e0_dp * t115 * t188) 00672 00673 e_ra(ii) = e_ra(ii) + 0.20e1_dp * lambda * (-0.40e1_dp * t124 * t125 + t132 - t138 - t145 & 00674 - t151 + t157 - t18 * t24 * t202) + t73 * t307 00675 00676 t316 = -t186 + t3 * t48 - t187 00677 t325 = my_rhoa * t58 + t25 * (0.33864555775321588795e2_dp * Cf * t32 * my_rhob & 00678 + t174 - t179 - t183 - 0.1111111111e0_dp * t51 * t316) - t194 + t200 & 00679 * t48 + (t199 - 0.6666666667e0_dp * my_rhob) * t47 00680 t348 = 0.40e1_dp * t1 * t80 - t212 + t220 - t222 + t231 + t237 + t83 * & 00681 t86 * t10 * t325 + t246 - t250 + t18 * t90 * t325 * d - t259 - t266 & 00682 + t270 - t18 * t24 * my_rhoa * t118 - t95 * t96 * my_rhob * (t288 - t294 & 00683 - t300 - 0.1111111111e0_dp * t115 * t316) 00684 00685 e_rb(ii) = e_rb(ii) + 0.20e1_dp * lambda * (-0.40e1_dp * t1 * t125 + t132 - t138 - t145 - & 00686 t151 + t157 - t18 * t24 * t325) + t73 * t348 00687 00688 t351 = lambda * a * b 00689 t355 = t3 * my_ndrhoa 00690 t362 = t25 * (-REAL(2 * t46 * my_ndrhoa,dp) - 0.2222222222e0_dp * t51 * my_rhoa 00691 * t355) + REAL(2 * t65 * my_ndrhoa,dp) 00692 00693 e_ndra(ii) = e_ndra(ii) -0.20e1_dp * t351 * t94 * t23 * t362 + t73 * (t83 * t86 * t10 * & 00694 t362 + t18 * t90 * t362 * d - t95 * t96 * my_rhob * (-REAL(2 * t113 * & my_ndrhoa,dp) - 0.2222222222e0_dp * t115 * my_rhoa * t355)) 00695 00696 t387 = t3 * my_ndrhob 00697 t394 = t25 * (-REAL(2 * t46 * my_ndrhob,dp) - 0.2222222222e0_dp * t51 * my_rhob 00698 * t387) + REAL(2 * t63 * my_ndrhob,dp) 00699 00700 e_ndrb(ii) = e_ndrb(ii) -0.20e1_dp * t351 * t94 * t23 * t394 + t73 * (t83 * t86 * t10 * & 00701 t394 + t18 * t90 * t394 * d - t95 * t96 * my_rhob * (-REAL(2 * t113 * & my_ndrhob,dp) - 0.2222222222e0_dp * t115 * my_rhob * t387)) 00702 00703 t421 = REAL(2 * t25 * t41 * my_ndrho,dp) - 0.1333333333e1_dp * REAL(t19,dp) * REAL(my_ndrho,dp) 00704 00705 e_ndr(ii) = e_ndr(ii) -0.20e1_dp * t351 * t94 * t23 * t421 + t73 * (t83 * t86 * t10 * t421 & 00706 + t18 * t90 * t421 * d - REAL(2 * t95 * t96 * my_rhob * t108 * my_ndrho,dp)) 00707 00708 END IF 00709 END IF 00710 END DO 00711 00712 !$omp end do 00713 00714 END SUBROUTINE lyp_adiabatic_lsd_calc 00715 00716 END MODULE xc_lyp_adiabatic 00717
1.7.3