|
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 ! ***************************************************************************** 00014 MODULE xc_tfw 00015 USE cp_array_r_utils, ONLY: cp_3d_r_p_type 00016 USE f77_blas 00017 USE kinds, ONLY: dp 00018 USE timings, ONLY: timeset,& 00019 timestop 00020 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00021 xc_dset_get_derivative 00022 USE xc_derivative_types, ONLY: xc_derivative_get,& 00023 xc_derivative_type 00024 USE xc_functionals_utilities, ONLY: set_util 00025 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00026 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 00027 xc_rho_set_type 00028 #include "cp_common_uses.h" 00029 00030 IMPLICIT NONE 00031 00032 PRIVATE 00033 00034 ! *** Global parameters *** 00035 00036 REAL(KIND=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp 00037 REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, 00038 f23 = 2.0_dp*f13, 00039 f43 = 4.0_dp*f13, 00040 f53 = 5.0_dp*f13 00041 00042 PUBLIC :: tfw_lda_info, tfw_lda_eval, tfw_lsd_info, tfw_lsd_eval 00043 00044 REAL(KIND=dp) :: cf, flda, flsd, fvw 00045 REAL(KIND=dp) :: eps_rho 00046 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_tfw' 00047 00048 CONTAINS 00049 00050 ! ***************************************************************************** 00051 SUBROUTINE tfw_init ( cutoff ) 00052 00053 REAL(KIND=dp), INTENT(IN) :: cutoff 00054 00055 eps_rho = cutoff 00056 CALL set_util ( cutoff ) 00057 00058 cf = 0.3_dp*(3.0_dp*pi*pi)**f23 00059 flda = cf 00060 flsd = flda * 2.0_dp**f23 00061 fvw = 1.0_dp/72.0_dp 00062 00063 END SUBROUTINE tfw_init 00064 00065 ! ***************************************************************************** 00066 SUBROUTINE tfw_lda_info ( reference, shortform, needs, max_deriv, error) 00067 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00068 TYPE(xc_rho_cflags_type), 00069 INTENT(inout), OPTIONAL :: needs 00070 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00071 TYPE(cp_error_type), INTENT(inout) :: error 00072 00073 IF ( PRESENT ( reference ) ) THEN 00074 reference = "Thomas-Fermi-Weizsaecker kinetic energy functional {LDA version}" 00075 END IF 00076 IF ( PRESENT ( shortform ) ) THEN 00077 shortform = "TF+vW kinetic energy functional {LDA}" 00078 END IF 00079 IF (PRESENT(needs)) THEN 00080 needs%rho=.TRUE. 00081 needs%rho_1_3=.TRUE. 00082 needs%norm_drho=.TRUE. 00083 END IF 00084 IF (PRESENT(max_deriv)) max_deriv=3 00085 00086 END SUBROUTINE tfw_lda_info 00087 00088 ! ***************************************************************************** 00089 SUBROUTINE tfw_lsd_info ( reference, shortform, needs, max_deriv, error) 00090 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00091 TYPE(xc_rho_cflags_type), 00092 INTENT(inout), OPTIONAL :: needs 00093 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00094 TYPE(cp_error_type), INTENT(inout) :: error 00095 00096 IF ( PRESENT ( reference ) ) THEN 00097 reference = "Thomas-Fermi-Weizsaecker kinetic energy functional" 00098 END IF 00099 IF ( PRESENT ( shortform ) ) THEN 00100 shortform = "TF+vW kinetic energy functional" 00101 END IF 00102 IF (PRESENT(needs)) THEN 00103 needs%rho_spin=.TRUE. 00104 needs%rho_spin_1_3=.TRUE. 00105 needs%norm_drho=.TRUE. 00106 END IF 00107 IF (PRESENT(max_deriv)) max_deriv=3 00108 00109 END SUBROUTINE tfw_lsd_info 00110 00111 ! ***************************************************************************** 00112 SUBROUTINE tfw_lda_eval (rho_set,deriv_set,order,error) 00113 TYPE(xc_rho_set_type), POINTER :: rho_set 00114 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00115 INTEGER, INTENT(in) :: order 00116 TYPE(cp_error_type), INTENT(inout) :: error 00117 00118 CHARACTER(len=*), PARAMETER :: routineN = 'tfw_lda_eval', 00119 routineP = moduleN//':'//routineN 00120 00121 INTEGER :: handle, npoints, stat 00122 INTEGER, DIMENSION(:, :), POINTER :: bo 00123 LOGICAL :: failure 00124 REAL(KIND=dp) :: epsilon_rho 00125 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s 00126 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, 00127 e_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, 00128 e_rho_rho_ndrho, e_rho_rho_rho, grho, r13, rho 00129 TYPE(xc_derivative_type), POINTER :: deriv 00130 00131 CALL timeset(routineN,handle) 00132 failure=.FALSE. 00133 NULLIFY(bo) 00134 00135 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00136 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00137 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00138 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00139 IF (.NOT. failure) THEN 00140 CALL xc_rho_set_get(rho_set,rho_1_3=r13,rho=rho,& 00141 norm_drho=grho,local_bounds=bo,rho_cutoff=epsilon_rho,& 00142 error=error) 00143 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00144 CALL tfw_init(epsilon_rho) 00145 00146 ALLOCATE ( s(npoints), STAT=stat ) 00147 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00148 CALL calc_s(rho,grho,s, npoints) 00149 00150 IF ( order>=0 ) THEN 00151 deriv => xc_dset_get_derivative(deriv_set,"",& 00152 allocate_deriv=.TRUE., error=error) 00153 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00154 00155 CALL tfw_u_0 ( rho, r13, s, e_0, npoints, error ) 00156 END IF 00157 IF ( order>=1.OR.order==-1 ) THEN 00158 deriv => xc_dset_get_derivative(deriv_set,"(rho)",& 00159 allocate_deriv=.TRUE.,error=error) 00160 CALL xc_derivative_get(deriv,deriv_data=e_rho,error=error) 00161 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",& 00162 allocate_deriv=.TRUE.,error=error) 00163 CALL xc_derivative_get(deriv,deriv_data=e_ndrho,error=error) 00164 00165 CALL tfw_u_1 ( rho, grho, r13, s, e_rho, e_ndrho, npoints, error ) 00166 END IF 00167 IF ( order>=2.OR.order==-2 ) THEN 00168 deriv => xc_dset_get_derivative(deriv_set,"(rho)(rho)",& 00169 allocate_deriv=.TRUE.,error=error) 00170 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho,error=error) 00171 deriv => xc_dset_get_derivative(deriv_set,"(rho)(norm_drho)",& 00172 allocate_deriv=.TRUE.,error=error) 00173 CALL xc_derivative_get(deriv,deriv_data=e_rho_ndrho,error=error) 00174 deriv => xc_dset_get_derivative(deriv_set,& 00175 "(norm_drho)(norm_drho)", allocate_deriv=.TRUE.,error=error) 00176 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho,error=error) 00177 00178 CALL tfw_u_2 ( rho, grho, r13, s, e_rho_rho, e_rho_ndrho,& 00179 e_ndrho_ndrho, npoints, error ) 00180 END IF 00181 IF ( order>=3.OR.order==-3 ) THEN 00182 deriv => xc_dset_get_derivative(deriv_set,"(rho)(rho)(rho)",& 00183 allocate_deriv=.TRUE.,error=error) 00184 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho_rho,error=error) 00185 deriv => xc_dset_get_derivative(deriv_set,& 00186 "(rho)(rho)(norm_drho)",allocate_deriv=.TRUE.,error=error) 00187 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho_ndrho,error=error) 00188 deriv => xc_dset_get_derivative(deriv_set,& 00189 "(rho)(norm_drho)(norm_drho)",allocate_deriv=.TRUE.,error=error) 00190 CALL xc_derivative_get(deriv,deriv_data=e_rho_ndrho_ndrho,error=error) 00191 00192 CALL tfw_u_3 ( rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho,& 00193 e_rho_ndrho_ndrho, npoints, error ) 00194 END IF 00195 IF ( order>3.OR.order<-3) THEN 00196 CALL cp_unimplemented_error(fromWhere=routineP, & 00197 message="derivatives bigger than 3 not implemented", & 00198 error=error, error_level=cp_failure_level) 00199 END IF 00200 00201 DEALLOCATE ( s, STAT=stat ) 00202 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00203 END IF 00204 CALL timestop(handle) 00205 END SUBROUTINE tfw_lda_eval 00206 00207 ! ***************************************************************************** 00208 SUBROUTINE calc_s(rho,grho,s, npoints) 00209 REAL(KIND=dp), DIMENSION(*), INTENT(in) :: rho, grho 00210 REAL(KIND=dp), DIMENSION(*), INTENT(out) :: s 00211 INTEGER, INTENT(in) :: npoints 00212 00213 INTEGER :: ip 00214 00215 !$omp parallel do private(ip) 00216 DO ip = 1, npoints 00217 IF ( rho(ip) < eps_rho ) THEN 00218 s(ip) = 0.0_dp 00219 ELSE 00220 s(ip) = grho(ip)*grho(ip) / rho(ip) 00221 END IF 00222 END DO 00223 END SUBROUTINE calc_s 00224 00225 ! ***************************************************************************** 00226 SUBROUTINE tfw_lsd_eval(rho_set,deriv_set,order,error) 00227 TYPE(xc_rho_set_type), POINTER :: rho_set 00228 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00229 INTEGER, INTENT(in) :: order 00230 TYPE(cp_error_type), INTENT(inout) :: error 00231 00232 CHARACTER(len=*), PARAMETER :: routineN = 'tfw_lsd_eval', 00233 routineP = moduleN//':'//routineN 00234 00235 CHARACTER(len=12), DIMENSION(2) :: 00236 norm_drho_spin_name = (/"(norm_drhoa)","(norm_drhob)"/) 00237 CHARACTER(len=6), DIMENSION(2) :: rho_spin_name = (/"(rhoa)","(rhob)"/) 00238 INTEGER :: handle, i, ispin, npoints, 00239 stat 00240 INTEGER, DIMENSION(:, :), POINTER :: bo 00241 LOGICAL :: failure 00242 REAL(KIND=dp) :: epsilon_rho 00243 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s 00244 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, 00245 e_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, 00246 e_rho_rho_ndrho, e_rho_rho_rho 00247 TYPE(cp_3d_r_p_type), DIMENSION(2) :: norm_drho, rho, rho_1_3 00248 TYPE(xc_derivative_type), POINTER :: deriv 00249 00250 CALL timeset(routineN,handle) 00251 failure=.FALSE. 00252 NULLIFY(deriv, bo) 00253 DO i=1,2 00254 NULLIFY(norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array) 00255 END DO 00256 00257 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00258 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00259 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00260 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00261 IF (.NOT. failure) THEN 00262 CALL xc_rho_set_get(rho_set,rhoa_1_3=rho_1_3(1)%array,& 00263 rhob_1_3=rho_1_3(2)%array,rhoa=rho(1)%array,& 00264 rhob=rho(2)%array,norm_drhoa=norm_drho(1)%array, & 00265 norm_drhob=norm_drho(2)%array,rho_cutoff=epsilon_rho,& 00266 local_bounds=bo, error=error) 00267 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00268 CALL tfw_init(epsilon_rho) 00269 00270 ALLOCATE ( s(npoints), STAT=stat ) 00271 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00272 00273 DO ispin=1,2 00274 CALL calc_s(rho(ispin)%array, norm_drho(ispin)%array, s, npoints) 00275 00276 IF ( order>=0 ) THEN 00277 deriv => xc_dset_get_derivative(deriv_set,"",& 00278 allocate_deriv=.TRUE., error=error) 00279 CALL xc_derivative_get(deriv, deriv_data=e_0,error=error) 00280 00281 CALL tfw_p_0 ( rho(ispin)%array, norm_drho(ispin)%array, & 00282 rho_1_3(ispin)%array, s, e_0, npoints,error ) 00283 END IF 00284 IF ( order>=1.OR.order==-1 ) THEN 00285 deriv => xc_dset_get_derivative(deriv_set,rho_spin_name(ispin),& 00286 allocate_deriv=.TRUE.,error=error) 00287 CALL xc_derivative_get(deriv,deriv_data=e_rho,error=error) 00288 deriv => xc_dset_get_derivative(deriv_set,norm_drho_spin_name(ispin),& 00289 allocate_deriv=.TRUE.,error=error) 00290 CALL xc_derivative_get(deriv,deriv_data=e_ndrho,error=error) 00291 00292 CALL tfw_p_1 ( rho(ispin)%array, norm_drho(ispin)%array, & 00293 rho_1_3(ispin)%array, s, e_rho, e_ndrho, npoints,error ) 00294 END IF 00295 IF ( order>=2.OR.order==-2 ) THEN 00296 deriv => xc_dset_get_derivative(deriv_set,rho_spin_name(ispin)//& 00297 rho_spin_name(ispin),allocate_deriv=.TRUE.,error=error) 00298 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho,error=error) 00299 deriv => xc_dset_get_derivative(deriv_set,rho_spin_name(ispin)//& 00300 norm_drho_spin_name(ispin),allocate_deriv=.TRUE.,error=error) 00301 CALL xc_derivative_get(deriv,deriv_data=e_rho_ndrho,error=error) 00302 deriv => xc_dset_get_derivative(deriv_set,norm_drho_spin_name(ispin)//& 00303 norm_drho_spin_name(ispin), allocate_deriv=.TRUE.,error=error) 00304 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho,error=error) 00305 00306 CALL tfw_p_2 ( rho(ispin)%array, norm_drho(ispin)%array, & 00307 rho_1_3(ispin)%array, s, e_rho_rho, e_rho_ndrho,& 00308 e_ndrho_ndrho, npoints,error ) 00309 END IF 00310 IF ( order>=3 .OR. order==-3 ) THEN 00311 deriv => xc_dset_get_derivative(deriv_set,rho_spin_name(ispin)//& 00312 rho_spin_name(ispin)//rho_spin_name(ispin),& 00313 allocate_deriv=.TRUE.,error=error) 00314 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho_rho,error=error) 00315 deriv => xc_dset_get_derivative(deriv_set,rho_spin_name(ispin)//& 00316 rho_spin_name(ispin)//norm_drho_spin_name(ispin),& 00317 allocate_deriv=.TRUE.,error=error) 00318 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho_ndrho,error=error) 00319 deriv => xc_dset_get_derivative(deriv_set,rho_spin_name(ispin)//& 00320 norm_drho_spin_name(ispin)//norm_drho_spin_name(ispin), & 00321 allocate_deriv=.TRUE.,error=error) 00322 CALL xc_derivative_get(deriv,deriv_data=e_rho_ndrho_ndrho,error=error) 00323 00324 CALL tfw_p_3 ( rho(ispin)%array, norm_drho(ispin)%array, & 00325 rho_1_3(ispin)%array, s, e_rho_rho_rho, e_rho_rho_ndrho,& 00326 e_rho_ndrho_ndrho, npoints,error ) 00327 END IF 00328 IF ( order>3.OR.order<-3) THEN 00329 CALL cp_unimplemented_error(fromWhere=routineP, & 00330 message="derivatives bigger than 3 not implemented", & 00331 error=error, error_level=cp_failure_level) 00332 END IF 00333 END DO 00334 00335 DEALLOCATE ( s, STAT=stat ) 00336 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00337 END IF 00338 CALL timestop(handle) 00339 END SUBROUTINE tfw_lsd_eval 00340 00341 ! ***************************************************************************** 00342 SUBROUTINE tfw_u_0 ( rho, r13, s, e_0, npoints, error ) 00343 00344 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s 00345 REAL(KIND=dp), DIMENSION(*), 00346 INTENT(INOUT) :: e_0 00347 INTEGER, INTENT(in) :: npoints 00348 TYPE(cp_error_type), INTENT(inout) :: error 00349 00350 INTEGER :: ip 00351 00352 !$omp parallel do private(ip) 00353 DO ip = 1, npoints 00354 00355 IF ( rho(ip) > eps_rho ) THEN 00356 00357 e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip) + fvw * s(ip) 00358 00359 END IF 00360 00361 END DO 00362 00363 END SUBROUTINE tfw_u_0 00364 00365 ! ***************************************************************************** 00366 SUBROUTINE tfw_u_1 ( rho, grho, r13, s, e_rho, e_ndrho, npoints, error ) 00367 00368 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s 00369 REAL(KIND=dp), DIMENSION(*), 00370 INTENT(INOUT) :: e_rho, e_ndrho 00371 INTEGER, INTENT(in) :: npoints 00372 TYPE(cp_error_type), INTENT(inout) :: error 00373 00374 INTEGER :: ip 00375 REAL(KIND=dp) :: f 00376 00377 f = f53 * flda 00378 00379 !$omp parallel do private(ip) 00380 DO ip = 1, npoints 00381 00382 IF ( rho(ip) > eps_rho ) THEN 00383 00384 e_rho(ip) = e_rho(ip) + f * r13(ip)*r13(ip) - fvw * s(ip)/rho(ip) 00385 e_ndrho(ip) = e_ndrho(ip) + 2.0_dp * fvw * grho(ip)/rho(ip) 00386 00387 END IF 00388 00389 END DO 00390 00391 END SUBROUTINE tfw_u_1 00392 00393 ! ***************************************************************************** 00394 SUBROUTINE tfw_u_2 ( rho, grho, r13, s, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho,& 00395 npoints, error) 00396 00397 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s 00398 REAL(KIND=dp), DIMENSION(*), 00399 INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, 00400 e_ndrho_ndrho 00401 INTEGER, INTENT(in) :: npoints 00402 TYPE(cp_error_type), INTENT(inout) :: error 00403 00404 INTEGER :: ip 00405 REAL(KIND=dp) :: f 00406 00407 f = f23 * f53 * flda 00408 00409 !$omp parallel do private(ip) 00410 DO ip = 1, npoints 00411 00412 IF ( rho(ip) > eps_rho ) THEN 00413 00414 e_rho_rho(ip) = e_rho_rho(ip) + f / r13(ip) + 2.0_dp * fvw * s(ip)/(rho(ip)*rho(ip)) 00415 e_rho_ndrho(ip) = e_rho_ndrho(ip) - 2.0_dp * fvw * grho(ip)/(rho(ip)*rho(ip)) 00416 e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp * fvw / rho(ip) 00417 00418 END IF 00419 00420 END DO 00421 00422 END SUBROUTINE tfw_u_2 00423 00424 ! ***************************************************************************** 00425 SUBROUTINE tfw_u_3 ( rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho,& 00426 e_rho_ndrho_ndrho, npoints, error) 00427 00428 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s 00429 REAL(KIND=dp), DIMENSION(*), 00430 INTENT(INOUT) :: e_rho_rho_rho, 00431 e_rho_rho_ndrho, 00432 e_rho_ndrho_ndrho 00433 INTEGER, INTENT(in) :: npoints 00434 TYPE(cp_error_type), INTENT(inout) :: error 00435 00436 INTEGER :: ip 00437 REAL(KIND=dp) :: f 00438 00439 f = -f13 * f23 * f53 * flda 00440 00441 !$omp parallel do private(ip) 00442 DO ip = 1, npoints 00443 00444 IF ( rho(ip) > eps_rho ) THEN 00445 00446 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f / ( r13(ip) * rho(ip) ) & 00447 - 6.0_dp * fvw * s(ip)/(rho(ip)*rho(ip)*rho(ip)) 00448 e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip)& 00449 + 4.0_dp * fvw * grho(ip)/(rho(ip)*rho(ip)*rho(ip)) 00450 e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip)& 00451 - 2.0_dp * fvw / (rho(ip)*rho(ip)) 00452 END IF 00453 00454 END DO 00455 00456 END SUBROUTINE tfw_u_3 00457 00458 ! ***************************************************************************** 00459 SUBROUTINE tfw_p_0 ( rhoa, grhoa, r13a, sa, e_0, npoints, error ) 00460 00461 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa 00462 REAL(KIND=dp), DIMENSION(*), 00463 INTENT(INOUT) :: e_0 00464 INTEGER, INTENT(in) :: npoints 00465 TYPE(cp_error_type), INTENT(inout) :: error 00466 00467 INTEGER :: ip 00468 00469 !$omp parallel do private(ip) 00470 DO ip = 1, npoints 00471 00472 IF ( rhoa(ip) > eps_rho ) THEN 00473 e_0(ip) = e_0(ip) + flsd * r13a(ip) * r13a(ip) * rhoa(ip) + fvw * sa(ip) 00474 END IF 00475 00476 END DO 00477 00478 END SUBROUTINE tfw_p_0 00479 00480 ! ***************************************************************************** 00481 SUBROUTINE tfw_p_1 ( rhoa, grhoa, r13a, sa, e_rho, e_ndrho, npoints,& 00482 error ) 00483 00484 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa 00485 REAL(KIND=dp), DIMENSION(*), 00486 INTENT(INOUT) :: e_rho, e_ndrho 00487 INTEGER, INTENT(in) :: npoints 00488 TYPE(cp_error_type), INTENT(inout) :: error 00489 00490 INTEGER :: ip 00491 REAL(KIND=dp) :: f 00492 00493 f = f53 * flsd 00494 00495 !$omp parallel do private(ip) 00496 DO ip = 1, npoints 00497 00498 IF ( rhoa(ip) > eps_rho ) THEN 00499 e_rho(ip) = e_rho(ip) + f * r13a(ip)*r13a(ip) - fvw * sa(ip)/rhoa(ip) 00500 e_ndrho(ip) = e_ndrho(ip) + 2.0_dp * fvw * grhoa(ip)/rhoa(ip) 00501 END IF 00502 00503 END DO 00504 00505 END SUBROUTINE tfw_p_1 00506 00507 ! ***************************************************************************** 00508 SUBROUTINE tfw_p_2 ( rhoa, grhoa, r13a, sa, e_rho_rho, e_rho_ndrho,& 00509 e_ndrho_ndrho, npoints, error ) 00510 00511 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa 00512 REAL(KIND=dp), DIMENSION(*), 00513 INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, 00514 e_ndrho_ndrho 00515 INTEGER, INTENT(in) :: npoints 00516 TYPE(cp_error_type), INTENT(inout) :: error 00517 00518 INTEGER :: ip 00519 REAL(KIND=dp) :: f 00520 00521 f = f23 * f53 * flsd 00522 00523 !$omp parallel do private(ip) 00524 DO ip = 1, npoints 00525 00526 IF ( rhoa(ip) > eps_rho ) THEN 00527 e_rho_rho(ip) = e_rho_rho(ip)& 00528 + f / r13a(ip) + 2.0_dp * fvw * sa(ip)/(rhoa(ip)*rhoa(ip)) 00529 e_rho_ndrho(ip) = e_rho_ndrho(ip)& 00530 - 2.0_dp * fvw * grhoa(ip)/(rhoa(ip)*rhoa(ip)) 00531 e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp * fvw / rhoa(ip) 00532 END IF 00533 00534 END DO 00535 00536 END SUBROUTINE tfw_p_2 00537 00538 ! ***************************************************************************** 00539 SUBROUTINE tfw_p_3 ( rhoa, grhoa, r13a, sa, e_rho_rho_rho, e_rho_rho_ndrho,& 00540 e_rho_ndrho_ndrho, npoints, error ) 00541 00542 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa 00543 REAL(KIND=dp), DIMENSION(*), 00544 INTENT(INOUT) :: e_rho_rho_rho, 00545 e_rho_rho_ndrho, 00546 e_rho_ndrho_ndrho 00547 INTEGER, INTENT(in) :: npoints 00548 TYPE(cp_error_type), INTENT(inout) :: error 00549 00550 INTEGER :: ip 00551 REAL(KIND=dp) :: f 00552 00553 f = -f13 * f23 * f53 * flsd 00554 00555 !$omp parallel do private(ip) 00556 DO ip = 1, npoints 00557 00558 IF ( rhoa(ip) > eps_rho ) THEN 00559 e_rho_rho_rho(ip) = e_rho_rho_rho(ip)& 00560 + f / ( r13a(ip) * rhoa(ip) ) & 00561 - 6.0_dp * fvw * sa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip)) 00562 e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip)& 00563 + 4.0_dp * fvw * grhoa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip)) 00564 e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip)& 00565 - 2.0_dp * fvw / (rhoa(ip)*rhoa(ip)) 00566 END IF 00567 00568 END DO 00569 00570 END SUBROUTINE tfw_p_3 00571 00572 END MODULE xc_tfw 00573
1.7.3