|
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 ! ***************************************************************************** 00026 00027 MODULE xc_xbr_pbe_lda_hole_t_c_lr 00028 USE f77_blas 00029 USE input_section_types, ONLY: section_vals_type,& 00030 section_vals_val_get 00031 USE kinds, ONLY: dp 00032 USE mathconstants, ONLY: pi 00033 USE timings, ONLY: timeset,& 00034 timestop 00035 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00036 xc_dset_get_derivative 00037 USE xc_derivative_types, ONLY: xc_derivative_get,& 00038 xc_derivative_type 00039 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00040 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 00041 xc_rho_set_type 00042 USE xc_xbecke_roussel, ONLY: x_br_lsd_y_gt_0,& 00043 x_br_lsd_y_gt_0_cutoff,& 00044 x_br_lsd_y_lte_0,& 00045 x_br_lsd_y_lte_0_cutoff 00046 USE xc_xlda_hole_t_c_lr, ONLY: xlda_hole_t_c_lr_lda_calc_0 00047 USE xc_xpbe_hole_t_c_lr, ONLY: xpbe_hole_t_c_lr_lda_calc_1,& 00048 xpbe_hole_t_c_lr_lda_calc_2 00049 #include "cp_common_uses.h" 00050 00051 IMPLICIT NONE 00052 PRIVATE 00053 00054 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00055 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbr_pbe_lda_hole_t_c_lr' 00056 00057 REAL(dp), PARAMETER, PRIVATE :: br_a1 = 1.5255251812009530_dp, 00058 br_a2 = 0.4576575543602858_dp, 00059 br_a3 = 0.4292036732051034_dp, 00060 br_c0 = 0.7566445420735584_dp, 00061 br_c1 = -2.6363977871370960_dp, 00062 br_c2 = 5.4745159964232880_dp, 00063 br_c3 = -12.657308127108290_dp, 00064 br_c4 = 4.1250584725121360_dp, 00065 br_c5 = -30.425133957163840_dp, 00066 br_b0 = 0.4771976183772063_dp, 00067 br_b1 = -1.7799813494556270_dp, 00068 br_b2 = 3.8433841862302150_dp, 00069 br_b3 = -9.5912050880518490_dp, 00070 br_b4 = 2.1730180285916720_dp, 00071 br_b5 = -30.425133851603660_dp, 00072 br_d0 = 0.00004435009886795587_dp, 00073 br_d1 = 0.58128653604457910_dp, 00074 br_d2 = 66.742764515940610_dp, 00075 br_d3 = 434.26780897229770_dp, 00076 br_d4 = 824.7765766052239000_dp, 00077 br_d5 = 1657.9652731582120_dp, 00078 br_e0 = 0.00003347285060926091_dp, 00079 br_e1 = 0.47917931023971350_dp, 00080 br_e2 = 62.392268338574240_dp, 00081 br_e3 = 463.14816427938120_dp, 00082 br_e4 = 785.2360350104029000_dp, 00083 br_e5 = 1657.962968223273000000_dp, 00084 br_BB = 2.085749716493756_dp 00085 00086 REAL(dp), PARAMETER, PRIVATE :: smax = 8.572844_dp, 00087 scutoff = 8.3_dp, 00088 sconst = 18.79622316_dp, 00089 gcutoff = 0.08_dp 00090 00091 REAL(dp), PARAMETER, PRIVATE :: alpha = 0.3956891_dp, 00092 N = -0.0009800242_dp, 00093 mu = 0.00118684_dp 00094 00095 PUBLIC :: xbr_pbe_lda_hole_tc_lr_lda_info, & 00096 xbr_pbe_lda_hole_tc_lr_lsd_info, & 00097 xbr_pbe_lda_hole_tc_lr_lda_eval, & 00098 xbr_pbe_lda_hole_tc_lr_lsd_eval 00099 CONTAINS 00100 00101 ! ***************************************************************************** 00111 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_info(reference,shortform, needs, max_deriv,& 00112 error) 00113 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00114 TYPE(xc_rho_cflags_type), 00115 INTENT(inout), OPTIONAL :: needs 00116 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00117 TYPE(cp_error_type), INTENT(inout) :: error 00118 00119 CHARACTER(len=*), PARAMETER :: 00120 routineN = 'xbr_pbe_lda_hole_tc_lr_lda_info', 00121 routineP = moduleN//':'//routineN 00122 00123 LOGICAL :: failure 00124 00125 failure=.FALSE. 00126 00127 IF ( PRESENT ( reference ) ) THEN 00128 reference = "{LDA version}" 00129 END IF 00130 IF ( PRESENT ( shortform ) ) THEN 00131 shortform = "{LDA}" 00132 END IF 00133 00134 IF (PRESENT(needs)) THEN 00135 needs%rho=.TRUE. 00136 needs%norm_drho=.TRUE. 00137 needs%tau = .TRUE. 00138 needs%laplace_rho = .TRUE. 00139 END IF 00140 00141 IF (PRESENT(max_deriv)) max_deriv=1 00142 00143 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_info 00144 00145 ! ***************************************************************************** 00155 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_info(reference,shortform, needs, max_deriv,& 00156 error) 00157 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00158 TYPE(xc_rho_cflags_type), 00159 INTENT(inout), OPTIONAL :: needs 00160 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00161 TYPE(cp_error_type), INTENT(inout) :: error 00162 00163 CHARACTER(len=*), PARAMETER :: 00164 routineN = 'xbr_pbe_lda_hole_tc_lr_lsd_info', 00165 routineP = moduleN//':'//routineN 00166 00167 LOGICAL :: failure 00168 00169 failure=.FALSE. 00170 00171 IF ( PRESENT ( reference ) ) THEN 00172 reference = "{LDA version}" 00173 END IF 00174 IF ( PRESENT ( shortform ) ) THEN 00175 shortform = "{LDA}" 00176 END IF 00177 00178 IF (PRESENT(needs)) THEN 00179 needs%rho_spin =.TRUE. 00180 needs%norm_drho_spin =.TRUE. 00181 needs%tau_spin = .TRUE. 00182 needs%laplace_rho_spin = .TRUE. 00183 END IF 00184 IF (PRESENT(max_deriv)) max_deriv=1 00185 00186 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_info 00187 00188 ! ***************************************************************************** 00201 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_eval(rho_set,deriv_set,grad_deriv,params,error) 00202 TYPE(xc_rho_set_type), POINTER :: rho_set 00203 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00204 INTEGER, INTENT(in) :: grad_deriv 00205 TYPE(section_vals_type), POINTER :: params 00206 TYPE(cp_error_type), INTENT(inout) :: error 00207 00208 CHARACTER(len=*), PARAMETER :: 00209 routineN = 'xbr_pbe_lda_hole_tc_lr_lda_eval', 00210 routineP = moduleN//':'//routineN 00211 00212 INTEGER :: handle, npoints, stat 00213 INTEGER, DIMENSION(:, :), POINTER :: bo 00214 LOGICAL :: failure 00215 REAL(dp) :: gamma, R, sx 00216 REAL(kind=dp) :: epsilon_norm_drho, epsilon_rho 00217 REAL(kind=dp), DIMENSION(:, :, :), 00218 POINTER :: dummy, e_0, e_laplace_rho, 00219 e_ndrho, e_rho, e_tau, 00220 laplace_rho, norm_drho, rho, 00221 tau 00222 TYPE(xc_derivative_type), POINTER :: deriv 00223 00224 CALL timeset(routineN,handle) 00225 00226 failure=.FALSE. 00227 NULLIFY(bo) 00228 00229 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00230 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00231 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00232 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00233 IF (.NOT. failure) THEN 00234 CALL xc_rho_set_get(rho_set, rho=rho, norm_drho=norm_drho, & 00235 tau=tau,laplace_rho=laplace_rho,local_bounds=bo,& 00236 rho_cutoff=epsilon_rho,drho_cutoff=epsilon_norm_drho,error=error) 00237 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00238 00239 ! meaningful default for the arrays we don't need: let us make compiler 00240 ! and debugger happy... 00241 IF (cp_debug) THEN 00242 ALLOCATE(dummy(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 00243 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00244 ELSE 00245 dummy=> rho 00246 END IF 00247 00248 e_0 => dummy 00249 e_rho => dummy 00250 e_ndrho => dummy 00251 e_tau => dummy 00252 e_laplace_rho => dummy 00253 00254 IF (grad_deriv>=0) THEN 00255 deriv => xc_dset_get_derivative(deriv_set,"",& 00256 allocate_deriv=.TRUE., error=error) 00257 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00258 END IF 00259 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 00260 deriv => xc_dset_get_derivative(deriv_set,"(rho)",& 00261 allocate_deriv=.TRUE.,error=error) 00262 CALL xc_derivative_get(deriv,deriv_data=e_rho,error=error) 00263 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",& 00264 allocate_deriv=.TRUE.,error=error) 00265 CALL xc_derivative_get(deriv,deriv_data=e_ndrho,error=error) 00266 deriv => xc_dset_get_derivative(deriv_set,"(tau)",& 00267 allocate_deriv=.TRUE.,error=error) 00268 CALL xc_derivative_get(deriv,deriv_data=e_tau,error=error) 00269 deriv => xc_dset_get_derivative(deriv_set,"(laplace_rho)",& 00270 allocate_deriv=.TRUE.,error=error) 00271 CALL xc_derivative_get(deriv,deriv_data=e_laplace_rho,error=error) 00272 END IF 00273 IF (grad_deriv>1.OR.grad_deriv<-1) THEN 00274 CALL cp_unimplemented_error(fromWhere=routineP, & 00275 message="derivatives bigger than 1 not implemented", & 00276 error=error, error_level=cp_failure_level) 00277 END IF 00278 00279 CALL section_vals_val_get(params,"scale_x",r_val=sx,error=error) 00280 CALL section_vals_val_get(params,"CUTOFF_RADIUS",r_val=R,error=error) 00281 CALL section_vals_val_get(params,"GAMMA",r_val=gamma,error=error) 00282 00283 IF ( R == 0.0_dp ) THEN 00284 CALL cp_unimplemented_error(fromWhere=routineP, & 00285 message="Cutoff_Radius 0.0 not implemented", & 00286 error=error, error_level=cp_failure_level) 00287 END IF 00288 00289 !$omp parallel default(none) & 00290 !$omp shared(rho, norm_drho, laplace_rho, tau, e_0, e_rho) & 00291 !$omp shared(e_ndrho, e_tau, e_laplace_rho, grad_deriv) & 00292 !$omp shared(npoints, epsilon_rho) & 00293 !$omp shared(epsilon_norm_drho, sx, r, gamma, error) 00294 00295 CALL xbr_pbe_lda_hole_tc_lr_lda_calc(rho=rho, norm_drho=norm_drho,& 00296 laplace_rho=laplace_rho,tau=tau,e_0=e_0,e_rho=e_rho,e_ndrho=e_ndrho,& 00297 e_tau=e_tau,e_laplace_rho=e_laplace_rho,grad_deriv=grad_deriv, & 00298 npoints=npoints,epsilon_rho=epsilon_rho,& 00299 epsilon_norm_drho=epsilon_norm_drho,sx=sx,R=R,gamma=gamma,error=error) 00300 00301 !$omp end parallel 00302 00303 IF (cp_debug) THEN 00304 DEALLOCATE(dummy,stat=stat) 00305 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00306 ELSE 00307 NULLIFY(dummy) 00308 END IF 00309 END IF 00310 CALL timestop(handle) 00311 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_eval 00312 00313 ! ***************************************************************************** 00328 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_calc(rho,norm_drho,laplace_rho,tau,e_0,e_rho,& 00329 e_ndrho,e_tau,e_laplace_rho,grad_deriv,npoints,& 00330 epsilon_rho,epsilon_norm_drho,sx,R,gamma,error) 00331 00332 INTEGER, INTENT(in) :: npoints, grad_deriv 00333 REAL(kind=dp), DIMENSION(1:npoints), 00334 INTENT(inout) :: e_laplace_rho, e_tau, 00335 e_ndrho, e_rho, e_0 00336 REAL(kind=dp), DIMENSION(1:npoints), 00337 INTENT(in) :: tau, laplace_rho, norm_drho, 00338 rho 00339 REAL(kind=dp), INTENT(in) :: epsilon_rho, 00340 epsilon_norm_drho, sx, R, 00341 gamma 00342 TYPE(cp_error_type), INTENT(inout) :: error 00343 00344 CHARACTER(len=*), PARAMETER :: 00345 routineN = 'xbr_pbe_lda_hole_tc_lr_lda_calc', 00346 routineP = moduleN//':'//routineN 00347 00348 INTEGER :: ip 00349 REAL(dp) :: dFermi_dlaplace_rho, dFermi_dndrho, dFermi_drho, dFermi_dtau, 00350 e_0_br, e_0_lda, e_0_pbe, e_laplace_rho_br, e_ndrho_br, e_ndrho_pbe, 00351 e_rho_br, e_rho_lda, e_rho_pbe, e_tau_br, Fermi, Fx, my_laplace_rho, 00352 my_ndrho, my_rho, my_tau, ss, ss2, sscale, t1, t15, t16, t2, t3, t4, 00353 t5, t6, t7, t8, t9, yval 00354 00355 !$omp do 00356 00357 DO ip = 1,npoints 00358 my_rho = 0.5_dp*MAX(rho(ip),0.0_dp) 00359 IF(my_rho > epsilon_rho) THEN 00360 my_ndrho = 0.5_dp*MAX(norm_drho(ip),EPSILON(0.0_dp)*1.e4_dp) 00361 my_tau = 0.5_dp*MAX(EPSILON(0.0_dp)*1.e4_dp,tau(ip)) 00362 my_laplace_rho = 0.5_dp*laplace_rho(ip) 00363 00364 ! ** We calculate first the Becke-Roussel part, saving everything in local variables 00365 t1 = pi ** (0.1e1_dp / 0.3e1_dp) 00366 t2 = t1 ** 2 00367 t3 = my_rho ** (0.1e1_dp / 0.3e1_dp) 00368 t4 = t3 ** 2 00369 t5 = t4 * my_rho 00370 t8 = my_ndrho ** 2 00371 t9 = 0.1e1_dp / my_rho 00372 t15 = my_laplace_rho / 0.6e1_dp - gamma * (2.0_dp*my_tau - t8 * t9 / 0.4e1_dp) / 0.3e1_dp 00373 t16 = 0.1e1_dp / t15 00374 yval = 0.2e1_dp / 0.3e1_dp * t2 * t5 * t16 00375 00376 e_0_br = 0.0_dp 00377 e_rho_br = 0.0_dp 00378 e_ndrho_br = 0.0_dp 00379 e_tau_br = 0.0_dp 00380 e_laplace_rho_br = 0.0_dp 00381 00382 IF( R == 0.0_dp ) THEN 00383 IF( yval <= 0.0_dp) THEN 00384 CALL x_br_lsd_y_lte_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, & 00385 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br,& 00386 sx, R, gamma, grad_deriv, error) 00387 ! VERY UGLY HACK e_0 has to multiplied by the factor 2 00388 e_0_br = 2.0_dp * e_0_br 00389 ELSE 00390 CALL x_br_lsd_y_gt_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, & 00391 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br,& 00392 sx, R, gamma, grad_deriv, error) 00393 ! VERY UGLY HACK e_0 has to multiplied by the factor 2 00394 e_0_br = 2.0_dp * e_0_br 00395 END IF 00396 ELSE 00397 IF( yval <= 0.0_dp) THEN 00398 CALL x_br_lsd_y_lte_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, & 00399 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br,& 00400 sx, R, gamma, grad_deriv, error) 00401 ! VERY UGLY HACK e_0 has to multiplied by the factor 2 00402 e_0_br = 2.0_dp * e_0_br 00403 ELSE 00404 CALL x_br_lsd_y_gt_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, & 00405 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br,& 00406 sx, R, gamma, grad_deriv, error) 00407 ! VERY UGLY HACK e_0 has to multiplied by the factor 2 00408 e_0_br = 2.0_dp * e_0_br 00409 END IF 00410 END IF 00411 00412 00413 ! ** Now we calculate the pbe cutoff part 00414 ! ** Attention we need to scale rho, ndrho first 00415 my_rho = my_rho * 2.0_dp 00416 my_ndrho = my_ndrho * 2.0_dp 00417 00418 ! ** Do some precalculation in order to catch the correct branch afterwards 00419 sscale = 1.0_dp 00420 t1 = pi ** 2 00421 t2 = t1 * my_rho 00422 t3 = t2 ** (0.1e1_dp / 0.3e1_dp) 00423 t4 = 0.1e1_dp / t3 00424 t6 = my_ndrho * t4 00425 t7 = 0.1e1_dp / my_rho 00426 t8 = t7 * sscale 00427 ss = 0.3466806371753173524216762e0_dp * t6 * t8 00428 IF( ss > scutoff) THEN 00429 ss2 = ss*ss 00430 sscale = (smax*ss2-sconst)/(ss2*ss) 00431 END IF 00432 e_0_pbe = 0.0_dp 00433 e_rho_pbe = 0.0_dp 00434 e_ndrho_pbe = 0.0_dp 00435 IF(ss*sscale>gcutoff) THEN 00436 CALL xpbe_hole_t_c_lr_lda_calc_1(e_0_pbe, e_rho_pbe, e_ndrho_pbe,& 00437 my_rho,& 00438 my_ndrho, sscale, sx, R, grad_deriv) 00439 ELSE 00440 CALL xpbe_hole_t_c_lr_lda_calc_2(e_0_pbe, e_rho_pbe, e_ndrho_pbe,& 00441 my_rho,& 00442 my_ndrho, sscale, sx, R, grad_deriv) 00443 END IF 00444 00445 ! ** Finally we get the LDA part 00446 00447 e_0_lda = 0.0_dp 00448 e_rho_lda = 0.0_dp 00449 CALL xlda_hole_t_c_lr_lda_calc_0(grad_deriv, my_rho, e_0_lda, e_rho_lda,& 00450 sx, R, error) 00451 00452 00453 Fx = e_0_br/e_0_lda 00454 00455 Fermi = alpha/(EXP( (Fx-mu)/N ) + 1.0_dp) 00456 00457 dFermi_drho = -Fermi**2/alpha/N*(e_rho_br/e_0_lda-e_0_br*e_rho_lda/e_0_lda**2)*EXP((Fx-mu)/N) 00458 dFermi_dndrho = -Fermi**2/alpha/N*(e_ndrho_br/e_0_lda)*EXP((Fx-mu)/N) 00459 dFermi_dtau = -Fermi**2/alpha/N*(e_tau_br/e_0_lda)*EXP((Fx-mu)/N) 00460 dFermi_dlaplace_rho = -Fermi**2/alpha/N*(e_laplace_rho_br/e_0_lda)*EXP((Fx-mu)/N) 00461 00462 00463 e_0(ip) = e_0(ip) + ( Fermi * e_0_pbe + (1.0_dp-Fermi) * e_0_br ) * sx 00464 00465 IF( grad_deriv >=1 .OR. grad_deriv == -1) THEN 00466 00467 e_rho(ip) = e_rho(ip) + ( Fermi * e_rho_pbe + dFermi_drho * e_0_pbe + & 00468 (1.0_dp - Fermi) * e_rho_br - dFermi_drho*e_0_br ) * sx 00469 00470 e_ndrho(ip) = e_ndrho(ip) + ( Fermi * e_ndrho_pbe + dFermi_dndrho * e_0_pbe + & 00471 (1.0_dp - Fermi) * e_ndrho_br - dFermi_dndrho*e_0_br ) * sx 00472 00473 e_tau(ip) = e_tau(ip) + ( dFermi_dtau * e_0_pbe + & 00474 (1.0_dp - Fermi) * e_tau_br - dFermi_dtau*e_0_br ) * sx 00475 00476 e_laplace_rho(ip) = e_laplace_rho(ip) + ( dFermi_dlaplace_rho * e_0_pbe + & 00477 (1.0_dp - Fermi) * e_laplace_rho_br - dFermi_dlaplace_rho*e_0_br ) * sx 00478 END IF 00479 00480 00481 END IF 00482 END DO 00483 00484 !$omp end do 00485 00486 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_calc 00487 00488 ! ***************************************************************************** 00501 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_eval(rho_set,deriv_set,grad_deriv,params,error) 00502 TYPE(xc_rho_set_type), POINTER :: rho_set 00503 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00504 INTEGER, INTENT(in) :: grad_deriv 00505 TYPE(section_vals_type), POINTER :: params 00506 TYPE(cp_error_type), INTENT(inout) :: error 00507 00508 CHARACTER(len=*), PARAMETER :: 00509 routineN = 'xbr_pbe_lda_hole_tc_lr_lsd_eval', 00510 routineP = moduleN//':'//routineN 00511 00512 INTEGER :: handle, npoints, stat 00513 INTEGER, DIMENSION(:, :), POINTER :: bo 00514 LOGICAL :: failure 00515 REAL(dp) :: gamma, R, sx 00516 REAL(kind=dp) :: epsilon_norm_drho, epsilon_rho 00517 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, 00518 e_laplace_rhob, e_ndrhoa, e_ndrhob, e_rhoa, e_rhob, e_tau_a, e_tau_b, 00519 laplace_rhoa, laplace_rhob, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, 00520 tau_b 00521 TYPE(xc_derivative_type), POINTER :: deriv 00522 00523 CALL timeset(routineN,handle) 00524 00525 failure=.FALSE. 00526 NULLIFY(bo) 00527 00528 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00529 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00530 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00531 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00532 IF (.NOT. failure) THEN 00533 CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, & 00534 norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, laplace_rhoa=laplace_rhoa,& 00535 laplace_rhob=laplace_rhob,local_bounds=bo,& 00536 rho_cutoff=epsilon_rho,drho_cutoff=epsilon_norm_drho,error=error) 00537 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00538 00539 ! meaningful default for the arrays we don't need: let us make compiler 00540 ! and debugger happy... 00541 IF (cp_debug) THEN 00542 ALLOCATE(dummy(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 00543 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00544 ELSE 00545 dummy=> rhoa 00546 END IF 00547 00548 e_0 => dummy 00549 e_rhoa => dummy 00550 e_rhob => dummy 00551 e_ndrhoa => dummy 00552 e_ndrhob => dummy 00553 e_tau_a => dummy 00554 e_tau_b => dummy 00555 e_laplace_rhoa => dummy 00556 e_laplace_rhob => dummy 00557 00558 IF (grad_deriv>=0) THEN 00559 deriv => xc_dset_get_derivative(deriv_set,"",& 00560 allocate_deriv=.TRUE., error=error) 00561 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00562 END IF 00563 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 00564 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)",& 00565 allocate_deriv=.TRUE.,error=error) 00566 CALL xc_derivative_get(deriv,deriv_data=e_rhoa,error=error) 00567 deriv => xc_dset_get_derivative(deriv_set,"(rhob)",& 00568 allocate_deriv=.TRUE.,error=error) 00569 CALL xc_derivative_get(deriv,deriv_data=e_rhob,error=error) 00570 00571 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)",& 00572 allocate_deriv=.TRUE.,error=error) 00573 CALL xc_derivative_get(deriv,deriv_data=e_ndrhoa,error=error) 00574 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)",& 00575 allocate_deriv=.TRUE.,error=error) 00576 CALL xc_derivative_get(deriv,deriv_data=e_ndrhob,error=error) 00577 00578 deriv => xc_dset_get_derivative(deriv_set,"(tau_a)",& 00579 allocate_deriv=.TRUE.,error=error) 00580 CALL xc_derivative_get(deriv,deriv_data=e_tau_a,error=error) 00581 deriv => xc_dset_get_derivative(deriv_set,"(tau_b)",& 00582 allocate_deriv=.TRUE.,error=error) 00583 CALL xc_derivative_get(deriv,deriv_data=e_tau_b,error=error) 00584 00585 deriv => xc_dset_get_derivative(deriv_set,"(laplace_rhoa)",& 00586 allocate_deriv=.TRUE.,error=error) 00587 CALL xc_derivative_get(deriv,deriv_data=e_laplace_rhoa,error=error) 00588 deriv => xc_dset_get_derivative(deriv_set,"(laplace_rhob)",& 00589 allocate_deriv=.TRUE.,error=error) 00590 CALL xc_derivative_get(deriv,deriv_data=e_laplace_rhob,error=error) 00591 END IF 00592 IF (grad_deriv>1.OR.grad_deriv<-1) THEN 00593 CALL cp_unimplemented_error(fromWhere=routineP, & 00594 message="derivatives bigger than 1 not implemented", & 00595 error=error, error_level=cp_failure_level) 00596 END IF 00597 00598 CALL section_vals_val_get(params,"scale_x",r_val=sx,error=error) 00599 CALL section_vals_val_get(params,"CUTOFF_RADIUS",r_val=R,error=error) 00600 CALL section_vals_val_get(params,"GAMMA",r_val=gamma,error=error) 00601 00602 IF ( R == 0.0_dp ) THEN 00603 CALL cp_unimplemented_error(fromWhere=routineP, & 00604 message="Cutoff_Radius 0.0 not implemented", & 00605 error=error, error_level=cp_failure_level) 00606 END IF 00607 00608 !$omp parallel default(none) & 00609 !$omp shared(rhoa, norm_drhoa, laplace_rhoa, tau_a, e_0) & 00610 !$omp shared(e_rhoa, e_ndrhoa, e_tau_a, e_laplace_rhoa) & 00611 !$omp shared(grad_deriv, npoints, epsilon_rho) & 00612 !$omp shared(epsilon_norm_drho, sx, r, gamma, error) & 00613 !$omp shared(rhob, norm_drhob, laplace_rhob, tau_b, e_rhob) & 00614 !$omp shared(e_ndrhob, e_tau_b, e_laplace_rhob) 00615 00616 CALL xbr_pbe_lda_hole_tc_lr_lsd_calc(rho=rhoa, norm_drho=norm_drhoa,& 00617 laplace_rho=laplace_rhoa,tau=tau_a,e_0=e_0,e_rho=e_rhoa,e_ndrho=e_ndrhoa,& 00618 e_tau=e_tau_a,e_laplace_rho=e_laplace_rhoa,grad_deriv=grad_deriv, & 00619 npoints=npoints,epsilon_rho=epsilon_rho,& 00620 epsilon_norm_drho=epsilon_norm_drho,sx=sx,R=R,gamma=gamma,error=error) 00621 00622 CALL xbr_pbe_lda_hole_tc_lr_lsd_calc(rho=rhob, norm_drho=norm_drhob,& 00623 laplace_rho=laplace_rhob,tau=tau_b,e_0=e_0,e_rho=e_rhob,e_ndrho=e_ndrhob,& 00624 e_tau=e_tau_b,e_laplace_rho=e_laplace_rhob,grad_deriv=grad_deriv, & 00625 npoints=npoints,epsilon_rho=epsilon_rho,& 00626 epsilon_norm_drho=epsilon_norm_drho,sx=sx,R=R,gamma=gamma,error=error) 00627 00628 !$omp end parallel 00629 00630 IF (cp_debug) THEN 00631 DEALLOCATE(dummy,stat=stat) 00632 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00633 ELSE 00634 NULLIFY(dummy) 00635 END IF 00636 END IF 00637 CALL timestop(handle) 00638 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_eval 00639 ! ***************************************************************************** 00654 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_calc(rho,norm_drho,laplace_rho,tau,e_0,e_rho,& 00655 e_ndrho,e_tau,e_laplace_rho,grad_deriv,npoints,& 00656 epsilon_rho,epsilon_norm_drho,sx,R,gamma,error) 00657 00658 INTEGER, INTENT(in) :: npoints, grad_deriv 00659 REAL(kind=dp), DIMENSION(1:npoints), 00660 INTENT(inout) :: e_laplace_rho, e_tau, 00661 e_ndrho, e_rho, e_0 00662 REAL(kind=dp), DIMENSION(1:npoints), 00663 INTENT(in) :: tau, laplace_rho, norm_drho, 00664 rho 00665 REAL(kind=dp), INTENT(in) :: epsilon_rho, 00666 epsilon_norm_drho, sx, R, 00667 gamma 00668 TYPE(cp_error_type), INTENT(inout) :: error 00669 00670 CHARACTER(len=*), PARAMETER :: 00671 routineN = 'xbr_pbe_lda_hole_tc_lr_lsd_calc', 00672 routineP = moduleN//':'//routineN 00673 00674 INTEGER :: ip 00675 REAL(dp) :: dFermi_dlaplace_rho, dFermi_dndrho, dFermi_drho, dFermi_dtau, 00676 e_0_br, e_0_lda, e_0_pbe, e_laplace_rho_br, e_ndrho_br, e_ndrho_pbe, 00677 e_rho_br, e_rho_lda, e_rho_pbe, e_tau_br, Fermi, Fx, my_laplace_rho, 00678 my_ndrho, my_rho, my_tau, ss, ss2, sscale, t1, t15, t16, t2, t3, t4, 00679 t5, t6, t7, t8, t9, yval 00680 00681 !$omp do 00682 00683 DO ip = 1,npoints 00684 my_rho = MAX(rho(ip),0.0_dp) 00685 IF(my_rho > epsilon_rho) THEN 00686 my_ndrho = MAX(norm_drho(ip),EPSILON(0.0_dp)*1.e4_dp) 00687 my_tau=1.0_dp*MAX(EPSILON(0.0_dp)*1.e4_dp,tau(ip)) 00688 my_laplace_rho = 1.0_dp*laplace_rho(ip) 00689 00690 t1 = pi ** (0.1e1_dp / 0.3e1_dp) 00691 t2 = t1 ** 2 00692 t3 = my_rho ** (0.1e1_dp / 0.3e1_dp) 00693 t4 = t3 ** 2 00694 t5 = t4 * my_rho 00695 t8 = my_ndrho ** 2 00696 t9 = 0.1e1_dp / my_rho 00697 t15 = my_laplace_rho / 0.6e1_dp - gamma * (2.0_dp*my_tau - t8 * t9 / 0.4e1_dp) / 0.3e1_dp 00698 t16 = 0.1e1_dp / t15 00699 yval = 0.2e1_dp / 0.3e1_dp * t2 * t5 * t16 00700 00701 e_0_br = 0.0_dp 00702 e_rho_br = 0.0_dp 00703 e_ndrho_br = 0.0_dp 00704 e_tau_br = 0.0_dp 00705 e_laplace_rho_br = 0.0_dp 00706 00707 IF( R == 0.0_dp ) THEN 00708 IF( yval <= 0.0_dp) THEN 00709 CALL x_br_lsd_y_lte_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, & 00710 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br,& 00711 sx, R, gamma, grad_deriv, error) 00712 ELSE 00713 CALL x_br_lsd_y_gt_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, & 00714 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br,& 00715 sx, R, gamma, grad_deriv, error) 00716 END IF 00717 ELSE 00718 IF( yval <= 0.0_dp) THEN 00719 CALL x_br_lsd_y_lte_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, & 00720 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br,& 00721 sx, R, gamma, grad_deriv, error) 00722 ELSE 00723 CALL x_br_lsd_y_gt_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, & 00724 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br,& 00725 sx, R, gamma, grad_deriv, error) 00726 END IF 00727 END IF 00728 00729 ! ** Now we calculate the pbe cutoff part 00730 ! ** Attention we need to scale rho, ndrho first 00731 my_rho = my_rho * 2.0_dp 00732 my_ndrho = my_ndrho * 2.0_dp 00733 00734 ! ** Do some precalculation in order to catch the correct branch afterwards 00735 sscale = 1.0_dp 00736 t1 = pi ** 2 00737 t2 = t1 * my_rho 00738 t3 = t2 ** (0.1e1_dp / 0.3e1_dp) 00739 t4 = 0.1e1_dp / t3 00740 t6 = my_ndrho * t4 00741 t7 = 0.1e1_dp / my_rho 00742 t8 = t7 * sscale 00743 ss = 0.3466806371753173524216762e0_dp * t6 * t8 00744 IF( ss > scutoff) THEN 00745 ss2 = ss*ss 00746 sscale = (smax*ss2-sconst)/(ss2*ss) 00747 END IF 00748 e_0_pbe = 0.0_dp 00749 e_rho_pbe = 0.0_dp 00750 e_ndrho_pbe = 0.0_dp 00751 IF(ss*sscale>gcutoff) THEN 00752 CALL xpbe_hole_t_c_lr_lda_calc_1(e_0_pbe, e_rho_pbe, e_ndrho_pbe,& 00753 my_rho,& 00754 my_ndrho, sscale, sx, R, grad_deriv) 00755 ELSE 00756 CALL xpbe_hole_t_c_lr_lda_calc_2(e_0_pbe, e_rho_pbe, e_ndrho_pbe,& 00757 my_rho,& 00758 my_ndrho, sscale, sx, R, grad_deriv) 00759 END IF 00760 00761 e_0_pbe = 0.5_dp * e_0_pbe 00762 00763 00764 ! ** Finally we get the LDA part 00765 00766 e_0_lda = 0.0_dp 00767 e_rho_lda = 0.0_dp 00768 CALL xlda_hole_t_c_lr_lda_calc_0(grad_deriv, my_rho, e_0_lda, e_rho_lda,& 00769 sx, R, error) 00770 e_0_lda = 0.5_dp * e_0_lda 00771 00772 Fx = e_0_br/e_0_lda 00773 00774 Fermi = alpha/(EXP( (Fx-mu)/N ) + 1.0_dp) 00775 00776 dFermi_drho = -Fermi**2/alpha/N*(e_rho_br/e_0_lda-e_0_br*e_rho_lda/e_0_lda**2)*EXP((Fx-mu)/N) 00777 dFermi_dndrho = -Fermi**2/alpha/N*(e_ndrho_br/e_0_lda)*EXP((Fx-mu)/N) 00778 dFermi_dtau = -Fermi**2/alpha/N*(e_tau_br/e_0_lda)*EXP((Fx-mu)/N) 00779 dFermi_dlaplace_rho = -Fermi**2/alpha/N*(e_laplace_rho_br/e_0_lda)*EXP((Fx-mu)/N) 00780 00781 00782 e_0(ip) = e_0(ip) + ( Fermi * e_0_pbe + (1.0_dp-Fermi) * e_0_br ) * sx 00783 00784 IF( grad_deriv >=1 .OR. grad_deriv == -1) THEN 00785 00786 e_rho(ip) = e_rho(ip) + ( Fermi * e_rho_pbe + dFermi_drho * e_0_pbe + & 00787 (1.0_dp - Fermi) * e_rho_br - dFermi_drho*e_0_br ) * sx 00788 00789 e_ndrho(ip) = e_ndrho(ip) + ( Fermi * e_ndrho_pbe + dFermi_dndrho * e_0_pbe + & 00790 (1.0_dp - Fermi) * e_ndrho_br - dFermi_dndrho*e_0_br ) * sx 00791 00792 e_tau(ip) = e_tau(ip) + ( dFermi_dtau * e_0_pbe + & 00793 (1.0_dp - Fermi) * e_tau_br - dFermi_dtau*e_0_br ) * sx 00794 00795 e_laplace_rho(ip) = e_laplace_rho(ip) + ( dFermi_dlaplace_rho * e_0_pbe + & 00796 (1.0_dp - Fermi) * e_laplace_rho_br - dFermi_dlaplace_rho*e_0_br ) * sx 00797 END IF 00798 00799 END IF 00800 END DO 00801 00802 !$omp end do 00803 00804 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_calc 00805 00806 END MODULE xc_xbr_pbe_lda_hole_t_c_lr
1.7.3