|
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 ! ***************************************************************************** 00012 MODULE xc_optx 00013 USE cp_array_r_utils, ONLY: cp_3d_r_p_type 00014 USE f77_blas 00015 USE input_section_types, ONLY: section_vals_type,& 00016 section_vals_val_get 00017 USE kinds, ONLY: dp 00018 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00019 xc_dset_get_derivative 00020 USE xc_derivative_types, ONLY: xc_derivative_get,& 00021 xc_derivative_type 00022 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00023 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 00024 xc_rho_set_type 00025 #include "cp_common_uses.h" 00026 00027 IMPLICIT NONE 00028 PRIVATE 00029 00030 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_optx' 00031 00032 PUBLIC :: optx_lda_info, optx_lda_eval, optx_lsd_info, optx_lsd_eval 00033 CONTAINS 00034 00035 ! ***************************************************************************** 00046 SUBROUTINE optx_lda_info(reference,shortform, needs, max_deriv, error) 00047 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00048 TYPE(xc_rho_cflags_type), 00049 INTENT(inout), OPTIONAL :: needs 00050 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00051 TYPE(cp_error_type), INTENT(inout) :: error 00052 00053 CHARACTER(len=*), PARAMETER :: routineN = 'optx_lda_info', 00054 routineP = moduleN//':'//routineN 00055 00056 IF ( PRESENT ( reference ) ) THEN 00057 reference = "OPTX, Handy NC and Cohen AJ, JCP 116, p. 5411 (2002) (LDA)" 00058 END IF 00059 IF ( PRESENT ( shortform ) ) THEN 00060 shortform = "OPTX exchange (LDA)" 00061 END IF 00062 IF (PRESENT(needs)) THEN 00063 needs%rho=.TRUE. 00064 needs%norm_drho=.TRUE. 00065 END IF 00066 IF (PRESENT(max_deriv)) max_deriv=1 00067 END SUBROUTINE optx_lda_info 00068 00069 ! ***************************************************************************** 00080 SUBROUTINE optx_lsd_info(reference,shortform, needs, max_deriv, error) 00081 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00082 TYPE(xc_rho_cflags_type), 00083 INTENT(inout), OPTIONAL :: needs 00084 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00085 TYPE(cp_error_type), INTENT(inout) :: error 00086 00087 CHARACTER(len=*), PARAMETER :: routineN = 'optx_lsd_info', 00088 routineP = moduleN//':'//routineN 00089 00090 IF ( PRESENT ( reference ) ) THEN 00091 reference = "OPTX, Handy NC and Cohen AJ, JCP 116, p. 5411 (2002), (LSD) " 00092 END IF 00093 IF ( PRESENT ( shortform ) ) THEN 00094 shortform = "OPTX exchange (LSD)" 00095 END IF 00096 IF (PRESENT(needs)) THEN 00097 needs%rho_spin=.TRUE. 00098 needs%norm_drho_spin=.TRUE. 00099 END IF 00100 IF (PRESENT(max_deriv)) max_deriv=1 00101 END SUBROUTINE optx_lsd_info 00102 00103 ! ***************************************************************************** 00118 SUBROUTINE optx_lda_eval(rho_set,deriv_set,grad_deriv,optx_params,error) 00119 TYPE(xc_rho_set_type), POINTER :: rho_set 00120 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00121 INTEGER, INTENT(in) :: grad_deriv 00122 TYPE(section_vals_type), POINTER :: optx_params 00123 TYPE(cp_error_type), INTENT(inout) :: error 00124 00125 CHARACTER(len=*), PARAMETER :: routineN = 'optx_lda_eval', 00126 routineP = moduleN//':'//routineN 00127 00128 INTEGER :: npoints 00129 INTEGER, DIMENSION(:, :), POINTER :: bo 00130 LOGICAL :: failure 00131 REAL(kind=dp) :: epsilon_drho, epsilon_rho, sx 00132 REAL(kind=dp), DIMENSION(:, :, :), 00133 POINTER :: e_0, e_ndrho, e_rho, 00134 norm_drho, rho 00135 TYPE(xc_derivative_type), POINTER :: deriv 00136 00137 failure=.FALSE. 00138 NULLIFY(bo,e_0, e_ndrho, e_rho, norm_drho, rho) 00139 00140 CALL section_vals_val_get(optx_params,"scale_x",r_val=sx,error=error) 00141 00142 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00143 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00144 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00145 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00146 IF (.NOT. failure) THEN 00147 CALL xc_rho_set_get(rho_set,rho=rho,& 00148 norm_drho=norm_drho,local_bounds=bo,rho_cutoff=epsilon_rho,& 00149 drho_cutoff=epsilon_drho,error=error) 00150 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00151 00152 deriv => xc_dset_get_derivative(deriv_set,"",& 00153 allocate_deriv=.TRUE., error=error) 00154 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00155 deriv => xc_dset_get_derivative(deriv_set,"(rho)",& 00156 allocate_deriv=.TRUE.,error=error) 00157 CALL xc_derivative_get(deriv,deriv_data=e_rho,error=error) 00158 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",& 00159 allocate_deriv=.TRUE.,error=error) 00160 CALL xc_derivative_get(deriv,deriv_data=e_ndrho,error=error) 00161 IF (grad_deriv>1.OR.grad_deriv<-1) THEN 00162 CALL cp_unimplemented_error(fromWhere=routineP, & 00163 message="derivatives bigger than 1 not implemented", & 00164 error=error, error_level=cp_failure_level) 00165 END IF 00166 00167 CALL optx_lda_calc(rho=rho, norm_drho=norm_drho,& 00168 e_0=e_0,e_rho=e_rho,e_ndrho=e_ndrho,& 00169 npoints=npoints,epsilon_rho=epsilon_rho,& 00170 epsilon_drho=epsilon_drho,sx=sx,error=error) 00171 END IF 00172 END SUBROUTINE optx_lda_eval 00173 00174 ! ***************************************************************************** 00189 SUBROUTINE optx_lsd_eval(rho_set,deriv_set,grad_deriv,optx_params,error) 00190 TYPE(xc_rho_set_type), POINTER :: rho_set 00191 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00192 INTEGER, INTENT(in) :: grad_deriv 00193 TYPE(section_vals_type), POINTER :: optx_params 00194 TYPE(cp_error_type), INTENT(inout) :: error 00195 00196 CHARACTER(len=*), PARAMETER :: routineN = 'optx_lsd_eval', 00197 routineP = moduleN//':'//routineN 00198 00199 INTEGER :: ispin, npoints 00200 INTEGER, DIMENSION(:, :), POINTER :: bo 00201 LOGICAL :: failure 00202 REAL(kind=dp) :: epsilon_drho, epsilon_rho, sx 00203 REAL(kind=dp), DIMENSION(:, :, :), 00204 POINTER :: e_0 00205 TYPE(cp_3d_r_p_type), DIMENSION(2) :: e_ndrho, e_rho, ndrho, rho 00206 TYPE(xc_derivative_type), POINTER :: deriv 00207 00208 failure=.FALSE. 00209 NULLIFY(bo,e_0) 00210 DO ispin=1,2 00211 NULLIFY(e_rho(ispin)%array) 00212 NULLIFY(e_ndrho(ispin)%array) 00213 NULLIFY(rho(ispin)%array) 00214 NULLIFY(ndrho(ispin)%array) 00215 ENDDO 00216 00217 CALL section_vals_val_get(optx_params,"scale_x",r_val=sx,error=error) 00218 00219 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00220 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00221 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00222 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00223 IF (.NOT. failure) THEN 00224 CALL xc_rho_set_get(rho_set,rhoa=rho(1)%array,rhob=rho(2)%array,& 00225 norm_drhoa=ndrho(1)%array, & 00226 norm_drhob=ndrho(2)%array,rho_cutoff=epsilon_rho,& 00227 drho_cutoff=epsilon_drho, local_bounds=bo, error=error) 00228 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00229 00230 deriv => xc_dset_get_derivative(deriv_set,"",& 00231 allocate_deriv=.TRUE., error=error) 00232 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00233 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)",& 00234 allocate_deriv=.TRUE.,error=error) 00235 CALL xc_derivative_get(deriv,deriv_data=e_rho(1)%array,error=error) 00236 deriv => xc_dset_get_derivative(deriv_set,"(rhob)",& 00237 allocate_deriv=.TRUE.,error=error) 00238 CALL xc_derivative_get(deriv,deriv_data=e_rho(2)%array,error=error) 00239 00240 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)",& 00241 allocate_deriv=.TRUE.,error=error) 00242 CALL xc_derivative_get(deriv,deriv_data=e_ndrho(1)%array,error=error) 00243 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)",& 00244 allocate_deriv=.TRUE.,error=error) 00245 CALL xc_derivative_get(deriv,deriv_data=e_ndrho(2)%array,error=error) 00246 00247 IF (grad_deriv>1.OR.grad_deriv<-1) THEN 00248 CALL cp_unimplemented_error(fromWhere=routineP, & 00249 message="derivatives bigger than 1 not implemented", & 00250 error=error, error_level=cp_failure_level) 00251 END IF 00252 DO ispin=1,2 00253 CALL optx_lsd_calc(rho=rho(ispin)%array, norm_drho=ndrho(ispin)%array,& 00254 e_0=e_0,e_rho=e_rho(ispin)%array,e_ndrho=e_ndrho(ispin)%array,& 00255 npoints=npoints,epsilon_rho=epsilon_rho,& 00256 epsilon_drho=epsilon_drho,sx=sx,error=error) 00257 ENDDO 00258 END IF 00259 END SUBROUTINE optx_lsd_eval 00260 00261 ! ***************************************************************************** 00276 SUBROUTINE optx_lda_calc(rho,norm_drho,e_0,e_rho,e_ndrho,& 00277 epsilon_rho,epsilon_drho,npoints,sx,error) 00278 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho 00279 REAL(KIND=dp), DIMENSION(*), 00280 INTENT(INOUT) :: e_0, e_rho, e_ndrho 00281 REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_drho 00282 INTEGER, INTENT(in) :: npoints 00283 REAL(kind=dp), INTENT(in) :: sx 00284 TYPE(cp_error_type), INTENT(inout) :: error 00285 00286 REAL(KIND=dp), PARAMETER :: a1cx = 0.9784571170284421_dp, 00287 a2 = 1.43169_dp, gam = 0.006_dp, o43 = 4.0_dp/3.0_dp 00288 00289 INTEGER :: ii 00290 REAL(KIND=dp) :: denom, ex, gamxsxs, myndrho, 00291 myrho, rho43, tmp, xs 00292 00293 !$omp parallel do default (none) & 00294 !$omp shared(rho, norm_drho, e_0, e_rho, e_ndrho) & 00295 !$omp shared(epsilon_rho, epsilon_drho, sx, npoints) & 00296 !$omp private(ii, myrho, myndrho, rho43, xs, gamxsxs) & 00297 !$omp private(denom, ex, tmp) 00298 00299 DO ii=1,npoints 00300 ! we get the full density and need spin parts -> 0.5 00301 myrho = 0.5_dp * rho(ii) 00302 myndrho = 0.5_dp * MAX(norm_drho(ii),epsilon_drho) 00303 IF (myrho> 0.5_dp * epsilon_rho) THEN 00304 rho43 = myrho**o43 00305 xs = (myndrho / rho43) 00306 gamxsxs = gam * xs * xs 00307 denom = 1.0_dp / (1.0_dp + gamxsxs ) 00308 ex = rho43*(a1cx+a2*(gamxsxs*denom)**2) 00309 ! 2.0 for both spins 00310 e_0(ii) = e_0(ii) - ( 2.0_dp * ex ) * sx 00311 tmp = rho43 * 2.0_dp * a2 * gamxsxs * denom**2 * ( 1.0_dp - gamxsxs * denom) 00312 ! derive e_0 wrt to rho (full) and ndrho (also full) 00313 e_rho(ii) = e_rho(ii) - ( ( o43 * ex + tmp * gamxsxs * (-2.0_dp * o43 ) ) / myrho ) * sx 00314 e_ndrho(ii)= e_ndrho(ii) - ( ( tmp * gam * 2.0_dp * myndrho / rho43**2 ) ) * sx 00315 END IF 00316 END DO 00317 00318 !$omp end parallel do 00319 00320 END SUBROUTINE optx_lda_calc 00321 00322 ! ***************************************************************************** 00337 SUBROUTINE optx_lsd_calc(rho,norm_drho,e_0,e_rho,e_ndrho,& 00338 epsilon_rho,epsilon_drho,npoints,sx,error) 00339 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho 00340 REAL(KIND=dp), DIMENSION(*), 00341 INTENT(INOUT) :: e_0, e_rho, e_ndrho 00342 REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_drho 00343 INTEGER, INTENT(in) :: npoints 00344 REAL(kind=dp), INTENT(in) :: sx 00345 TYPE(cp_error_type), INTENT(inout) :: error 00346 00347 REAL(KIND=dp), PARAMETER :: a1cx = 0.9784571170284421_dp, 00348 a2 = 1.43169_dp, gam = 0.006_dp, o43 = 4.0_dp/3.0_dp 00349 00350 INTEGER :: ii 00351 REAL(KIND=dp) :: denom, ex, gamxsxs, myndrho, 00352 myrho, rho43, tmp, xs 00353 00354 !$omp parallel do default(none) & 00355 !$omp shared(rho, norm_drho, e_0, e_rho, e_ndrho) & 00356 !$omp shared(epsilon_rho, epsilon_drho, npoints, sx) & 00357 !$omp private(ii, denom, ex, gamxsxs, myndrho, myrho) & 00358 !$omp private(rho43, tmp, xs) 00359 00360 DO ii=1,npoints 00361 ! we do have the spin density already 00362 myrho = rho(ii) 00363 myndrho = MAX(norm_drho(ii),epsilon_drho) 00364 IF (myrho> epsilon_rho) THEN 00365 rho43 = myrho**o43 00366 xs = (myndrho / rho43) 00367 gamxsxs = gam * xs * xs 00368 denom = 1.0_dp / (1.0_dp + gamxsxs ) 00369 ex = rho43*(a1cx+a2*(gamxsxs*denom)**2) 00370 ! for a single spin 00371 e_0(ii) = e_0(ii) - ex * sx 00372 tmp = rho43 * 2.0_dp * a2 * gamxsxs * denom**2 * ( 1.0_dp - gamxsxs * denom) 00373 ! derive e_0 wrt to rho and ndrho 00374 e_rho(ii) = e_rho(ii) - ( ( o43 * ex + tmp * gamxsxs * (-2.0_dp * o43 ) ) / myrho ) * sx 00375 e_ndrho(ii)= e_ndrho(ii) - ( ( tmp * gam * 2.0_dp * myndrho / rho43**2 ) ) * sx 00376 END IF 00377 END DO 00378 00379 !$omp end parallel do 00380 00381 END SUBROUTINE optx_lsd_calc 00382 00383 END MODULE xc_optx
1.7.3