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