|
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 ! ***************************************************************************** 00016 MODULE xc_thomas_fermi 00017 USE cp_array_r_utils, ONLY: cp_3d_r_p_type 00018 USE f77_blas 00019 USE kinds, ONLY: dp 00020 USE timings, ONLY: timeset,& 00021 timestop 00022 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00023 xc_dset_get_derivative 00024 USE xc_derivative_types, ONLY: xc_derivative_get,& 00025 xc_derivative_type 00026 USE xc_functionals_utilities, ONLY: set_util 00027 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00028 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 00029 xc_rho_set_type 00030 #include "cp_common_uses.h" 00031 00032 IMPLICIT NONE 00033 00034 PRIVATE 00035 00036 ! *** Global parameters *** 00037 00038 REAL(KIND=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp 00039 REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, 00040 f23 = 2.0_dp*f13, 00041 f43 = 4.0_dp*f13, 00042 f53 = 5.0_dp*f13 00043 00044 PUBLIC :: thomas_fermi_info, thomas_fermi_lda_eval, thomas_fermi_lsd_eval 00045 00046 REAL(KIND=dp) :: cf, flda, flsd 00047 REAL(KIND=dp) :: eps_rho 00048 00049 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_thomas_fermi' 00050 00051 CONTAINS 00052 00053 ! ***************************************************************************** 00054 SUBROUTINE thomas_fermi_init ( cutoff ) 00055 00056 REAL(KIND=dp), INTENT(IN) :: cutoff 00057 00058 eps_rho = cutoff 00059 CALL set_util ( cutoff ) 00060 00061 cf = 0.3_dp*(3.0_dp*pi*pi)**f23 00062 flda = cf 00063 flsd = flda * 2.0_dp**f23 00064 00065 END SUBROUTINE thomas_fermi_init 00066 00067 ! ***************************************************************************** 00068 SUBROUTINE thomas_fermi_info ( lsd, reference, shortform, needs, max_deriv, error) 00069 LOGICAL, INTENT(in) :: lsd 00070 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00071 TYPE(xc_rho_cflags_type), 00072 INTENT(inout), OPTIONAL :: needs 00073 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00074 TYPE(cp_error_type), INTENT(inout) :: error 00075 00076 IF ( PRESENT ( reference ) ) THEN 00077 reference = "Thomas-Fermi kinetic energy functional: see Parr and Yang" 00078 IF (.not.lsd) THEN 00079 IF (LEN_TRIM(reference)+6<LEN(reference)) THEN 00080 reference(LEN_TRIM(reference):LEN_TRIM(reference)+6)=' {LDA}' 00081 END IF 00082 END IF 00083 END IF 00084 IF ( PRESENT ( shortform ) ) THEN 00085 shortform = "Thomas-Fermi kinetic energy functional" 00086 IF (.not.lsd) THEN 00087 IF (LEN_TRIM(shortform)+6<LEN(shortform)) THEN 00088 shortform(LEN_TRIM(shortform):LEN_TRIM(shortform)+6)=' {LDA}' 00089 END IF 00090 END IF 00091 END IF 00092 IF (PRESENT(needs)) THEN 00093 IF (lsd) THEN 00094 needs%rho_spin=.TRUE. 00095 needs%rho_spin_1_3=.TRUE. 00096 ELSE 00097 needs%rho=.TRUE. 00098 needs%rho_1_3=.TRUE. 00099 END IF 00100 END IF 00101 IF (PRESENT(max_deriv)) max_deriv=3 00102 00103 END SUBROUTINE thomas_fermi_info 00104 00105 ! ***************************************************************************** 00106 SUBROUTINE thomas_fermi_lda_eval(rho_set,deriv_set,order,error) 00107 TYPE(xc_rho_set_type), POINTER :: rho_set 00108 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00109 INTEGER, INTENT(in) :: order 00110 TYPE(cp_error_type), INTENT(inout) :: error 00111 00112 CHARACTER(len=*), PARAMETER :: routineN = 'thomas_fermi_lda_eval', 00113 routineP = moduleN//':'//routineN 00114 00115 INTEGER :: handle, npoints 00116 INTEGER, DIMENSION(:, :), POINTER :: bo 00117 LOGICAL :: failure 00118 REAL(KIND=dp) :: epsilon_rho 00119 REAL(KIND=dp), DIMENSION(:, :, :), 00120 POINTER :: e_0, e_rho, e_rho_rho, 00121 e_rho_rho_rho, r13, rho 00122 TYPE(xc_derivative_type), POINTER :: deriv 00123 00124 CALL timeset(routineN,handle) 00125 failure=.FALSE. 00126 NULLIFY(bo) 00127 00128 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00129 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00130 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00131 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00132 IF (.NOT. failure) THEN 00133 CALL xc_rho_set_get(rho_set,rho_1_3=r13,rho=rho,& 00134 local_bounds=bo,rho_cutoff=epsilon_rho,& 00135 error=error) 00136 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00137 CALL thomas_fermi_init(epsilon_rho) 00138 00139 IF ( order>=0 ) THEN 00140 deriv => xc_dset_get_derivative(deriv_set,"",& 00141 allocate_deriv=.TRUE., error=error) 00142 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00143 00144 CALL thomas_fermi_lda_0 ( rho, r13, e_0,npoints,error ) 00145 END IF 00146 IF ( order>=1.OR.order==-1 ) THEN 00147 deriv => xc_dset_get_derivative(deriv_set,"(rho)",& 00148 allocate_deriv=.TRUE.,error=error) 00149 CALL xc_derivative_get(deriv,deriv_data=e_rho,error=error) 00150 00151 CALL thomas_fermi_lda_1 ( rho, r13, e_rho, npoints, error ) 00152 END IF 00153 IF ( order>=2.OR.order==-2 ) THEN 00154 deriv => xc_dset_get_derivative(deriv_set,"(rho)(rho)",& 00155 allocate_deriv=.TRUE.,error=error) 00156 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho,error=error) 00157 00158 CALL thomas_fermi_lda_2 ( rho, r13, e_rho_rho, npoints, error ) 00159 END IF 00160 IF ( order>=3.OR.order==-3 ) THEN 00161 deriv => xc_dset_get_derivative(deriv_set,"(rho)(rho)(rho)",& 00162 allocate_deriv=.TRUE.,error=error) 00163 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho_rho,error=error) 00164 00165 CALL thomas_fermi_lda_3 ( rho, r13, e_rho_rho_rho, npoints,& 00166 error) 00167 END IF 00168 IF ( order>3.OR.order<-3) THEN 00169 CALL cp_unimplemented_error(fromWhere=routineP, & 00170 message="derivatives bigger than 3 not implemented", & 00171 error=error, error_level=cp_failure_level) 00172 END IF 00173 END IF 00174 CALL timestop(handle) 00175 END SUBROUTINE thomas_fermi_lda_eval 00176 00177 ! ***************************************************************************** 00178 SUBROUTINE thomas_fermi_lsd_eval(rho_set,deriv_set,order,error) 00179 TYPE(xc_rho_set_type), POINTER :: rho_set 00180 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00181 INTEGER, INTENT(in) :: order 00182 TYPE(cp_error_type), INTENT(inout) :: error 00183 00184 CHARACTER(len=*), PARAMETER :: routineN = 'thomas_fermi_lsd_eval', 00185 routineP = moduleN//':'//routineN 00186 00187 CHARACTER(len=6), DIMENSION(2) :: rho_spin_name = (/"(rhoa)","(rhob)"/) 00188 INTEGER :: handle, i, ispin, npoints 00189 INTEGER, DIMENSION(:, :), POINTER :: bo 00190 LOGICAL :: failure 00191 REAL(KIND=dp) :: epsilon_rho 00192 REAL(KIND=dp), DIMENSION(:, :, :), 00193 POINTER :: e_0, e_rho, e_rho_rho, 00194 e_rho_rho_rho 00195 TYPE(cp_3d_r_p_type), DIMENSION(2) :: rho, rho_1_3 00196 TYPE(xc_derivative_type), POINTER :: deriv 00197 00198 CALL timeset(routineN,handle) 00199 failure=.FALSE. 00200 NULLIFY(deriv, bo) 00201 DO i=1,2 00202 NULLIFY(rho(i)%array, rho_1_3(i)%array) 00203 END DO 00204 00205 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00206 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00207 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00208 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00209 IF (.NOT. failure) THEN 00210 00211 CALL xc_rho_set_get(rho_set,rhoa_1_3=rho_1_3(1)%array,& 00212 rhob_1_3=rho_1_3(2)%array,rhoa=rho(1)%array,& 00213 rhob=rho(2)%array,& 00214 rho_cutoff=epsilon_rho,& 00215 local_bounds=bo, error=error) 00216 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00217 CALL thomas_fermi_init(epsilon_rho) 00218 00219 DO ispin=1,2 00220 IF ( order>=0 ) THEN 00221 deriv => xc_dset_get_derivative(deriv_set,"",& 00222 allocate_deriv=.TRUE., error=error) 00223 CALL xc_derivative_get(deriv, deriv_data=e_0,error=error) 00224 00225 CALL thomas_fermi_lsd_0 ( rho(ispin)%array, rho_1_3(ispin)%array,& 00226 e_0, npoints, error ) 00227 END IF 00228 IF ( order>=1.OR.order==-1 ) THEN 00229 deriv => xc_dset_get_derivative(deriv_set,rho_spin_name(ispin),& 00230 allocate_deriv=.TRUE.,error=error) 00231 CALL xc_derivative_get(deriv,deriv_data=e_rho,error=error) 00232 00233 CALL thomas_fermi_lsd_1( rho(ispin)%array, rho_1_3(ispin)%array,& 00234 e_rho, npoints, error ) 00235 END IF 00236 IF ( order>=2.OR.order==-2 ) THEN 00237 deriv => xc_dset_get_derivative(deriv_set,rho_spin_name(ispin)//& 00238 rho_spin_name(ispin),allocate_deriv=.TRUE.,error=error) 00239 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho,error=error) 00240 00241 CALL thomas_fermi_lsd_2( rho(ispin)%array, rho_1_3(ispin)%array,& 00242 e_rho_rho, npoints, error ) 00243 END IF 00244 IF ( order>=3 .OR. order==-3 ) THEN 00245 deriv => xc_dset_get_derivative(deriv_set,rho_spin_name(ispin)//& 00246 rho_spin_name(ispin)//rho_spin_name(ispin),& 00247 allocate_deriv=.TRUE.,error=error) 00248 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho_rho,error=error) 00249 00250 CALL thomas_fermi_lsd_3( rho(ispin)%array, rho_1_3(ispin)%array,& 00251 e_rho_rho_rho, npoints, error ) 00252 END IF 00253 IF ( order>3.OR.order<-3) THEN 00254 CALL cp_unimplemented_error(fromWhere=routineP, & 00255 message="derivatives bigger than 3 not implemented", & 00256 error=error, error_level=cp_failure_level) 00257 END IF 00258 END DO 00259 END IF 00260 CALL timestop(handle) 00261 END SUBROUTINE thomas_fermi_lsd_eval 00262 00263 ! ***************************************************************************** 00264 SUBROUTINE thomas_fermi_lda_0 ( rho, r13, e_0, npoints, error ) 00265 00266 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13 00267 REAL(KIND=dp), DIMENSION(*), 00268 INTENT(INOUT) :: e_0 00269 INTEGER, INTENT(in) :: npoints 00270 TYPE(cp_error_type), INTENT(inout) :: error 00271 00272 INTEGER :: ip 00273 00274 !$omp parallel do private(ip) 00275 DO ip = 1, npoints 00276 00277 IF ( rho(ip) > eps_rho ) THEN 00278 00279 e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip) 00280 00281 END IF 00282 00283 END DO 00284 00285 END SUBROUTINE thomas_fermi_lda_0 00286 00287 ! ***************************************************************************** 00288 SUBROUTINE thomas_fermi_lda_1 ( rho, r13, e_rho, npoints, error ) 00289 00290 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13 00291 REAL(KIND=dp), DIMENSION(*), 00292 INTENT(INOUT) :: e_rho 00293 INTEGER, INTENT(in) :: npoints 00294 TYPE(cp_error_type), INTENT(inout) :: error 00295 00296 INTEGER :: ip 00297 REAL(KIND=dp) :: f 00298 00299 f = f53 * flda 00300 00301 !$omp parallel do private(ip) 00302 DO ip = 1, npoints 00303 00304 IF ( rho(ip) > eps_rho ) THEN 00305 00306 e_rho(ip) = e_rho(ip) + f * r13(ip)*r13(ip) 00307 00308 END IF 00309 00310 END DO 00311 00312 END SUBROUTINE thomas_fermi_lda_1 00313 00314 ! ***************************************************************************** 00315 SUBROUTINE thomas_fermi_lda_2 ( rho, r13, e_rho_rho, npoints, error ) 00316 00317 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13 00318 REAL(KIND=dp), DIMENSION(*), 00319 INTENT(INOUT) :: e_rho_rho 00320 INTEGER, INTENT(in) :: npoints 00321 TYPE(cp_error_type), INTENT(inout) :: error 00322 00323 INTEGER :: ip 00324 REAL(KIND=dp) :: f 00325 00326 f = f23 * f53 * flda 00327 00328 !$omp parallel do private(ip) 00329 DO ip = 1, npoints 00330 00331 IF ( rho(ip) > eps_rho ) THEN 00332 00333 e_rho_rho(ip) = e_rho_rho(ip) + f / r13(ip) 00334 00335 END IF 00336 00337 END DO 00338 00339 END SUBROUTINE thomas_fermi_lda_2 00340 00341 ! ***************************************************************************** 00342 SUBROUTINE thomas_fermi_lda_3 ( rho, r13, e_rho_rho_rho, npoints, error ) 00343 00344 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13 00345 REAL(KIND=dp), DIMENSION(*), 00346 INTENT(INOUT) :: e_rho_rho_rho 00347 INTEGER, INTENT(in) :: npoints 00348 TYPE(cp_error_type), INTENT(inout) :: error 00349 00350 INTEGER :: ip 00351 REAL(KIND=dp) :: f 00352 00353 f = -f13 * f23 * f53 * flda 00354 00355 !$omp parallel do private(ip) 00356 DO ip = 1, npoints 00357 00358 IF ( rho(ip) > eps_rho ) THEN 00359 00360 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f / ( r13(ip) * rho(ip) ) 00361 00362 END IF 00363 00364 END DO 00365 00366 END SUBROUTINE thomas_fermi_lda_3 00367 00368 ! ***************************************************************************** 00369 SUBROUTINE thomas_fermi_lsd_0 ( rhoa, r13a, e_0, npoints, error ) 00370 00371 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a 00372 REAL(KIND=dp), DIMENSION(*), 00373 INTENT(INOUT) :: e_0 00374 INTEGER, INTENT(in) :: npoints 00375 TYPE(cp_error_type), INTENT(inout) :: error 00376 00377 INTEGER :: ip 00378 00379 !$omp parallel do private(ip) 00380 DO ip = 1, npoints 00381 00382 IF ( rhoa(ip) > eps_rho ) THEN 00383 e_0(ip) = e_0(ip) + flsd * r13a(ip) * r13a(ip) * rhoa(ip) 00384 END IF 00385 00386 END DO 00387 00388 END SUBROUTINE thomas_fermi_lsd_0 00389 00390 ! ***************************************************************************** 00391 SUBROUTINE thomas_fermi_lsd_1 ( rhoa, r13a, e_rho, npoints, error ) 00392 00393 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a 00394 REAL(KIND=dp), DIMENSION(*), 00395 INTENT(INOUT) :: e_rho 00396 INTEGER, INTENT(in) :: npoints 00397 TYPE(cp_error_type), INTENT(inout) :: error 00398 00399 INTEGER :: ip 00400 REAL(KIND=dp) :: f 00401 00402 f = f53 * flsd 00403 00404 !$omp parallel do private(ip) 00405 DO ip = 1, npoints 00406 00407 IF ( rhoa(ip) > eps_rho ) THEN 00408 e_rho(ip) = e_rho(ip) + f * r13a(ip) * r13a(ip) 00409 END IF 00410 00411 END DO 00412 00413 END SUBROUTINE thomas_fermi_lsd_1 00414 00415 ! ***************************************************************************** 00416 SUBROUTINE thomas_fermi_lsd_2 ( rhoa, r13a, e_rho_rho, npoints, error ) 00417 00418 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a 00419 REAL(KIND=dp), DIMENSION(*), 00420 INTENT(INOUT) :: e_rho_rho 00421 INTEGER, INTENT(in) :: npoints 00422 TYPE(cp_error_type), INTENT(inout) :: error 00423 00424 INTEGER :: ip 00425 REAL(KIND=dp) :: f 00426 00427 f = f23 * f53 * flsd 00428 00429 !$omp parallel do private(ip) 00430 DO ip = 1, npoints 00431 00432 IF ( rhoa(ip) > eps_rho ) THEN 00433 e_rho_rho(ip) = e_rho_rho(ip) + f / r13a(ip) 00434 END IF 00435 00436 END DO 00437 00438 END SUBROUTINE thomas_fermi_lsd_2 00439 00440 ! ***************************************************************************** 00441 SUBROUTINE thomas_fermi_lsd_3 ( rhoa, r13a, e_rho_rho_rho, npoints, error ) 00442 00443 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a 00444 REAL(KIND=dp), DIMENSION(*), 00445 INTENT(INOUT) :: e_rho_rho_rho 00446 INTEGER, INTENT(in) :: npoints 00447 TYPE(cp_error_type), INTENT(inout) :: error 00448 00449 INTEGER :: ip 00450 REAL(KIND=dp) :: f 00451 00452 f = -f13 * f23 * f53 * flsd 00453 00454 !$omp parallel do private(ip) 00455 DO ip = 1, npoints 00456 00457 IF ( rhoa(ip) > eps_rho ) THEN 00458 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f / ( r13a(ip)*rhoa(ip)) 00459 END IF 00460 00461 END DO 00462 00463 END SUBROUTINE thomas_fermi_lsd_3 00464 00465 END MODULE xc_thomas_fermi 00466
1.7.3