|
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 ! ***************************************************************************** 00019 MODULE xc_pbe 00020 USE bibliography, ONLY: Perdew1996,& 00021 Perdew2008,& 00022 Zhang1998,& 00023 cite_reference 00024 USE f77_blas 00025 USE input_constants, ONLY: xc_pbe_orig,& 00026 xc_pbe_rev,& 00027 xc_pbe_sol 00028 USE input_section_types, ONLY: section_vals_type,& 00029 section_vals_val_get 00030 USE kinds, ONLY: dp 00031 USE mathconstants, ONLY: pi 00032 USE timings, ONLY: timeset,& 00033 timestop 00034 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00035 xc_dset_get_derivative 00036 USE xc_derivative_types, ONLY: xc_derivative_get,& 00037 xc_derivative_type 00038 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00039 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 00040 xc_rho_set_type 00041 #include "cp_common_uses.h" 00042 00043 IMPLICIT NONE 00044 PRIVATE 00045 00046 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00047 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_pbe' 00048 REAL(kind=dp), PARAMETER, PRIVATE :: a=0.04918_dp, b=0.132_dp, 00049 c=0.2533_dp,d=0.349_dp 00050 00051 PUBLIC :: pbe_lda_info, pbe_lsd_info, pbe_lda_eval, pbe_lsd_eval 00052 CONTAINS 00053 00054 ! ***************************************************************************** 00065 SUBROUTINE pbe_lda_info(pbe_params,reference,shortform, needs, max_deriv,& 00066 error) 00067 TYPE(section_vals_type), POINTER :: pbe_params 00068 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00069 TYPE(xc_rho_cflags_type), 00070 INTENT(inout), OPTIONAL :: needs 00071 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00072 TYPE(cp_error_type), INTENT(inout) :: error 00073 00074 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lda_info', 00075 routineP = moduleN//':'//routineN 00076 00077 INTEGER :: param 00078 LOGICAL :: failure 00079 REAL(kind=dp) :: sc, sx 00080 00081 failure=.FALSE. 00082 CALL section_vals_val_get(pbe_params,"parametrization",i_val=param,error=error) 00083 CALL section_vals_val_get(pbe_params,"scale_x",r_val=sx,error=error) 00084 CALL section_vals_val_get(pbe_params,"scale_c",r_val=sc,error=error) 00085 00086 SELECT CASE(param) 00087 CASE(xc_pbe_orig) 00088 CALL cite_reference(Perdew1996) 00089 IF (sx==1._dp .AND.sc==1._dp) THEN 00090 IF ( PRESENT ( reference ) ) THEN 00091 reference = "J.P.Perdew, K.Burke, M.Ernzerhof, "//& 00092 "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)"//& 00093 "{spin unpolarized}" 00094 END IF 00095 IF ( PRESENT ( shortform ) ) THEN 00096 shortform = "PBE Perdew-Burke-Ernzerhof xc functional (unpolarized)" 00097 END IF 00098 ELSE 00099 IF ( PRESENT ( reference ) ) THEN 00100 WRITE (reference,"(a,a,'sx=',f5.3,'sc=',f5.3,a)")& 00101 "J.P.Perdew, K.Burke, M.Ernzerhof, ",& 00102 "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)",& 00103 sx,sc,"{spin unpolarized}" 00104 END IF 00105 IF ( PRESENT ( shortform ) ) THEN 00106 WRITE (shortform,"(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')")& 00107 "PBE, Perdew-Burke-Ernzerhof xc functional (unpolarized)",sx,sc 00108 END IF 00109 END IF 00110 CASE(xc_pbe_rev) 00111 CALL cite_reference(Perdew1996) 00112 CALL cite_reference(Zhang1998) 00113 IF (sx==1._dp .AND.sc==1._dp) THEN 00114 IF ( PRESENT ( reference ) ) THEN 00115 reference = "revPBE, Yingkay Zhang and Weitao Yang,"//& 00116 " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)"//& 00117 "{spin unpolarized}" 00118 END IF 00119 IF ( PRESENT ( shortform ) ) THEN 00120 shortform = "revPBE, revised PBE xc functional (unpolarized)" 00121 END IF 00122 ELSE 00123 IF ( PRESENT ( reference ) ) THEN 00124 WRITE (reference,"(a,a,'sx=',f5.3,'sc=',f5.3,a)")& 00125 "revPBE, Yingkay Zhang and Weitao Yang,",& 00126 " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)",& 00127 sx,sc,"{spin unpolarized}" 00128 END IF 00129 IF ( PRESENT ( shortform ) ) THEN 00130 WRITE (shortform,"(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')")& 00131 "revPBE, revised PBE xc functional (unpolarized)",sx,sc 00132 END IF 00133 END IF 00134 CASE(xc_pbe_sol) 00135 CALL cite_reference(Perdew1996) 00136 CALL cite_reference(Perdew2008) 00137 IF (sx==1._dp .AND.sc==1._dp) THEN 00138 IF ( PRESENT ( reference ) ) THEN 00139 reference = "PBEsol, J.P. Perdew et al., "//& 00140 "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) "//& 00141 "{spin unpolarized}" 00142 END IF 00143 IF ( PRESENT ( shortform ) ) THEN 00144 shortform = "PBEsol, PBE for solids and surfaces xc functional (unpolarized)" 00145 END IF 00146 ELSE 00147 IF ( PRESENT ( reference ) ) THEN 00148 WRITE (reference,"(a,a,'sx=',f5.3,'sc=',f5.3,a)")& 00149 "PBEsol, J.P. Perdew et al., ",& 00150 "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) ",& 00151 sx,sc,"{spin unpolarized}" 00152 END IF 00153 IF ( PRESENT ( shortform ) ) THEN 00154 WRITE (shortform,"(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')")& 00155 "PBEsol, PBE for solids and surfaces xc functional (unpolarized)",sx,sc 00156 END IF 00157 END IF 00158 CASE default 00159 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00160 END SELECT 00161 IF (PRESENT(needs)) THEN 00162 needs%rho=.TRUE. 00163 needs%rho_1_3=.TRUE. 00164 needs%norm_drho=.TRUE. 00165 END IF 00166 IF (PRESENT(max_deriv)) max_deriv=3 00167 00168 END SUBROUTINE pbe_lda_info 00169 00170 ! ***************************************************************************** 00181 SUBROUTINE pbe_lsd_info(pbe_params,reference,shortform, needs, max_deriv,& 00182 error) 00183 TYPE(section_vals_type), POINTER :: pbe_params 00184 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00185 TYPE(xc_rho_cflags_type), 00186 INTENT(inout), OPTIONAL :: needs 00187 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00188 TYPE(cp_error_type), INTENT(inout) :: error 00189 00190 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lsd_info', 00191 routineP = moduleN//':'//routineN 00192 00193 INTEGER :: param 00194 LOGICAL :: failure 00195 REAL(kind=dp) :: sc, sx 00196 00197 failure=.FALSE. 00198 CALL section_vals_val_get(pbe_params,"parametrization",i_val=param,error=error) 00199 CALL section_vals_val_get(pbe_params,"scale_x",r_val=sx,error=error) 00200 CALL section_vals_val_get(pbe_params,"scale_c",r_val=sc,error=error) 00201 00202 SELECT CASE(param) 00203 CASE(xc_pbe_orig) 00204 CALL cite_reference(Perdew1996) 00205 IF (sx==1._dp .AND.sc==1._dp) THEN 00206 IF ( PRESENT ( reference ) ) THEN 00207 reference = "J.P.Perdew, K.Burke, M.Ernzerhof, "//& 00208 "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)"//& 00209 "{spin polarized}" 00210 END IF 00211 IF ( PRESENT ( shortform ) ) THEN 00212 shortform = "PBE xc functional (spin polarized)" 00213 END IF 00214 ELSE 00215 IF ( PRESENT ( reference ) ) THEN 00216 WRITE (reference,"(a,a,'sx=',f5.3,'sc=',f5.3,a)")& 00217 "J.P.Perdew, K.Burke, M.Ernzerhof, ",& 00218 "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)",& 00219 sx,sc,"{spin polarized}" 00220 END IF 00221 IF ( PRESENT ( shortform ) ) THEN 00222 WRITE (shortform,"(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')")& 00223 "PBE xc functional (spin polarized)",sx,sc 00224 END IF 00225 END IF 00226 CASE(xc_pbe_rev) 00227 CALL cite_reference(Perdew1996) 00228 CALL cite_reference(Zhang1998) 00229 IF (sx==1._dp .AND.sc==1._dp) THEN 00230 IF ( PRESENT ( reference ) ) THEN 00231 reference="revPBE, Yingkay Zhang and Weitao Yang,"//& 00232 " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)"//& 00233 "{spin polarized}" 00234 END IF 00235 IF ( PRESENT ( shortform ) ) THEN 00236 shortform = "revPBE, revised PBE xc functional (spin polarized)" 00237 END IF 00238 ELSE 00239 IF ( PRESENT ( reference ) ) THEN 00240 WRITE (reference,"(a,a,'sx=',f5.3,'sc=',f5.3,a)")& 00241 "revPBE, Yingkay Zhang and Weitao Yang,",& 00242 " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)",& 00243 sx,sc,"{spin polarized}" 00244 END IF 00245 IF ( PRESENT ( shortform ) ) THEN 00246 WRITE (shortform,"(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')")& 00247 "revPBE, revised PBE xc functional (spin polarized)",sx,sc 00248 END IF 00249 END IF 00250 CASE(xc_pbe_sol) 00251 CALL cite_reference(Perdew1996) 00252 CALL cite_reference(Perdew2008) 00253 IF (sx==1._dp .AND.sc==1._dp) THEN 00254 IF ( PRESENT ( reference ) ) THEN 00255 reference = "PBEsol, J.P. Perdew et al., "//& 00256 "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) "//& 00257 "{spin polarized}" 00258 END IF 00259 IF ( PRESENT ( shortform ) ) THEN 00260 shortform = "PBEsol, PBE for solids and surfaces xc functional (spin polarized)" 00261 END IF 00262 ELSE 00263 IF ( PRESENT ( reference ) ) THEN 00264 WRITE (reference,"(a,a,'sx=',f5.3,'sc=',f5.3,a)")& 00265 "PBEsol, J.P. Perdew et al., ",& 00266 "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) ",& 00267 sx,sc,"{spin unpolarized}" 00268 END IF 00269 IF ( PRESENT ( shortform ) ) THEN 00270 WRITE (shortform,"(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')")& 00271 "PBEsol, PBE for solids and surfaces xc functional (spin polarized)",sx,sc 00272 END IF 00273 END IF 00274 CASE default 00275 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00276 END SELECT 00277 IF (PRESENT(needs)) THEN 00278 needs%rho_spin=.TRUE. 00279 needs%norm_drho_spin=.TRUE. 00280 needs%norm_drho=.TRUE. 00281 END IF 00282 IF (PRESENT(max_deriv)) max_deriv=2 00283 00284 END SUBROUTINE pbe_lsd_info 00285 00286 ! ***************************************************************************** 00298 SUBROUTINE pbe_lda_eval(rho_set,deriv_set,grad_deriv,pbe_params,error) 00299 TYPE(xc_rho_set_type), POINTER :: rho_set 00300 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00301 INTEGER, INTENT(in) :: grad_deriv 00302 TYPE(section_vals_type), POINTER :: pbe_params 00303 TYPE(cp_error_type), INTENT(inout) :: error 00304 00305 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lda_eval', 00306 routineP = moduleN//':'//routineN 00307 00308 INTEGER :: handle, npoints, param, stat 00309 INTEGER, DIMENSION(:, :), POINTER :: bo 00310 LOGICAL :: failure 00311 REAL(kind=dp) :: epsilon_norm_drho, 00312 epsilon_rho, scale_ec, 00313 scale_ex 00314 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, 00315 e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, 00316 e_ndrho_rho_rho, e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho, 00317 rho_1_3 00318 TYPE(xc_derivative_type), POINTER :: deriv 00319 00320 CALL timeset(routineN,handle) 00321 00322 failure=.FALSE. 00323 NULLIFY(bo) 00324 00325 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00326 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00327 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00328 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00329 IF (.NOT. failure) THEN 00330 CALL xc_rho_set_get(rho_set,rho_1_3=rho_1_3,rho=rho,& 00331 norm_drho=norm_drho,local_bounds=bo,rho_cutoff=epsilon_rho,& 00332 drho_cutoff=epsilon_norm_drho,error=error) 00333 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00334 00335 ! meaningful default for the arrays we don't need: let us make compiler 00336 ! and debugger happy... 00337 IF (cp_debug) THEN 00338 ALLOCATE(dummy(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 00339 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00340 ELSE 00341 dummy=> rho 00342 END IF 00343 00344 e_0 => dummy 00345 e_rho => dummy 00346 e_ndrho => dummy 00347 e_rho_rho => dummy 00348 e_ndrho_rho => dummy 00349 e_ndrho_ndrho => dummy 00350 e_rho_rho_rho => dummy 00351 e_ndrho_rho_rho => dummy 00352 e_ndrho_ndrho_rho => dummy 00353 e_ndrho_ndrho_ndrho => dummy 00354 00355 IF (grad_deriv>=0) THEN 00356 deriv => xc_dset_get_derivative(deriv_set,"",& 00357 allocate_deriv=.TRUE., error=error) 00358 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00359 END IF 00360 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 00361 deriv => xc_dset_get_derivative(deriv_set,"(rho)",& 00362 allocate_deriv=.TRUE.,error=error) 00363 CALL xc_derivative_get(deriv,deriv_data=e_rho,error=error) 00364 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",& 00365 allocate_deriv=.TRUE.,error=error) 00366 CALL xc_derivative_get(deriv,deriv_data=e_ndrho,error=error) 00367 END IF 00368 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 00369 deriv => xc_dset_get_derivative(deriv_set,"(rho)(rho)",& 00370 allocate_deriv=.TRUE.,error=error) 00371 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho,error=error) 00372 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(rho)",& 00373 allocate_deriv=.TRUE.,error=error) 00374 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_rho,error=error) 00375 deriv => xc_dset_get_derivative(deriv_set,& 00376 "(norm_drho)(norm_drho)", allocate_deriv=.TRUE.,error=error) 00377 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho,error=error) 00378 END IF 00379 IF (grad_deriv>=3.OR.grad_deriv==-3) THEN 00380 deriv => xc_dset_get_derivative(deriv_set,"(rho)(rho)(rho)",& 00381 allocate_deriv=.TRUE.,error=error) 00382 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho_rho,error=error) 00383 deriv => xc_dset_get_derivative(deriv_set,& 00384 "(norm_drho)(rho)(rho)",allocate_deriv=.TRUE.,error=error) 00385 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_rho_rho,error=error) 00386 deriv => xc_dset_get_derivative(deriv_set,& 00387 "(norm_drho)(norm_drho)(rho)",allocate_deriv=.TRUE.,error=error) 00388 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho_rho,error=error) 00389 deriv => xc_dset_get_derivative(deriv_set,& 00390 "(norm_drho)(norm_drho)(norm_drho)", allocate_deriv=.TRUE.,& 00391 error=error) 00392 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho_ndrho,error=error) 00393 END IF 00394 IF (grad_deriv>3.OR.grad_deriv<-3) THEN 00395 CALL cp_unimplemented_error(fromWhere=routineP, & 00396 message="derivatives bigger than 3 not implemented", & 00397 error=error, error_level=cp_failure_level) 00398 END IF 00399 00400 00401 CALL section_vals_val_get(pbe_params,"scale_c",r_val=scale_ec,error=error) 00402 CALL section_vals_val_get(pbe_params,"scale_x",r_val=scale_ex,error=error) 00403 CALL section_vals_val_get(pbe_params,"parametrization",i_val=param,error=error) 00404 00405 !$omp parallel default(none), & 00406 !$omp shared(rho,rho_1_3,norm_drho,e_0,e_rho,e_ndrho,e_rho_rho,e_ndrho_rho),& 00407 !$omp shared(e_ndrho_ndrho,e_rho_rho_rho,e_ndrho_rho_rho,e_ndrho_ndrho_rho),& 00408 !$omp shared(e_ndrho_ndrho_ndrho,grad_deriv,npoints,epsilon_rho),& 00409 !$omp shared(epsilon_norm_drho,pbe_params),& 00410 !$omp shared(error,param,scale_ec,scale_ex) 00411 CALL pbe_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho,& 00412 e_0=e_0,e_rho=e_rho,e_ndrho=e_ndrho,e_rho_rho=e_rho_rho,& 00413 e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, & 00414 e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho,& 00415 e_ndrho_ndrho_rho=e_ndrho_ndrho_rho,e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho,& 00416 grad_deriv=grad_deriv,& 00417 npoints=npoints,epsilon_rho=epsilon_rho,epsilon_norm_drho=epsilon_norm_drho,& 00418 param=param,scale_ec=scale_ec,scale_ex=scale_ex,error=error) 00419 !$omp end parallel 00420 00421 IF (cp_debug) THEN 00422 DEALLOCATE(dummy,stat=stat) 00423 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00424 ELSE 00425 NULLIFY(dummy) 00426 END IF 00427 END IF 00428 00429 CALL timestop(handle) 00430 END SUBROUTINE pbe_lda_eval 00431 00432 ! ***************************************************************************** 00445 SUBROUTINE pbe_lda_calc(rho, rho_1_3, norm_drho,& 00446 e_0,e_rho,e_ndrho,e_rho_rho,e_ndrho_rho,& 00447 e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho,& 00448 e_ndrho_ndrho_ndrho,grad_deriv,npoints,epsilon_rho, epsilon_norm_drho,& 00449 ! pbe_params,error) 00450 param,scale_ec,scale_ex,error) 00451 INTEGER, INTENT(in) :: npoints, grad_deriv 00452 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: 00453 e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho_rho, e_rho_rho_rho, 00454 e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, e_ndrho, e_rho, e_0 00455 REAL(kind=dp), DIMENSION(1:npoints), 00456 INTENT(in) :: norm_drho, rho_1_3, rho 00457 REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_norm_drho 00458 INTEGER, INTENT(in) :: param 00459 REAL(kind=dp), INTENT(in) :: scale_ec, scale_ex 00460 TYPE(cp_error_type), INTENT(inout) :: error 00461 00462 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lda_calc', 00463 routineP = moduleN//':'//routineN 00464 00465 INTEGER :: ii 00466 LOGICAL :: failure 00467 REAL(kind=dp) :: A, A1rho, A1rhorho, A2rho, A_1, alpha_1_1, Arho, 00468 Arho1rho, Arhorho, Arhorhorho, beta, beta_1_1, beta_2_1, beta_3_1, 00469 beta_4_1, e_c_u_0, e_c_u_01rho, e_c_u_01rhorho, e_c_u_02rho, 00470 e_c_u_0rho, e_c_u_0rho1rho, e_c_u_0rhorho, e_c_u_0rhorhorho, 00471 epsilon_cGGA, epsilon_cGGArho, epsilon_cGGArhorho, ex_ldarhorhorho, 00472 ex_unif, ex_unif1rho, ex_unif1rhorho, ex_unif2rho, ex_unifrho, 00473 ex_unifrho1rho, ex_unifrhorho, Fx, Fx1rho, Fx1rhorho, Fx2rho, 00474 Fxnorm_drho, Fxnorm_drho1rho, Fxnorm_drhonorm_drho, Fxnorm_drhorho, 00475 Fxrho, Fxrho1rho, Fxrhorho, gamma_var, Hnorm_drho, Hnorm_drhonorm_drho, 00476 Hnorm_drhorho, k_f, k_f2rho, k_frho 00477 REAL(kind=dp) :: k_frhorho, k_frhorhorho, k_s, k_s1rho, k_s1rhorho, 00478 k_s2rho, k_srho, k_srho1rho, k_srhorho, kappa, kf, kf2rho, kfrho, 00479 kfrhorho, kfrhorhorho, mu, my_norm_drho, my_rho, p_1, p_2, rs, rs2rho, 00480 rsrho, rsrhorho, rsrhorhorho, s, s1rho, s1rhorho, s2norm_drho, s2rho, 00481 snorm_drho, snorm_drho1rho, snorm_drhorho, srho, srho1rho, srhorho, t, 00482 t1, t1001, t1004, t1005, t1006, t1008, t101, t1011, t1012, t1013, 00483 t1014, t1016, t1017, t1018, t1019, t1022, t1024, t1026, t1028, t103, 00484 t1030, t1031, t1035, t1042, t1046, t1048, t105, t1050, t1052, t1054, 00485 t1055, t1060, t1061, t1062, t1067, t108, t1093, t1097, t11, t1103, t1104 00486 REAL(kind=dp) :: t1106, t1109, t111, t1115, t1118, t1121, t1122, t1126, 00487 t1129, t113, t114, t1148, t115, t1152, t1157, t1167, t1187, t119, 00488 t1196, t1203, t1218, t1226, t123, t124, t1249, t125, t126, t1262, 00489 t1263, t127, t1291, t1292, t1295, t13, t131, t1329, t1342, t135, t138, 00490 t1380, t1385, t1386, t1387, t1389, t139, t1393, t14, t140, t142, t146, 00491 t1468, t148, t1482, t149, t1492, t1493, t150, t1504, t1505, t151, 00492 t1511, t1515, t1521, t1525, t1528, t153, t1532, t1545, t155, t1565, 00493 t157, t1573, t158, t1584, t159, t1608, t1612, t162, t1628, t1629, t163, 00494 t1633, t164, t1646, t1652, t1660, t167, t1672, t1676, t168, t17, t170 00495 REAL(kind=dp) :: t171, t1715, t1722, t175, t1758, t1759, t176, t1761, 00496 t177, t178, t179, t1797, t182, t183, t1838, t1851, t186, t187, t1878, 00497 t1889, t189, t19, t190, t1907, t191, t1922, t1927, t1933, t1937, t195, 00498 t1952, t196, t1964, t1965, t1968, t1969, t197, t1972, t198, t1990, 00499 t1rho, t1rhorho, t2, t20, t200, t202, t2020, t2024, t2028, t2031, t204, 00500 t2041, t208, t21, t214, t217, t218, t22, t226, t229, t230, t238, t239, 00501 t240, t241, t246, t252, t253, t255, t259, t26, t260, t261, t262, t265, 00502 t266, t267, t27, t271, t272, t273, t277, t278, t280, t281, t286, t290, 00503 t291, t293, t294, t295, t296, t297, t299, t2norm_drho 00504 REAL(kind=dp) :: t2rho, t3, t305, t309, t310, t315, t317, t318, t321, 00505 t323, t324, t327, t329, t330, t331, t336, t338, t339, t340, t341, t348, 00506 t349, t354, t357, t359, t360, t361, t362, t369, t370, t371, t374, t377, 00507 t378, t381, t382, t384, t385, t387, t388, t390, t392, t393, t397, t4, 00508 t400, t401, t404, t408, t409, t410, t414, t417, t418, t423, t426, t427, 00509 t432, t435, t436, t438, t439, t440, t448, t449, t451, t456, t457, t458, 00510 t459, t461, t462, t463, t465, t466, t469, t470, t471, t472, t476, t487, 00511 t491, t496, t498, t5, t500, t503, t506, t507, t510, t513, t517, t521, 00512 t525, t529, t530, t535, t541, t548, t553, t556 00513 REAL(kind=dp) :: t557, t559, t562, t566, t581, t586, t590, t591, t594, 00514 t598, t6, t605, t609, t611, t612, t614, t627, t645, t65, t654, t656, 00515 t661, t664, t669, t670, t671, t673, t675, t685, t69, t693, t7, t70, 00516 t701, t702, t71, t714, t717, t72, t720, t73, t74, t740, t743, t748, 00517 t75, t758, t77, t776, t78, t8, t80, t801, t809, t81, t812, t821, t825, 00518 t83, t831, t84, t85, t86, t868, t87, t877, t879, t88, t880, t885, t90, 00519 t91, t94, t940, t942, t943, t945, t946, t948, t95, t950, t951, t954, 00520 t967, t976, t98, t980, t982, t984, t989, t99, t990, t994, t995, t998, 00521 t999, tnorm_drho, tnorm_drho1rho, tnorm_drhorho, tnorm_drhorhorho 00522 REAL(kind=dp) :: trho, trho1rho, trhorho, trhorhorho 00523 00524 !TYPE(section_vals_type), POINTER :: pbe_params 00525 !, param 00526 ! scale_ec, scale_ex, snorm_drho, snorm_drho1rho, snorm_drhorho, srho, & 00527 00528 failure=.FALSE. 00529 00530 ! Parametrization of PBE 00531 SELECT CASE(param) 00532 ! Original PBE 00533 CASE (xc_pbe_orig) 00534 kappa = 0.804e0_dp 00535 beta = 0.66725e-1_dp 00536 mu = beta * pi ** 2 / 0.3e1_dp 00537 ! Revised PBE (revPBE) 00538 CASE (xc_pbe_rev) 00539 kappa = 0.1245e1_dp 00540 beta = 0.66725e-1_dp 00541 mu = beta * pi ** 2 / 0.3e1_dp 00542 ! PBE for solids (PBEsol) 00543 CASE (xc_pbe_sol) 00544 kappa = 0.804e0_dp 00545 beta = 0.46e-1_dp 00546 mu = 0.1e2_dp / 0.81e2_dp 00547 CASE default 00548 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00549 END SELECT 00550 00551 gamma_var = (0.1e1_dp - LOG(0.2e1_dp)) / pi ** 2 00552 p_1 = 0.10e1_dp 00553 A_1 = 0.31091e-1_dp 00554 alpha_1_1 = 0.21370e0_dp 00555 beta_1_1 = 0.75957e1_dp 00556 beta_2_1 = 0.35876e1_dp 00557 beta_3_1 = 0.16382e1_dp 00558 beta_4_1 = 0.49294e0_dp 00559 p_2 = 0.10e1_dp 00560 t1 = 3 ** (0.1e1_dp / 0.3e1_dp) 00561 t2 = 4 ** (0.1e1_dp / 0.3e1_dp) 00562 t3 = t2 ** 2 00563 t4 = t1 * t3 00564 t5 = 0.1e1_dp / pi 00565 t69 = pi ** 2 00566 t627 = 0.1e1_dp / t69 / pi 00567 00568 SELECT CASE(grad_deriv) 00569 CASE default 00570 !$omp do 00571 DO ii=1,npoints 00572 my_rho = rho(ii) 00573 IF (my_rho>epsilon_rho) THEN 00574 my_norm_drho = norm_drho(ii) 00575 00576 t6 = 0.1e1_dp / my_rho 00577 t7 = t5 * t6 00578 t8 = t7 ** (0.1e1_dp / 0.3e1_dp) 00579 rs = t4 * t8 / 0.4e1_dp 00580 t11 = 0.1e1_dp + alpha_1_1 * rs 00581 t13 = 0.1e1_dp / A_1 00582 t14 = SQRT(rs) 00583 t17 = t14 * rs 00584 t19 = p_1 + 0.1e1_dp 00585 t20 = rs ** t19 00586 t21 = beta_4_1 * t20 00587 t22 = beta_1_1 * t14 + beta_2_1 * rs + beta_3_1 * t17 + t21 00588 t26 = 0.1e1_dp + t13 / t22 / 0.2e1_dp 00589 t27 = LOG(t26) 00590 e_c_u_0 = -0.2e1_dp * A_1 * t11 * t27 00591 t65 = 2 ** (0.1e1_dp / 0.3e1_dp) 00592 t70 = t69 * my_rho 00593 t71 = t70 ** (0.1e1_dp / 0.3e1_dp) 00594 k_f = t1 * t71 00595 t72 = k_f * t5 00596 t73 = SQRT(t72) 00597 k_s = 0.2e1_dp * t73 00598 t74 = 0.1e1_dp / k_s 00599 t75 = my_norm_drho * t74 00600 t = t75 * t6 / 0.2e1_dp 00601 t77 = 0.1e1_dp / gamma_var 00602 t78 = beta * t77 00603 t80 = EXP(-e_c_u_0 * t77) 00604 t81 = -0.1e1_dp + t80 00605 A = t78 / t81 00606 t83 = t ** 2 00607 t84 = A * t83 00608 t85 = 0.1e1_dp + t84 00609 t86 = t83 * t85 00610 t87 = A ** 2 00611 t88 = t83 ** 2 00612 t90 = 0.1e1_dp + t84 + t87 * t88 00613 t91 = 0.1e1_dp / t90 00614 t94 = 0.1e1_dp + t78 * t86 * t91 00615 t95 = LOG(t94) 00616 epsilon_cGGA = e_c_u_0 + gamma_var * t95 00617 kf = k_f 00618 ex_unif = -0.3e1_dp / 0.4e1_dp * t5 * kf 00619 t98 = 0.1e1_dp / kf 00620 t99 = my_norm_drho * t98 00621 s = t99 * t6 / 0.2e1_dp 00622 t101 = s ** 2 00623 t103 = 0.1e1_dp / kappa 00624 t105 = 0.1e1_dp + mu * t101 * t103 00625 Fx = 0.1e1_dp + kappa - kappa / t105 00626 t108 = my_rho * ex_unif 00627 00628 IF (grad_deriv>=0) THEN 00629 e_0(ii) = e_0(ii)+& 00630 scale_ex * t108 * Fx + scale_ec * my_rho * epsilon_cGGA 00631 END IF 00632 00633 t111 = t8 ** 2 00634 t113 = 0.1e1_dp / t111 * t5 00635 t114 = my_rho ** 2 00636 t115 = 0.1e1_dp / t114 00637 rsrho = -t4 * t113 * t115 / 0.12e2_dp 00638 t119 = A_1 * alpha_1_1 00639 t123 = t22 ** 2 00640 t124 = 0.1e1_dp / t123 00641 t125 = t11 * t124 00642 t126 = 0.1e1_dp / t14 00643 t127 = beta_1_1 * t126 00644 t131 = beta_3_1 * t14 00645 t135 = 0.1e1_dp / rs 00646 t138 = t127 * rsrho / 0.2e1_dp + beta_2_1 * rsrho + 0.3e1_dp / & 00647 0.2e1_dp * t131 * rsrho + t21 * t19 * rsrho * t135 00648 t139 = 0.1e1_dp / t26 00649 t140 = t138 * t139 00650 e_c_u_0rho = -0.2e1_dp * t119 * rsrho * t27 + t125 * t140 00651 t142 = t71 ** 2 00652 k_frho = t1 / t142 * t69 / 0.3e1_dp 00653 t146 = 0.1e1_dp / t73 00654 k_srho = t146 * k_frho * t5 00655 t148 = k_s ** 2 00656 t149 = 0.1e1_dp / t148 00657 t150 = my_norm_drho * t149 00658 t151 = t6 * k_srho 00659 t153 = t75 * t115 00660 trho = -t150 * t151 / 0.2e1_dp - t153 / 0.2e1_dp 00661 t155 = gamma_var ** 2 00662 t157 = beta / t155 00663 t158 = t81 ** 2 00664 t159 = 0.1e1_dp / t158 00665 Arho = t157 * t159 * e_c_u_0rho * t80 00666 t162 = t78 * t 00667 t163 = t85 * t91 00668 t164 = t163 * trho 00669 t167 = Arho * t83 00670 t168 = A * t 00671 t170 = 0.2e1_dp * t168 * trho 00672 t171 = t167 + t170 00673 t175 = t78 * t83 00674 t176 = t90 ** 2 00675 t177 = 0.1e1_dp / t176 00676 t178 = t85 * t177 00677 t179 = A * t88 00678 t182 = t83 * t 00679 t183 = t87 * t182 00680 t186 = t167 + t170 + 0.2e1_dp * t179 * Arho + 0.4e1_dp * t183 * trho 00681 t187 = t178 * t186 00682 t189 = 0.2e1_dp * t162 * t164 + t78 * t83 * t171 * t91 - t175 * t187 00683 t190 = gamma_var * t189 00684 t191 = 0.1e1_dp / t94 00685 epsilon_cGGArho = e_c_u_0rho + t190 * t191 00686 kfrho = k_frho 00687 ex_unifrho = -0.3e1_dp / 0.4e1_dp * t5 * kfrho 00688 t195 = kf ** 2 00689 t196 = 0.1e1_dp / t195 00690 t197 = my_norm_drho * t196 00691 t198 = t6 * kfrho 00692 t200 = t99 * t115 00693 srho = -t197 * t198 / 0.2e1_dp - t200 / 0.2e1_dp 00694 t202 = t105 ** 2 00695 t204 = 0.1e1_dp / t202 * mu 00696 Fxrho = 0.2e1_dp * t204 * s * srho 00697 t208 = my_rho * ex_unifrho 00698 00699 IF (grad_deriv>=1 .OR. grad_deriv==-1) THEN 00700 e_rho(ii) = e_rho(ii)+& 00701 scale_ex * (ex_unif * Fx + t208 * Fx + t108 * Fxrho) + & 00702 scale_ec * (epsilon_cGGA + my_rho * epsilon_cGGArho) 00703 END IF 00704 00705 tnorm_drho = t74 * t6 / 0.2e1_dp 00706 t214 = t163 * tnorm_drho 00707 t217 = t78 * t182 00708 t218 = A * tnorm_drho 00709 t226 = 0.2e1_dp * t168 * tnorm_drho + 0.4e1_dp * t183 * tnorm_drho 00710 t229 = 0.2e1_dp * t162 * t214 + 0.2e1_dp * t217 * t218 * t91 - & 00711 t175 * t178 * t226 00712 t230 = gamma_var * t229 00713 Hnorm_drho = t230 * t191 00714 snorm_drho = t98 * t6 / 0.2e1_dp 00715 Fxnorm_drho = 0.2e1_dp * t204 * s * snorm_drho 00716 00717 IF (grad_deriv>=1 .OR. grad_deriv==-1) THEN 00718 e_ndrho(ii) = e_ndrho(ii)+& 00719 scale_ex * t108 * Fxnorm_drho + scale_ec * my_rho * & 00720 Hnorm_drho 00721 END IF 00722 00723 t238 = 0.1e1_dp / t69 00724 t239 = 0.1e1_dp / t111 / t7 * t238 00725 t240 = t114 ** 2 00726 t241 = 0.1e1_dp / t240 00727 t246 = 0.1e1_dp / t114 / my_rho 00728 rsrhorho = -t4 * t239 * t241 / 0.18e2_dp + t4 * t113 *& 00729 t246 / 0.6e1_dp 00730 t252 = 0.2e1_dp * t119 * rsrhorho * t27 00731 t253 = alpha_1_1 * rsrho 00732 t255 = t124 * t138 * t139 00733 t259 = 0.1e1_dp / t123 / t22 00734 t260 = t11 * t259 00735 t261 = t138 ** 2 00736 t262 = t261 * t139 00737 t265 = 0.1e1_dp / t17 00738 t266 = beta_1_1 * t265 00739 t267 = rsrho ** 2 00740 t271 = t127 * rsrhorho / 0.2e1_dp 00741 t272 = beta_2_1 * rsrhorho 00742 t273 = beta_3_1 * t126 00743 t277 = 0.3e1_dp / 0.2e1_dp * t131 * rsrhorho 00744 t278 = t19 ** 2 00745 t280 = rs ** 2 00746 t281 = 0.1e1_dp / t280 00747 t286 = t21 * t19 * rsrhorho * t135 00748 t290 = -t266 * t267 / 0.4e1_dp + t271 + t272 + 0.3e1_dp / 0.4e1_dp& 00749 * t273 * t267 + t277 + t21 * t278 * t267 * t281 + t286 - t21 * t19 & 00750 * t267 * t281 00751 t291 = t290 * t139 00752 t293 = t123 ** 2 00753 t294 = 0.1e1_dp / t293 00754 t295 = t11 * t294 00755 t296 = t26 ** 2 00756 t297 = 0.1e1_dp / t296 00757 t299 = t261 * t297 * t13 00758 e_c_u_0rhorho = -t252 + 0.2e1_dp * t253 * t255 - 0.2e1_dp * t260 *& 00759 t262 + t125 * t291 + t295 * t299 / 0.2e1_dp 00760 e_c_u_01rho = e_c_u_0rho 00761 t305 = t69 ** 2 00762 k_frhorho = -0.2e1_dp / 0.9e1_dp * t1 / t142 / t70 * t305 00763 t309 = 0.1e1_dp / t73 / t72 00764 t310 = k_frho ** 2 00765 t315 = t146 * k_frhorho * t5 00766 k_srhorho = -t309 * t310 * t238 / 0.2e1_dp + t315 00767 k_s1rho = k_srho 00768 t317 = 0.1e1_dp / t148 / k_s 00769 t318 = my_norm_drho * t317 00770 t321 = t115 * k_srho 00771 t323 = t150 * t321 / 0.2e1_dp 00772 t324 = t6 * k_srhorho 00773 t327 = t115 * k_s1rho 00774 t329 = t150 * t327 / 0.2e1_dp 00775 t330 = t75 * t246 00776 trhorho = t318 * t151 * k_s1rho + t323 - t150 * t324 / 0.2e1_dp + & 00777 t329 + t330 00778 t331 = t6 * k_s1rho 00779 t1rho = -t150 * t331 / 0.2e1_dp - t153 / 0.2e1_dp 00780 t336 = beta / t155 / gamma_var 00781 t338 = 0.1e1_dp / t158 / t81 00782 t339 = t336 * t338 00783 t340 = t80 ** 2 00784 t341 = e_c_u_0rho * t340 00785 t348 = t336 * t159 00786 t349 = e_c_u_0rho * e_c_u_01rho 00787 Arhorho = 0.2e1_dp * t339 * t341 * e_c_u_01rho + t157 * t159 * & 00788 e_c_u_0rhorho * t80 - t348 * t349 * t80 00789 A1rho = t157 * t159 * e_c_u_01rho * t80 00790 t354 = t78 * t1rho 00791 t357 = A1rho * t83 00792 t359 = 0.2e1_dp * t168 * t1rho 00793 t360 = t357 + t359 00794 t361 = t360 * t91 00795 t362 = t361 * trho 00796 t369 = t357 + t359 + 0.2e1_dp * t179 * A1rho + 0.4e1_dp * t183 * t1rho 00797 t370 = trho * t369 00798 t371 = t178 * t370 00799 t374 = t163 * trhorho 00800 t377 = t171 * t91 00801 t378 = t377 * t1rho 00802 t381 = Arhorho * t83 00803 t382 = Arho * t 00804 t384 = 0.2e1_dp * t382 * t1rho 00805 t385 = A1rho * t 00806 t387 = 0.2e1_dp * t385 * trho 00807 t388 = A * t1rho 00808 t390 = 0.2e1_dp * t388 * trho 00809 t392 = 0.2e1_dp * t168 * trhorho 00810 t393 = t381 + t384 + t387 + t390 + t392 00811 t397 = t171 * t177 00812 t400 = t186 * t1rho 00813 t401 = t178 * t400 00814 t404 = t360 * t177 00815 t408 = 0.1e1_dp / t176 / t90 00816 t409 = t85 * t408 00817 t410 = t186 * t369 00818 t414 = A1rho * t88 00819 t417 = A * t182 00820 t418 = Arho * t1rho 00821 t423 = trho * A1rho 00822 t426 = t87 * t83 00823 t427 = trho * t1rho 00824 t432 = t381 + t384 + t387 + t390 + t392 + 0.2e1_dp * t414 * Arho +& 00825 0.8e1_dp * t417 * t418 + 0.2e1_dp * t179 * Arhorho + 0.8e1_dp * & 00826 t417 * t423 + 0.12e2_dp * t426 * t427 + 0.4e1_dp * t183 * trhorho 00827 t435 = 0.2e1_dp * t354 * t164 + 0.2e1_dp * t162 * t362 - 0.2e1_dp & 00828 * t162 * t371 + 0.2e1_dp * t162 * t374 + 0.2e1_dp * t162 * t378 + & 00829 t78 * t83 * t393 * t91 - t175 * t397 * t369 - 0.2e1_dp * t162 * t401& 00830 - t175 * t404 * t186 + 0.2e1_dp * t175 * t409 * t410 - t175 * t178 & 00831 * t432 00832 t436 = gamma_var * t435 00833 t438 = t94 ** 2 00834 t439 = 0.1e1_dp / t438 00835 t440 = t163 * t1rho 00836 t448 = 0.2e1_dp * t162 * t440 + t78 * t83 * t360 * t91 - t175 * & 00837 t178 * t369 00838 t449 = t439 * t448 00839 t451 = gamma_var * t448 00840 epsilon_cGGArhorho = e_c_u_0rhorho + t436 * t191 - t190 * t449 00841 kfrhorho = k_frhorho 00842 ex_unifrhorho = -0.3e1_dp / 0.4e1_dp * t5 * kfrhorho 00843 ex_unif1rho = ex_unifrho 00844 t456 = 0.1e1_dp / t195 / kf 00845 t457 = my_norm_drho * t456 00846 t458 = kfrho ** 2 00847 t459 = t6 * t458 00848 t461 = t115 * kfrho 00849 t462 = t197 * t461 00850 t463 = t6 * kfrhorho 00851 t465 = t197 * t463 / 0.2e1_dp 00852 t466 = t99 * t246 00853 srhorho = t457 * t459 + t462 - t465 + t466 00854 s1rho = srho 00855 t469 = mu ** 2 00856 t470 = 0.1e1_dp / t202 / t105 * t469 00857 t471 = t470 * t101 00858 t472 = srho * t103 00859 t476 = s1rho * srho 00860 Fxrhorho = -0.8e1_dp * t471 * t472 * s1rho + 0.2e1_dp * t204 * & 00861 t476 + 0.2e1_dp * t204 * s * srhorho 00862 Fx1rho = 0.2e1_dp * t204 * s * s1rho 00863 t487 = my_rho * ex_unifrhorho 00864 t491 = my_rho * ex_unif1rho 00865 00866 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 00867 e_rho_rho(ii) = e_rho_rho(ii)+& 00868 scale_ex * (ex_unif1rho * Fx + ex_unif * Fx1rho + & 00869 ex_unifrho * Fx + t487 * Fx + t208 * Fx1rho + ex_unif * Fxrho + t491& 00870 * Fxrho + t108 * Fxrhorho) + scale_ec * (e_c_u_01rho + t451 * t191 & 00871 + epsilon_cGGArho + my_rho * epsilon_cGGArhorho) 00872 END IF 00873 00874 t496 = t149 * t6 00875 t498 = t74 * t115 00876 tnorm_drhorho = -t496 * k_srho / 0.2e1_dp - t498 / 0.2e1_dp 00877 t500 = t78 * trho 00878 t503 = t377 * tnorm_drho 00879 t506 = tnorm_drho * t186 00880 t507 = t178 * t506 00881 t510 = t163 * tnorm_drhorho 00882 t513 = t91 * trho 00883 t517 = Arho * tnorm_drho 00884 t521 = A * tnorm_drhorho 00885 t525 = t177 * t186 00886 t529 = t226 * trho 00887 t530 = t178 * t529 00888 t535 = t226 * t186 00889 t541 = A * trho 00890 t548 = tnorm_drho * trho 00891 t553 = 0.2e1_dp * t382 * tnorm_drho + 0.2e1_dp * t541 * tnorm_drho& 00892 + 0.2e1_dp * t168 * tnorm_drhorho + 0.8e1_dp * t417 * t517 + & 00893 0.12e2_dp * t426 * t548 + 0.4e1_dp * t183 * tnorm_drhorho 00894 t556 = 0.2e1_dp * t500 * t214 + 0.2e1_dp * t162 * t503 - 0.2e1_dp & 00895 * t162 * t507 + 0.2e1_dp * t162 * t510 + 0.6e1_dp * t175 * t218 * & 00896 t513 + 0.2e1_dp * t217 * t517 * t91 + 0.2e1_dp * t217 * t521 * t91 -& 00897 0.2e1_dp * t217 * t218 * t525 - 0.2e1_dp * t162 * t530 - t175 * & 00898 t397 * t226 + 0.2e1_dp * t175 * t409 * t535 - t175 * t178 * t553 00899 t557 = gamma_var * t556 00900 t559 = t439 * t189 00901 Hnorm_drhorho = t557 * t191 - t230 * t559 00902 t562 = t196 * t6 00903 snorm_drhorho = -t562 * kfrho / 0.2e1_dp - t98 * t115 / 0.2e1_dp 00904 t566 = snorm_drho * t103 00905 Fxnorm_drhorho = -0.8e1_dp * t471 * t566 * srho + 0.2e1_dp * t204 & 00906 * srho * snorm_drho + 0.2e1_dp * t204 * s * snorm_drhorho 00907 00908 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 00909 e_ndrho_rho(ii) = e_ndrho_rho(ii)+& 00910 scale_ex * (ex_unif * Fxnorm_drho + t208 * & 00911 Fxnorm_drho + t108 * Fxnorm_drhorho) + scale_ec * (Hnorm_drho + my_rho & 00912 * Hnorm_drhorho) 00913 END IF 00914 00915 t581 = tnorm_drho ** 2 00916 t586 = A * t581 00917 t590 = tnorm_drho * t226 00918 t591 = t178 * t590 00919 t594 = t177 * t226 00920 t598 = t226 ** 2 00921 t605 = 0.2e1_dp * t586 + 0.12e2_dp * t426 * t581 00922 t609 = gamma_var * (0.2e1_dp * t78 * t581 * t85 * t91 + 0.10e2_dp & 00923 * t175 * t586 * t91 - 0.4e1_dp * t162 * t591 - 0.4e1_dp * t217 * & 00924 t218 * t594 + 0.2e1_dp * t175 * t409 * t598 - t175 * t178 * t605) 00925 t611 = t229 ** 2 00926 t612 = gamma_var * t611 00927 Hnorm_drhonorm_drho = t609 * t191 - t612 * t439 00928 t614 = snorm_drho ** 2 00929 Fxnorm_drhonorm_drho = -0.8e1_dp * t470 * t101 * t614 * t103 + & 00930 0.2e1_dp * t204 * t614 00931 00932 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 00933 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii)+& 00934 scale_ex * t108 * Fxnorm_drhonorm_drho +& 00935 scale_ec * my_rho * Hnorm_drhonorm_drho 00936 END IF 00937 00938 IF (grad_deriv>=3.OR.grad_deriv==-3) THEN 00939 rsrhorhorho = -0.5e1_dp / 0.54e2_dp * t4 / t111 / t238 / & 00940 t115 * t627 / t240 / t114 + t4 * t239 / t240 / my_rho / 0.3e1_dp & 00941 - t4 * t113 * t241 / 0.2e1_dp 00942 rs2rho = rsrho 00943 t645 = alpha_1_1 * rsrhorho 00944 t654 = t127 * rs2rho / 0.2e1_dp + beta_2_1 * rs2rho + 0.3e1_dp / & 00945 0.2e1_dp * t131 * rs2rho + t21 * t19 * rs2rho * t135 00946 t656 = t124 * t654 * t139 00947 t661 = t140 * t654 00948 t664 = rsrho * rs2rho 00949 t669 = t21 * t278 00950 t670 = rs2rho * t281 00951 t671 = t670 * rsrho 00952 t673 = t21 * t19 00953 t675 = -t266 * t664 / 0.4e1_dp + t271 + t272 + 0.3e1_dp / 0.4e1_dp& 00954 * t273 * t664 + t277 + t669 * t671 + t286 - t673 * t671 00955 t685 = alpha_1_1 * rs2rho 00956 t693 = t675 * t139 00957 t701 = t297 * t13 00958 t702 = t701 * t654 00959 t714 = t267 * rs2rho 00960 t717 = rsrhorho * rsrho 00961 t720 = rsrhorho * rs2rho 00962 t740 = rs2rho / t280 / rs * t267 00963 t743 = rsrhorho * t281 * rsrho 00964 t748 = t670 * rsrhorho 00965 t758 = 0.3e1_dp / 0.8e1_dp * beta_1_1 / t14 / t280 * t714 - t266 *& 00966 t717 / 0.2e1_dp - t266 * t720 / 0.4e1_dp + t127 * rsrhorhorho / & 00967 0.2e1_dp + beta_2_1 * rsrhorhorho - 0.3e1_dp / 0.8e1_dp * beta_3_1 *& 00968 t265 * t714 + 0.3e1_dp / 0.2e1_dp * t273 * t717 + 0.3e1_dp / & 00969 0.4e1_dp * t273 * t720 + 0.3e1_dp / 0.2e1_dp * t131 * rsrhorhorho + & 00970 t21 * t278 * t19 * t740 + 0.2e1_dp * t669 * t743 - 0.3e1_dp * t669 *& 00971 t740 + t669 * t748 + t21 * t19 * rsrhorhorho * t135 - t673 * t748 -& 00972 0.2e1_dp * t673 * t743 + 0.2e1_dp * t673 * t740 00973 t776 = A_1 ** 2 00974 e_c_u_0rhorhorho = -0.2e1_dp * t119 * rsrhorhorho * t27 + t645 * & 00975 t656 + 0.2e1_dp * t645 * t255 - 0.4e1_dp * t253 * t259 * t661 + & 00976 0.2e1_dp * t253 * t124 * t675 * t139 + t253 * t294 * t138 * t297 * & 00977 t13 * t654 - 0.2e1_dp * t685 * t259 * t261 * t139 + 0.6e1_dp * t295 & 00978 * t262 * t654 - 0.4e1_dp * t260 * t693 * t138 - 0.3e1_dp * t11 / & 00979 t293 / t22 * t261 * t702 + t685 * t124 * t290 * t139 - 0.2e1_dp * & 00980 t260 * t291 * t654 + t125 * t758 * t139 + t295 * t290 * t702 / & 00981 0.2e1_dp + t685 * t294 * t299 / 0.2e1_dp + t295 * t675 * t701 * t138& 00982 + t11 / t293 / t123 * t261 / t296 / t26 / t776 * t654 / 0.2e1_dp 00983 e_c_u_0rho1rho = -t252 + t253 * t656 + t685 * t255 - 0.2e1_dp * & 00984 t260 * t661 + t125 * t693 + t295 * t138 * t702 / 0.2e1_dp 00985 e_c_u_01rhorho = e_c_u_0rho1rho 00986 e_c_u_02rho = -0.2e1_dp * t119 * rs2rho * t27 + t125 * t654 * t139 00987 k_frhorhorho = 0.10e2_dp / 0.27e2_dp * t1 / t142 / t114 * t69 00988 k_f2rho = kfrho 00989 t801 = k_f ** 2 00990 t809 = t309 * k_frhorho 00991 t812 = t238 * k_f2rho 00992 k_srho1rho = -t309 * k_frho * t812 / 0.2e1_dp + t315 00993 k_s1rhorho = k_srho1rho 00994 k_s2rho = t146 * k_f2rho * t5 00995 t821 = t148 ** 2 00996 t825 = k_srho * k_s1rho 00997 t831 = t6 * k_srho1rho 00998 trhorhorho = -0.3e1_dp * my_norm_drho / t821 * t6 * t825 * k_s2rho - & 00999 t318 * t321 * k_s1rho + t318 * t831 * k_s1rho + t318 * t151 * & 01000 k_s1rhorho - t318 * t321 * k_s2rho - t150 * t246 * k_srho + t150 * & 01001 t115 * k_srho1rho / 0.2e1_dp + t318 * t324 * k_s2rho + t150 * t115 *& 01002 k_srhorho / 0.2e1_dp - t150 * t6 * (0.3e1_dp / 0.4e1_dp / t73 / & 01003 t801 / t238 * t310 * t627 * k_f2rho - t809 * t238 * k_frho - t809 * & 01004 t812 / 0.2e1_dp + t146 * k_frhorhorho * t5) / 0.2e1_dp - t318 * t327& 01005 * k_s2rho - t150 * t246 * k_s1rho + t150 * t115 * k_s1rhorho / & 01006 0.2e1_dp - t150 * t246 * k_s2rho - 0.3e1_dp * t75 * t241 01007 t868 = t150 * t115 * k_s2rho / 0.2e1_dp 01008 trho1rho = t318 * t151 * k_s2rho + t323 - t150 * t831 / 0.2e1_dp +& 01009 t868 + t330 01010 t1rhorho = t318 * t331 * k_s2rho + t329 - t150 * t6 * k_s1rhorho /& 01011 0.2e1_dp + t868 + t330 01012 t2rho = -t150 * t6 * k_s2rho / 0.2e1_dp - t153 / 0.2e1_dp 01013 t877 = t155 ** 2 01014 t879 = beta / t877 01015 t880 = t158 ** 2 01016 t885 = e_c_u_01rho * e_c_u_02rho 01017 Arhorhorho = 0.6e1_dp * t879 / t880 * e_c_u_0rho * t340 * t80 * & 01018 t885 + 0.2e1_dp * t339 * e_c_u_0rho1rho * t340 * e_c_u_01rho - & 01019 0.6e1_dp * t879 * t338 * t341 * t885 + 0.2e1_dp * t339 * t341 * & 01020 e_c_u_01rhorho + 0.2e1_dp * t339 * e_c_u_0rhorho * t340 * & 01021 e_c_u_02rho + t157 * t159 * e_c_u_0rhorhorho * t80 - t348 * & 01022 e_c_u_0rhorho * e_c_u_02rho * t80 - t348 * e_c_u_0rho1rho * & 01023 e_c_u_01rho * t80 - t348 * e_c_u_0rho * e_c_u_01rhorho * t80 + t879 & 01024 * t159 * t349 * e_c_u_02rho * t80 01025 Arho1rho = 0.2e1_dp * t339 * t341 * e_c_u_02rho + t157 * t159 * & 01026 e_c_u_0rho1rho * t80 - t348 * e_c_u_0rho * e_c_u_02rho * t80 01027 A1rhorho = 0.2e1_dp * t339 * e_c_u_01rho * t340 * e_c_u_02rho + & 01028 t157 * t159 * e_c_u_01rhorho * t80 - t348 * t885 * t80 01029 A2rho = t157 * t159 * e_c_u_02rho * t80 01030 t940 = Arho1rho * t83 01031 t942 = 0.2e1_dp * t382 * t2rho 01032 t943 = A2rho * t 01033 t945 = 0.2e1_dp * t943 * trho 01034 t946 = A * t2rho 01035 t948 = 0.2e1_dp * t946 * trho 01036 t950 = 0.2e1_dp * t168 * trho1rho 01037 t951 = A2rho * t88 01038 t954 = Arho * t2rho 01039 t967 = t940 + t942 + t945 + t948 + t950 + 0.2e1_dp * t951 * Arho +& 01040 0.8e1_dp * t417 * t954 + 0.2e1_dp * t179 * Arho1rho + 0.8e1_dp * & 01041 t417 * trho * A2rho + 0.12e2_dp * t426 * trho * t2rho + 0.4e1_dp * & 01042 t183 * trho1rho 01043 t976 = t78 * t2rho 01044 t980 = t78 * t * t85 01045 t982 = A2rho * t83 01046 t984 = 0.2e1_dp * t168 * t2rho 01047 t989 = t982 + t984 + 0.2e1_dp * t179 * A2rho + 0.4e1_dp * t183 * t2rho 01048 t990 = t369 * t989 01049 t994 = t982 + t984 01050 t995 = t994 * t177 01051 t998 = Arhorhorho * t83 01052 t999 = Arhorho * t 01053 t1001 = 0.2e1_dp * t999 * t2rho 01054 t1004 = 0.2e1_dp * Arho1rho * t * t1rho 01055 t1005 = t954 * t1rho 01056 t1006 = 0.2e1_dp * t1005 01057 t1008 = 0.2e1_dp * t382 * t1rhorho 01058 t1011 = 0.2e1_dp * A1rhorho * t * trho 01059 t1012 = A1rho * t2rho 01060 t1013 = t1012 * trho 01061 t1014 = 0.2e1_dp * t1013 01062 t1016 = 0.2e1_dp * t385 * trho1rho 01063 t1017 = A2rho * t1rho 01064 t1018 = t1017 * trho 01065 t1019 = 0.2e1_dp * t1018 01066 t1022 = 0.2e1_dp * A * t1rhorho * trho 01067 t1024 = 0.2e1_dp * t388 * trho1rho 01068 t1026 = 0.2e1_dp * t943 * trhorho 01069 t1028 = 0.2e1_dp * t946 * trhorho 01070 t1030 = 0.2e1_dp * t168 * trhorhorho 01071 t1031 = t998 + t1001 + t1004 + t1006 + t1008 + t1011 + t1014 + & 01072 t1016 + t1019 + t1022 + t1024 + t1026 + t1028 + t1030 01073 t1035 = t940 + t942 + t945 + t948 + t950 01074 t1042 = t369 * t2rho 01075 t1046 = A1rhorho * t83 01076 t1048 = 0.2e1_dp * t385 * t2rho 01077 t1050 = 0.2e1_dp * t943 * t1rho 01078 t1052 = 0.2e1_dp * t946 * t1rho 01079 t1054 = 0.2e1_dp * t168 * t1rhorho 01080 t1055 = t1046 + t1048 + t1050 + t1052 + t1054 01081 t1060 = t78 * t86 01082 t1061 = t176 ** 2 01083 t1062 = 0.1e1_dp / t1061 01084 t1067 = -0.2e1_dp * t162 * t178 * t967 * t1rho + 0.2e1_dp * t175 *& 01085 t409 * t967 * t369 + 0.2e1_dp * t976 * t374 + 0.4e1_dp * t980 * & 01086 t408 * trho * t990 - t175 * t995 * t432 + t78 * t83 * t1031 * t91 - & 01087 t175 * t1035 * t177 * t369 - 0.2e1_dp * t162 * t995 * t370 - & 01088 0.2e1_dp * t162 * t397 * t1042 + 0.2e1_dp * t162 * t1055 * t91 * & 01089 trho - 0.6e1_dp * t1060 * t1062 * t186 * t990 01090 t1093 = t1046 + t1048 + t1050 + t1052 + t1054 + 0.2e1_dp * t951 * & 01091 A1rho + 0.8e1_dp * t417 * t1012 + 0.2e1_dp * t179 * A1rhorho + & 01092 0.8e1_dp * t417 * t1017 + 0.12e2_dp * t426 * t1rho * t2rho + & 01093 0.4e1_dp * t183 * t1rhorho 01094 t1097 = t1rho * t989 01095 t1103 = trho * t989 01096 t1104 = t178 * t1103 01097 t1106 = t994 * t91 01098 t1109 = t171 * t408 01099 t1115 = t976 * t362 + t162 * t361 * trho1rho - t162 * t178 * t432 & 01100 * t2rho + t175 * t994 * t408 * t410 - t162 * t178 * trho1rho * t369 & 01101 - t162 * t178 * trho * t1093 - t162 * t397 * t1097 + t175 * t409 * & 01102 t186 * t1093 - t354 * t1104 + t162 * t1106 * trhorho + t175 * t1109 & 01103 * t990 + t175 * t409 * t432 * t989 01104 t1118 = t1106 * trho 01105 t1121 = t360 * t408 01106 t1122 = t186 * t989 01107 t1126 = t393 * t177 01108 t1129 = t408 * t186 01109 t1148 = t186 * t2rho 01110 t1152 = 0.2e1_dp * t354 * t1118 + 0.2e1_dp * t175 * t1121 * t1122 & 01111 - t175 * t1126 * t989 + 0.4e1_dp * t980 * t1129 * t1042 + 0.2e1_dp *& 01112 t976 * t378 - t175 * t404 * t967 + 0.2e1_dp * t162 * t377 * & 01113 t1rhorho - 0.2e1_dp * t976 * t401 + 0.2e1_dp * t78 * t1rhorho * t164& 01114 - 0.2e1_dp * t162 * t404 * t1103 - 0.2e1_dp * t162 * t404 * t1148 01115 t1157 = t393 * t91 01116 t1167 = A2rho * t182 01117 t1187 = A1rho * t182 01118 t1196 = t87 * t 01119 t1203 = t1008 + 0.8e1_dp * t417 * trhorho * A2rho + 0.12e2_dp * & 01120 t426 * trhorho * t2rho + 0.8e1_dp * t1167 * t423 + 0.8e1_dp * t417 *& 01121 trho * A1rhorho + 0.8e1_dp * t417 * Arho * t1rhorho + 0.8e1_dp * & 01122 t1167 * t418 + 0.8e1_dp * t417 * Arho1rho * t1rho + 0.12e2_dp * t426& 01123 * trho1rho * t1rho + 0.12e2_dp * t426 * trho * t1rhorho + t998 + & 01124 0.8e1_dp * t1187 * t954 + 0.8e1_dp * t417 * trho1rho * A1rho + & 01125 0.8e1_dp * t417 * Arhorho * t2rho + 0.24e2_dp * t1196 * t427 * t2rho& 01126 + 0.2e1_dp * A1rhorho * t88 * Arho + t1014 01127 t1218 = t1016 + t1001 + 0.2e1_dp * t951 * Arhorho + t1026 + t1028 & 01128 + t1022 + t1011 + t1030 + 0.2e1_dp * t179 * Arhorhorho + 0.4e1_dp * & 01129 t183 * trhorhorho + 0.2e1_dp * t414 * Arho1rho + t1004 + t1006 + & 01130 t1019 + t1024 + 0.24e2_dp * t84 * t1018 + 0.24e2_dp * t84 * t1005 + & 01131 0.24e2_dp * t84 * t1013 01132 t1226 = t163 * trho1rho 01133 t1249 = -0.2e1_dp * t162 * t178 * t186 * t1rhorho + 0.2e1_dp * & 01134 t162 * t1157 * t2rho - t175 * t178 * (t1203 + t1218) + 0.2e1_dp * & 01135 t162 * t1035 * t91 * t1rho + 0.2e1_dp * t354 * t1226 - 0.2e1_dp * & 01136 t162 * t995 * t400 + 0.2e1_dp * t162 * t163 * trhorhorho - 0.2e1_dp & 01137 * t162 * t178 * trhorho * t989 - t175 * t1055 * t177 * t186 - & 01138 0.2e1_dp * t976 * t371 - t175 * t397 * t1093 + 0.4e1_dp * t980 * & 01139 t1129 * t1097 01140 t1262 = 0.2e1_dp * t162 * t163 * t2rho + t78 * t83 * t994 * t91 - & 01141 t175 * t178 * t989 01142 t1263 = t439 * t1262 01143 t1291 = 0.2e1_dp * t976 * t164 + 0.2e1_dp * t162 * t1118 - & 01144 0.2e1_dp * t162 * t1104 + 0.2e1_dp * t162 * t1226 + 0.2e1_dp * t162 & 01145 * t377 * t2rho + t78 * t83 * t1035 * t91 - t175 * t397 * t989 - & 01146 0.2e1_dp * t162 * t178 * t1148 - t175 * t995 * t186 + 0.2e1_dp * & 01147 t175 * t409 * t1122 - t175 * t178 * t967 01148 t1292 = gamma_var * t1291 01149 t1295 = 0.1e1_dp / t438 / t94 01150 t1329 = 0.2e1_dp * t976 * t440 + 0.2e1_dp * t162 * t1106 * t1rho -& 01151 0.2e1_dp * t162 * t178 * t1097 + 0.2e1_dp * t162 * t163 * t1rhorho & 01152 + 0.2e1_dp * t162 * t361 * t2rho + t78 * t83 * t1055 * t91 - t175 * & 01153 t404 * t989 - 0.2e1_dp * t162 * t178 * t1042 - t175 * t995 * t369 + & 01154 0.2e1_dp * t175 * t409 * t990 - t175 * t178 * t1093 01155 kfrhorhorho = k_frhorhorho 01156 kf2rho = k_f2rho 01157 ex_unifrho1rho = ex_unifrhorho 01158 ex_unif1rhorho = ex_unifrho1rho 01159 ex_unif2rho = -0.3e1_dp / 0.4e1_dp * t5 * kf2rho 01160 t1342 = t195 ** 2 01161 srho1rho = t457 * t198 * kf2rho + t462 / 0.2e1_dp - t465 + t197 * & 01162 t115 * kf2rho / 0.2e1_dp + t466 01163 s1rhorho = srho1rho 01164 s2rho = -t197 * t6 * kf2rho / 0.2e1_dp - t200 / 0.2e1_dp 01165 t1380 = t202 ** 2 01166 t1385 = 0.1e1_dp / t1380 * t469 * mu * t101 * s 01167 t1386 = kappa ** 2 01168 t1387 = 0.1e1_dp / t1386 01169 t1389 = s1rho * s2rho 01170 t1393 = t470 * s 01171 Fxrho1rho = -0.8e1_dp * t471 * t472 * s2rho + 0.2e1_dp * t204 * & 01172 s2rho * srho + 0.2e1_dp * t204 * s * srho1rho 01173 Fx1rhorho = -0.8e1_dp * t471 * s1rho * t103 * s2rho + 0.2e1_dp * & 01174 t204 * t1389 + 0.2e1_dp * t204 * s * s1rhorho 01175 Fx2rho = 0.2e1_dp * t204 * s * s2rho 01176 ex_ldarhorhorho = ex_unif1rhorho * Fx + ex_unif1rho * Fx2rho + & 01177 ex_unif2rho * Fx1rho + ex_unif * Fx1rhorho + ex_unifrho1rho * Fx + & 01178 ex_unifrho * Fx2rho + ex_unifrhorho * Fx - 0.3e1_dp / 0.4e1_dp * my_rho& 01179 * t5 * kfrhorhorho * Fx + t487 * Fx2rho + ex_unifrho * Fx1rho + my_rho& 01180 * ex_unifrho1rho * Fx1rho + t208 * Fx1rhorho + ex_unif2rho * Fxrho & 01181 + ex_unif * Fxrho1rho + ex_unif1rho * Fxrho + my_rho * ex_unif1rhorho *& 01182 Fxrho + t491 * Fxrho1rho + ex_unif * Fxrhorho + my_rho * ex_unif2rho *& 01183 Fxrhorho + t108 * (0.48e2_dp * t1385 * srho * t1387 * t1389 - & 01184 0.24e2_dp * t1393 * t472 * t1389 - 0.8e1_dp * t471 * srho1rho * t103& 01185 * s1rho - 0.8e1_dp * t471 * t472 * s1rhorho + 0.2e1_dp * t204 * & 01186 s1rhorho * srho + 0.2e1_dp * t204 * s1rho * srho1rho - 0.8e1_dp * & 01187 t471 * srhorho * t103 * s2rho + 0.2e1_dp * t204 * s2rho * srhorho + & 01188 0.2e1_dp * t204 * s * (-0.3e1_dp * my_norm_drho / t1342 * t459 * kf2rho& 01189 - t457 * t115 * t458 + 0.2e1_dp * t457 * t463 * kfrho - 0.2e1_dp * & 01190 t457 * t461 * kf2rho - 0.2e1_dp * t197 * t246 * kfrho + 0.3e1_dp / & 01191 0.2e1_dp * t197 * t115 * kfrhorho + t457 * t463 * kf2rho - t197 * t6& 01192 * kfrhorhorho / 0.2e1_dp - t197 * t246 * kf2rho - 0.3e1_dp * t99 * & 01193 t241)) 01194 01195 e_rho_rho_rho(ii) = e_rho_rho_rho(ii)+& 01196 scale_ex * ex_ldarhorhorho + scale_ec * (& 01197 e_c_u_01rhorho + gamma_var * t1329 * t191 - t451 * t1263 + & 01198 e_c_u_0rho1rho + t1292 * t191 - t190 * t1263 + epsilon_cGGArhorho + & 01199 my_rho * (e_c_u_0rhorhorho + gamma_var * (t1067 + 0.2e1_dp * t1115 + & 01200 t1152 + t1249) * t191 - t436 * t1263 - t1292 * t449 + 0.2e1_dp * & 01201 t190 * t1295 * t448 * t1262 - t190 * t439 * t1329)) 01202 01203 t1468 = t149 * t115 01204 tnorm_drhorhorho = t317 * t6 * t825 + t1468 * k_srho / 0.2e1_dp - & 01205 t496 * k_srhorho / 0.2e1_dp + t1468 * k_s1rho / 0.2e1_dp + t74 * & 01206 t246 01207 tnorm_drho1rho = -t496 * k_s1rho / 0.2e1_dp - t498 / 0.2e1_dp 01208 t1482 = A1rho * tnorm_drhorho 01209 t1492 = t78 * t84 01210 t1493 = tnorm_drho * t177 01211 t1504 = t408 * tnorm_drho 01212 t1505 = t1504 * t410 01213 t1511 = A1rho * tnorm_drho 01214 t1515 = t226 * t1rho 01215 t1521 = A * tnorm_drho1rho 01216 t1525 = 0.2e1_dp * t175 * t409 * t553 * t369 + 0.2e1_dp * t217 * & 01217 t1482 * t91 - 0.2e1_dp * t162 * t178 * tnorm_drhorho * t369 + & 01218 0.2e1_dp * t354 * t510 - 0.6e1_dp * t1492 * t1493 * t400 + 0.6e1_dp & 01219 * t175 * t218 * t91 * trhorho + 0.2e1_dp * t162 * t1157 * tnorm_drho& 01220 + 0.4e1_dp * t980 * t1505 - 0.6e1_dp * t1492 * t1493 * t370 + & 01221 0.6e1_dp * t175 * t1511 * t513 - 0.2e1_dp * t162 * t397 * t1515 - & 01222 0.2e1_dp * t354 * t507 + 0.6e1_dp * t175 * t1521 * t513 01223 t1528 = Arhorho * tnorm_drho 01224 t1532 = t177 * t369 01225 t1545 = t78 * t417 01226 t1565 = 0.2e1_dp * t385 * tnorm_drho + 0.2e1_dp * t388 * & 01227 tnorm_drho + 0.2e1_dp * t168 * tnorm_drho1rho + 0.8e1_dp * t417 * & 01228 t1511 + 0.12e2_dp * t426 * tnorm_drho * t1rho + 0.4e1_dp * t183 * & 01229 tnorm_drho1rho 01230 t1573 = t91 * t1rho 01231 t1584 = -0.2e1_dp * t354 * t530 + 0.2e1_dp * t217 * t1528 * t91 - & 01232 0.2e1_dp * t217 * t517 * t1532 - 0.2e1_dp * t217 * t1511 * t525 + & 01233 0.2e1_dp * t162 * t377 * tnorm_drho1rho + 0.2e1_dp * t162 * t361 * & 01234 tnorm_drhorho + 0.4e1_dp * t1545 * t1505 + 0.2e1_dp * t217 * A * & 01235 tnorm_drhorhorho * t91 + 0.2e1_dp * t175 * t409 * t1565 * t186 - & 01236 0.2e1_dp * t217 * t521 * t1532 + 0.6e1_dp * t175 * t521 * t1573 - & 01237 0.6e1_dp * t1060 * t1062 * t226 * t410 + 0.6e1_dp * t175 * t517 * & 01238 t1573 01239 t1608 = t408 * t226 01240 t1612 = t163 * tnorm_drho1rho 01241 t1628 = -0.2e1_dp * t162 * t404 * t506 + 0.2e1_dp * t354 * t503 - & 01242 0.2e1_dp * t162 * t178 * tnorm_drho * t432 - 0.2e1_dp * t162 * t178 & 01243 * t553 * t1rho - 0.2e1_dp * t162 * t404 * t529 - 0.2e1_dp * t162 * & 01244 t178 * t1565 * trho - t175 * t1126 * t226 + 0.4e1_dp * t980 * t1608 & 01245 * t370 + 0.2e1_dp * t500 * t1612 + 0.12e2_dp * t78 * t168 * & 01246 tnorm_drho * t91 * t427 + 0.2e1_dp * t162 * t163 * tnorm_drhorhorho & 01247 - t175 * t404 * t553 + 0.2e1_dp * t175 * t1121 * t535 01248 t1629 = Arho * tnorm_drho1rho 01249 t1633 = tnorm_drho * t369 01250 t1646 = t361 * tnorm_drho 01251 t1652 = t226 * t369 01252 t1660 = t178 * t1633 01253 t1672 = t418 * tnorm_drho 01254 t1676 = t423 * tnorm_drho 01255 t1715 = 0.2e1_dp * t999 * tnorm_drho + 0.2e1_dp * t1672 + 0.2e1_dp& 01256 * t382 * tnorm_drho1rho + 0.2e1_dp * t1676 + 0.2e1_dp * A * trhorho& 01257 * tnorm_drho + 0.2e1_dp * t541 * tnorm_drho1rho + 0.2e1_dp * t385 *& 01258 tnorm_drhorho + 0.2e1_dp * t388 * tnorm_drhorho + 0.2e1_dp * t168 *& 01259 tnorm_drhorhorho + 0.8e1_dp * t1187 * t517 + 0.24e2_dp * t84 * & 01260 t1672 + 0.8e1_dp * t417 * t1629 + 0.8e1_dp * t417 * t1528 + & 01261 0.24e2_dp * t84 * t1676 + 0.24e2_dp * t1196 * t548 * t1rho + & 01262 0.12e2_dp * t426 * tnorm_drho1rho * trho + 0.12e2_dp * t426 * & 01263 tnorm_drho * trhorho + 0.8e1_dp * t417 * t1482 + 0.12e2_dp * t426 * & 01264 tnorm_drhorho * t1rho + 0.4e1_dp * t183 * tnorm_drhorhorho 01265 t1722 = 0.2e1_dp * t217 * t1629 * t91 - 0.2e1_dp * t162 * t397 * & 01266 t1633 - t175 * t397 * t1565 - 0.2e1_dp * t162 * t178 * t226 * & 01267 trhorho + 0.2e1_dp * t78 * trhorho * t214 + 0.2e1_dp * t500 * t1646 & 01268 - 0.2e1_dp * t217 * t1521 * t525 + 0.2e1_dp * t175 * t1109 * t1652 -& 01269 0.2e1_dp * t162 * t178 * tnorm_drho1rho * t186 - 0.2e1_dp * t500 * & 01270 t1660 + 0.4e1_dp * t980 * t1608 * t400 + 0.2e1_dp * t175 * t409 * & 01271 t226 * t432 - t175 * t178 * t1715 - 0.2e1_dp * t217 * t218 * t177 * & 01272 t432 01273 t1758 = 0.2e1_dp * t354 * t214 + 0.2e1_dp * t162 * t1646 - & 01274 0.2e1_dp * t162 * t1660 + 0.2e1_dp * t162 * t1612 + 0.6e1_dp * t175 & 01275 * t218 * t1573 + 0.2e1_dp * t217 * t1511 * t91 + 0.2e1_dp * t217 * & 01276 t1521 * t91 - 0.2e1_dp * t217 * t218 * t1532 - 0.2e1_dp * t162 * & 01277 t178 * t1515 - t175 * t404 * t226 + 0.2e1_dp * t175 * t409 * t1652 -& 01278 t175 * t178 * t1565 01279 t1759 = gamma_var * t1758 01280 t1761 = t1295 * t189 01281 snorm_drho1rho = snorm_drhorho 01282 t1797 = snorm_drhorho * t103 01283 Fxnorm_drho1rho = -0.8e1_dp * t471 * t566 * s1rho + 0.2e1_dp * & 01284 t204 * s1rho * snorm_drho + 0.2e1_dp * t204 * s * snorm_drho1rho 01285 01286 e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii)+& 01287 scale_ex * (ex_unif1rho * Fxnorm_drho + & 01288 ex_unif * Fxnorm_drho1rho + ex_unifrho * Fxnorm_drho + t487 * & 01289 Fxnorm_drho + t208 * Fxnorm_drho1rho + ex_unif * Fxnorm_drhorho + & 01290 t491 * Fxnorm_drhorho + t108 * (0.48e2_dp * t1385 * snorm_drho * & 01291 t1387 * t476 - 0.24e2_dp * t1393 * t566 * t476 - 0.8e1_dp * t471 * & 01292 snorm_drho1rho * t103 * srho - 0.8e1_dp * t471 * t566 * srhorho + & 01293 0.2e1_dp * t204 * srhorho * snorm_drho + 0.2e1_dp * t204 * srho * & 01294 snorm_drho1rho - 0.8e1_dp * t471 * t1797 * s1rho + 0.2e1_dp * t204 *& 01295 s1rho * snorm_drhorho + 0.2e1_dp * t204 * s * (t456 * t6 * t458 + & 01296 t196 * t115 * kfrho - t562 * kfrhorho / 0.2e1_dp + t98 * t246))) + & 01297 scale_ec * (t1759 * t191 - t230 * t449 + Hnorm_drhorho + my_rho * (& 01298 gamma_var * (t1525 + t1584 + t1628 + t1722) * t191 - t557 * t449 - & 01299 t1759 * t559 + 0.2e1_dp * t230 * t1761 * t448 - t230 * t439 * t435)) 01300 01301 t1838 = t1504 * t535 01302 t1851 = Arho * t581 01303 t1878 = 0.4e1_dp * t175 * t409 * t226 * t553 - 0.2e1_dp * t162 * & 01304 t178 * t605 * trho - 0.4e1_dp * t217 * t218 * t177 * t553 + 0.8e1_dp& 01305 * t1545 * t1838 + 0.2e1_dp * t78 * t581 * t171 * t91 - 0.4e1_dp * & 01306 t162 * t397 * t590 - 0.10e2_dp * t175 * t586 * t525 - t175 * t178 * & 01307 (0.2e1_dp * t1851 + 0.4e1_dp * t521 * tnorm_drho + 0.24e2_dp * t84 *& 01308 t1851 + 0.24e2_dp * t1196 * t581 * trho + 0.24e2_dp * t426 * & 01309 tnorm_drhorho * tnorm_drho) - 0.12e2_dp * t1492 * t1493 * t529 + & 01310 0.8e1_dp * t980 * t1838 - 0.4e1_dp * t162 * t178 * tnorm_drhorho * & 01311 t226 + 0.20e2_dp * t162 * t586 * t513 01312 t1889 = t78 * t581 01313 t1907 = t85 * t1062 01314 t1922 = -0.4e1_dp * t500 * t591 + 0.2e1_dp * t175 * t409 * t605 * & 01315 t186 + 0.4e1_dp * t162 * t409 * t598 * trho - 0.2e1_dp * t1889 * & 01316 t187 + 0.2e1_dp * t175 * t1109 * t598 + 0.20e2_dp * t175 * t218 * & 01317 t91 * tnorm_drhorho + 0.10e2_dp * t175 * t1851 * t91 - 0.4e1_dp * & 01318 t217 * t517 * t594 - t175 * t397 * t605 - 0.6e1_dp * t175 * t1907 * & 01319 t598 * t186 - 0.4e1_dp * t162 * t178 * tnorm_drho * t553 + 0.4e1_dp & 01320 * t78 * tnorm_drhorho * t214 - 0.4e1_dp * t217 * t521 * t594 01321 t1927 = t439 * t229 01322 t1933 = t614 * t1387 01323 t1937 = t614 * t103 01324 01325 e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii)+& 01326 scale_ex * (ex_unif * & 01327 Fxnorm_drhonorm_drho + t208 * Fxnorm_drhonorm_drho + t108 * (& 01328 0.48e2_dp * t1385 * t1933 * srho - 0.24e2_dp * t1393 * t1937 * srho & 01329 - 0.16e2_dp * t471 * t1797 * snorm_drho + 0.4e1_dp * t204 * & 01330 snorm_drhorho * snorm_drho)) + scale_ec * (Hnorm_drhonorm_drho + my_rho& 01331 * (gamma_var * (t1878 + t1922) * t191 - t609 * t559 - 0.2e1_dp * & 01332 t557 * t1927 + 0.2e1_dp * t612 * t1761)) 01333 01334 t2norm_drho = tnorm_drho 01335 t1952 = t226 * t2norm_drho 01336 t1964 = 0.2e1_dp * t168 * t2norm_drho + 0.4e1_dp * t183 * t2norm_drho 01337 t1965 = t178 * t1964 01338 t1968 = t226 * t1964 01339 t1969 = t1504 * t1968 01340 t1972 = A * t2norm_drho 01341 t1990 = 0.2e1_dp * t1972 * tnorm_drho + 0.12e2_dp * t426 * & 01342 tnorm_drho * t2norm_drho 01343 t2020 = t177 * t1964 01344 t2024 = t91 * t2norm_drho 01345 t2028 = t78 * t2norm_drho 01346 t2031 = -0.20e2_dp * t1492 * t1493 * t1952 + 0.4e1_dp * t162 * & 01347 t409 * t598 * t2norm_drho - 0.2e1_dp * t1889 * t1965 + 0.8e1_dp * & 01348 t1545 * t1969 + 0.4e1_dp * t217 * t1972 * t408 * t598 - 0.6e1_dp * & 01349 t175 * t1907 * t598 * t1964 - 0.2e1_dp * t217 * t1972 * t177 * t605 & 01350 - 0.4e1_dp * t217 * t218 * t177 * t1990 + 0.4e1_dp * t175 * t409 * & 01351 t1990 * t226 - 0.24e2_dp * t78 * t182 * t85 * t177 * t87 * t581 * & 01352 t2norm_drho + 0.2e1_dp * t175 * t409 * t605 * t1964 + 0.8e1_dp * & 01353 t980 * t1969 - 0.2e1_dp * t162 * t178 * t605 * t2norm_drho - & 01354 0.4e1_dp * t162 * t178 * t1990 * tnorm_drho - 0.10e2_dp * t175 * & 01355 t586 * t2020 + 0.24e2_dp * t162 * t586 * t2024 - 0.4e1_dp * t2028 * & 01356 t591 01357 t2041 = 0.2e1_dp * t162 * t163 * t2norm_drho + 0.2e1_dp * t217 * & 01358 t1972 * t91 - t175 * t1965 01359 s2norm_drho = snorm_drho 01360 01361 e_ndrho_ndrho_ndrho(ii) = e_ndrho_ndrho_ndrho(ii)+& 01362 scale_ex * t108 * (0.48e2_dp *& 01363 t1385 * t1933 * s2norm_drho - 0.24e2_dp * t1393 * t1937 * & 01364 s2norm_drho) + scale_ec * my_rho * (gamma_var * t2031 * t191 - t609 * & 01365 t439 * t2041 - 0.2e1_dp * gamma_var * (0.2e1_dp * t2028 * t214 + & 01366 0.10e2_dp * t175 * t218 * t2024 - 0.2e1_dp * t162 * t178 * & 01367 tnorm_drho * t1964 - 0.2e1_dp * t217 * t218 * t2020 - 0.2e1_dp * & 01368 t162 * t178 * t1952 - 0.2e1_dp * t217 * t1972 * t594 + 0.2e1_dp * & 01369 t175 * t409 * t1968 - t175 * t178 * t1990) * t1927 + 0.2e1_dp * t612& 01370 * t1295 * t2041) 01371 END IF 01372 END IF 01373 END DO 01374 !$omp end do 01375 END SELECT 01376 01377 END SUBROUTINE pbe_lda_calc 01378 01379 ! ***************************************************************************** 01385 SUBROUTINE pbe_lsd_eval(rho_set,deriv_set,grad_deriv,pbe_params,error) 01386 TYPE(xc_rho_set_type), POINTER :: rho_set 01387 TYPE(xc_derivative_set_type), POINTER :: deriv_set 01388 INTEGER, INTENT(in) :: grad_deriv 01389 TYPE(section_vals_type), POINTER :: pbe_params 01390 TYPE(cp_error_type), INTENT(inout) :: error 01391 01392 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lsd_eval', 01393 routineP = moduleN//':'//routineN 01394 01395 INTEGER :: handle, npoints, param, stat 01396 INTEGER, DIMENSION(:, :), POINTER :: bo 01397 LOGICAL :: failure 01398 REAL(kind=dp) :: epsilon_drho, epsilon_rho, 01399 scale_ec, scale_ex 01400 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, 01401 e_ndr_ndr, e_ndr_ra, e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, e_ndrb, 01402 e_ndrb_ndrb, e_ndrb_rb, e_ra, e_ra_ra, e_ra_rb, e_rb, e_rb_rb, 01403 norm_drho, norm_drhoa, norm_drhob, rhoa, rhob 01404 TYPE(xc_derivative_type), POINTER :: deriv 01405 01406 CALL timeset(routineN,handle) 01407 failure=.FALSE. 01408 NULLIFY(deriv, bo) 01409 01410 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 01411 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 01412 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 01413 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 01414 IF (.NOT. failure) THEN 01415 CALL xc_rho_set_get(rho_set,& 01416 rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, & 01417 norm_drhob=norm_drhob, norm_drho=norm_drho, & 01418 rho_cutoff=epsilon_rho,& 01419 drho_cutoff=epsilon_drho, local_bounds=bo, error=error) 01420 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 01421 01422 ! meaningful default for the arrays we don't need: let us make compiler 01423 ! and debugger happy... 01424 IF (cp_debug) THEN 01425 ALLOCATE(dummy(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 01426 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01427 ELSE 01428 dummy=> rhoa 01429 END IF 01430 e_0 => dummy 01431 e_ra => dummy 01432 e_rb => dummy 01433 e_ndra_ra => dummy 01434 e_ndrb_rb => dummy 01435 e_ndr_ndr => dummy 01436 e_ndra_ndra => dummy 01437 e_ndrb_ndrb => dummy 01438 e_ndr => dummy 01439 e_ndra => dummy 01440 e_ndrb => dummy 01441 e_ra_ra => dummy 01442 e_ra_rb => dummy 01443 e_rb_rb => dummy 01444 e_ndr_ra => dummy 01445 e_ndr_rb => dummy 01446 01447 IF (grad_deriv>=0) THEN 01448 deriv => xc_dset_get_derivative(deriv_set,"",& 01449 allocate_deriv=.TRUE., error=error) 01450 CALL xc_derivative_get(deriv, deriv_data=e_0,error=error) 01451 END IF 01452 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 01453 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)",& 01454 allocate_deriv=.TRUE.,error=error) 01455 CALL xc_derivative_get(deriv,deriv_data=e_ra,error=error) 01456 deriv => xc_dset_get_derivative(deriv_set,"(rhob)",& 01457 allocate_deriv=.TRUE.,error=error) 01458 CALL xc_derivative_get(deriv,deriv_data=e_rb,error=error) 01459 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",& 01460 allocate_deriv=.TRUE.,error=error) 01461 CALL xc_derivative_get(deriv,deriv_data=e_ndr,error=error) 01462 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)",& 01463 allocate_deriv=.TRUE.,error=error) 01464 CALL xc_derivative_get(deriv,deriv_data=e_ndra,error=error) 01465 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)",& 01466 allocate_deriv=.TRUE.,error=error) 01467 CALL xc_derivative_get(deriv,deriv_data=e_ndrb,error=error) 01468 END IF 01469 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 01470 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)(rhoa)",& 01471 allocate_deriv=.TRUE.,error=error) 01472 CALL xc_derivative_get(deriv,deriv_data=e_ra_ra,error=error) 01473 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)(rhob)",& 01474 allocate_deriv=.TRUE.,error=error) 01475 CALL xc_derivative_get(deriv,deriv_data=e_ra_rb,error=error) 01476 deriv => xc_dset_get_derivative(deriv_set,"(rhob)(rhob)",& 01477 allocate_deriv=.TRUE.,error=error) 01478 CALL xc_derivative_get(deriv,deriv_data=e_rb_rb,error=error) 01479 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(rhoa)",& 01480 allocate_deriv=.TRUE.,error=error) 01481 CALL xc_derivative_get(deriv,deriv_data=e_ndr_ra,error=error) 01482 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(rhob)",& 01483 allocate_deriv=.TRUE.,error=error) 01484 CALL xc_derivative_get(deriv,deriv_data=e_ndr_rb,error=error) 01485 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(rhoa)",& 01486 allocate_deriv=.TRUE.,error=error) 01487 CALL xc_derivative_get(deriv,deriv_data=e_ndra_ra,error=error) 01488 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)(rhob)",& 01489 allocate_deriv=.TRUE.,error=error) 01490 CALL xc_derivative_get(deriv,deriv_data=e_ndrb_rb,error=error) 01491 deriv => xc_dset_get_derivative(deriv_set,& 01492 "(norm_drho)(norm_drho)", allocate_deriv=.TRUE.,error=error) 01493 CALL xc_derivative_get(deriv,deriv_data=e_ndr_ndr,error=error) 01494 deriv => xc_dset_get_derivative(deriv_set,& 01495 "(norm_drhoa)(norm_drhoa)", allocate_deriv=.TRUE.,error=error) 01496 CALL xc_derivative_get(deriv,deriv_data=e_ndra_ndra,error=error) 01497 deriv => xc_dset_get_derivative(deriv_set,& 01498 "(norm_drhob)(norm_drhob)", allocate_deriv=.TRUE.,error=error) 01499 CALL xc_derivative_get(deriv,deriv_data=e_ndrb_ndrb,error=error) 01500 END IF 01501 IF (grad_deriv>=3.OR.grad_deriv==-3) THEN 01502 ! to do 01503 END IF 01504 01505 CALL section_vals_val_get(pbe_params,"scale_c",r_val=scale_ec,error=error) 01506 CALL section_vals_val_get(pbe_params,"scale_x",r_val=scale_ex,error=error) 01507 CALL section_vals_val_get(pbe_params,"parametrization",i_val=param,error=error) 01508 01509 !$omp parallel default(none),& 01510 !$omp shared(rhoa,rhob,norm_drho,norm_drhoa,norm_drhob,e_0,e_ra,e_rb,e_ndra_ra),& 01511 !$omp shared(e_ndrb_rb,e_ndr_ndr,e_ndra_ndra,e_ndrb_ndrb,e_ndr,e_ndra,e_ndrb),& 01512 !$omp shared(e_ra_ra,e_ra_rb,e_rb_rb,e_ndr_ra,e_ndr_rb,grad_deriv,npoints),& 01513 !$omp shared(epsilon_rho,epsilon_drho,error,param,scale_ec,scale_ex) 01514 CALL pbe_lsd_calc(& 01515 rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa,& 01516 norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb,& 01517 e_ra_ndra=e_ndra_ra,& 01518 e_rb_ndrb=e_ndrb_rb, e_ndr_ndr=e_ndr_ndr,& 01519 e_ndra_ndra=e_ndra_ndra, e_ndrb_ndrb=e_ndrb_ndrb, e_ndr=e_ndr,& 01520 e_ndra=e_ndra, e_ndrb=e_ndrb, e_ra_ra=e_ra_ra, & 01521 e_ra_rb=e_ra_rb, e_rb_rb=e_rb_rb, e_ra_ndr=e_ndr_ra,& 01522 e_rb_ndr=e_ndr_rb,& 01523 grad_deriv=grad_deriv, npoints=npoints, & 01524 epsilon_rho=epsilon_rho,epsilon_drho=epsilon_drho,& 01525 param=param,scale_ec=scale_ec,scale_ex=scale_ex,error=error) 01526 !$omp end parallel 01527 01528 IF (cp_debug) THEN 01529 DEALLOCATE(dummy,stat=stat) 01530 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01531 ELSE 01532 NULLIFY(dummy) 01533 END IF 01534 END IF 01535 CALL timestop(handle) 01536 END SUBROUTINE pbe_lsd_eval 01537 01538 ! ***************************************************************************** 01550 SUBROUTINE pbe_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob,& 01551 e_0, e_ra, e_rb, e_ra_ndra, e_rb_ndrb, e_ndr_ndr,& 01552 e_ndra_ndra, e_ndrb_ndrb, e_ndr,& 01553 e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, e_ra_ndr, e_rb_ndr,& 01554 grad_deriv,npoints,epsilon_rho,epsilon_drho,param,scale_ec,scale_ex,error) 01555 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob, norm_drho, 01556 norm_drhoa, norm_drhob 01557 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_ra, e_rb, e_ra_ndra, 01558 e_rb_ndrb, e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb, e_ndr, e_ndra, e_ndrb, 01559 e_ra_ra, e_ra_rb, e_rb_rb, e_ra_ndr, e_rb_ndr 01560 INTEGER, INTENT(in) :: grad_deriv, npoints 01561 REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_drho 01562 INTEGER, INTENT(in) :: param 01563 REAL(kind=dp), INTENT(in) :: scale_ec, scale_ex 01564 TYPE(cp_error_type), INTENT(inout) :: error 01565 01566 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lsd_calc', 01567 routineP = moduleN//':'//routineN 01568 REAL(kind=dp), PARAMETER :: small = 1.0e-20, t_1_9 = 1.0_dp/9.0_dp, 01569 t_2_3 = 2.0_dp/3.0_dp, t_4_27 = 4.0_dp/27.0_dp, t_4_3 = 4.0_dp/3.0_dp, 01570 t_4_9 = 4.0_dp/9.0_dp, t_8_27 = 8.0_dp/27.0_dp 01571 01572 INTEGER :: ii 01573 LOGICAL :: failure 01574 REAL(kind=dp) :: A, A1rhoa, A1rhob, A_1, A_2, A_3, alpha_1_1, alpha_1_2, 01575 alpha_1_3, alpha_c, alpha_c1rhoa, alpha_c1rhob, alpha_crhoa, 01576 alpha_crhob, Arhoa, Arhoarhoa, Arhoarhob, Arhob, Arhobrhob, beta, 01577 beta_1_1, beta_1_2, beta_1_3, beta_2_1, beta_2_2, beta_2_3, beta_3_1, 01578 beta_3_2, beta_3_3, beta_4_1, beta_4_2, beta_4_3, chi, chirhoa, 01579 chirhoarhoa, chirhoarhob, chirhob, chirhobrhob, e_c_u_0, e_c_u_01rhoa, 01580 e_c_u_01rhob, e_c_u_0rhoa, e_c_u_0rhoarhoa, e_c_u_0rhoarhob, 01581 e_c_u_0rhob, e_c_u_0rhobrhob, e_c_u_1rhoa, e_c_u_1rhob, epsilon_c_unif, 01582 epsilon_c_unif1rhoa, epsilon_c_unif1rhob, epsilon_c_unifrhoa, 01583 epsilon_c_unifrhoarhoa 01584 REAL(kind=dp) :: epsilon_c_unifrhoarhob, epsilon_c_unifrhob, 01585 epsilon_c_unifrhobrhob, epsilon_cGGA, epsilon_cGGArhoa, 01586 epsilon_cGGArhob, ex_unif_a, ex_unif_a1rhoa, ex_unif_arhoa, ex_unif_b, 01587 ex_unif_b1rhob, ex_unif_brhob, f, f1rhoa, f1rhob, frhoa, frhoarhoa, 01588 frhoarhob, frhob, frhobrhob, Fx_a, Fx_a1rhoa, Fx_anorm_drhoa, Fx_arhoa, 01589 Fx_b, Fx_b1rhob, Fx_bnorm_drhob, Fx_brhob, gamma_var, Hnorm_drho, 01590 k_frhoa, k_frhoarhoa, k_frhoarhob, k_frhob, k_s, k_s1rhoa, k_s1rhob, 01591 k_srhoa, k_srhob, kappa, kf_a, kf_arhoa, kf_arhoarhoa, kf_b, kf_brhob, 01592 kf_brhobrhob, mu, my_norm_drho, my_norm_drhoa, my_norm_drhob, my_rho, 01593 my_rhoa, my_rhob 01594 REAL(kind=dp) :: p_1, p_2, p_3, phi, phi1rhoa, phi1rhob, phirhoa, 01595 phirhoarhoa, phirhoarhob, phirhob, phirhobrhob, rs, rsrhoa, rsrhoarhoa, 01596 rsrhoarhob, rsrhob, rsrhobrhob, s_a, s_a1rhoa, s_anorm_drhoa, s_arhoa, 01597 s_b, s_b1rhob, s_bnorm_drhob, s_brhob, t, t1, t100, t1000, t1001, t101, 01598 t102, t103, t1031, t1033, t1038, t104, t105, t1050, t1069, t107, t1071, 01599 t108, t110, t1104, t1106, t111, t112, t113, t115, t116, t1164, t118, 01600 t119, t1193, t12, t122, t1228, t123, t1231, t124, t125, t126, t1269, 01601 t1283, t1285, t1286, t1288, t129, t1291, t1293, t130, t1304, t131, 01602 t1327, t133, t1330, t1342, t1346, t135, t137, t1377, t14, t140 01603 REAL(kind=dp) :: t1411, t142, t143, t1440, t146, t1462, t1465, t147, 01604 t148, t1482, t1489, t1492, t15, t150, t1501, t1514, t1520, t1529, t153, 01605 t1550, t1552, t1555, t156, t1588, t1590, t1591, t1599, t1602, t1603, 01606 t162, t163, t1630, t1632, t1635, t1638, t164, t1640, t165, t167, t1674, 01607 t1677, t1680, t1698, t171, t1711, t1712, t1713, t1739, t1741, t1743, 01608 t1748, t175, t176, t1765, t1766, t1767, t177, t178, t179, t18, t1801, 01609 t1829, t183, t1830, t1831, t1865, t187, t1871, t1876, t1888, t190, 01610 t1901, t191, t192, t1922, t194, t1949, t198, t199, t1rhoa, t1rhob, t2, 01611 t20, t200, t201, t205, t21, t211, t212, t213, t215, t219, t22 01612 REAL(kind=dp) :: t220, t221, t222, t226, t23, t232, t233, t234, t240, 01613 t242, t244, t245, t246, t248, t249, t250, t252, t254, t256, t257, t259, 01614 t262, t266, t268, t269, t27, t272, t273, t274, t277, t278, t28, t280, 01615 t281, t282, t284, t285, t286, t289, t293, t294, t297, t298, t299, t3, 01616 t302, t303, t305, t306, t310, t311, t312, t313, t314, t317, t318, t32, 01617 t321, t324, t325, t326, t329, t335, t336, t337, t34, t340, t341, t344, 01618 t346, t350, t368, t38, t382, t39, t396, t4, t40, t403, t405, t407, 01619 t409, t41, t410, t411, t413, t415, t417, t427, t429, t432, t436, t439, 01620 t442, t444, t445, t45, t453, t456, t457, t46, t460, t466 01621 REAL(kind=dp) :: t467, t468, t471, t472, t475, t477, t481, t488, t493, 01622 t494, t5, t50, t502, t505, t506, t518, t519, t52, t523, t525, t536, 01623 t538, t543, t544, t548, t549, t550, t556, t56, t561, t563, t564, t57, 01624 t576, t578, t579, t58, t580, t588, t59, t590, t595, t596, t6, t600, 01625 t606, t611, t624, t626, t627, t628, t63, t636, t638, t64, t643, t644, 01626 t648, t654, t659, t66, t672, t674, t675, t676, t681, t682, t687, t69, 01627 t7, t70, t705, t708, t71, t711, t72, t726, t73, t733, t736, t74, t745, 01628 t75, t750, t755, t763, t767, t768, t77, t775, t776, t779, t78, t785, 01629 t789, t79, t795, t798, t8, t80, t801, t812, t82, t820, t821 01630 REAL(kind=dp) :: t822, t823, t825, t828, t839, t84, t840, t85, t851, 01631 t858, t865, t867, t868, t87, t876, t879, t88, t880, t9, t90, t904, 01632 t908, t909, t91, t911, t914, t917, t919, t92, t924, t93, t936, t94, 01633 t944, t95, t953, t959, t96, t962, t965, t966, t967, t97, t98, t985, 01634 t998, t999, tnorm_drho, trhoa, trhoanorm_drho, trhoarhoa, trhoarhob, 01635 trhob, trhobnorm_drho, trhobrhob 01636 01637 failure=.FALSE. 01638 01639 ! Parametrization of PBE 01640 SELECT CASE(param) 01641 ! Original PBE 01642 CASE (xc_pbe_orig) 01643 kappa = 0.804e0_dp 01644 beta = 0.66725e-1_dp 01645 mu = beta * pi ** 2 / 0.3e1_dp 01646 ! Revised PBE (revPBE) 01647 CASE (xc_pbe_rev) 01648 kappa = 0.1245e1_dp 01649 beta = 0.66725e-1_dp 01650 mu = beta * pi ** 2 / 0.3e1_dp 01651 ! PBE for solids (PBEsol) 01652 CASE (xc_pbe_sol) 01653 kappa = 0.804e0_dp 01654 beta = 0.46e-1_dp 01655 mu = 0.1e2_dp / 0.81e2_dp 01656 CASE default 01657 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01658 END SELECT 01659 01660 gamma_var = (0.1e1_dp - LOG(0.2e1_dp)) / pi ** 2 01661 p_1 = 0.10e1_dp 01662 A_1 = 0.31091e-1_dp 01663 alpha_1_1 = 0.21370e0_dp 01664 beta_1_1 = 0.75957e1_dp 01665 beta_2_1 = 0.35876e1_dp 01666 beta_3_1 = 0.16382e1_dp 01667 beta_4_1 = 0.49294e0_dp 01668 p_2 = 0.10e1_dp 01669 A_2 = 0.15545e-1_dp 01670 alpha_1_2 = 0.20548e0_dp 01671 beta_1_2 = 0.141189e2_dp 01672 beta_2_2 = 0.61977e1_dp 01673 beta_3_2 = 0.33662e1_dp 01674 beta_4_2 = 0.62517e0_dp 01675 p_3 = 0.10e1_dp 01676 A_3 = 0.16887e-1_dp 01677 alpha_1_3 = 0.11125e0_dp 01678 beta_1_3 = 0.10357e2_dp 01679 beta_2_3 = 0.36231e1_dp 01680 beta_3_3 = 0.88026e0_dp 01681 beta_4_3 = 0.49671e0_dp 01682 t3 = 3 ** (0.1e1_dp / 0.3e1_dp) 01683 t4 = 4 ** (0.1e1_dp / 0.3e1_dp) 01684 t5 = t4 ** 2 01685 t6 = t3 * t5 01686 t7 = 0.1e1_dp / pi 01687 t90 = pi ** 2 01688 01689 SELECT CASE(grad_deriv) 01690 CASE default 01691 !$omp do 01692 DO ii=1,npoints 01693 my_rhoa=MAX(rhoa(ii),0.0_dp) 01694 my_rhob=MAX(rhob(ii),0.0_dp) 01695 my_rho=my_rhoa+my_rhob 01696 IF (my_rho>epsilon_rho) THEN 01697 my_rhoa=MAX(EPSILON(0.0_dp)*1.e4_dp,my_rhoa) 01698 my_rhob=MAX(EPSILON(0.0_dp)*1.e4_dp,my_rhob) 01699 my_rho = my_rhoa + my_rhob 01700 my_norm_drho=norm_drho(ii) 01701 my_norm_drhoa=norm_drhoa(ii) 01702 my_norm_drhob=norm_drhob(ii) 01703 01704 t1 = my_rhoa - my_rhob 01705 t2 = 0.1e1_dp / my_rho 01706 chi = t1 * t2 01707 t8 = t7 * t2 01708 t9 = t8 ** (0.1e1_dp / 0.3e1_dp) 01709 rs = t6 * t9 / 0.4e1_dp 01710 t12 = 0.1e1_dp + alpha_1_1 * rs 01711 t14 = 0.1e1_dp / A_1 01712 t15 = SQRT(rs) 01713 t18 = t15 * rs 01714 t20 = p_1 + 0.1e1_dp 01715 t21 = rs ** t20 01716 t22 = beta_4_1 * t21 01717 t23 = beta_1_1 * t15 + beta_2_1 * rs + beta_3_1 * t18 + t22 01718 t27 = 0.1e1_dp + t14 / t23 / 0.2e1_dp 01719 t28 = LOG(t27) 01720 e_c_u_0 = -0.2e1_dp * A_1 * t12 * t28 01721 t32 = 0.1e1_dp + alpha_1_2 * rs 01722 t34 = 0.1e1_dp / A_2 01723 t38 = p_2 + 0.1e1_dp 01724 t39 = rs ** t38 01725 t40 = beta_4_2 * t39 01726 t41 = beta_1_2 * t15 + beta_2_2 * rs + beta_3_2 * t18 + t40 01727 t45 = 0.1e1_dp + t34 / t41 / 0.2e1_dp 01728 t46 = LOG(t45) 01729 t50 = 0.1e1_dp + alpha_1_3 * rs 01730 t52 = 0.1e1_dp / A_3 01731 t56 = p_3 + 0.1e1_dp 01732 t57 = rs ** t56 01733 t58 = beta_4_3 * t57 01734 t59 = beta_1_3 * t15 + beta_2_3 * rs + beta_3_3 * t18 + t58 01735 t63 = 0.1e1_dp + t52 / t59 / 0.2e1_dp 01736 t64 = LOG(t63) 01737 alpha_c = 0.2e1_dp * A_3 * t50 * t64 01738 t66 = 2 ** (0.1e1_dp / 0.3e1_dp) 01739 t69 = 1 / (2 * t66 - 2) 01740 t70 = 0.1e1_dp + chi 01741 t71 = t70 ** (0.1e1_dp / 0.3e1_dp) 01742 t72 = t71 * t70 01743 t73 = 0.1e1_dp - chi 01744 t74 = t73 ** (0.1e1_dp / 0.3e1_dp) 01745 t75 = t74 * t73 01746 f = (t72 + t75 - 0.2e1_dp) * t69 01747 t77 = alpha_c * f 01748 t78 = 0.9e1_dp / 0.8e1_dp / t69 01749 t79 = chi ** 2 01750 t80 = t79 ** 2 01751 t82 = t78 * (0.1e1_dp - t80) 01752 t84 = -0.2e1_dp * A_2 * t32 * t46 - e_c_u_0 01753 t85 = t84 * f 01754 epsilon_c_unif = e_c_u_0 + t77 * t82 + t85 * t80 01755 t87 = t71 ** 2 01756 t88 = t74 ** 2 01757 phi = t87 / 0.2e1_dp + t88 / 0.2e1_dp 01758 t91 = t90 * my_rho 01759 t92 = t91 ** (0.1e1_dp / 0.3e1_dp) 01760 t93 = t3 * t92 * t7 01761 t94 = SQRT(t93) 01762 k_s = 0.2e1_dp * t94 01763 t95 = 0.1e1_dp / phi 01764 t96 = my_norm_drho * t95 01765 t97 = 0.1e1_dp / k_s 01766 t98 = t97 * t2 01767 t = t96 * t98 / 0.2e1_dp 01768 t100 = 0.1e1_dp / gamma_var 01769 t101 = beta * t100 01770 t102 = epsilon_c_unif * t100 01771 t103 = phi ** 2 01772 t104 = t103 * phi 01773 t105 = 0.1e1_dp / t104 01774 t107 = EXP(-t102 * t105) 01775 t108 = t107 - 0.1e1_dp 01776 A = t101 / t108 01777 t110 = gamma_var * t104 01778 t111 = t ** 2 01779 t112 = A * t111 01780 t113 = 0.1e1_dp + t112 01781 t115 = A ** 2 01782 t116 = t111 ** 2 01783 t118 = 0.1e1_dp + t112 + t115 * t116 01784 t119 = 0.1e1_dp / t118 01785 t122 = 0.1e1_dp + t101 * t111 * t113 * t119 01786 t123 = LOG(t122) 01787 epsilon_cGGA = epsilon_c_unif + t110 * t123 01788 t124 = t3 * t66 01789 t125 = t90 * my_rhoa 01790 t126 = t125 ** (0.1e1_dp / 0.3e1_dp) 01791 kf_a = t124 * t126 01792 ex_unif_a = -0.3e1_dp / 0.4e1_dp * t7 * kf_a 01793 t129 = 0.1e1_dp / kf_a 01794 t130 = my_norm_drhoa * t129 01795 t131 = 0.1e1_dp / my_rhoa 01796 s_a = t130 * t131 / 0.2e1_dp 01797 t133 = s_a ** 2 01798 t135 = 0.1e1_dp / kappa 01799 t137 = 0.1e1_dp + mu * t133 * t135 01800 Fx_a = 0.1e1_dp + kappa - kappa / t137 01801 t140 = my_rhoa * ex_unif_a 01802 t142 = t90 * my_rhob 01803 t143 = t142 ** (0.1e1_dp / 0.3e1_dp) 01804 kf_b = t124 * t143 01805 ex_unif_b = -0.3e1_dp / 0.4e1_dp * t7 * kf_b 01806 t146 = 0.1e1_dp / kf_b 01807 t147 = my_norm_drhob * t146 01808 t148 = 0.1e1_dp / my_rhob 01809 s_b = t147 * t148 / 0.2e1_dp 01810 t150 = s_b ** 2 01811 t153 = 0.1e1_dp + mu * t150 * t135 01812 Fx_b = 0.1e1_dp + kappa - kappa / t153 01813 t156 = my_rhob * ex_unif_b 01814 01815 IF (grad_deriv>=0) THEN 01816 e_0(ii) = e_0(ii)+& 01817 scale_ex * (0.2e1_dp * t140 * Fx_a + 0.2e1_dp * t156 * Fx_b)& 01818 / 0.2e1_dp + scale_ec * my_rho * epsilon_cGGA 01819 END IF 01820 01821 t162 = my_rho ** 2 01822 t163 = 0.1e1_dp / t162 01823 t164 = t1 * t163 01824 chirhoa = t2 - t164 01825 t165 = t9 ** 2 01826 t167 = 0.1e1_dp / t165 * t7 01827 rsrhoa = -t6 * t167 * t163 / 0.12e2_dp 01828 t171 = A_1 * alpha_1_1 01829 t175 = t23 ** 2 01830 t176 = 0.1e1_dp / t175 01831 t177 = t12 * t176 01832 t178 = 0.1e1_dp / t15 01833 t179 = beta_1_1 * t178 01834 t183 = beta_3_1 * t15 01835 t187 = 0.1e1_dp / rs 01836 t190 = t179 * rsrhoa / 0.2e1_dp + beta_2_1 * rsrhoa + 0.3e1_dp / & 01837 0.2e1_dp * t183 * rsrhoa + t22 * t20 * rsrhoa * t187 01838 t191 = 0.1e1_dp / t27 01839 t192 = t190 * t191 01840 e_c_u_0rhoa = -0.2e1_dp * t171 * rsrhoa * t28 + t177 * t192 01841 t194 = A_2 * alpha_1_2 01842 t198 = t41 ** 2 01843 t199 = 0.1e1_dp / t198 01844 t200 = t32 * t199 01845 t201 = beta_1_2 * t178 01846 t205 = beta_3_2 * t15 01847 t211 = t201 * rsrhoa / 0.2e1_dp + beta_2_2 * rsrhoa + 0.3e1_dp / & 01848 0.2e1_dp * t205 * rsrhoa + t40 * t38 * rsrhoa * t187 01849 t212 = 0.1e1_dp / t45 01850 t213 = t211 * t212 01851 e_c_u_1rhoa = -0.2e1_dp * t194 * rsrhoa * t46 + t200 * t213 01852 t215 = A_3 * alpha_1_3 01853 t219 = t59 ** 2 01854 t220 = 0.1e1_dp / t219 01855 t221 = t50 * t220 01856 t222 = beta_1_3 * t178 01857 t226 = beta_3_3 * t15 01858 t232 = t222 * rsrhoa / 0.2e1_dp + beta_2_3 * rsrhoa + 0.3e1_dp / & 01859 0.2e1_dp * t226 * rsrhoa + t58 * t56 * rsrhoa * t187 01860 t233 = 0.1e1_dp / t63 01861 t234 = t232 * t233 01862 alpha_crhoa = 0.2e1_dp * t215 * rsrhoa * t64 - t221 * t234 01863 frhoa = (0.4e1_dp / 0.3e1_dp * t71 * chirhoa - 0.4e1_dp / 0.3e1_dp& 01864 * t74 * chirhoa) * t69 01865 t240 = alpha_crhoa * f 01866 t242 = alpha_c * frhoa 01867 t244 = t79 * chi 01868 t245 = t78 * t244 01869 t246 = t245 * chirhoa 01870 t248 = 0.4e1_dp * t77 * t246 01871 t249 = e_c_u_1rhoa - e_c_u_0rhoa 01872 t250 = t249 * f 01873 t252 = t84 * frhoa 01874 t254 = t244 * chirhoa 01875 t256 = 0.4e1_dp * t85 * t254 01876 epsilon_c_unifrhoa = e_c_u_0rhoa + t240 * t82 + t242 * t82 - t248 & 01877 + t250 * t80 + t252 * t80 + t256 01878 t257 = 0.1e1_dp / t71 01879 t259 = 0.1e1_dp / t74 01880 phirhoa = t257 * chirhoa / 0.3e1_dp - t259 * chirhoa / 0.3e1_dp 01881 t262 = t92 ** 2 01882 k_frhoa = t3 / t262 * t90 / 0.3e1_dp 01883 t266 = 0.1e1_dp / t94 01884 k_srhoa = t266 * k_frhoa * t7 01885 t268 = 0.1e1_dp / t103 01886 t269 = my_norm_drho * t268 01887 t272 = k_s ** 2 01888 t273 = 0.1e1_dp / t272 01889 t274 = t273 * t2 01890 t277 = t97 * t163 01891 t278 = t96 * t277 01892 trhoa = -t269 * t98 * phirhoa / 0.2e1_dp - t96 * t274 * k_srhoa / & 01893 0.2e1_dp - t278 / 0.2e1_dp 01894 t280 = t108 ** 2 01895 t281 = 0.1e1_dp / t280 01896 t282 = epsilon_c_unifrhoa * t100 01897 t284 = t103 ** 2 01898 t285 = 0.1e1_dp / t284 01899 t286 = t285 * phirhoa 01900 t289 = -t282 * t105 + 0.3e1_dp * t102 * t286 01901 Arhoa = -t101 * t281 * t289 * t107 01902 t293 = gamma_var * t103 01903 t294 = t123 * phirhoa 01904 t297 = t101 * t 01905 t298 = t113 * t119 01906 t299 = t298 * trhoa 01907 t302 = Arhoa * t111 01908 t303 = A * t 01909 t305 = 0.2e1_dp * t303 * trhoa 01910 t306 = t302 + t305 01911 t310 = t101 * t111 01912 t311 = t118 ** 2 01913 t312 = 0.1e1_dp / t311 01914 t313 = t113 * t312 01915 t314 = A * t116 01916 t317 = t111 * t 01917 t318 = t115 * t317 01918 t321 = t302 + t305 + 0.2e1_dp * t314 * Arhoa + 0.4e1_dp * t318 * trhoa 01919 t324 = 0.2e1_dp * t297 * t299 + t101 * t111 * t306 * t119 - t310 *& 01920 t313 * t321 01921 t325 = 0.1e1_dp / t122 01922 t326 = t324 * t325 01923 epsilon_cGGArhoa = epsilon_c_unifrhoa + 0.3e1_dp * t293 * t294 + & 01924 t110 * t326 01925 t329 = t126 ** 2 01926 kf_arhoa = t124 / t329 * t90 / 0.3e1_dp 01927 ex_unif_arhoa = -0.3e1_dp / 0.4e1_dp * t7 * kf_arhoa 01928 t335 = kf_a ** 2 01929 t336 = 0.1e1_dp / t335 01930 t337 = my_norm_drhoa * t336 01931 t340 = my_rhoa ** 2 01932 t341 = 0.1e1_dp / t340 01933 s_arhoa = -t337 * t131 * kf_arhoa / 0.2e1_dp - t130 * t341 / 0.2e1_dp 01934 t344 = t137 ** 2 01935 t346 = 0.1e1_dp / t344 * mu 01936 Fx_arhoa = 0.2e1_dp * t346 * s_a * s_arhoa 01937 t350 = my_rhoa * ex_unif_arhoa 01938 01939 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 01940 e_ra(ii) = e_ra(ii)+& 01941 scale_ex * (0.2e1_dp * ex_unif_a * Fx_a + 0.2e1_dp * & 01942 t350 * Fx_a + 0.2e1_dp * t140 * Fx_arhoa) / 0.2e1_dp + scale_ec * (& 01943 epsilon_cGGA + my_rho * epsilon_cGGArhoa) 01944 END IF 01945 01946 chirhob = -t2 - t164 01947 rsrhob = rsrhoa 01948 t368 = t179 * rsrhob / 0.2e1_dp + beta_2_1 * rsrhob + 0.3e1_dp / & 01949 0.2e1_dp * t183 * rsrhob + t22 * t20 * rsrhob * t187 01950 e_c_u_0rhob = -0.2e1_dp * t171 * rsrhob * t28 + t177 * t368 * t191 01951 t382 = t201 * rsrhob / 0.2e1_dp + beta_2_2 * rsrhob + 0.3e1_dp / & 01952 0.2e1_dp * t205 * rsrhob + t40 * t38 * rsrhob * t187 01953 e_c_u_1rhob = -0.2e1_dp * t194 * rsrhob * t46 + t200 * t382 * t212 01954 t396 = t222 * rsrhob / 0.2e1_dp + beta_2_3 * rsrhob + 0.3e1_dp / & 01955 0.2e1_dp * t226 * rsrhob + t58 * t56 * rsrhob * t187 01956 alpha_crhob = 0.2e1_dp * t215 * rsrhob * t64 - t221 * t396 * t233 01957 frhob = (0.4e1_dp / 0.3e1_dp * t71 * chirhob - 0.4e1_dp / 0.3e1_dp& 01958 * t74 * chirhob) * t69 01959 t403 = alpha_crhob * f 01960 t405 = alpha_c * frhob 01961 t407 = t245 * chirhob 01962 t409 = 0.4e1_dp * t77 * t407 01963 t410 = e_c_u_1rhob - e_c_u_0rhob 01964 t411 = t410 * f 01965 t413 = t84 * frhob 01966 t415 = t244 * chirhob 01967 t417 = 0.4e1_dp * t85 * t415 01968 epsilon_c_unifrhob = e_c_u_0rhob + t403 * t82 + t405 * t82 - t409 & 01969 + t411 * t80 + t413 * t80 + t417 01970 phirhob = t257 * chirhob / 0.3e1_dp - t259 * chirhob / 0.3e1_dp 01971 k_frhob = k_frhoa 01972 k_srhob = t266 * k_frhob * t7 01973 trhob = -t269 * t98 * phirhob / 0.2e1_dp - t96 * t274 * k_srhob / & 01974 0.2e1_dp - t278 / 0.2e1_dp 01975 t427 = epsilon_c_unifrhob * t100 01976 t429 = t285 * phirhob 01977 t432 = -t427 * t105 + 0.3e1_dp * t102 * t429 01978 Arhob = -t101 * t281 * t432 * t107 01979 t436 = t123 * phirhob 01980 t439 = t298 * trhob 01981 t442 = Arhob * t111 01982 t444 = 0.2e1_dp * t303 * trhob 01983 t445 = t442 + t444 01984 t453 = t442 + t444 + 0.2e1_dp * t314 * Arhob + 0.4e1_dp * t318 * trhob 01985 t456 = 0.2e1_dp * t297 * t439 + t101 * t111 * t445 * t119 - t310 *& 01986 t313 * t453 01987 t457 = t456 * t325 01988 epsilon_cGGArhob = epsilon_c_unifrhob + 0.3e1_dp * t293 * t436 + & 01989 t110 * t457 01990 t460 = t143 ** 2 01991 kf_brhob = t124 / t460 * t90 / 0.3e1_dp 01992 ex_unif_brhob = -0.3e1_dp / 0.4e1_dp * t7 * kf_brhob 01993 t466 = kf_b ** 2 01994 t467 = 0.1e1_dp / t466 01995 t468 = my_norm_drhob * t467 01996 t471 = my_rhob ** 2 01997 t472 = 0.1e1_dp / t471 01998 s_brhob = -t468 * t148 * kf_brhob / 0.2e1_dp - t147 * t472 / 0.2e1_dp 01999 t475 = t153 ** 2 02000 t477 = 0.1e1_dp / t475 * mu 02001 Fx_brhob = 0.2e1_dp * t477 * s_b * s_brhob 02002 t481 = my_rhob * ex_unif_brhob 02003 02004 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 02005 e_rb(ii) = e_rb(ii)+& 02006 scale_ex * (0.2e1_dp * ex_unif_b * Fx_b + 0.2e1_dp * & 02007 t481 * Fx_b + 0.2e1_dp * t156 * Fx_brhob) / 0.2e1_dp + scale_ec * (& 02008 epsilon_cGGA + my_rho * epsilon_cGGArhob) 02009 END IF 02010 02011 t488 = t95 * t97 02012 tnorm_drho = t488 * t2 / 0.2e1_dp 02013 t493 = t101 * t317 02014 t494 = A * tnorm_drho 02015 t502 = 0.2e1_dp * t303 * tnorm_drho + 0.4e1_dp * t318 * tnorm_drho 02016 t505 = 0.2e1_dp * t297 * t298 * tnorm_drho + 0.2e1_dp * t493 * & 02017 t494 * t119 - t310 * t313 * t502 02018 t506 = t505 * t325 02019 Hnorm_drho = t110 * t506 02020 02021 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 02022 e_ndr(ii) = e_ndr(ii)+& 02023 scale_ec * my_rho * Hnorm_drho 02024 END IF 02025 02026 s_anorm_drhoa = t129 * t131 / 0.2e1_dp 02027 Fx_anorm_drhoa = 0.2e1_dp * t346 * s_a * s_anorm_drhoa 02028 02029 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 02030 e_ndra(ii) = e_ndra(ii)+& 02031 scale_ex * t140 * Fx_anorm_drhoa 02032 END IF 02033 02034 s_bnorm_drhob = t146 * t148 / 0.2e1_dp 02035 Fx_bnorm_drhob = 0.2e1_dp * t477 * s_b * s_bnorm_drhob 02036 02037 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 02038 e_ndrb(ii) = e_ndrb(ii)+& 02039 scale_ex * t156 * Fx_bnorm_drhob 02040 END IF 02041 02042 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 02043 t518 = 0.1e1_dp / t162 / my_rho 02044 t519 = t1 * t518 02045 chirhoarhoa = -0.2e1_dp * t163 + 0.2e1_dp * t519 02046 t523 = 0.1e1_dp / t90 02047 t525 = t162 ** 2 02048 rsrhoarhoa = -t6 / t165 / t8 * t523 / t525 / 0.18e2_dp + & 02049 t6 * t167 * t518 / 0.6e1_dp 02050 t536 = alpha_1_1 * rsrhoa 02051 t538 = t176 * t190 * t191 02052 t543 = t12 / t175 / t23 02053 t544 = t190 ** 2 02054 t548 = 0.1e1_dp / t18 02055 t549 = beta_1_1 * t548 02056 t550 = rsrhoa ** 2 02057 t556 = beta_3_1 * t178 02058 t561 = t20 ** 2 02059 t563 = rs ** 2 02060 t564 = 0.1e1_dp / t563 02061 t576 = t175 ** 2 02062 t578 = t12 / t576 02063 t579 = t27 ** 2 02064 t580 = 0.1e1_dp / t579 02065 e_c_u_0rhoarhoa = -0.2e1_dp * t171 * rsrhoarhoa * t28 + 0.2e1_dp *& 02066 t536 * t538 - 0.2e1_dp * t543 * t544 * t191 + t177 * (-t549 * t550 & 02067 / 0.4e1_dp + t179 * rsrhoarhoa / 0.2e1_dp + beta_2_1 * rsrhoarhoa + & 02068 0.3e1_dp / 0.4e1_dp * t556 * t550 + 0.3e1_dp / 0.2e1_dp * t183 * & 02069 rsrhoarhoa + t22 * t561 * t550 * t564 + t22 * t20 * rsrhoarhoa * & 02070 t187 - t22 * t20 * t550 * t564) * t191 + t578 * t544 * t580 * t14 / & 02071 0.2e1_dp 02072 e_c_u_01rhoa = e_c_u_0rhoa 02073 t588 = alpha_1_2 * rsrhoa 02074 t590 = t199 * t211 * t212 02075 t595 = t32 / t198 / t41 02076 t596 = t211 ** 2 02077 t600 = beta_1_2 * t548 02078 t606 = beta_3_2 * t178 02079 t611 = t38 ** 2 02080 t624 = t198 ** 2 02081 t626 = t32 / t624 02082 t627 = t45 ** 2 02083 t628 = 0.1e1_dp / t627 02084 t636 = alpha_1_3 * rsrhoa 02085 t638 = t220 * t232 * t233 02086 t643 = t50 / t219 / t59 02087 t644 = t232 ** 2 02088 t648 = beta_1_3 * t548 02089 t654 = beta_3_3 * t178 02090 t659 = t56 ** 2 02091 t672 = t219 ** 2 02092 t674 = t50 / t672 02093 t675 = t63 ** 2 02094 t676 = 0.1e1_dp / t675 02095 alpha_c1rhoa = alpha_crhoa 02096 t681 = 0.1e1_dp / t87 02097 t682 = chirhoa ** 2 02098 t687 = 0.1e1_dp / t88 02099 frhoarhoa = (0.4e1_dp / 0.9e1_dp * t681 * t682 + 0.4e1_dp / & 02100 0.3e1_dp * t71 * chirhoarhoa + 0.4e1_dp / 0.9e1_dp * t687 * t682 - & 02101 0.4e1_dp / 0.3e1_dp * t74 * chirhoarhoa) * t69 02102 f1rhoa = frhoa 02103 t705 = alpha_c1rhoa * f 02104 t708 = alpha_c * f1rhoa 02105 t711 = t78 * t79 02106 t726 = e_c_u_1rhoa - e_c_u_01rhoa 02107 t733 = t726 * f 02108 t736 = t84 * f1rhoa 02109 t745 = -0.4e1_dp * t77 * t245 * chirhoarhoa + (-0.2e1_dp * t194 * & 02110 rsrhoarhoa * t46 + 0.2e1_dp * t588 * t590 - 0.2e1_dp * t595 * t596 *& 02111 t212 + t200 * (-t600 * t550 / 0.4e1_dp + t201 * rsrhoarhoa / & 02112 0.2e1_dp + beta_2_2 * rsrhoarhoa + 0.3e1_dp / 0.4e1_dp * t606 * t550& 02113 + 0.3e1_dp / 0.2e1_dp * t205 * rsrhoarhoa + t40 * t611 * t550 * & 02114 t564 + t40 * t38 * rsrhoarhoa * t187 - t40 * t38 * t550 * t564) * & 02115 t212 + t626 * t596 * t628 * t34 / 0.2e1_dp - e_c_u_0rhoarhoa) * f * & 02116 t80 + t249 * f1rhoa * t80 + 0.4e1_dp * t250 * t254 + t726 * frhoa * & 02117 t80 + t84 * frhoarhoa * t80 + 0.4e1_dp * t252 * t254 + 0.4e1_dp * & 02118 t733 * t254 + 0.4e1_dp * t736 * t254 + 0.12e2_dp * t85 * t79 * t682 & 02119 + 0.4e1_dp * t85 * t244 * chirhoarhoa 02120 epsilon_c_unifrhoarhoa = e_c_u_0rhoarhoa + (0.2e1_dp * t215 * & 02121 rsrhoarhoa * t64 - 0.2e1_dp * t636 * t638 + 0.2e1_dp * t643 * t644 *& 02122 t233 - t221 * (-t648 * t550 / 0.4e1_dp + t222 * rsrhoarhoa / & 02123 0.2e1_dp + beta_2_3 * rsrhoarhoa + 0.3e1_dp / 0.4e1_dp * t654 * t550& 02124 + 0.3e1_dp / 0.2e1_dp * t226 * rsrhoarhoa + t58 * t659 * t550 * & 02125 t564 + t58 * t56 * rsrhoarhoa * t187 - t58 * t56 * t550 * t564) * & 02126 t233 - t674 * t644 * t676 * t52 / 0.2e1_dp) * f * t82 + alpha_crhoa & 02127 * f1rhoa * t82 - 0.4e1_dp * t240 * t246 + alpha_c1rhoa * frhoa * t82& 02128 + alpha_c * frhoarhoa * t82 - 0.4e1_dp * t242 * t246 - 0.4e1_dp * & 02129 t705 * t246 - 0.4e1_dp * t708 * t246 - 0.12e2_dp * t77 * t711 * t682& 02130 + t745 02131 epsilon_c_unif1rhoa = e_c_u_01rhoa + t705 * t82 + t708 * t82 - & 02132 t248 + t733 * t80 + t736 * t80 + t256 02133 t750 = 0.1e1_dp / t72 02134 t755 = 0.1e1_dp / t75 02135 phirhoarhoa = -t750 * t682 / 0.9e1_dp + t257 * chirhoarhoa / & 02136 0.3e1_dp - t755 * t682 / 0.9e1_dp - t259 * chirhoarhoa / 0.3e1_dp 02137 phi1rhoa = phirhoa 02138 t763 = t90 ** 2 02139 k_frhoarhoa = -0.2e1_dp / 0.9e1_dp * t3 / t262 / t91 * t763 02140 t767 = 0.1e1_dp / t94 / t93 02141 t768 = k_frhoa ** 2 02142 k_s1rhoa = k_srhoa 02143 t775 = my_norm_drho * t105 * t97 02144 t776 = t2 * phirhoa 02145 t779 = t269 * t273 02146 t785 = t269 * t277 * phirhoa / 0.2e1_dp 02147 t789 = t2 * k_srhoa 02148 t795 = t96 / t272 / k_s 02149 t798 = t273 * t163 02150 t801 = t96 * t798 * k_srhoa / 0.2e1_dp 02151 t812 = t96 * t97 * t518 02152 trhoarhoa = t775 * t776 * phi1rhoa + t779 * t776 * k_s1rhoa / & 02153 0.2e1_dp + t785 - t269 * t98 * phirhoarhoa / 0.2e1_dp + t779 * t789 & 02154 * phi1rhoa / 0.2e1_dp + t795 * t789 * k_s1rhoa + t801 - t96 * t274 *& 02155 (-t767 * t768 * t523 / 0.2e1_dp + t266 * k_frhoarhoa * t7) / & 02156 0.2e1_dp + t269 * t277 * phi1rhoa / 0.2e1_dp + t96 * t798 * k_s1rhoa& 02157 / 0.2e1_dp + t812 02158 t1rhoa = -t269 * t98 * phi1rhoa / 0.2e1_dp - t96 * t274 * k_s1rhoa& 02159 / 0.2e1_dp - t278 / 0.2e1_dp 02160 t820 = t101 / t280 / t108 02161 t821 = t107 ** 2 02162 t822 = t289 * t821 02163 t823 = epsilon_c_unif1rhoa * t100 02164 t825 = t285 * phi1rhoa 02165 t828 = -t823 * t105 + 0.3e1_dp * t102 * t825 02166 t839 = 0.1e1_dp / t284 / phi 02167 t840 = t839 * phirhoa 02168 t851 = t101 * t281 02169 Arhoarhoa = 0.2e1_dp * t820 * t822 * t828 - t101 * t281 * (& 02170 -epsilon_c_unifrhoarhoa * t100 * t105 + 0.3e1_dp * t282 * t825 + & 02171 0.3e1_dp * t823 * t286 - 0.12e2_dp * t102 * t840 * phi1rhoa + & 02172 0.3e1_dp * t102 * t285 * phirhoarhoa) * t107 - t851 * t289 * t828 * & 02173 t107 02174 A1rhoa = -t101 * t281 * t828 * t107 02175 t858 = gamma_var * phi 02176 t865 = A1rhoa * t111 02177 t867 = 0.2e1_dp * t303 * t1rhoa 02178 t868 = t865 + t867 02179 t876 = t865 + t867 + 0.2e1_dp * t314 * A1rhoa + 0.4e1_dp * t318 * t1rhoa 02180 t879 = 0.2e1_dp * t297 * t298 * t1rhoa + t101 * t111 * t868 * t119& 02181 - t310 * t313 * t876 02182 t880 = t879 * t325 02183 t904 = t306 * t119 02184 t908 = Arhoarhoa * t111 02185 t909 = Arhoa * t 02186 t911 = 0.2e1_dp * t909 * t1rhoa 02187 t914 = 0.2e1_dp * A1rhoa * t * trhoa 02188 t917 = 0.2e1_dp * A * t1rhoa * trhoa 02189 t919 = 0.2e1_dp * t303 * trhoarhoa 02190 t924 = t306 * t312 02191 t936 = t113 / t311 / t118 02192 t944 = A * t317 02193 t953 = t115 * t111 02194 t959 = t908 + t911 + t914 + t917 + t919 + 0.2e1_dp * A1rhoa * t116& 02195 * Arhoa + 0.8e1_dp * t944 * Arhoa * t1rhoa + 0.2e1_dp * t314 * & 02196 Arhoarhoa + 0.8e1_dp * t944 * trhoa * A1rhoa + 0.12e2_dp * t953 * & 02197 trhoa * t1rhoa + 0.4e1_dp * t318 * trhoarhoa 02198 t962 = 0.2e1_dp * t101 * t1rhoa * t299 + 0.2e1_dp * t297 * t868 * & 02199 t119 * trhoa - 0.2e1_dp * t297 * t313 * trhoa * t876 + 0.2e1_dp * & 02200 t297 * t298 * trhoarhoa + 0.2e1_dp * t297 * t904 * t1rhoa + t101 * & 02201 t111 * (t908 + t911 + t914 + t917 + t919) * t119 - t310 * t924 * & 02202 t876 - 0.2e1_dp * t297 * t313 * t321 * t1rhoa - t310 * t868 * t312 *& 02203 t321 + 0.2e1_dp * t310 * t936 * t321 * t876 - t310 * t313 * t959 02204 t965 = t122 ** 2 02205 t966 = 0.1e1_dp / t965 02206 t967 = t324 * t966 02207 kf_arhoarhoa = -0.2e1_dp / 0.9e1_dp * t124 / t329 / t125 * t763 02208 ex_unif_a1rhoa = ex_unif_arhoa 02209 t985 = kf_arhoa ** 2 02210 s_a1rhoa = s_arhoa 02211 t998 = mu ** 2 02212 t999 = 0.1e1_dp / t344 / t137 * t998 02213 t1000 = t999 * t133 02214 t1001 = s_arhoa * t135 02215 Fx_a1rhoa = 0.2e1_dp * t346 * s_a * s_a1rhoa 02216 02217 e_ra_ra(ii) = e_ra_ra(ii)+& 02218 scale_ex * (0.2e1_dp * ex_unif_a1rhoa * Fx_a + & 02219 0.2e1_dp * ex_unif_a * Fx_a1rhoa + 0.2e1_dp * ex_unif_arhoa * Fx_a -& 02220 0.3e1_dp / 0.2e1_dp * my_rhoa * t7 * kf_arhoarhoa * Fx_a + 0.2e1_dp * & 02221 t350 * Fx_a1rhoa + 0.2e1_dp * ex_unif_a * Fx_arhoa + 0.2e1_dp * my_rhoa& 02222 * ex_unif_a1rhoa * Fx_arhoa + 0.2e1_dp * t140 * (-0.8e1_dp * t1000 & 02223 * t1001 * s_a1rhoa + 0.2e1_dp * t346 * s_a1rhoa * s_arhoa + 0.2e1_dp& 02224 * t346 * s_a * (my_norm_drhoa / t335 / kf_a * t131 * t985 + t337 * & 02225 t341 * kf_arhoa - t337 * t131 * kf_arhoarhoa / 0.2e1_dp + t130 / & 02226 t340 / my_rhoa))) / 0.2e1_dp + scale_ec * (epsilon_c_unif1rhoa + & 02227 0.3e1_dp * t293 * t123 * phi1rhoa + t110 * t880 + epsilon_cGGArhoa +& 02228 my_rho * (epsilon_c_unifrhoarhoa + 0.6e1_dp * t858 * t294 * phi1rhoa +& 02229 0.3e1_dp * t293 * t880 * phirhoa + 0.3e1_dp * t293 * t123 * & 02230 phirhoarhoa + 0.3e1_dp * t293 * t326 * phi1rhoa + t110 * t962 * t325& 02231 - t110 * t967 * t879)) 02232 02233 chirhoarhob = 0.2e1_dp * t519 02234 rsrhoarhob = rsrhoarhoa 02235 t1031 = t176 * t368 * t191 02236 t1033 = alpha_1_1 * rsrhob 02237 t1038 = rsrhoa * rsrhob 02238 t1050 = rsrhob * t564 * rsrhoa 02239 e_c_u_0rhoarhob = -0.2e1_dp * t171 * rsrhoarhob * t28 + t536 * & 02240 t1031 + t1033 * t538 - 0.2e1_dp * t543 * t192 * t368 + t177 * (-t549& 02241 * t1038 / 0.4e1_dp + t179 * rsrhoarhob / 0.2e1_dp + beta_2_1 * & 02242 rsrhoarhob + 0.3e1_dp / 0.4e1_dp * t556 * t1038 + 0.3e1_dp / & 02243 0.2e1_dp * t183 * rsrhoarhob + t22 * t561 * t1050 + t22 * t20 * & 02244 rsrhoarhob * t187 - t22 * t20 * t1050) * t191 + t578 * t190 * t580 *& 02245 t14 * t368 / 0.2e1_dp 02246 t1069 = t199 * t382 * t212 02247 t1071 = alpha_1_2 * rsrhob 02248 t1104 = t220 * t396 * t233 02249 t1106 = alpha_1_3 * rsrhob 02250 frhoarhob = (0.4e1_dp / 0.9e1_dp * t681 * chirhoa * chirhob + & 02251 0.4e1_dp / 0.3e1_dp * t71 * chirhoarhob + 0.4e1_dp / 0.9e1_dp * t687& 02252 * chirhoa * chirhob - 0.4e1_dp / 0.3e1_dp * t74 * chirhoarhob) * & 02253 t69 02254 t1164 = t79 * chirhoa * chirhob 02255 t1193 = -0.4e1_dp * t77 * t245 * chirhoarhob + (-0.2e1_dp * t194 *& 02256 rsrhoarhob * t46 + t588 * t1069 + t1071 * t590 - 0.2e1_dp * t595 * & 02257 t213 * t382 + t200 * (-t600 * t1038 / 0.4e1_dp + t201 * rsrhoarhob /& 02258 0.2e1_dp + beta_2_2 * rsrhoarhob + 0.3e1_dp / 0.4e1_dp * t606 * & 02259 t1038 + 0.3e1_dp / 0.2e1_dp * t205 * rsrhoarhob + t40 * t611 * t1050& 02260 + t40 * t38 * rsrhoarhob * t187 - t40 * t38 * t1050) * t212 + t626 & 02261 * t211 * t628 * t34 * t382 / 0.2e1_dp - e_c_u_0rhoarhob) * f * t80 +& 02262 t249 * frhob * t80 + 0.4e1_dp * t250 * t415 + t410 * frhoa * t80 + & 02263 t84 * frhoarhob * t80 + 0.4e1_dp * t252 * t415 + 0.4e1_dp * t411 * & 02264 t254 + 0.4e1_dp * t413 * t254 + 0.12e2_dp * t85 * t1164 + 0.4e1_dp *& 02265 t85 * t244 * chirhoarhob 02266 epsilon_c_unifrhoarhob = e_c_u_0rhoarhob + (0.2e1_dp * t215 * & 02267 rsrhoarhob * t64 - t636 * t1104 - t1106 * t638 + 0.2e1_dp * t643 * & 02268 t234 * t396 - t221 * (-t648 * t1038 / 0.4e1_dp + t222 * rsrhoarhob /& 02269 0.2e1_dp + beta_2_3 * rsrhoarhob + 0.3e1_dp / 0.4e1_dp * t654 * & 02270 t1038 + 0.3e1_dp / 0.2e1_dp * t226 * rsrhoarhob + t58 * t659 * t1050& 02271 + t58 * t56 * rsrhoarhob * t187 - t58 * t56 * t1050) * t233 - t674 & 02272 * t232 * t676 * t52 * t396 / 0.2e1_dp) * f * t82 + alpha_crhoa * & 02273 frhob * t82 - 0.4e1_dp * t240 * t407 + alpha_crhob * frhoa * t82 + & 02274 alpha_c * frhoarhob * t82 - 0.4e1_dp * t242 * t407 - 0.4e1_dp * t403& 02275 * t246 - 0.4e1_dp * t405 * t246 - 0.12e2_dp * t77 * t78 * t1164 + & 02276 t1193 02277 phirhoarhob = -t750 * chirhoa * chirhob / 0.9e1_dp + t257 * & 02278 chirhoarhob / 0.3e1_dp - t755 * chirhoa * chirhob / 0.9e1_dp - t259 & 02279 * chirhoarhob / 0.3e1_dp 02280 k_frhoarhob = k_frhoarhoa 02281 t1228 = t269 * t277 * phirhob / 0.2e1_dp 02282 t1231 = t96 * t798 * k_srhob / 0.2e1_dp 02283 trhoarhob = t775 * t776 * phirhob + t779 * t776 * k_srhob / & 02284 0.2e1_dp + t785 - t269 * t98 * phirhoarhob / 0.2e1_dp + t779 * t789 & 02285 * phirhob / 0.2e1_dp + t795 * t789 * k_srhob + t801 - t96 * t274 * (& 02286 -t767 * k_frhoa * t523 * k_frhob / 0.2e1_dp + t266 * k_frhoarhob * & 02287 t7) / 0.2e1_dp + t1228 + t1231 + t812 02288 Arhoarhob = 0.2e1_dp * t820 * t822 * t432 - t101 * t281 * (& 02289 -epsilon_c_unifrhoarhob * t100 * t105 + 0.3e1_dp * t282 * t429 + & 02290 0.3e1_dp * t427 * t286 - 0.12e2_dp * t102 * t840 * phirhob + & 02291 0.3e1_dp * t102 * t285 * phirhoarhob) * t107 - t851 * t289 * t432 * & 02292 t107 02293 t1269 = t445 * t119 02294 t1283 = Arhoarhob * t111 02295 t1285 = 0.2e1_dp * t909 * trhob 02296 t1286 = Arhob * t 02297 t1288 = 0.2e1_dp * t1286 * trhoa 02298 t1291 = 0.2e1_dp * A * trhob * trhoa 02299 t1293 = 0.2e1_dp * t303 * trhoarhob 02300 t1304 = t445 * t312 02301 t1327 = t1283 + t1285 + t1288 + t1291 + t1293 + 0.2e1_dp * Arhob *& 02302 t116 * Arhoa + 0.8e1_dp * t944 * Arhoa * trhob + 0.2e1_dp * t314 * & 02303 Arhoarhob + 0.8e1_dp * t944 * trhoa * Arhob + 0.12e2_dp * t953 * & 02304 trhoa * trhob + 0.4e1_dp * t318 * trhoarhob 02305 t1330 = 0.2e1_dp * t101 * trhob * t299 + 0.2e1_dp * t297 * t1269 *& 02306 trhoa - 0.2e1_dp * t297 * t313 * trhoa * t453 + 0.2e1_dp * t297 * & 02307 t298 * trhoarhob + 0.2e1_dp * t297 * t904 * trhob + t101 * t111 * (& 02308 t1283 + t1285 + t1288 + t1291 + t1293) * t119 - t310 * t924 * t453 -& 02309 0.2e1_dp * t297 * t313 * t321 * trhob - t310 * t1304 * t321 + & 02310 0.2e1_dp * t310 * t936 * t321 * t453 - t310 * t313 * t1327 02311 02312 e_ra_rb(ii) = e_ra_rb(ii)+& 02313 scale_ec * (epsilon_cGGArhob + epsilon_cGGArhoa + & 02314 my_rho * (epsilon_c_unifrhoarhob + 0.6e1_dp * t858 * t294 * phirhob + & 02315 0.3e1_dp * t293 * t457 * phirhoa + 0.3e1_dp * t293 * t123 * & 02316 phirhoarhob + 0.3e1_dp * t293 * t326 * phirhob + t110 * t1330 * t325& 02317 - t110 * t967 * t456)) 02318 02319 chirhobrhob = 0.2e1_dp * t163 + 0.2e1_dp * t519 02320 rsrhobrhob = rsrhoarhob 02321 t1342 = t368 ** 2 02322 t1346 = rsrhob ** 2 02323 e_c_u_0rhobrhob = -0.2e1_dp * t171 * rsrhobrhob * t28 + 0.2e1_dp *& 02324 t1033 * t1031 - 0.2e1_dp * t543 * t1342 * t191 + t177 * (-t549 * & 02325 t1346 / 0.4e1_dp + t179 * rsrhobrhob / 0.2e1_dp + beta_2_1 * & 02326 rsrhobrhob + 0.3e1_dp / 0.4e1_dp * t556 * t1346 + 0.3e1_dp / & 02327 0.2e1_dp * t183 * rsrhobrhob + t22 * t561 * t1346 * t564 + t22 * t20& 02328 * rsrhobrhob * t187 - t22 * t20 * t1346 * t564) * t191 + t578 * & 02329 t1342 * t580 * t14 / 0.2e1_dp 02330 e_c_u_01rhob = e_c_u_0rhob 02331 t1377 = t382 ** 2 02332 t1411 = t396 ** 2 02333 alpha_c1rhob = alpha_crhob 02334 t1440 = chirhob ** 2 02335 frhobrhob = (0.4e1_dp / 0.9e1_dp * t681 * t1440 + 0.4e1_dp / & 02336 0.3e1_dp * t71 * chirhobrhob + 0.4e1_dp / 0.9e1_dp * t687 * t1440 - & 02337 0.4e1_dp / 0.3e1_dp * t74 * chirhobrhob) * t69 02338 f1rhob = frhob 02339 t1462 = alpha_c1rhob * f 02340 t1465 = alpha_c * f1rhob 02341 t1482 = e_c_u_1rhob - e_c_u_01rhob 02342 t1489 = t1482 * f 02343 t1492 = t84 * f1rhob 02344 t1501 = -0.4e1_dp * t77 * t245 * chirhobrhob + (-0.2e1_dp * t194 *& 02345 rsrhobrhob * t46 + 0.2e1_dp * t1071 * t1069 - 0.2e1_dp * t595 * & 02346 t1377 * t212 + t200 * (-t600 * t1346 / 0.4e1_dp + t201 * rsrhobrhob & 02347 / 0.2e1_dp + beta_2_2 * rsrhobrhob + 0.3e1_dp / 0.4e1_dp * t606 * & 02348 t1346 + 0.3e1_dp / 0.2e1_dp * t205 * rsrhobrhob + t40 * t611 * t1346& 02349 * t564 + t40 * t38 * rsrhobrhob * t187 - t40 * t38 * t1346 * t564) & 02350 * t212 + t626 * t1377 * t628 * t34 / 0.2e1_dp - e_c_u_0rhobrhob) * f& 02351 * t80 + t410 * f1rhob * t80 + 0.4e1_dp * t411 * t415 + t1482 * & 02352 frhob * t80 + t84 * frhobrhob * t80 + 0.4e1_dp * t413 * t415 + & 02353 0.4e1_dp * t1489 * t415 + 0.4e1_dp * t1492 * t415 + 0.12e2_dp * t85 & 02354 * t79 * t1440 + 0.4e1_dp * t85 * t244 * chirhobrhob 02355 epsilon_c_unifrhobrhob = e_c_u_0rhobrhob + (0.2e1_dp * t215 * & 02356 rsrhobrhob * t64 - 0.2e1_dp * t1106 * t1104 + 0.2e1_dp * t643 * & 02357 t1411 * t233 - t221 * (-t648 * t1346 / 0.4e1_dp + t222 * rsrhobrhob & 02358 / 0.2e1_dp + beta_2_3 * rsrhobrhob + 0.3e1_dp / 0.4e1_dp * t654 * & 02359 t1346 + 0.3e1_dp / 0.2e1_dp * t226 * rsrhobrhob + t58 * t659 * t1346& 02360 * t564 + t58 * t56 * rsrhobrhob * t187 - t58 * t56 * t1346 * t564) & 02361 * t233 - t674 * t1411 * t676 * t52 / 0.2e1_dp) * f * t82 + & 02362 alpha_crhob * f1rhob * t82 - 0.4e1_dp * t403 * t407 + alpha_c1rhob *& 02363 frhob * t82 + alpha_c * frhobrhob * t82 - 0.4e1_dp * t405 * t407 - & 02364 0.4e1_dp * t1462 * t407 - 0.4e1_dp * t1465 * t407 - 0.12e2_dp * t77 & 02365 * t711 * t1440 + t1501 02366 epsilon_c_unif1rhob = e_c_u_01rhob + t1462 * t82 + t1465 * t82 - & 02367 t409 + t1489 * t80 + t1492 * t80 + t417 02368 phirhobrhob = -t750 * t1440 / 0.9e1_dp + t257 * chirhobrhob / & 02369 0.3e1_dp - t755 * t1440 / 0.9e1_dp - t259 * chirhobrhob / 0.3e1_dp 02370 phi1rhob = phirhob 02371 t1514 = k_frhob ** 2 02372 k_s1rhob = k_srhob 02373 t1520 = t2 * phirhob 02374 t1529 = t2 * k_srhob 02375 trhobrhob = t775 * t1520 * phi1rhob + t779 * t1520 * k_s1rhob / & 02376 0.2e1_dp + t1228 - t269 * t98 * phirhobrhob / 0.2e1_dp + t779 * & 02377 t1529 * phi1rhob / 0.2e1_dp + t795 * t1529 * k_s1rhob + t1231 - t96 & 02378 * t274 * (-t767 * t1514 * t523 / 0.2e1_dp + t266 * k_frhoarhob * t7)& 02379 / 0.2e1_dp + t269 * t277 * phi1rhob / 0.2e1_dp + t96 * t798 * & 02380 k_s1rhob / 0.2e1_dp + t812 02381 t1rhob = -t269 * t98 * phi1rhob / 0.2e1_dp - t96 * t274 * k_s1rhob& 02382 / 0.2e1_dp - t278 / 0.2e1_dp 02383 t1550 = epsilon_c_unif1rhob * t100 02384 t1552 = t285 * phi1rhob 02385 t1555 = -t1550 * t105 + 0.3e1_dp * t102 * t1552 02386 Arhobrhob = 0.2e1_dp * t820 * t432 * t821 * t1555 - t101 * t281 * & 02387 (-epsilon_c_unifrhobrhob * t100 * t105 + 0.3e1_dp * t427 * t1552 + & 02388 0.3e1_dp * t1550 * t429 - 0.12e2_dp * t102 * t839 * phirhob * & 02389 phi1rhob + 0.3e1_dp * t102 * t285 * phirhobrhob) * t107 - t851 * & 02390 t432 * t1555 * t107 02391 A1rhob = -t101 * t281 * t1555 * t107 02392 t1588 = A1rhob * t111 02393 t1590 = 0.2e1_dp * t303 * t1rhob 02394 t1591 = t1588 + t1590 02395 t1599 = t1588 + t1590 + 0.2e1_dp * t314 * A1rhob + 0.4e1_dp * t318& 02396 * t1rhob 02397 t1602 = 0.2e1_dp * t297 * t298 * t1rhob + t101 * t111 * t1591 * & 02398 t119 - t310 * t313 * t1599 02399 t1603 = t1602 * t325 02400 t1630 = Arhobrhob * t111 02401 t1632 = 0.2e1_dp * t1286 * t1rhob 02402 t1635 = 0.2e1_dp * A1rhob * t * trhob 02403 t1638 = 0.2e1_dp * A * t1rhob * trhob 02404 t1640 = 0.2e1_dp * t303 * trhobrhob 02405 t1674 = t1630 + t1632 + t1635 + t1638 + t1640 + 0.2e1_dp * A1rhob & 02406 * t116 * Arhob + 0.8e1_dp * t944 * Arhob * t1rhob + 0.2e1_dp * t314 & 02407 * Arhobrhob + 0.8e1_dp * t944 * trhob * A1rhob + 0.12e2_dp * t953 * & 02408 trhob * t1rhob + 0.4e1_dp * t318 * trhobrhob 02409 t1677 = 0.2e1_dp * t101 * t1rhob * t439 + 0.2e1_dp * t297 * t1591 & 02410 * t119 * trhob - 0.2e1_dp * t297 * t313 * trhob * t1599 + 0.2e1_dp *& 02411 t297 * t298 * trhobrhob + 0.2e1_dp * t297 * t1269 * t1rhob + t101 *& 02412 t111 * (t1630 + t1632 + t1635 + t1638 + t1640) * t119 - t310 * & 02413 t1304 * t1599 - 0.2e1_dp * t297 * t313 * t453 * t1rhob - t310 * & 02414 t1591 * t312 * t453 + 0.2e1_dp * t310 * t936 * t453 * t1599 - t310 *& 02415 t313 * t1674 02416 t1680 = t456 * t966 02417 kf_brhobrhob = -0.2e1_dp / 0.9e1_dp * t124 / t460 / t142 * t763 02418 ex_unif_b1rhob = ex_unif_brhob 02419 t1698 = kf_brhob ** 2 02420 s_b1rhob = s_brhob 02421 t1711 = 0.1e1_dp / t475 / t153 * t998 02422 t1712 = t1711 * t150 02423 t1713 = s_brhob * t135 02424 Fx_b1rhob = 0.2e1_dp * t477 * s_b * s_b1rhob 02425 02426 e_rb_rb(ii) = e_rb_rb(ii)+& 02427 scale_ex * (0.2e1_dp * ex_unif_b1rhob * Fx_b + & 02428 0.2e1_dp * ex_unif_b * Fx_b1rhob + 0.2e1_dp * ex_unif_brhob * Fx_b -& 02429 0.3e1_dp / 0.2e1_dp * my_rhob * t7 * kf_brhobrhob * Fx_b + 0.2e1_dp * & 02430 t481 * Fx_b1rhob + 0.2e1_dp * ex_unif_b * Fx_brhob + 0.2e1_dp * my_rhob& 02431 * ex_unif_b1rhob * Fx_brhob + 0.2e1_dp * t156 * (-0.8e1_dp * t1712 & 02432 * t1713 * s_b1rhob + 0.2e1_dp * t477 * s_b1rhob * s_brhob + 0.2e1_dp& 02433 * t477 * s_b * (my_norm_drhob / t466 / kf_b * t148 * t1698 + t468 * & 02434 t472 * kf_brhob - t468 * t148 * kf_brhobrhob / 0.2e1_dp + t147 / & 02435 t471 / my_rhob))) / 0.2e1_dp + scale_ec * (epsilon_c_unif1rhob + & 02436 0.3e1_dp * t293 * t123 * phi1rhob + t110 * t1603 + epsilon_cGGArhob & 02437 + my_rho * (epsilon_c_unifrhobrhob + 0.6e1_dp * t858 * t436 * phi1rhob & 02438 + 0.3e1_dp * t293 * t1603 * phirhob + 0.3e1_dp * t293 * t123 * & 02439 phirhobrhob + 0.3e1_dp * t293 * t457 * phi1rhob + t110 * t1677 * & 02440 t325 - t110 * t1680 * t1602)) 02441 t1739 = t268 * t97 02442 t1741 = t95 * t273 02443 t1743 = t488 * t163 02444 trhoanorm_drho = -t1739 * t776 / 0.2e1_dp - t1741 * t789 / & 02445 0.2e1_dp - t1743 / 0.2e1_dp 02446 t1748 = t101 * tnorm_drho 02447 t1765 = t909 * tnorm_drho 02448 t1766 = t494 * trhoa 02449 t1767 = t303 * trhoanorm_drho 02450 t1801 = 0.2e1_dp * t1748 * t299 + 0.4e1_dp * t310 * t494 * t119 * & 02451 trhoa - 0.2e1_dp * t297 * t313 * trhoa * t502 + 0.2e1_dp * t297 * & 02452 t298 * trhoanorm_drho + 0.2e1_dp * t297 * t904 * tnorm_drho + t101 *& 02453 t111 * (0.2e1_dp * t1765 + 0.2e1_dp * t1766 + 0.2e1_dp * t1767) * & 02454 t119 - t310 * t924 * t502 - 0.2e1_dp * t297 * t313 * t321 * & 02455 tnorm_drho - 0.2e1_dp * t493 * t494 * t312 * t321 + 0.2e1_dp * t310 & 02456 * t936 * t321 * t502 - t310 * t313 * (0.2e1_dp * t1765 + 0.2e1_dp * & 02457 t1766 + 0.2e1_dp * t1767 + 0.8e1_dp * t944 * Arhoa * tnorm_drho + & 02458 0.12e2_dp * t953 * trhoa * tnorm_drho + 0.4e1_dp * t318 * & 02459 trhoanorm_drho) 02460 02461 e_ra_ndr(ii) = e_ra_ndr(ii)+& 02462 scale_ec * (Hnorm_drho + my_rho * (0.3e1_dp * & 02463 t293 * t506 * phirhoa + t110 * t1801 * t325 - t110 * t967 * t505)) 02464 02465 trhobnorm_drho = -t1739 * t1520 / 0.2e1_dp - t1741 * t1529 / & 02466 0.2e1_dp - t1743 / 0.2e1_dp 02467 t1829 = t1286 * tnorm_drho 02468 t1830 = t494 * trhob 02469 t1831 = t303 * trhobnorm_drho 02470 t1865 = 0.2e1_dp * t1748 * t439 + 0.4e1_dp * t310 * t494 * t119 * & 02471 trhob - 0.2e1_dp * t297 * t313 * trhob * t502 + 0.2e1_dp * t297 * & 02472 t298 * trhobnorm_drho + 0.2e1_dp * t297 * t1269 * tnorm_drho + t101 & 02473 * t111 * (0.2e1_dp * t1829 + 0.2e1_dp * t1830 + 0.2e1_dp * t1831) * & 02474 t119 - t310 * t1304 * t502 - 0.2e1_dp * t297 * t313 * t453 * & 02475 tnorm_drho - 0.2e1_dp * t493 * t494 * t312 * t453 + 0.2e1_dp * t310 & 02476 * t936 * t453 * t502 - t310 * t313 * (0.2e1_dp * t1829 + 0.2e1_dp * & 02477 t1830 + 0.2e1_dp * t1831 + 0.8e1_dp * t944 * Arhob * tnorm_drho + & 02478 0.12e2_dp * t953 * trhob * tnorm_drho + 0.4e1_dp * t318 * & 02479 trhobnorm_drho) 02480 02481 e_rb_ndr(ii) = e_rb_ndr(ii)+& 02482 scale_ec * (Hnorm_drho + my_rho * (0.3e1_dp * & 02483 t293 * t506 * phirhob + t110 * t1865 * t325 - t110 * t1680 * t505)) 02484 02485 t1871 = tnorm_drho ** 2 02486 t1876 = A * t1871 02487 t1888 = t502 ** 2 02488 t1901 = t505 ** 2 02489 02490 e_ndr_ndr(ii) = e_ndr_ndr(ii)+& 02491 scale_ec * my_rho * (t110 * (0.2e1_dp * & 02492 t101 * t1871 * t113 * t119 + 0.10e2_dp * t310 * t1876 * t119 - & 02493 0.4e1_dp * t297 * t313 * tnorm_drho * t502 - 0.4e1_dp * t493 * t494 & 02494 * t312 * t502 + 0.2e1_dp * t310 * t936 * t1888 - t310 * t313 * (& 02495 0.2e1_dp * t1876 + 0.12e2_dp * t953 * t1871)) * t325 - t110 * t1901 & 02496 * t966) 02497 02498 e_ra_ndra(ii) = e_ra_ndra(ii)+& 02499 scale_ex * (0.2e1_dp * ex_unif_a * & 02500 Fx_anorm_drhoa + 0.2e1_dp * t350 * Fx_anorm_drhoa + 0.2e1_dp * t140 & 02501 * (-0.8e1_dp * t1000 * t1001 * s_anorm_drhoa + 0.2e1_dp * t346 * & 02502 s_anorm_drhoa * s_arhoa + 0.2e1_dp * t346 * s_a * (-t336 * t131 * & 02503 kf_arhoa / 0.2e1_dp - t129 * t341 / 0.2e1_dp))) / 0.2e1_dp 02504 02505 t1922 = s_anorm_drhoa ** 2 02506 02507 e_ndra_ndra(ii) = e_ndra_ndra(ii)+& 02508 scale_ex * t140 * (-0.8e1_dp * t999 * & 02509 t133 * t1922 * t135 + 0.2e1_dp * t346 * t1922) 02510 02511 e_rb_ndrb(ii) = e_rb_ndrb(ii)+& 02512 scale_ex * (0.2e1_dp * ex_unif_b * & 02513 Fx_bnorm_drhob + 0.2e1_dp * t481 * Fx_bnorm_drhob + 0.2e1_dp * t156 & 02514 * (-0.8e1_dp * t1712 * t1713 * s_bnorm_drhob + 0.2e1_dp * t477 * & 02515 s_bnorm_drhob * s_brhob + 0.2e1_dp * t477 * s_b * (-t467 * t148 * & 02516 kf_brhob / 0.2e1_dp - t146 * t472 / 0.2e1_dp))) / 0.2e1_dp 02517 02518 t1949 = s_bnorm_drhob ** 2 02519 e_ndrb_ndrb(ii) = e_ndrb_ndrb(ii)+& 02520 scale_ex * t156 * (-0.8e1_dp * t1711 *& 02521 t150 * t1949 * t135 + 0.2e1_dp * t477 * t1949) 02522 END IF 02523 END IF 02524 END DO 02525 !$omp end do 02526 END SELECT 02527 END SUBROUTINE pbe_lsd_calc 02528 02529 END MODULE xc_pbe
1.7.3