1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief calculates the pbe correlation functional 8!> \note 9!> This was generated with the help of a maple worksheet that you can 10!> find under doc/pbe.mw . 11!> I did not add 3. derivatives in the polazied (lsd) case because it 12!> would have added 2500 lines of code. If you need them the worksheet 13!> is already prepared for them, and by uncommenting a couple of lines 14!> you should be able to generate them. 15!> \par History 16!> 09.2004 created [fawzi] 17!> \author fawzi 18! ************************************************************************************************** 19MODULE xc_pbe 20 USE bibliography, ONLY: Perdew1996,& 21 Perdew2008,& 22 Zhang1998,& 23 cite_reference 24 USE input_section_types, ONLY: section_vals_type,& 25 section_vals_val_get 26 USE kinds, ONLY: dp 27 USE mathconstants, ONLY: pi 28 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 29 xc_dset_get_derivative 30 USE xc_derivative_types, ONLY: xc_derivative_get,& 31 xc_derivative_type 32 USE xc_input_constants, ONLY: xc_pbe_orig,& 33 xc_pbe_rev,& 34 xc_pbe_sol 35 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 36 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 37 xc_rho_set_type 38#include "../base/base_uses.f90" 39 40 IMPLICIT NONE 41 PRIVATE 42 43 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 44 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_pbe' 45 REAL(kind=dp), PARAMETER, PRIVATE :: a = 0.04918_dp, b = 0.132_dp, & 46 c = 0.2533_dp, d = 0.349_dp 47 48 PUBLIC :: pbe_lda_info, pbe_lsd_info, pbe_lda_eval, pbe_lsd_eval 49CONTAINS 50 51! ************************************************************************************************** 52!> \brief return various information on the functional 53!> \param pbe_params section selecting the varius parameters for the functional 54!> \param reference string with the reference of the actual functional 55!> \param shortform string with the shortform of the functional name 56!> \param needs the components needed by this functional are set to 57!> true (does not set the unneeded components to false) 58!> \param max_deriv ... 59!> \author fawzi 60! ************************************************************************************************** 61 SUBROUTINE pbe_lda_info(pbe_params, reference, shortform, needs, max_deriv) 62 TYPE(section_vals_type), POINTER :: pbe_params 63 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 64 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs 65 INTEGER, INTENT(out), OPTIONAL :: max_deriv 66 67 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lda_info', routineP = moduleN//':'//routineN 68 69 INTEGER :: param 70 REAL(kind=dp) :: sc, sx 71 72 CALL section_vals_val_get(pbe_params, "parametrization", i_val=param) 73 CALL section_vals_val_get(pbe_params, "scale_x", r_val=sx) 74 CALL section_vals_val_get(pbe_params, "scale_c", r_val=sc) 75 76 SELECT CASE (param) 77 CASE (xc_pbe_orig) 78 CALL cite_reference(Perdew1996) 79 IF (sx == 1._dp .AND. sc == 1._dp) THEN 80 IF (PRESENT(reference)) THEN 81 reference = "J.P.Perdew, K.Burke, M.Ernzerhof, "// & 82 "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)"// & 83 "{spin unpolarized}" 84 END IF 85 IF (PRESENT(shortform)) THEN 86 shortform = "PBE Perdew-Burke-Ernzerhof xc functional (unpolarized)" 87 END IF 88 ELSE 89 IF (PRESENT(reference)) THEN 90 WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") & 91 "J.P.Perdew, K.Burke, M.Ernzerhof, ", & 92 "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)", & 93 sx, sc, "{spin unpolarized}" 94 END IF 95 IF (PRESENT(shortform)) THEN 96 WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')") & 97 "PBE, Perdew-Burke-Ernzerhof xc functional (unpolarized)", sx, sc 98 END IF 99 END IF 100 CASE (xc_pbe_rev) 101 CALL cite_reference(Perdew1996) 102 CALL cite_reference(Zhang1998) 103 IF (sx == 1._dp .AND. sc == 1._dp) THEN 104 IF (PRESENT(reference)) THEN 105 reference = "revPBE, Yingkay Zhang and Weitao Yang,"// & 106 " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)"// & 107 "{spin unpolarized}" 108 END IF 109 IF (PRESENT(shortform)) THEN 110 shortform = "revPBE, revised PBE xc functional (unpolarized)" 111 END IF 112 ELSE 113 IF (PRESENT(reference)) THEN 114 WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") & 115 "revPBE, Yingkay Zhang and Weitao Yang,", & 116 " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)", & 117 sx, sc, "{spin unpolarized}" 118 END IF 119 IF (PRESENT(shortform)) THEN 120 WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')") & 121 "revPBE, revised PBE xc functional (unpolarized)", sx, sc 122 END IF 123 END IF 124 CASE (xc_pbe_sol) 125 CALL cite_reference(Perdew1996) 126 CALL cite_reference(Perdew2008) 127 IF (sx == 1._dp .AND. sc == 1._dp) THEN 128 IF (PRESENT(reference)) THEN 129 reference = "PBEsol, J.P. Perdew et al., "// & 130 "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) "// & 131 "{spin unpolarized}" 132 END IF 133 IF (PRESENT(shortform)) THEN 134 shortform = "PBEsol, PBE for solids and surfaces xc functional (unpolarized)" 135 END IF 136 ELSE 137 IF (PRESENT(reference)) THEN 138 WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") & 139 "PBEsol, J.P. Perdew et al., ", & 140 "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) ", & 141 sx, sc, "{spin unpolarized}" 142 END IF 143 IF (PRESENT(shortform)) THEN 144 WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LDA)')") & 145 "PBEsol, PBE for solids and surfaces xc functional (unpolarized)", sx, sc 146 END IF 147 END IF 148 CASE default 149 CPABORT("") 150 END SELECT 151 IF (PRESENT(needs)) THEN 152 needs%rho = .TRUE. 153 needs%norm_drho = .TRUE. 154 END IF 155 IF (PRESENT(max_deriv)) max_deriv = 3 156 157 END SUBROUTINE pbe_lda_info 158 159! ************************************************************************************************** 160!> \brief return various information on the functional 161!> \param pbe_params section selecting the varius parameters for the functional 162!> \param reference string with the reference of the actual functional 163!> \param shortform string with the shortform of the functional name 164!> \param needs the components needed by this functional are set to 165!> true (does not set the unneeded components to false) 166!> \param max_deriv ... 167!> \author fawzi 168! ************************************************************************************************** 169 SUBROUTINE pbe_lsd_info(pbe_params, reference, shortform, needs, max_deriv) 170 TYPE(section_vals_type), POINTER :: pbe_params 171 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 172 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs 173 INTEGER, INTENT(out), OPTIONAL :: max_deriv 174 175 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lsd_info', routineP = moduleN//':'//routineN 176 177 INTEGER :: param 178 REAL(kind=dp) :: sc, sx 179 180 CALL section_vals_val_get(pbe_params, "parametrization", i_val=param) 181 CALL section_vals_val_get(pbe_params, "scale_x", r_val=sx) 182 CALL section_vals_val_get(pbe_params, "scale_c", r_val=sc) 183 184 SELECT CASE (param) 185 CASE (xc_pbe_orig) 186 CALL cite_reference(Perdew1996) 187 IF (sx == 1._dp .AND. sc == 1._dp) THEN 188 IF (PRESENT(reference)) THEN 189 reference = "J.P.Perdew, K.Burke, M.Ernzerhof, "// & 190 "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)"// & 191 "{spin polarized}" 192 END IF 193 IF (PRESENT(shortform)) THEN 194 shortform = "PBE xc functional (spin polarized)" 195 END IF 196 ELSE 197 IF (PRESENT(reference)) THEN 198 WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") & 199 "J.P.Perdew, K.Burke, M.Ernzerhof, ", & 200 "Phys. Rev. Letter, vol. 77, n 18, pp. 3865-3868, (1996)", & 201 sx, sc, "{spin polarized}" 202 END IF 203 IF (PRESENT(shortform)) THEN 204 WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')") & 205 "PBE xc functional (spin polarized)", sx, sc 206 END IF 207 END IF 208 CASE (xc_pbe_rev) 209 CALL cite_reference(Perdew1996) 210 CALL cite_reference(Zhang1998) 211 IF (sx == 1._dp .AND. sc == 1._dp) THEN 212 IF (PRESENT(reference)) THEN 213 reference = "revPBE, Yingkay Zhang and Weitao Yang,"// & 214 " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)"// & 215 "{spin polarized}" 216 END IF 217 IF (PRESENT(shortform)) THEN 218 shortform = "revPBE, revised PBE xc functional (spin polarized)" 219 END IF 220 ELSE 221 IF (PRESENT(reference)) THEN 222 WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") & 223 "revPBE, Yingkay Zhang and Weitao Yang,", & 224 " Phys. Rev. Letter, vol 80,n 4, p. 890, (1998)", & 225 sx, sc, "{spin polarized}" 226 END IF 227 IF (PRESENT(shortform)) THEN 228 WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')") & 229 "revPBE, revised PBE xc functional (spin polarized)", sx, sc 230 END IF 231 END IF 232 CASE (xc_pbe_sol) 233 CALL cite_reference(Perdew1996) 234 CALL cite_reference(Perdew2008) 235 IF (sx == 1._dp .AND. sc == 1._dp) THEN 236 IF (PRESENT(reference)) THEN 237 reference = "PBEsol, J.P. Perdew et al., "// & 238 "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) "// & 239 "{spin polarized}" 240 END IF 241 IF (PRESENT(shortform)) THEN 242 shortform = "PBEsol, PBE for solids and surfaces xc functional (spin polarized)" 243 END IF 244 ELSE 245 IF (PRESENT(reference)) THEN 246 WRITE (reference, "(a,a,'sx=',f5.3,'sc=',f5.3,a)") & 247 "PBEsol, J.P. Perdew et al., ", & 248 "Phys. Rev. Letter, vol 100,n 13, p. 136406, (2008) ", & 249 sx, sc, "{spin unpolarized}" 250 END IF 251 IF (PRESENT(shortform)) THEN 252 WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,'(LSD)')") & 253 "PBEsol, PBE for solids and surfaces xc functional (spin polarized)", sx, sc 254 END IF 255 END IF 256 CASE default 257 CPABORT("") 258 END SELECT 259 IF (PRESENT(needs)) THEN 260 needs%rho_spin = .TRUE. 261 needs%norm_drho_spin = .TRUE. 262 needs%norm_drho = .TRUE. 263 END IF 264 IF (PRESENT(max_deriv)) max_deriv = 2 265 266 END SUBROUTINE pbe_lsd_info 267 268! ************************************************************************************************** 269!> \brief evaluates the pbe correlation functional for lda 270!> \param rho_set the density where you want to evaluate the functional 271!> \param deriv_set place where to store the functional derivatives (they are 272!> added to the derivatives) 273!> \param grad_deriv degree of the derivative that should be evalated, 274!> if positive all the derivatives up to the given degree are evaluated, 275!> if negative only the given degree is calculated 276!> \param pbe_params ... 277!> \author fawzi 278! ************************************************************************************************** 279 SUBROUTINE pbe_lda_eval(rho_set, deriv_set, grad_deriv, pbe_params) 280 TYPE(xc_rho_set_type), POINTER :: rho_set 281 TYPE(xc_derivative_set_type), POINTER :: deriv_set 282 INTEGER, INTENT(in) :: grad_deriv 283 TYPE(section_vals_type), POINTER :: pbe_params 284 285 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lda_eval', routineP = moduleN//':'//routineN 286 287 INTEGER :: handle, npoints, param 288 INTEGER, DIMENSION(:, :), POINTER :: bo 289 REAL(kind=dp) :: epsilon_rho, scale_ec, scale_ex 290 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, e_ndrho_ndrho, & 291 e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, & 292 e_rho_rho_rho, norm_drho, rho 293 TYPE(xc_derivative_type), POINTER :: deriv 294 295 CALL timeset(routineN, handle) 296 297 NULLIFY (bo) 298 299 CPASSERT(ASSOCIATED(rho_set)) 300 CPASSERT(rho_set%ref_count > 0) 301 CPASSERT(ASSOCIATED(deriv_set)) 302 CPASSERT(deriv_set%ref_count > 0) 303 CALL xc_rho_set_get(rho_set, rho=rho, & 304 norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho) 305 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) 306 307 dummy => rho 308 309 e_0 => dummy 310 e_rho => dummy 311 e_ndrho => dummy 312 e_rho_rho => dummy 313 e_ndrho_rho => dummy 314 e_ndrho_ndrho => dummy 315 e_rho_rho_rho => dummy 316 e_ndrho_rho_rho => dummy 317 e_ndrho_ndrho_rho => dummy 318 e_ndrho_ndrho_ndrho => dummy 319 320 IF (grad_deriv >= 0) THEN 321 deriv => xc_dset_get_derivative(deriv_set, "", & 322 allocate_deriv=.TRUE.) 323 CALL xc_derivative_get(deriv, deriv_data=e_0) 324 END IF 325 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 326 deriv => xc_dset_get_derivative(deriv_set, "(rho)", & 327 allocate_deriv=.TRUE.) 328 CALL xc_derivative_get(deriv, deriv_data=e_rho) 329 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", & 330 allocate_deriv=.TRUE.) 331 CALL xc_derivative_get(deriv, deriv_data=e_ndrho) 332 END IF 333 IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN 334 deriv => xc_dset_get_derivative(deriv_set, "(rho)(rho)", & 335 allocate_deriv=.TRUE.) 336 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho) 337 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)(rho)", & 338 allocate_deriv=.TRUE.) 339 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho) 340 deriv => xc_dset_get_derivative(deriv_set, & 341 "(norm_drho)(norm_drho)", allocate_deriv=.TRUE.) 342 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho) 343 END IF 344 IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN 345 deriv => xc_dset_get_derivative(deriv_set, "(rho)(rho)(rho)", & 346 allocate_deriv=.TRUE.) 347 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho) 348 deriv => xc_dset_get_derivative(deriv_set, & 349 "(norm_drho)(rho)(rho)", allocate_deriv=.TRUE.) 350 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho) 351 deriv => xc_dset_get_derivative(deriv_set, & 352 "(norm_drho)(norm_drho)(rho)", allocate_deriv=.TRUE.) 353 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho) 354 deriv => xc_dset_get_derivative(deriv_set, & 355 "(norm_drho)(norm_drho)(norm_drho)", allocate_deriv=.TRUE.) 356 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho) 357 END IF 358 IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN 359 CPABORT("derivatives bigger than 3 not implemented") 360 END IF 361 362 CALL section_vals_val_get(pbe_params, "scale_c", r_val=scale_ec) 363 CALL section_vals_val_get(pbe_params, "scale_x", r_val=scale_ex) 364 CALL section_vals_val_get(pbe_params, "parametrization", i_val=param) 365 366!$OMP PARALLEL DEFAULT(NONE), & 367!$OMP SHARED(rho,norm_drho,e_0,e_rho,e_ndrho,e_rho_rho,e_ndrho_rho),& 368!$OMP SHARED(e_ndrho_ndrho,e_rho_rho_rho,e_ndrho_rho_rho,e_ndrho_ndrho_rho),& 369!$OMP SHARED(e_ndrho_ndrho_ndrho,grad_deriv,npoints,epsilon_rho),& 370!$OMP SHARED(pbe_params),& 371!$OMP SHARED(param,scale_ec,scale_ex) 372 CALL pbe_lda_calc(rho=rho, norm_drho=norm_drho, & 373 e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, & 374 e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, & 375 e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, & 376 e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, & 377 grad_deriv=grad_deriv, & 378 npoints=npoints, epsilon_rho=epsilon_rho, & 379 param=param, scale_ec=scale_ec, scale_ex=scale_ex) 380!$OMP END PARALLEL 381 382 CALL timestop(handle) 383 END SUBROUTINE pbe_lda_eval 384 385! ************************************************************************************************** 386!> \brief evaluates the pbe correlation functional for lda 387!> \param rho the density where you want to evaluate the functional 388!> \param norm_drho ... 389!> \param e_0 ... 390!> \param e_rho ... 391!> \param e_ndrho ... 392!> \param e_rho_rho ... 393!> \param e_ndrho_rho ... 394!> \param e_ndrho_ndrho ... 395!> \param e_rho_rho_rho ... 396!> \param e_ndrho_rho_rho ... 397!> \param e_ndrho_ndrho_rho ... 398!> \param e_ndrho_ndrho_ndrho ... 399!> \param grad_deriv degree of the derivative that should be evalated, 400!> if positive all the derivatives up to the given degree are evaluated, 401!> if negative only the given degree is calculated 402!> \param npoints ... 403!> \param epsilon_rho ... 404!> \param ! ... 405!> \param pbe_params parameters for the pbe functional 406!> \author fawzi 407! ************************************************************************************************** 408 SUBROUTINE pbe_lda_calc(rho, norm_drho, & 409 e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, & 410 e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, & 411 e_ndrho_ndrho_ndrho, grad_deriv, npoints, epsilon_rho, & 412 ! pbe_params) 413 param, scale_ec, scale_ex) 414 INTEGER, INTENT(in) :: npoints, grad_deriv 415 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_ndrho, & 416 e_ndrho_ndrho_rho, e_ndrho_rho_rho, e_rho_rho_rho, e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, & 417 e_ndrho, e_rho, e_0 418 REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho 419 REAL(kind=dp), INTENT(in) :: epsilon_rho 420 INTEGER, INTENT(in) :: param 421 REAL(kind=dp), INTENT(in) :: scale_ec, scale_ex 422 423 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lda_calc', routineP = moduleN//':'//routineN 424 425 INTEGER :: ii 426 REAL(kind=dp) :: A, A1rho, A1rhorho, A2rho, A_1, alpha_1_1, Arho, Arho1rho, Arhorho, & 427 Arhorhorho, beta, beta_1_1, beta_2_1, beta_3_1, beta_4_1, e_c_u_0, e_c_u_01rho, & 428 e_c_u_01rhorho, e_c_u_02rho, e_c_u_0rho, e_c_u_0rho1rho, e_c_u_0rhorho, e_c_u_0rhorhorho, & 429 epsilon_cGGA, epsilon_cGGArho, epsilon_cGGArhorho, ex_ldarhorhorho, ex_unif, ex_unif1rho, & 430 ex_unif1rhorho, ex_unif2rho, ex_unifrho, ex_unifrho1rho, ex_unifrhorho, Fx, Fx1rho, & 431 Fx1rhorho, Fx2rho, Fxnorm_drho, Fxnorm_drho1rho, Fxnorm_drhonorm_drho, Fxnorm_drhorho, & 432 Fxrho, Fxrho1rho, Fxrhorho, gamma_var, Hnorm_drho, Hnorm_drhonorm_drho 433 REAL(kind=dp) :: Hnorm_drhorho, k_f, k_f2rho, k_frho, k_frhorho, k_frhorhorho, k_s, k_s1rho, & 434 k_s1rhorho, k_s2rho, k_srho, k_srho1rho, k_srhorho, kappa, kf, kf2rho, kfrho, kfrhorho, & 435 kfrhorhorho, mu, my_norm_drho, my_rho, p_1, p_2, rs, rs2rho, rsrho, rsrhorho, & 436 rsrhorhorho, s, s1rho, s1rhorho, s2norm_drho, s2rho, snorm_drho, snorm_drho1rho, & 437 snorm_drhorho, srho, srho1rho, srhorho, t, t1, t1001, t1004, t1005, t1006, t1008, t101, & 438 t1011, t1012, t1013, t1014, t1016, t1017, t1018, t1019, t1022, t1024, t1026, t1028, t103, & 439 t1030, t1031, t1035, t1042, t1046, t1048, t105, t1050, t1052, t1054, t1055 440 REAL(kind=dp) :: t1060, t1061, t1062, t1067, t108, t1093, t1097, t11, t1103, t1104, t1106, & 441 t1109, t111, t1115, t1118, t1121, t1122, t1126, t1129, t113, t114, t1148, t115, t1152, & 442 t1157, t1167, t1187, t119, t1196, t1203, t1218, t1226, t123, t124, t1249, t125, t126, & 443 t1262, t1263, t127, t1291, t1292, t1295, t13, t131, t1329, t1342, t135, t138, t1380, & 444 t1385, t1386, t1387, t1389, t139, t1393, t14, t140, t142, t146, t1468, t148, t1482, t149, & 445 t1492, t1493, t150, t1504, t1505, t151, t1511, t1515, t1521, t1525, t1528, t153, t1532, & 446 t1545, t155, t1565, t157, t1573, t158, t1584, t159, t1608, t1612 447 REAL(kind=dp) :: t162, t1628, t1629, t163, t1633, t164, t1646, t1652, t1660, t167, t1672, & 448 t1676, t168, t17, t170, t171, t1715, t1722, t175, t1758, t1759, t176, t1761, t177, t178, & 449 t179, t1797, t182, t183, t1838, t1851, t186, t187, t1878, t1889, t189, t19, t190, t1907, & 450 t191, t1922, t1927, t1933, t1937, t195, t1952, t196, t1964, t1965, t1968, t1969, t197, & 451 t1972, t198, t1990, t1rho, t1rhorho, t2, t20, t200, t202, t2020, t2024, t2028, t2031, & 452 t204, t2041, t208, t21, t214, t217, t218, t22, t226, t229, t230, t238, t239, t240, t241, & 453 t246, t252, t253, t255, t259, t26, t260, t261, t262, t265, t266 454 REAL(kind=dp) :: t267, t27, t271, t272, t273, t277, t278, t280, t281, t286, t290, t291, & 455 t293, t294, t295, t296, t297, t299, t2norm_drho, t2rho, t3, t305, t309, t310, t315, t317, & 456 t318, t321, t323, t324, t327, t329, t330, t331, t336, t338, t339, t340, t341, t348, t349, & 457 t354, t357, t359, t360, t361, t362, t369, t370, t371, t374, t377, t378, t381, t382, t384, & 458 t385, t387, t388, t390, t392, t393, t397, t4, t400, t401, t404, t408, t409, t410, t414, & 459 t417, t418, t423, t426, t427, t432, t435, t436, t438, t439, t440, t448, t449, t451, t456, & 460 t457, t458, t459, t461, t462, t463, t465, t466, t469, t470 461 REAL(kind=dp) :: t471, t472, t476, t487, t491, t496, t498, t5, t500, t503, t506, t507, t510, & 462 t513, t517, t521, t525, t529, t530, t535, t541, t548, t553, t556, t557, t559, t562, t566, & 463 t581, t586, t590, t591, t594, t598, t6, t605, t609, t611, t612, t614, t627, t645, t65, & 464 t654, t656, t661, t664, t669, t670, t671, t673, t675, t685, t69, t693, t7, t70, t701, & 465 t702, t71, t714, t717, t72, t720, t73, t74, t740, t743, t748, t75, t758, t77, t776, t78, & 466 t8, t80, t801, t809, t81, t812, t821, t825, t83, t831, t84, t85, t86, t868, t87, t877, & 467 t879, t88, t880, t885, t90, t91, t94, t940, t942, t943, t945 468 REAL(kind=dp) :: t946, t948, t95, t950, t951, t954, t967, t976, t98, t980, t982, t984, t989, & 469 t99, t990, t994, t995, t998, t999, tnorm_drho, tnorm_drho1rho, tnorm_drhorho, & 470 tnorm_drhorhorho, trho, trho1rho, trhorho, trhorhorho 471 472!TYPE(section_vals_type), POINTER :: pbe_params 473!, param 474! scale_ec, scale_ex, snorm_drho, snorm_drho1rho, snorm_drhorho, srho, & 475! Parametrization of PBE 476 477 SELECT CASE (param) 478 ! Original PBE 479 CASE (xc_pbe_orig) 480 kappa = 0.804e0_dp 481 beta = 0.66725e-1_dp 482 mu = beta*pi**2/0.3e1_dp 483 ! Revised PBE (revPBE) 484 CASE (xc_pbe_rev) 485 kappa = 0.1245e1_dp 486 beta = 0.66725e-1_dp 487 mu = beta*pi**2/0.3e1_dp 488 ! PBE for solids (PBEsol) 489 CASE (xc_pbe_sol) 490 kappa = 0.804e0_dp 491 beta = 0.46e-1_dp 492 mu = 0.1e2_dp/0.81e2_dp 493 CASE default 494 CPABORT("") 495 END SELECT 496 497 gamma_var = (0.1e1_dp - LOG(0.2e1_dp))/pi**2 498 p_1 = 0.10e1_dp 499 A_1 = 0.31091e-1_dp 500 alpha_1_1 = 0.21370e0_dp 501 beta_1_1 = 0.75957e1_dp 502 beta_2_1 = 0.35876e1_dp 503 beta_3_1 = 0.16382e1_dp 504 beta_4_1 = 0.49294e0_dp 505 p_2 = 0.10e1_dp 506 t1 = 3**(0.1e1_dp/0.3e1_dp) 507 t2 = 4**(0.1e1_dp/0.3e1_dp) 508 t3 = t2**2 509 t4 = t1*t3 510 t5 = 0.1e1_dp/pi 511 t69 = pi**2 512 t627 = 0.1e1_dp/t69/pi 513 514 SELECT CASE (grad_deriv) 515 CASE default 516!$OMP DO 517 DO ii = 1, npoints 518 my_rho = rho(ii) 519 IF (my_rho > epsilon_rho) THEN 520 my_norm_drho = norm_drho(ii) 521 522 t6 = 0.1e1_dp/my_rho 523 t7 = t5*t6 524 t8 = t7**(0.1e1_dp/0.3e1_dp) 525 rs = t4*t8/0.4e1_dp 526 t11 = 0.1e1_dp + alpha_1_1*rs 527 t13 = 0.1e1_dp/A_1 528 t14 = SQRT(rs) 529 t17 = t14*rs 530 t19 = p_1 + 0.1e1_dp 531 t20 = rs**t19 532 t21 = beta_4_1*t20 533 t22 = beta_1_1*t14 + beta_2_1*rs + beta_3_1*t17 + t21 534 t26 = 0.1e1_dp + t13/t22/0.2e1_dp 535 t27 = LOG(t26) 536 e_c_u_0 = -0.2e1_dp*A_1*t11*t27 537 t65 = 2**(0.1e1_dp/0.3e1_dp) 538 t70 = t69*my_rho 539 t71 = t70**(0.1e1_dp/0.3e1_dp) 540 k_f = t1*t71 541 t72 = k_f*t5 542 t73 = SQRT(t72) 543 k_s = 0.2e1_dp*t73 544 t74 = 0.1e1_dp/k_s 545 t75 = my_norm_drho*t74 546 t = t75*t6/0.2e1_dp 547 t77 = 0.1e1_dp/gamma_var 548 t78 = beta*t77 549 t80 = EXP(-e_c_u_0*t77) 550 t81 = -0.1e1_dp + t80 551 A = t78/t81 552 t83 = t**2 553 t84 = A*t83 554 t85 = 0.1e1_dp + t84 555 t86 = t83*t85 556 t87 = A**2 557 t88 = t83**2 558 t90 = 0.1e1_dp + t84 + t87*t88 559 t91 = 0.1e1_dp/t90 560 t94 = 0.1e1_dp + t78*t86*t91 561 t95 = LOG(t94) 562 epsilon_cGGA = e_c_u_0 + gamma_var*t95 563 kf = k_f 564 ex_unif = -0.3e1_dp/0.4e1_dp*t5*kf 565 t98 = 0.1e1_dp/kf 566 t99 = my_norm_drho*t98 567 s = t99*t6/0.2e1_dp 568 t101 = s**2 569 t103 = 0.1e1_dp/kappa 570 t105 = 0.1e1_dp + mu*t101*t103 571 Fx = 0.1e1_dp + kappa - kappa/t105 572 t108 = my_rho*ex_unif 573 574 IF (grad_deriv >= 0) THEN 575 e_0(ii) = e_0(ii) + & 576 scale_ex*t108*Fx + scale_ec*my_rho*epsilon_cGGA 577 END IF 578 579 t111 = t8**2 580 t113 = 0.1e1_dp/t111*t5 581 t114 = my_rho**2 582 t115 = 0.1e1_dp/t114 583 rsrho = -t4*t113*t115/0.12e2_dp 584 t119 = A_1*alpha_1_1 585 t123 = t22**2 586 t124 = 0.1e1_dp/t123 587 t125 = t11*t124 588 t126 = 0.1e1_dp/t14 589 t127 = beta_1_1*t126 590 t131 = beta_3_1*t14 591 t135 = 0.1e1_dp/rs 592 t138 = t127*rsrho/0.2e1_dp + beta_2_1*rsrho + 0.3e1_dp/ & 593 0.2e1_dp*t131*rsrho + t21*t19*rsrho*t135 594 t139 = 0.1e1_dp/t26 595 t140 = t138*t139 596 e_c_u_0rho = -0.2e1_dp*t119*rsrho*t27 + t125*t140 597 t142 = t71**2 598 k_frho = t1/t142*t69/0.3e1_dp 599 t146 = 0.1e1_dp/t73 600 k_srho = t146*k_frho*t5 601 t148 = k_s**2 602 t149 = 0.1e1_dp/t148 603 t150 = my_norm_drho*t149 604 t151 = t6*k_srho 605 t153 = t75*t115 606 trho = -t150*t151/0.2e1_dp - t153/0.2e1_dp 607 t155 = gamma_var**2 608 t157 = beta/t155 609 t158 = t81**2 610 t159 = 0.1e1_dp/t158 611 Arho = t157*t159*e_c_u_0rho*t80 612 t162 = t78*t 613 t163 = t85*t91 614 t164 = t163*trho 615 t167 = Arho*t83 616 t168 = A*t 617 t170 = 0.2e1_dp*t168*trho 618 t171 = t167 + t170 619 t175 = t78*t83 620 t176 = t90**2 621 t177 = 0.1e1_dp/t176 622 t178 = t85*t177 623 t179 = A*t88 624 t182 = t83*t 625 t183 = t87*t182 626 t186 = t167 + t170 + 0.2e1_dp*t179*Arho + 0.4e1_dp*t183*trho 627 t187 = t178*t186 628 t189 = 0.2e1_dp*t162*t164 + t78*t83*t171*t91 - t175*t187 629 t190 = gamma_var*t189 630 t191 = 0.1e1_dp/t94 631 epsilon_cGGArho = e_c_u_0rho + t190*t191 632 kfrho = k_frho 633 ex_unifrho = -0.3e1_dp/0.4e1_dp*t5*kfrho 634 t195 = kf**2 635 t196 = 0.1e1_dp/t195 636 t197 = my_norm_drho*t196 637 t198 = t6*kfrho 638 t200 = t99*t115 639 srho = -t197*t198/0.2e1_dp - t200/0.2e1_dp 640 t202 = t105**2 641 t204 = 0.1e1_dp/t202*mu 642 Fxrho = 0.2e1_dp*t204*s*srho 643 t208 = my_rho*ex_unifrho 644 645 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 646 e_rho(ii) = e_rho(ii) + & 647 scale_ex*(ex_unif*Fx + t208*Fx + t108*Fxrho) + & 648 scale_ec*(epsilon_cGGA + my_rho*epsilon_cGGArho) 649 END IF 650 651 tnorm_drho = t74*t6/0.2e1_dp 652 t214 = t163*tnorm_drho 653 t217 = t78*t182 654 t218 = A*tnorm_drho 655 t226 = 0.2e1_dp*t168*tnorm_drho + 0.4e1_dp*t183*tnorm_drho 656 t229 = 0.2e1_dp*t162*t214 + 0.2e1_dp*t217*t218*t91 - & 657 t175*t178*t226 658 t230 = gamma_var*t229 659 Hnorm_drho = t230*t191 660 snorm_drho = t98*t6/0.2e1_dp 661 Fxnorm_drho = 0.2e1_dp*t204*s*snorm_drho 662 663 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 664 e_ndrho(ii) = e_ndrho(ii) + & 665 scale_ex*t108*Fxnorm_drho + scale_ec*my_rho* & 666 Hnorm_drho 667 END IF 668 669 t238 = 0.1e1_dp/t69 670 t239 = 0.1e1_dp/t111/t7*t238 671 t240 = t114**2 672 t241 = 0.1e1_dp/t240 673 t246 = 0.1e1_dp/t114/my_rho 674 rsrhorho = -t4*t239*t241/0.18e2_dp + t4*t113* & 675 t246/0.6e1_dp 676 t252 = 0.2e1_dp*t119*rsrhorho*t27 677 t253 = alpha_1_1*rsrho 678 t255 = t124*t138*t139 679 t259 = 0.1e1_dp/t123/t22 680 t260 = t11*t259 681 t261 = t138**2 682 t262 = t261*t139 683 t265 = 0.1e1_dp/t17 684 t266 = beta_1_1*t265 685 t267 = rsrho**2 686 t271 = t127*rsrhorho/0.2e1_dp 687 t272 = beta_2_1*rsrhorho 688 t273 = beta_3_1*t126 689 t277 = 0.3e1_dp/0.2e1_dp*t131*rsrhorho 690 t278 = t19**2 691 t280 = rs**2 692 t281 = 0.1e1_dp/t280 693 t286 = t21*t19*rsrhorho*t135 694 t290 = -t266*t267/0.4e1_dp + t271 + t272 + 0.3e1_dp/0.4e1_dp & 695 *t273*t267 + t277 + t21*t278*t267*t281 + t286 - t21*t19 & 696 *t267*t281 697 t291 = t290*t139 698 t293 = t123**2 699 t294 = 0.1e1_dp/t293 700 t295 = t11*t294 701 t296 = t26**2 702 t297 = 0.1e1_dp/t296 703 t299 = t261*t297*t13 704 e_c_u_0rhorho = -t252 + 0.2e1_dp*t253*t255 - 0.2e1_dp*t260* & 705 t262 + t125*t291 + t295*t299/0.2e1_dp 706 e_c_u_01rho = e_c_u_0rho 707 t305 = t69**2 708 k_frhorho = -0.2e1_dp/0.9e1_dp*t1/t142/t70*t305 709 t309 = 0.1e1_dp/t73/t72 710 t310 = k_frho**2 711 t315 = t146*k_frhorho*t5 712 k_srhorho = -t309*t310*t238/0.2e1_dp + t315 713 k_s1rho = k_srho 714 t317 = 0.1e1_dp/t148/k_s 715 t318 = my_norm_drho*t317 716 t321 = t115*k_srho 717 t323 = t150*t321/0.2e1_dp 718 t324 = t6*k_srhorho 719 t327 = t115*k_s1rho 720 t329 = t150*t327/0.2e1_dp 721 t330 = t75*t246 722 trhorho = t318*t151*k_s1rho + t323 - t150*t324/0.2e1_dp + & 723 t329 + t330 724 t331 = t6*k_s1rho 725 t1rho = -t150*t331/0.2e1_dp - t153/0.2e1_dp 726 t336 = beta/t155/gamma_var 727 t338 = 0.1e1_dp/t158/t81 728 t339 = t336*t338 729 t340 = t80**2 730 t341 = e_c_u_0rho*t340 731 t348 = t336*t159 732 t349 = e_c_u_0rho*e_c_u_01rho 733 Arhorho = 0.2e1_dp*t339*t341*e_c_u_01rho + t157*t159* & 734 e_c_u_0rhorho*t80 - t348*t349*t80 735 A1rho = t157*t159*e_c_u_01rho*t80 736 t354 = t78*t1rho 737 t357 = A1rho*t83 738 t359 = 0.2e1_dp*t168*t1rho 739 t360 = t357 + t359 740 t361 = t360*t91 741 t362 = t361*trho 742 t369 = t357 + t359 + 0.2e1_dp*t179*A1rho + 0.4e1_dp*t183*t1rho 743 t370 = trho*t369 744 t371 = t178*t370 745 t374 = t163*trhorho 746 t377 = t171*t91 747 t378 = t377*t1rho 748 t381 = Arhorho*t83 749 t382 = Arho*t 750 t384 = 0.2e1_dp*t382*t1rho 751 t385 = A1rho*t 752 t387 = 0.2e1_dp*t385*trho 753 t388 = A*t1rho 754 t390 = 0.2e1_dp*t388*trho 755 t392 = 0.2e1_dp*t168*trhorho 756 t393 = t381 + t384 + t387 + t390 + t392 757 t397 = t171*t177 758 t400 = t186*t1rho 759 t401 = t178*t400 760 t404 = t360*t177 761 t408 = 0.1e1_dp/t176/t90 762 t409 = t85*t408 763 t410 = t186*t369 764 t414 = A1rho*t88 765 t417 = A*t182 766 t418 = Arho*t1rho 767 t423 = trho*A1rho 768 t426 = t87*t83 769 t427 = trho*t1rho 770 t432 = t381 + t384 + t387 + t390 + t392 + 0.2e1_dp*t414*Arho + & 771 0.8e1_dp*t417*t418 + 0.2e1_dp*t179*Arhorho + 0.8e1_dp* & 772 t417*t423 + 0.12e2_dp*t426*t427 + 0.4e1_dp*t183*trhorho 773 t435 = 0.2e1_dp*t354*t164 + 0.2e1_dp*t162*t362 - 0.2e1_dp & 774 *t162*t371 + 0.2e1_dp*t162*t374 + 0.2e1_dp*t162*t378 + & 775 t78*t83*t393*t91 - t175*t397*t369 - 0.2e1_dp*t162*t401 & 776 - t175*t404*t186 + 0.2e1_dp*t175*t409*t410 - t175*t178 & 777 *t432 778 t436 = gamma_var*t435 779 t438 = t94**2 780 t439 = 0.1e1_dp/t438 781 t440 = t163*t1rho 782 t448 = 0.2e1_dp*t162*t440 + t78*t83*t360*t91 - t175* & 783 t178*t369 784 t449 = t439*t448 785 t451 = gamma_var*t448 786 epsilon_cGGArhorho = e_c_u_0rhorho + t436*t191 - t190*t449 787 kfrhorho = k_frhorho 788 ex_unifrhorho = -0.3e1_dp/0.4e1_dp*t5*kfrhorho 789 ex_unif1rho = ex_unifrho 790 t456 = 0.1e1_dp/t195/kf 791 t457 = my_norm_drho*t456 792 t458 = kfrho**2 793 t459 = t6*t458 794 t461 = t115*kfrho 795 t462 = t197*t461 796 t463 = t6*kfrhorho 797 t465 = t197*t463/0.2e1_dp 798 t466 = t99*t246 799 srhorho = t457*t459 + t462 - t465 + t466 800 s1rho = srho 801 t469 = mu**2 802 t470 = 0.1e1_dp/t202/t105*t469 803 t471 = t470*t101 804 t472 = srho*t103 805 t476 = s1rho*srho 806 Fxrhorho = -0.8e1_dp*t471*t472*s1rho + 0.2e1_dp*t204* & 807 t476 + 0.2e1_dp*t204*s*srhorho 808 Fx1rho = 0.2e1_dp*t204*s*s1rho 809 t487 = my_rho*ex_unifrhorho 810 t491 = my_rho*ex_unif1rho 811 812 IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN 813 e_rho_rho(ii) = e_rho_rho(ii) + & 814 scale_ex*(ex_unif1rho*Fx + ex_unif*Fx1rho + & 815 ex_unifrho*Fx + t487*Fx + t208*Fx1rho + ex_unif*Fxrho + t491 & 816 *Fxrho + t108*Fxrhorho) + scale_ec*(e_c_u_01rho + t451*t191 & 817 + epsilon_cGGArho + my_rho*epsilon_cGGArhorho) 818 END IF 819 820 t496 = t149*t6 821 t498 = t74*t115 822 tnorm_drhorho = -t496*k_srho/0.2e1_dp - t498/0.2e1_dp 823 t500 = t78*trho 824 t503 = t377*tnorm_drho 825 t506 = tnorm_drho*t186 826 t507 = t178*t506 827 t510 = t163*tnorm_drhorho 828 t513 = t91*trho 829 t517 = Arho*tnorm_drho 830 t521 = A*tnorm_drhorho 831 t525 = t177*t186 832 t529 = t226*trho 833 t530 = t178*t529 834 t535 = t226*t186 835 t541 = A*trho 836 t548 = tnorm_drho*trho 837 t553 = 0.2e1_dp*t382*tnorm_drho + 0.2e1_dp*t541*tnorm_drho & 838 + 0.2e1_dp*t168*tnorm_drhorho + 0.8e1_dp*t417*t517 + & 839 0.12e2_dp*t426*t548 + 0.4e1_dp*t183*tnorm_drhorho 840 t556 = 0.2e1_dp*t500*t214 + 0.2e1_dp*t162*t503 - 0.2e1_dp & 841 *t162*t507 + 0.2e1_dp*t162*t510 + 0.6e1_dp*t175*t218* & 842 t513 + 0.2e1_dp*t217*t517*t91 + 0.2e1_dp*t217*t521*t91 - & 843 0.2e1_dp*t217*t218*t525 - 0.2e1_dp*t162*t530 - t175* & 844 t397*t226 + 0.2e1_dp*t175*t409*t535 - t175*t178*t553 845 t557 = gamma_var*t556 846 t559 = t439*t189 847 Hnorm_drhorho = t557*t191 - t230*t559 848 t562 = t196*t6 849 snorm_drhorho = -t562*kfrho/0.2e1_dp - t98*t115/0.2e1_dp 850 t566 = snorm_drho*t103 851 Fxnorm_drhorho = -0.8e1_dp*t471*t566*srho + 0.2e1_dp*t204 & 852 *srho*snorm_drho + 0.2e1_dp*t204*s*snorm_drhorho 853 854 IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN 855 e_ndrho_rho(ii) = e_ndrho_rho(ii) + & 856 scale_ex*(ex_unif*Fxnorm_drho + t208* & 857 Fxnorm_drho + t108*Fxnorm_drhorho) + scale_ec*(Hnorm_drho + my_rho & 858 *Hnorm_drhorho) 859 END IF 860 861 t581 = tnorm_drho**2 862 t586 = A*t581 863 t590 = tnorm_drho*t226 864 t591 = t178*t590 865 t594 = t177*t226 866 t598 = t226**2 867 t605 = 0.2e1_dp*t586 + 0.12e2_dp*t426*t581 868 t609 = gamma_var*(0.2e1_dp*t78*t581*t85*t91 + 0.10e2_dp & 869 *t175*t586*t91 - 0.4e1_dp*t162*t591 - 0.4e1_dp*t217* & 870 t218*t594 + 0.2e1_dp*t175*t409*t598 - t175*t178*t605) 871 t611 = t229**2 872 t612 = gamma_var*t611 873 Hnorm_drhonorm_drho = t609*t191 - t612*t439 874 t614 = snorm_drho**2 875 Fxnorm_drhonorm_drho = -0.8e1_dp*t470*t101*t614*t103 + & 876 0.2e1_dp*t204*t614 877 878 IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN 879 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + & 880 scale_ex*t108*Fxnorm_drhonorm_drho + & 881 scale_ec*my_rho*Hnorm_drhonorm_drho 882 END IF 883 884 IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN 885 rsrhorhorho = -0.5e1_dp/0.54e2_dp*t4/t111/t238/ & 886 t115*t627/t240/t114 + t4*t239/t240/my_rho/0.3e1_dp & 887 - t4*t113*t241/0.2e1_dp 888 rs2rho = rsrho 889 t645 = alpha_1_1*rsrhorho 890 t654 = t127*rs2rho/0.2e1_dp + beta_2_1*rs2rho + 0.3e1_dp/ & 891 0.2e1_dp*t131*rs2rho + t21*t19*rs2rho*t135 892 t656 = t124*t654*t139 893 t661 = t140*t654 894 t664 = rsrho*rs2rho 895 t669 = t21*t278 896 t670 = rs2rho*t281 897 t671 = t670*rsrho 898 t673 = t21*t19 899 t675 = -t266*t664/0.4e1_dp + t271 + t272 + 0.3e1_dp/0.4e1_dp & 900 *t273*t664 + t277 + t669*t671 + t286 - t673*t671 901 t685 = alpha_1_1*rs2rho 902 t693 = t675*t139 903 t701 = t297*t13 904 t702 = t701*t654 905 t714 = t267*rs2rho 906 t717 = rsrhorho*rsrho 907 t720 = rsrhorho*rs2rho 908 t740 = rs2rho/t280/rs*t267 909 t743 = rsrhorho*t281*rsrho 910 t748 = t670*rsrhorho 911 t758 = 0.3e1_dp/0.8e1_dp*beta_1_1/t14/t280*t714 - t266* & 912 t717/0.2e1_dp - t266*t720/0.4e1_dp + t127*rsrhorhorho/ & 913 0.2e1_dp + beta_2_1*rsrhorhorho - 0.3e1_dp/0.8e1_dp*beta_3_1* & 914 t265*t714 + 0.3e1_dp/0.2e1_dp*t273*t717 + 0.3e1_dp/ & 915 0.4e1_dp*t273*t720 + 0.3e1_dp/0.2e1_dp*t131*rsrhorhorho + & 916 t21*t278*t19*t740 + 0.2e1_dp*t669*t743 - 0.3e1_dp*t669* & 917 t740 + t669*t748 + t21*t19*rsrhorhorho*t135 - t673*t748 - & 918 0.2e1_dp*t673*t743 + 0.2e1_dp*t673*t740 919 t776 = A_1**2 920 e_c_u_0rhorhorho = -0.2e1_dp*t119*rsrhorhorho*t27 + t645* & 921 t656 + 0.2e1_dp*t645*t255 - 0.4e1_dp*t253*t259*t661 + & 922 0.2e1_dp*t253*t124*t675*t139 + t253*t294*t138*t297* & 923 t13*t654 - 0.2e1_dp*t685*t259*t261*t139 + 0.6e1_dp*t295 & 924 *t262*t654 - 0.4e1_dp*t260*t693*t138 - 0.3e1_dp*t11/ & 925 t293/t22*t261*t702 + t685*t124*t290*t139 - 0.2e1_dp* & 926 t260*t291*t654 + t125*t758*t139 + t295*t290*t702/ & 927 0.2e1_dp + t685*t294*t299/0.2e1_dp + t295*t675*t701*t138 & 928 + t11/t293/t123*t261/t296/t26/t776*t654/0.2e1_dp 929 e_c_u_0rho1rho = -t252 + t253*t656 + t685*t255 - 0.2e1_dp* & 930 t260*t661 + t125*t693 + t295*t138*t702/0.2e1_dp 931 e_c_u_01rhorho = e_c_u_0rho1rho 932 e_c_u_02rho = -0.2e1_dp*t119*rs2rho*t27 + t125*t654*t139 933 k_frhorhorho = 0.10e2_dp/0.27e2_dp*t1/t142/t114*t69 934 k_f2rho = kfrho 935 t801 = k_f**2 936 t809 = t309*k_frhorho 937 t812 = t238*k_f2rho 938 k_srho1rho = -t309*k_frho*t812/0.2e1_dp + t315 939 k_s1rhorho = k_srho1rho 940 k_s2rho = t146*k_f2rho*t5 941 t821 = t148**2 942 t825 = k_srho*k_s1rho 943 t831 = t6*k_srho1rho 944 trhorhorho = -0.3e1_dp*my_norm_drho/t821*t6*t825*k_s2rho - & 945 t318*t321*k_s1rho + t318*t831*k_s1rho + t318*t151* & 946 k_s1rhorho - t318*t321*k_s2rho - t150*t246*k_srho + t150* & 947 t115*k_srho1rho/0.2e1_dp + t318*t324*k_s2rho + t150*t115* & 948 k_srhorho/0.2e1_dp - t150*t6*(0.3e1_dp/0.4e1_dp/t73/ & 949 t801/t238*t310*t627*k_f2rho - t809*t238*k_frho - t809* & 950 t812/0.2e1_dp + t146*k_frhorhorho*t5)/0.2e1_dp - t318*t327 & 951 *k_s2rho - t150*t246*k_s1rho + t150*t115*k_s1rhorho/ & 952 0.2e1_dp - t150*t246*k_s2rho - 0.3e1_dp*t75*t241 953 t868 = t150*t115*k_s2rho/0.2e1_dp 954 trho1rho = t318*t151*k_s2rho + t323 - t150*t831/0.2e1_dp + & 955 t868 + t330 956 t1rhorho = t318*t331*k_s2rho + t329 - t150*t6*k_s1rhorho/ & 957 0.2e1_dp + t868 + t330 958 t2rho = -t150*t6*k_s2rho/0.2e1_dp - t153/0.2e1_dp 959 t877 = t155**2 960 t879 = beta/t877 961 t880 = t158**2 962 t885 = e_c_u_01rho*e_c_u_02rho 963 Arhorhorho = 0.6e1_dp*t879/t880*e_c_u_0rho*t340*t80* & 964 t885 + 0.2e1_dp*t339*e_c_u_0rho1rho*t340*e_c_u_01rho - & 965 0.6e1_dp*t879*t338*t341*t885 + 0.2e1_dp*t339*t341* & 966 e_c_u_01rhorho + 0.2e1_dp*t339*e_c_u_0rhorho*t340* & 967 e_c_u_02rho + t157*t159*e_c_u_0rhorhorho*t80 - t348* & 968 e_c_u_0rhorho*e_c_u_02rho*t80 - t348*e_c_u_0rho1rho* & 969 e_c_u_01rho*t80 - t348*e_c_u_0rho*e_c_u_01rhorho*t80 + t879 & 970 *t159*t349*e_c_u_02rho*t80 971 Arho1rho = 0.2e1_dp*t339*t341*e_c_u_02rho + t157*t159* & 972 e_c_u_0rho1rho*t80 - t348*e_c_u_0rho*e_c_u_02rho*t80 973 A1rhorho = 0.2e1_dp*t339*e_c_u_01rho*t340*e_c_u_02rho + & 974 t157*t159*e_c_u_01rhorho*t80 - t348*t885*t80 975 A2rho = t157*t159*e_c_u_02rho*t80 976 t940 = Arho1rho*t83 977 t942 = 0.2e1_dp*t382*t2rho 978 t943 = A2rho*t 979 t945 = 0.2e1_dp*t943*trho 980 t946 = A*t2rho 981 t948 = 0.2e1_dp*t946*trho 982 t950 = 0.2e1_dp*t168*trho1rho 983 t951 = A2rho*t88 984 t954 = Arho*t2rho 985 t967 = t940 + t942 + t945 + t948 + t950 + 0.2e1_dp*t951*Arho + & 986 0.8e1_dp*t417*t954 + 0.2e1_dp*t179*Arho1rho + 0.8e1_dp* & 987 t417*trho*A2rho + 0.12e2_dp*t426*trho*t2rho + 0.4e1_dp* & 988 t183*trho1rho 989 t976 = t78*t2rho 990 t980 = t78*t*t85 991 t982 = A2rho*t83 992 t984 = 0.2e1_dp*t168*t2rho 993 t989 = t982 + t984 + 0.2e1_dp*t179*A2rho + 0.4e1_dp*t183*t2rho 994 t990 = t369*t989 995 t994 = t982 + t984 996 t995 = t994*t177 997 t998 = Arhorhorho*t83 998 t999 = Arhorho*t 999 t1001 = 0.2e1_dp*t999*t2rho 1000 t1004 = 0.2e1_dp*Arho1rho*t*t1rho 1001 t1005 = t954*t1rho 1002 t1006 = 0.2e1_dp*t1005 1003 t1008 = 0.2e1_dp*t382*t1rhorho 1004 t1011 = 0.2e1_dp*A1rhorho*t*trho 1005 t1012 = A1rho*t2rho 1006 t1013 = t1012*trho 1007 t1014 = 0.2e1_dp*t1013 1008 t1016 = 0.2e1_dp*t385*trho1rho 1009 t1017 = A2rho*t1rho 1010 t1018 = t1017*trho 1011 t1019 = 0.2e1_dp*t1018 1012 t1022 = 0.2e1_dp*A*t1rhorho*trho 1013 t1024 = 0.2e1_dp*t388*trho1rho 1014 t1026 = 0.2e1_dp*t943*trhorho 1015 t1028 = 0.2e1_dp*t946*trhorho 1016 t1030 = 0.2e1_dp*t168*trhorhorho 1017 t1031 = t998 + t1001 + t1004 + t1006 + t1008 + t1011 + t1014 + & 1018 t1016 + t1019 + t1022 + t1024 + t1026 + t1028 + t1030 1019 t1035 = t940 + t942 + t945 + t948 + t950 1020 t1042 = t369*t2rho 1021 t1046 = A1rhorho*t83 1022 t1048 = 0.2e1_dp*t385*t2rho 1023 t1050 = 0.2e1_dp*t943*t1rho 1024 t1052 = 0.2e1_dp*t946*t1rho 1025 t1054 = 0.2e1_dp*t168*t1rhorho 1026 t1055 = t1046 + t1048 + t1050 + t1052 + t1054 1027 t1060 = t78*t86 1028 t1061 = t176**2 1029 t1062 = 0.1e1_dp/t1061 1030 t1067 = -0.2e1_dp*t162*t178*t967*t1rho + 0.2e1_dp*t175* & 1031 t409*t967*t369 + 0.2e1_dp*t976*t374 + 0.4e1_dp*t980* & 1032 t408*trho*t990 - t175*t995*t432 + t78*t83*t1031*t91 - & 1033 t175*t1035*t177*t369 - 0.2e1_dp*t162*t995*t370 - & 1034 0.2e1_dp*t162*t397*t1042 + 0.2e1_dp*t162*t1055*t91* & 1035 trho - 0.6e1_dp*t1060*t1062*t186*t990 1036 t1093 = t1046 + t1048 + t1050 + t1052 + t1054 + 0.2e1_dp*t951* & 1037 A1rho + 0.8e1_dp*t417*t1012 + 0.2e1_dp*t179*A1rhorho + & 1038 0.8e1_dp*t417*t1017 + 0.12e2_dp*t426*t1rho*t2rho + & 1039 0.4e1_dp*t183*t1rhorho 1040 t1097 = t1rho*t989 1041 t1103 = trho*t989 1042 t1104 = t178*t1103 1043 t1106 = t994*t91 1044 t1109 = t171*t408 1045 t1115 = t976*t362 + t162*t361*trho1rho - t162*t178*t432 & 1046 *t2rho + t175*t994*t408*t410 - t162*t178*trho1rho*t369 & 1047 - t162*t178*trho*t1093 - t162*t397*t1097 + t175*t409* & 1048 t186*t1093 - t354*t1104 + t162*t1106*trhorho + t175*t1109 & 1049 *t990 + t175*t409*t432*t989 1050 t1118 = t1106*trho 1051 t1121 = t360*t408 1052 t1122 = t186*t989 1053 t1126 = t393*t177 1054 t1129 = t408*t186 1055 t1148 = t186*t2rho 1056 t1152 = 0.2e1_dp*t354*t1118 + 0.2e1_dp*t175*t1121*t1122 & 1057 - t175*t1126*t989 + 0.4e1_dp*t980*t1129*t1042 + 0.2e1_dp* & 1058 t976*t378 - t175*t404*t967 + 0.2e1_dp*t162*t377* & 1059 t1rhorho - 0.2e1_dp*t976*t401 + 0.2e1_dp*t78*t1rhorho*t164 & 1060 - 0.2e1_dp*t162*t404*t1103 - 0.2e1_dp*t162*t404*t1148 1061 t1157 = t393*t91 1062 t1167 = A2rho*t182 1063 t1187 = A1rho*t182 1064 t1196 = t87*t 1065 t1203 = t1008 + 0.8e1_dp*t417*trhorho*A2rho + 0.12e2_dp* & 1066 t426*trhorho*t2rho + 0.8e1_dp*t1167*t423 + 0.8e1_dp*t417* & 1067 trho*A1rhorho + 0.8e1_dp*t417*Arho*t1rhorho + 0.8e1_dp* & 1068 t1167*t418 + 0.8e1_dp*t417*Arho1rho*t1rho + 0.12e2_dp*t426 & 1069 *trho1rho*t1rho + 0.12e2_dp*t426*trho*t1rhorho + t998 + & 1070 0.8e1_dp*t1187*t954 + 0.8e1_dp*t417*trho1rho*A1rho + & 1071 0.8e1_dp*t417*Arhorho*t2rho + 0.24e2_dp*t1196*t427*t2rho & 1072 + 0.2e1_dp*A1rhorho*t88*Arho + t1014 1073 t1218 = t1016 + t1001 + 0.2e1_dp*t951*Arhorho + t1026 + t1028 & 1074 + t1022 + t1011 + t1030 + 0.2e1_dp*t179*Arhorhorho + 0.4e1_dp* & 1075 t183*trhorhorho + 0.2e1_dp*t414*Arho1rho + t1004 + t1006 + & 1076 t1019 + t1024 + 0.24e2_dp*t84*t1018 + 0.24e2_dp*t84*t1005 + & 1077 0.24e2_dp*t84*t1013 1078 t1226 = t163*trho1rho 1079 t1249 = -0.2e1_dp*t162*t178*t186*t1rhorho + 0.2e1_dp* & 1080 t162*t1157*t2rho - t175*t178*(t1203 + t1218) + 0.2e1_dp* & 1081 t162*t1035*t91*t1rho + 0.2e1_dp*t354*t1226 - 0.2e1_dp* & 1082 t162*t995*t400 + 0.2e1_dp*t162*t163*trhorhorho - 0.2e1_dp & 1083 *t162*t178*trhorho*t989 - t175*t1055*t177*t186 - & 1084 0.2e1_dp*t976*t371 - t175*t397*t1093 + 0.4e1_dp*t980* & 1085 t1129*t1097 1086 t1262 = 0.2e1_dp*t162*t163*t2rho + t78*t83*t994*t91 - & 1087 t175*t178*t989 1088 t1263 = t439*t1262 1089 t1291 = 0.2e1_dp*t976*t164 + 0.2e1_dp*t162*t1118 - & 1090 0.2e1_dp*t162*t1104 + 0.2e1_dp*t162*t1226 + 0.2e1_dp*t162 & 1091 *t377*t2rho + t78*t83*t1035*t91 - t175*t397*t989 - & 1092 0.2e1_dp*t162*t178*t1148 - t175*t995*t186 + 0.2e1_dp* & 1093 t175*t409*t1122 - t175*t178*t967 1094 t1292 = gamma_var*t1291 1095 t1295 = 0.1e1_dp/t438/t94 1096 t1329 = 0.2e1_dp*t976*t440 + 0.2e1_dp*t162*t1106*t1rho - & 1097 0.2e1_dp*t162*t178*t1097 + 0.2e1_dp*t162*t163*t1rhorho & 1098 + 0.2e1_dp*t162*t361*t2rho + t78*t83*t1055*t91 - t175* & 1099 t404*t989 - 0.2e1_dp*t162*t178*t1042 - t175*t995*t369 + & 1100 0.2e1_dp*t175*t409*t990 - t175*t178*t1093 1101 kfrhorhorho = k_frhorhorho 1102 kf2rho = k_f2rho 1103 ex_unifrho1rho = ex_unifrhorho 1104 ex_unif1rhorho = ex_unifrho1rho 1105 ex_unif2rho = -0.3e1_dp/0.4e1_dp*t5*kf2rho 1106 t1342 = t195**2 1107 srho1rho = t457*t198*kf2rho + t462/0.2e1_dp - t465 + t197* & 1108 t115*kf2rho/0.2e1_dp + t466 1109 s1rhorho = srho1rho 1110 s2rho = -t197*t6*kf2rho/0.2e1_dp - t200/0.2e1_dp 1111 t1380 = t202**2 1112 t1385 = 0.1e1_dp/t1380*t469*mu*t101*s 1113 t1386 = kappa**2 1114 t1387 = 0.1e1_dp/t1386 1115 t1389 = s1rho*s2rho 1116 t1393 = t470*s 1117 Fxrho1rho = -0.8e1_dp*t471*t472*s2rho + 0.2e1_dp*t204* & 1118 s2rho*srho + 0.2e1_dp*t204*s*srho1rho 1119 Fx1rhorho = -0.8e1_dp*t471*s1rho*t103*s2rho + 0.2e1_dp* & 1120 t204*t1389 + 0.2e1_dp*t204*s*s1rhorho 1121 Fx2rho = 0.2e1_dp*t204*s*s2rho 1122 ex_ldarhorhorho = ex_unif1rhorho*Fx + ex_unif1rho*Fx2rho + & 1123 ex_unif2rho*Fx1rho + ex_unif*Fx1rhorho + ex_unifrho1rho*Fx + & 1124 ex_unifrho*Fx2rho + ex_unifrhorho*Fx - 0.3e1_dp/0.4e1_dp*my_rho & 1125 *t5*kfrhorhorho*Fx + t487*Fx2rho + ex_unifrho*Fx1rho + my_rho & 1126 *ex_unifrho1rho*Fx1rho + t208*Fx1rhorho + ex_unif2rho*Fxrho & 1127 + ex_unif*Fxrho1rho + ex_unif1rho*Fxrho + my_rho*ex_unif1rhorho* & 1128 Fxrho + t491*Fxrho1rho + ex_unif*Fxrhorho + my_rho*ex_unif2rho* & 1129 Fxrhorho + t108*(0.48e2_dp*t1385*srho*t1387*t1389 - & 1130 0.24e2_dp*t1393*t472*t1389 - 0.8e1_dp*t471*srho1rho*t103 & 1131 *s1rho - 0.8e1_dp*t471*t472*s1rhorho + 0.2e1_dp*t204* & 1132 s1rhorho*srho + 0.2e1_dp*t204*s1rho*srho1rho - 0.8e1_dp* & 1133 t471*srhorho*t103*s2rho + 0.2e1_dp*t204*s2rho*srhorho + & 1134 0.2e1_dp*t204*s*(-0.3e1_dp*my_norm_drho/t1342*t459*kf2rho & 1135 - t457*t115*t458 + 0.2e1_dp*t457*t463*kfrho - 0.2e1_dp* & 1136 t457*t461*kf2rho - 0.2e1_dp*t197*t246*kfrho + 0.3e1_dp/ & 1137 0.2e1_dp*t197*t115*kfrhorho + t457*t463*kf2rho - t197*t6 & 1138 *kfrhorhorho/0.2e1_dp - t197*t246*kf2rho - 0.3e1_dp*t99* & 1139 t241)) 1140 1141 e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + & 1142 scale_ex*ex_ldarhorhorho + scale_ec*( & 1143 e_c_u_01rhorho + gamma_var*t1329*t191 - t451*t1263 + & 1144 e_c_u_0rho1rho + t1292*t191 - t190*t1263 + epsilon_cGGArhorho + & 1145 my_rho*(e_c_u_0rhorhorho + gamma_var*(t1067 + 0.2e1_dp*t1115 + & 1146 t1152 + t1249)*t191 - t436*t1263 - t1292*t449 + 0.2e1_dp* & 1147 t190*t1295*t448*t1262 - t190*t439*t1329)) 1148 1149 t1468 = t149*t115 1150 tnorm_drhorhorho = t317*t6*t825 + t1468*k_srho/0.2e1_dp - & 1151 t496*k_srhorho/0.2e1_dp + t1468*k_s1rho/0.2e1_dp + t74* & 1152 t246 1153 tnorm_drho1rho = -t496*k_s1rho/0.2e1_dp - t498/0.2e1_dp 1154 t1482 = A1rho*tnorm_drhorho 1155 t1492 = t78*t84 1156 t1493 = tnorm_drho*t177 1157 t1504 = t408*tnorm_drho 1158 t1505 = t1504*t410 1159 t1511 = A1rho*tnorm_drho 1160 t1515 = t226*t1rho 1161 t1521 = A*tnorm_drho1rho 1162 t1525 = 0.2e1_dp*t175*t409*t553*t369 + 0.2e1_dp*t217* & 1163 t1482*t91 - 0.2e1_dp*t162*t178*tnorm_drhorho*t369 + & 1164 0.2e1_dp*t354*t510 - 0.6e1_dp*t1492*t1493*t400 + 0.6e1_dp & 1165 *t175*t218*t91*trhorho + 0.2e1_dp*t162*t1157*tnorm_drho & 1166 + 0.4e1_dp*t980*t1505 - 0.6e1_dp*t1492*t1493*t370 + & 1167 0.6e1_dp*t175*t1511*t513 - 0.2e1_dp*t162*t397*t1515 - & 1168 0.2e1_dp*t354*t507 + 0.6e1_dp*t175*t1521*t513 1169 t1528 = Arhorho*tnorm_drho 1170 t1532 = t177*t369 1171 t1545 = t78*t417 1172 t1565 = 0.2e1_dp*t385*tnorm_drho + 0.2e1_dp*t388* & 1173 tnorm_drho + 0.2e1_dp*t168*tnorm_drho1rho + 0.8e1_dp*t417* & 1174 t1511 + 0.12e2_dp*t426*tnorm_drho*t1rho + 0.4e1_dp*t183* & 1175 tnorm_drho1rho 1176 t1573 = t91*t1rho 1177 t1584 = -0.2e1_dp*t354*t530 + 0.2e1_dp*t217*t1528*t91 - & 1178 0.2e1_dp*t217*t517*t1532 - 0.2e1_dp*t217*t1511*t525 + & 1179 0.2e1_dp*t162*t377*tnorm_drho1rho + 0.2e1_dp*t162*t361* & 1180 tnorm_drhorho + 0.4e1_dp*t1545*t1505 + 0.2e1_dp*t217*A* & 1181 tnorm_drhorhorho*t91 + 0.2e1_dp*t175*t409*t1565*t186 - & 1182 0.2e1_dp*t217*t521*t1532 + 0.6e1_dp*t175*t521*t1573 - & 1183 0.6e1_dp*t1060*t1062*t226*t410 + 0.6e1_dp*t175*t517* & 1184 t1573 1185 t1608 = t408*t226 1186 t1612 = t163*tnorm_drho1rho 1187 t1628 = -0.2e1_dp*t162*t404*t506 + 0.2e1_dp*t354*t503 - & 1188 0.2e1_dp*t162*t178*tnorm_drho*t432 - 0.2e1_dp*t162*t178 & 1189 *t553*t1rho - 0.2e1_dp*t162*t404*t529 - 0.2e1_dp*t162* & 1190 t178*t1565*trho - t175*t1126*t226 + 0.4e1_dp*t980*t1608 & 1191 *t370 + 0.2e1_dp*t500*t1612 + 0.12e2_dp*t78*t168* & 1192 tnorm_drho*t91*t427 + 0.2e1_dp*t162*t163*tnorm_drhorhorho & 1193 - t175*t404*t553 + 0.2e1_dp*t175*t1121*t535 1194 t1629 = Arho*tnorm_drho1rho 1195 t1633 = tnorm_drho*t369 1196 t1646 = t361*tnorm_drho 1197 t1652 = t226*t369 1198 t1660 = t178*t1633 1199 t1672 = t418*tnorm_drho 1200 t1676 = t423*tnorm_drho 1201 t1715 = 0.2e1_dp*t999*tnorm_drho + 0.2e1_dp*t1672 + 0.2e1_dp & 1202 *t382*tnorm_drho1rho + 0.2e1_dp*t1676 + 0.2e1_dp*A*trhorho & 1203 *tnorm_drho + 0.2e1_dp*t541*tnorm_drho1rho + 0.2e1_dp*t385* & 1204 tnorm_drhorho + 0.2e1_dp*t388*tnorm_drhorho + 0.2e1_dp*t168* & 1205 tnorm_drhorhorho + 0.8e1_dp*t1187*t517 + 0.24e2_dp*t84* & 1206 t1672 + 0.8e1_dp*t417*t1629 + 0.8e1_dp*t417*t1528 + & 1207 0.24e2_dp*t84*t1676 + 0.24e2_dp*t1196*t548*t1rho + & 1208 0.12e2_dp*t426*tnorm_drho1rho*trho + 0.12e2_dp*t426* & 1209 tnorm_drho*trhorho + 0.8e1_dp*t417*t1482 + 0.12e2_dp*t426* & 1210 tnorm_drhorho*t1rho + 0.4e1_dp*t183*tnorm_drhorhorho 1211 t1722 = 0.2e1_dp*t217*t1629*t91 - 0.2e1_dp*t162*t397* & 1212 t1633 - t175*t397*t1565 - 0.2e1_dp*t162*t178*t226* & 1213 trhorho + 0.2e1_dp*t78*trhorho*t214 + 0.2e1_dp*t500*t1646 & 1214 - 0.2e1_dp*t217*t1521*t525 + 0.2e1_dp*t175*t1109*t1652 - & 1215 0.2e1_dp*t162*t178*tnorm_drho1rho*t186 - 0.2e1_dp*t500* & 1216 t1660 + 0.4e1_dp*t980*t1608*t400 + 0.2e1_dp*t175*t409* & 1217 t226*t432 - t175*t178*t1715 - 0.2e1_dp*t217*t218*t177* & 1218 t432 1219 t1758 = 0.2e1_dp*t354*t214 + 0.2e1_dp*t162*t1646 - & 1220 0.2e1_dp*t162*t1660 + 0.2e1_dp*t162*t1612 + 0.6e1_dp*t175 & 1221 *t218*t1573 + 0.2e1_dp*t217*t1511*t91 + 0.2e1_dp*t217* & 1222 t1521*t91 - 0.2e1_dp*t217*t218*t1532 - 0.2e1_dp*t162* & 1223 t178*t1515 - t175*t404*t226 + 0.2e1_dp*t175*t409*t1652 - & 1224 t175*t178*t1565 1225 t1759 = gamma_var*t1758 1226 t1761 = t1295*t189 1227 snorm_drho1rho = snorm_drhorho 1228 t1797 = snorm_drhorho*t103 1229 Fxnorm_drho1rho = -0.8e1_dp*t471*t566*s1rho + 0.2e1_dp* & 1230 t204*s1rho*snorm_drho + 0.2e1_dp*t204*s*snorm_drho1rho 1231 1232 e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) + & 1233 scale_ex*(ex_unif1rho*Fxnorm_drho + & 1234 ex_unif*Fxnorm_drho1rho + ex_unifrho*Fxnorm_drho + t487* & 1235 Fxnorm_drho + t208*Fxnorm_drho1rho + ex_unif*Fxnorm_drhorho + & 1236 t491*Fxnorm_drhorho + t108*(0.48e2_dp*t1385*snorm_drho* & 1237 t1387*t476 - 0.24e2_dp*t1393*t566*t476 - 0.8e1_dp*t471* & 1238 snorm_drho1rho*t103*srho - 0.8e1_dp*t471*t566*srhorho + & 1239 0.2e1_dp*t204*srhorho*snorm_drho + 0.2e1_dp*t204*srho* & 1240 snorm_drho1rho - 0.8e1_dp*t471*t1797*s1rho + 0.2e1_dp*t204* & 1241 s1rho*snorm_drhorho + 0.2e1_dp*t204*s*(t456*t6*t458 + & 1242 t196*t115*kfrho - t562*kfrhorho/0.2e1_dp + t98*t246))) + & 1243 scale_ec*(t1759*t191 - t230*t449 + Hnorm_drhorho + my_rho*( & 1244 gamma_var*(t1525 + t1584 + t1628 + t1722)*t191 - t557*t449 - & 1245 t1759*t559 + 0.2e1_dp*t230*t1761*t448 - t230*t439*t435)) 1246 1247 t1838 = t1504*t535 1248 t1851 = Arho*t581 1249 t1878 = 0.4e1_dp*t175*t409*t226*t553 - 0.2e1_dp*t162* & 1250 t178*t605*trho - 0.4e1_dp*t217*t218*t177*t553 + 0.8e1_dp & 1251 *t1545*t1838 + 0.2e1_dp*t78*t581*t171*t91 - 0.4e1_dp* & 1252 t162*t397*t590 - 0.10e2_dp*t175*t586*t525 - t175*t178* & 1253 (0.2e1_dp*t1851 + 0.4e1_dp*t521*tnorm_drho + 0.24e2_dp*t84* & 1254 t1851 + 0.24e2_dp*t1196*t581*trho + 0.24e2_dp*t426* & 1255 tnorm_drhorho*tnorm_drho) - 0.12e2_dp*t1492*t1493*t529 + & 1256 0.8e1_dp*t980*t1838 - 0.4e1_dp*t162*t178*tnorm_drhorho* & 1257 t226 + 0.20e2_dp*t162*t586*t513 1258 t1889 = t78*t581 1259 t1907 = t85*t1062 1260 t1922 = -0.4e1_dp*t500*t591 + 0.2e1_dp*t175*t409*t605* & 1261 t186 + 0.4e1_dp*t162*t409*t598*trho - 0.2e1_dp*t1889* & 1262 t187 + 0.2e1_dp*t175*t1109*t598 + 0.20e2_dp*t175*t218* & 1263 t91*tnorm_drhorho + 0.10e2_dp*t175*t1851*t91 - 0.4e1_dp* & 1264 t217*t517*t594 - t175*t397*t605 - 0.6e1_dp*t175*t1907* & 1265 t598*t186 - 0.4e1_dp*t162*t178*tnorm_drho*t553 + 0.4e1_dp & 1266 *t78*tnorm_drhorho*t214 - 0.4e1_dp*t217*t521*t594 1267 t1927 = t439*t229 1268 t1933 = t614*t1387 1269 t1937 = t614*t103 1270 1271 e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) + & 1272 scale_ex*(ex_unif* & 1273 Fxnorm_drhonorm_drho + t208*Fxnorm_drhonorm_drho + t108*( & 1274 0.48e2_dp*t1385*t1933*srho - 0.24e2_dp*t1393*t1937*srho & 1275 - 0.16e2_dp*t471*t1797*snorm_drho + 0.4e1_dp*t204* & 1276 snorm_drhorho*snorm_drho)) + scale_ec*(Hnorm_drhonorm_drho + my_rho & 1277 *(gamma_var*(t1878 + t1922)*t191 - t609*t559 - 0.2e1_dp* & 1278 t557*t1927 + 0.2e1_dp*t612*t1761)) 1279 1280 t2norm_drho = tnorm_drho 1281 t1952 = t226*t2norm_drho 1282 t1964 = 0.2e1_dp*t168*t2norm_drho + 0.4e1_dp*t183*t2norm_drho 1283 t1965 = t178*t1964 1284 t1968 = t226*t1964 1285 t1969 = t1504*t1968 1286 t1972 = A*t2norm_drho 1287 t1990 = 0.2e1_dp*t1972*tnorm_drho + 0.12e2_dp*t426* & 1288 tnorm_drho*t2norm_drho 1289 t2020 = t177*t1964 1290 t2024 = t91*t2norm_drho 1291 t2028 = t78*t2norm_drho 1292 t2031 = -0.20e2_dp*t1492*t1493*t1952 + 0.4e1_dp*t162* & 1293 t409*t598*t2norm_drho - 0.2e1_dp*t1889*t1965 + 0.8e1_dp* & 1294 t1545*t1969 + 0.4e1_dp*t217*t1972*t408*t598 - 0.6e1_dp* & 1295 t175*t1907*t598*t1964 - 0.2e1_dp*t217*t1972*t177*t605 & 1296 - 0.4e1_dp*t217*t218*t177*t1990 + 0.4e1_dp*t175*t409* & 1297 t1990*t226 - 0.24e2_dp*t78*t182*t85*t177*t87*t581* & 1298 t2norm_drho + 0.2e1_dp*t175*t409*t605*t1964 + 0.8e1_dp* & 1299 t980*t1969 - 0.2e1_dp*t162*t178*t605*t2norm_drho - & 1300 0.4e1_dp*t162*t178*t1990*tnorm_drho - 0.10e2_dp*t175* & 1301 t586*t2020 + 0.24e2_dp*t162*t586*t2024 - 0.4e1_dp*t2028* & 1302 t591 1303 t2041 = 0.2e1_dp*t162*t163*t2norm_drho + 0.2e1_dp*t217* & 1304 t1972*t91 - t175*t1965 1305 s2norm_drho = snorm_drho 1306 1307 e_ndrho_ndrho_ndrho(ii) = e_ndrho_ndrho_ndrho(ii) + & 1308 scale_ex*t108*(0.48e2_dp* & 1309 t1385*t1933*s2norm_drho - 0.24e2_dp*t1393*t1937* & 1310 s2norm_drho) + scale_ec*my_rho*(gamma_var*t2031*t191 - t609* & 1311 t439*t2041 - 0.2e1_dp*gamma_var*(0.2e1_dp*t2028*t214 + & 1312 0.10e2_dp*t175*t218*t2024 - 0.2e1_dp*t162*t178* & 1313 tnorm_drho*t1964 - 0.2e1_dp*t217*t218*t2020 - 0.2e1_dp* & 1314 t162*t178*t1952 - 0.2e1_dp*t217*t1972*t594 + 0.2e1_dp* & 1315 t175*t409*t1968 - t175*t178*t1990)*t1927 + 0.2e1_dp*t612 & 1316 *t1295*t2041) 1317 END IF 1318 END IF 1319 END DO 1320!$OMP END DO 1321 END SELECT 1322 1323 END SUBROUTINE pbe_lda_calc 1324 1325! ************************************************************************************************** 1326!> \brief evaluates the becke 88 exchange functional for lsd 1327!> \param rho_set ... 1328!> \param deriv_set ... 1329!> \param grad_deriv ... 1330!> \param pbe_params ... 1331!> \author fawzi 1332! ************************************************************************************************** 1333 SUBROUTINE pbe_lsd_eval(rho_set, deriv_set, grad_deriv, pbe_params) 1334 TYPE(xc_rho_set_type), POINTER :: rho_set 1335 TYPE(xc_derivative_set_type), POINTER :: deriv_set 1336 INTEGER, INTENT(in) :: grad_deriv 1337 TYPE(section_vals_type), POINTER :: pbe_params 1338 1339 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lsd_eval', routineP = moduleN//':'//routineN 1340 1341 INTEGER :: handle, npoints, param 1342 INTEGER, DIMENSION(:, :), POINTER :: bo 1343 REAL(kind=dp) :: epsilon_rho, scale_ec, scale_ex 1344 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, e_ndr_ndr, e_ndr_ra, & 1345 e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, e_ndrb, e_ndrb_ndrb, e_ndrb_rb, e_ra, e_ra_ra, & 1346 e_ra_rb, e_rb, e_rb_rb, norm_drho, norm_drhoa, norm_drhob, rhoa, rhob 1347 TYPE(xc_derivative_type), POINTER :: deriv 1348 1349 CALL timeset(routineN, handle) 1350 NULLIFY (deriv, bo) 1351 1352 CPASSERT(ASSOCIATED(rho_set)) 1353 CPASSERT(rho_set%ref_count > 0) 1354 CPASSERT(ASSOCIATED(deriv_set)) 1355 CPASSERT(deriv_set%ref_count > 0) 1356 CALL xc_rho_set_get(rho_set, & 1357 rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, & 1358 norm_drhob=norm_drhob, norm_drho=norm_drho, & 1359 rho_cutoff=epsilon_rho, & 1360 local_bounds=bo) 1361 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) 1362 1363 dummy => rhoa 1364 e_0 => dummy 1365 e_ra => dummy 1366 e_rb => dummy 1367 e_ndra_ra => dummy 1368 e_ndrb_rb => dummy 1369 e_ndr_ndr => dummy 1370 e_ndra_ndra => dummy 1371 e_ndrb_ndrb => dummy 1372 e_ndr => dummy 1373 e_ndra => dummy 1374 e_ndrb => dummy 1375 e_ra_ra => dummy 1376 e_ra_rb => dummy 1377 e_rb_rb => dummy 1378 e_ndr_ra => dummy 1379 e_ndr_rb => dummy 1380 1381 IF (grad_deriv >= 0) THEN 1382 deriv => xc_dset_get_derivative(deriv_set, "", & 1383 allocate_deriv=.TRUE.) 1384 CALL xc_derivative_get(deriv, deriv_data=e_0) 1385 END IF 1386 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 1387 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)", & 1388 allocate_deriv=.TRUE.) 1389 CALL xc_derivative_get(deriv, deriv_data=e_ra) 1390 deriv => xc_dset_get_derivative(deriv_set, "(rhob)", & 1391 allocate_deriv=.TRUE.) 1392 CALL xc_derivative_get(deriv, deriv_data=e_rb) 1393 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", & 1394 allocate_deriv=.TRUE.) 1395 CALL xc_derivative_get(deriv, deriv_data=e_ndr) 1396 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)", & 1397 allocate_deriv=.TRUE.) 1398 CALL xc_derivative_get(deriv, deriv_data=e_ndra) 1399 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhob)", & 1400 allocate_deriv=.TRUE.) 1401 CALL xc_derivative_get(deriv, deriv_data=e_ndrb) 1402 END IF 1403 IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN 1404 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)(rhoa)", & 1405 allocate_deriv=.TRUE.) 1406 CALL xc_derivative_get(deriv, deriv_data=e_ra_ra) 1407 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)(rhob)", & 1408 allocate_deriv=.TRUE.) 1409 CALL xc_derivative_get(deriv, deriv_data=e_ra_rb) 1410 deriv => xc_dset_get_derivative(deriv_set, "(rhob)(rhob)", & 1411 allocate_deriv=.TRUE.) 1412 CALL xc_derivative_get(deriv, deriv_data=e_rb_rb) 1413 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)(rhoa)", & 1414 allocate_deriv=.TRUE.) 1415 CALL xc_derivative_get(deriv, deriv_data=e_ndr_ra) 1416 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)(rhob)", & 1417 allocate_deriv=.TRUE.) 1418 CALL xc_derivative_get(deriv, deriv_data=e_ndr_rb) 1419 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)(rhoa)", & 1420 allocate_deriv=.TRUE.) 1421 CALL xc_derivative_get(deriv, deriv_data=e_ndra_ra) 1422 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhob)(rhob)", & 1423 allocate_deriv=.TRUE.) 1424 CALL xc_derivative_get(deriv, deriv_data=e_ndrb_rb) 1425 deriv => xc_dset_get_derivative(deriv_set, & 1426 "(norm_drho)(norm_drho)", allocate_deriv=.TRUE.) 1427 CALL xc_derivative_get(deriv, deriv_data=e_ndr_ndr) 1428 deriv => xc_dset_get_derivative(deriv_set, & 1429 "(norm_drhoa)(norm_drhoa)", allocate_deriv=.TRUE.) 1430 CALL xc_derivative_get(deriv, deriv_data=e_ndra_ndra) 1431 deriv => xc_dset_get_derivative(deriv_set, & 1432 "(norm_drhob)(norm_drhob)", allocate_deriv=.TRUE.) 1433 CALL xc_derivative_get(deriv, deriv_data=e_ndrb_ndrb) 1434 END IF 1435 IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN 1436 ! to do 1437 END IF 1438 1439 CALL section_vals_val_get(pbe_params, "scale_c", r_val=scale_ec) 1440 CALL section_vals_val_get(pbe_params, "scale_x", r_val=scale_ex) 1441 CALL section_vals_val_get(pbe_params, "parametrization", i_val=param) 1442 1443!$OMP PARALLEL DEFAULT(NONE),& 1444!$OMP SHARED(rhoa,rhob,norm_drho,norm_drhoa,norm_drhob,e_0,e_ra,e_rb,e_ndra_ra),& 1445!$OMP SHARED(e_ndrb_rb,e_ndr_ndr,e_ndra_ndra,e_ndrb_ndrb,e_ndr,e_ndra,e_ndrb),& 1446!$OMP SHARED(e_ra_ra,e_ra_rb,e_rb_rb,e_ndr_ra,e_ndr_rb,grad_deriv,npoints),& 1447!$OMP SHARED(epsilon_rho,param,scale_ec,scale_ex) 1448 CALL pbe_lsd_calc( & 1449 rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa, & 1450 norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb, & 1451 e_ra_ndra=e_ndra_ra, & 1452 e_rb_ndrb=e_ndrb_rb, e_ndr_ndr=e_ndr_ndr, & 1453 e_ndra_ndra=e_ndra_ndra, e_ndrb_ndrb=e_ndrb_ndrb, e_ndr=e_ndr, & 1454 e_ndra=e_ndra, e_ndrb=e_ndrb, e_ra_ra=e_ra_ra, & 1455 e_ra_rb=e_ra_rb, e_rb_rb=e_rb_rb, e_ra_ndr=e_ndr_ra, & 1456 e_rb_ndr=e_ndr_rb, & 1457 grad_deriv=grad_deriv, npoints=npoints, & 1458 epsilon_rho=epsilon_rho, & 1459 param=param, scale_ec=scale_ec, scale_ex=scale_ex) 1460!$OMP END PARALLEL 1461 1462 CALL timestop(handle) 1463 END SUBROUTINE pbe_lsd_eval 1464 1465! ************************************************************************************************** 1466!> \brief low level calculation of the pbe exchange-correlation functional for lsd 1467!> \param rhoa ,rhob: alpha or beta spin density 1468!> \param rhob ... 1469!> \param norm_drho ... 1470!> \param norm_drhoa ,norm_drhob,norm_drho: || grad rhoa |||,| grad rhoa ||, 1471!> || grad (rhoa+rhob) || 1472!> \param norm_drhob ... 1473!> \param e_0 adds to it the local value of the functional 1474!> \param e_ra derivative of the functional wrt. to the variables 1475!> named where the * is 1476!> \param e_rb ... 1477!> \param e_ra_ndra ... 1478!> \param e_rb_ndrb ... 1479!> \param e_ndr_ndr ... 1480!> \param e_ndra_ndra ... 1481!> \param e_ndrb_ndrb ... 1482!> \param e_ndr ... 1483!> \param e_ndra ... 1484!> \param e_ndrb ... 1485!> \param e_ra_ra ... 1486!> \param e_ra_rb ... 1487!> \param e_rb_rb ... 1488!> \param e_ra_ndr ... 1489!> \param e_rb_ndr ... 1490!> \param grad_deriv ... 1491!> \param npoints ... 1492!> \param epsilon_rho ... 1493!> \param param ... 1494!> \param scale_ec ... 1495!> \param scale_ex ... 1496!> \author fawzi 1497! ************************************************************************************************** 1498 SUBROUTINE pbe_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, & 1499 e_0, e_ra, e_rb, e_ra_ndra, e_rb_ndrb, e_ndr_ndr, & 1500 e_ndra_ndra, e_ndrb_ndrb, e_ndr, & 1501 e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, e_ra_ndr, e_rb_ndr, & 1502 grad_deriv, npoints, epsilon_rho, param, scale_ec, scale_ex) 1503 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob, norm_drho, norm_drhoa, & 1504 norm_drhob 1505 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_ra, e_rb, e_ra_ndra, e_rb_ndrb, & 1506 e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb, e_ndr, e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, & 1507 e_ra_ndr, e_rb_ndr 1508 INTEGER, INTENT(in) :: grad_deriv, npoints 1509 REAL(kind=dp), INTENT(in) :: epsilon_rho 1510 INTEGER, INTENT(in) :: param 1511 REAL(kind=dp), INTENT(in) :: scale_ec, scale_ex 1512 1513 CHARACTER(len=*), PARAMETER :: routineN = 'pbe_lsd_calc', routineP = moduleN//':'//routineN 1514 1515 INTEGER :: ii 1516 REAL(kind=dp) :: A, A1rhoa, A1rhob, A_1, A_2, A_3, alpha_1_1, alpha_1_2, alpha_1_3, alpha_c, & 1517 alpha_c1rhoa, alpha_c1rhob, alpha_crhoa, alpha_crhob, Arhoa, Arhoarhoa, Arhoarhob, Arhob, & 1518 Arhobrhob, beta, beta_1_1, beta_1_2, beta_1_3, beta_2_1, beta_2_2, beta_2_3, beta_3_1, & 1519 beta_3_2, beta_3_3, beta_4_1, beta_4_2, beta_4_3, chi, chirhoa, chirhoarhoa, chirhoarhob, & 1520 chirhob, chirhobrhob, e_c_u_0, e_c_u_01rhoa, e_c_u_01rhob, e_c_u_0rhoa, e_c_u_0rhoarhoa, & 1521 e_c_u_0rhoarhob, e_c_u_0rhob, e_c_u_0rhobrhob, e_c_u_1rhoa, e_c_u_1rhob, epsilon_c_unif, & 1522 epsilon_c_unif1rhoa, epsilon_c_unif1rhob 1523 REAL(kind=dp) :: epsilon_c_unifrhoa, epsilon_c_unifrhoarhoa, epsilon_c_unifrhoarhob, & 1524 epsilon_c_unifrhob, epsilon_c_unifrhobrhob, epsilon_cGGA, epsilon_cGGArhoa, & 1525 epsilon_cGGArhob, ex_unif_a, ex_unif_a1rhoa, ex_unif_arhoa, ex_unif_b, ex_unif_b1rhob, & 1526 ex_unif_brhob, f, f1rhoa, f1rhob, frhoa, frhoarhoa, frhoarhob, frhob, frhobrhob, Fx_a, & 1527 Fx_a1rhoa, Fx_anorm_drhoa, Fx_arhoa, Fx_b, Fx_b1rhob, Fx_bnorm_drhob, Fx_brhob, & 1528 gamma_var, Hnorm_drho, k_frhoa, k_frhoarhoa, k_frhoarhob, k_frhob, k_s, k_s1rhoa, & 1529 k_s1rhob, k_srhoa, k_srhob, kappa, kf_a, kf_arhoa, kf_arhoarhoa, kf_b, kf_brhob, & 1530 kf_brhobrhob, mu 1531 REAL(kind=dp) :: my_norm_drho, my_norm_drhoa, my_norm_drhob, my_rho, my_rhoa, my_rhob, p_1, & 1532 p_2, p_3, phi, phi1rhoa, phi1rhob, phirhoa, phirhoarhoa, phirhoarhob, phirhob, & 1533 phirhobrhob, rs, rsrhoa, rsrhoarhoa, rsrhoarhob, rsrhob, rsrhobrhob, s_a, s_a1rhoa, & 1534 s_anorm_drhoa, s_arhoa, s_b, s_b1rhob, s_bnorm_drhob, s_brhob, t, t1, t100, t1000, t1001, & 1535 t101, t102, t103, t1031, t1033, t1038, t104, t105, t1050, t1069, t107, t1071, t108, t110, & 1536 t1104, t1106, t111, t112, t113, t115, t116, t1164, t118, t119, t1193, t12, t122, t1228, & 1537 t123, t1231, t124, t125, t126, t1269, t1283, t1285, t1286, t1288, t129 1538 REAL(kind=dp) :: t1291, t1293, t130, t1304, t131, t1327, t133, t1330, t1342, t1346, t135, & 1539 t137, t1377, t14, t140, t1411, t142, t143, t1440, t146, t1462, t1465, t147, t148, t1482, & 1540 t1489, t1492, t15, t150, t1501, t1514, t1520, t1529, t153, t1550, t1552, t1555, t156, & 1541 t1588, t1590, t1591, t1599, t1602, t1603, t162, t163, t1630, t1632, t1635, t1638, t164, & 1542 t1640, t165, t167, t1674, t1677, t1680, t1698, t171, t1711, t1712, t1713, t1739, t1741, & 1543 t1743, t1748, t175, t176, t1765, t1766, t1767, t177, t178, t179, t18, t1801, t1829, t183, & 1544 t1830, t1831, t1865, t187, t1871, t1876, t1888, t190, t1901, t191 1545 REAL(kind=dp) :: t192, t1922, t194, t1949, t198, t199, t1rhoa, t1rhob, t2, t20, t200, t201, & 1546 t205, t21, t211, t212, t213, t215, t219, t22, t220, t221, t222, t226, t23, t232, t233, & 1547 t234, t240, t242, t244, t245, t246, t248, t249, t250, t252, t254, t256, t257, t259, t262, & 1548 t266, t268, t269, t27, t272, t273, t274, t277, t278, t28, t280, t281, t282, t284, t285, & 1549 t286, t289, t293, t294, t297, t298, t299, t3, t302, t303, t305, t306, t310, t311, t312, & 1550 t313, t314, t317, t318, t32, t321, t324, t325, t326, t329, t335, t336, t337, t34, t340, & 1551 t341, t344, t346, t350, t368, t38, t382, t39, t396, t4, t40 1552 REAL(kind=dp) :: t403, t405, t407, t409, t41, t410, t411, t413, t415, t417, t427, t429, & 1553 t432, t436, t439, t442, t444, t445, t45, t453, t456, t457, t46, t460, t466, t467, t468, & 1554 t471, t472, t475, t477, t481, t488, t493, t494, t5, t50, t502, t505, t506, t518, t519, & 1555 t52, t523, t525, t536, t538, t543, t544, t548, t549, t550, t556, t56, t561, t563, t564, & 1556 t57, t576, t578, t579, t58, t580, t588, t59, t590, t595, t596, t6, t600, t606, t611, & 1557 t624, t626, t627, t628, t63, t636, t638, t64, t643, t644, t648, t654, t659, t66, t672, & 1558 t674, t675, t676, t681, t682, t687, t69, t7, t70, t705, t708, t71, t711 1559 REAL(kind=dp) :: t72, t726, t73, t733, t736, t74, t745, t75, t750, t755, t763, t767, t768, & 1560 t77, t775, t776, t779, t78, t785, t789, t79, t795, t798, t8, t80, t801, t812, t82, t820, & 1561 t821, t822, t823, t825, t828, t839, t84, t840, t85, t851, t858, t865, t867, t868, t87, & 1562 t876, t879, t88, t880, t9, t90, t904, t908, t909, t91, t911, t914, t917, t919, t92, t924, & 1563 t93, t936, t94, t944, t95, t953, t959, t96, t962, t965, t966, t967, t97, t98, t985, t998, & 1564 t999, tnorm_drho, trhoa, trhoanorm_drho, trhoarhoa, trhoarhob, trhob, trhobnorm_drho, & 1565 trhobrhob 1566 1567! Parametrization of PBE 1568 1569 SELECT CASE (param) 1570 ! Original PBE 1571 CASE (xc_pbe_orig) 1572 kappa = 0.804e0_dp 1573 beta = 0.66725e-1_dp 1574 mu = beta*pi**2/0.3e1_dp 1575 ! Revised PBE (revPBE) 1576 CASE (xc_pbe_rev) 1577 kappa = 0.1245e1_dp 1578 beta = 0.66725e-1_dp 1579 mu = beta*pi**2/0.3e1_dp 1580 ! PBE for solids (PBEsol) 1581 CASE (xc_pbe_sol) 1582 kappa = 0.804e0_dp 1583 beta = 0.46e-1_dp 1584 mu = 0.1e2_dp/0.81e2_dp 1585 CASE default 1586 CPABORT("") 1587 END SELECT 1588 1589 gamma_var = (0.1e1_dp - LOG(0.2e1_dp))/pi**2 1590 p_1 = 0.10e1_dp 1591 A_1 = 0.31091e-1_dp 1592 alpha_1_1 = 0.21370e0_dp 1593 beta_1_1 = 0.75957e1_dp 1594 beta_2_1 = 0.35876e1_dp 1595 beta_3_1 = 0.16382e1_dp 1596 beta_4_1 = 0.49294e0_dp 1597 p_2 = 0.10e1_dp 1598 A_2 = 0.15545e-1_dp 1599 alpha_1_2 = 0.20548e0_dp 1600 beta_1_2 = 0.141189e2_dp 1601 beta_2_2 = 0.61977e1_dp 1602 beta_3_2 = 0.33662e1_dp 1603 beta_4_2 = 0.62517e0_dp 1604 p_3 = 0.10e1_dp 1605 A_3 = 0.16887e-1_dp 1606 alpha_1_3 = 0.11125e0_dp 1607 beta_1_3 = 0.10357e2_dp 1608 beta_2_3 = 0.36231e1_dp 1609 beta_3_3 = 0.88026e0_dp 1610 beta_4_3 = 0.49671e0_dp 1611 t3 = 3**(0.1e1_dp/0.3e1_dp) 1612 t4 = 4**(0.1e1_dp/0.3e1_dp) 1613 t5 = t4**2 1614 t6 = t3*t5 1615 t7 = 0.1e1_dp/pi 1616 t90 = pi**2 1617 1618 SELECT CASE (grad_deriv) 1619 CASE default 1620!$OMP DO 1621 DO ii = 1, npoints 1622 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1623 my_rhob = MAX(rhob(ii), 0.0_dp) 1624 my_rho = my_rhoa + my_rhob 1625 IF (my_rho > epsilon_rho) THEN 1626 my_rhoa = MAX(EPSILON(0.0_dp)*1.e4_dp, my_rhoa) 1627 my_rhob = MAX(EPSILON(0.0_dp)*1.e4_dp, my_rhob) 1628 my_rho = my_rhoa + my_rhob 1629 my_norm_drho = norm_drho(ii) 1630 my_norm_drhoa = norm_drhoa(ii) 1631 my_norm_drhob = norm_drhob(ii) 1632 1633 t1 = my_rhoa - my_rhob 1634 t2 = 0.1e1_dp/my_rho 1635 chi = t1*t2 1636 t8 = t7*t2 1637 t9 = t8**(0.1e1_dp/0.3e1_dp) 1638 rs = t6*t9/0.4e1_dp 1639 t12 = 0.1e1_dp + alpha_1_1*rs 1640 t14 = 0.1e1_dp/A_1 1641 t15 = SQRT(rs) 1642 t18 = t15*rs 1643 t20 = p_1 + 0.1e1_dp 1644 t21 = rs**t20 1645 t22 = beta_4_1*t21 1646 t23 = beta_1_1*t15 + beta_2_1*rs + beta_3_1*t18 + t22 1647 t27 = 0.1e1_dp + t14/t23/0.2e1_dp 1648 t28 = LOG(t27) 1649 e_c_u_0 = -0.2e1_dp*A_1*t12*t28 1650 t32 = 0.1e1_dp + alpha_1_2*rs 1651 t34 = 0.1e1_dp/A_2 1652 t38 = p_2 + 0.1e1_dp 1653 t39 = rs**t38 1654 t40 = beta_4_2*t39 1655 t41 = beta_1_2*t15 + beta_2_2*rs + beta_3_2*t18 + t40 1656 t45 = 0.1e1_dp + t34/t41/0.2e1_dp 1657 t46 = LOG(t45) 1658 t50 = 0.1e1_dp + alpha_1_3*rs 1659 t52 = 0.1e1_dp/A_3 1660 t56 = p_3 + 0.1e1_dp 1661 t57 = rs**t56 1662 t58 = beta_4_3*t57 1663 t59 = beta_1_3*t15 + beta_2_3*rs + beta_3_3*t18 + t58 1664 t63 = 0.1e1_dp + t52/t59/0.2e1_dp 1665 t64 = LOG(t63) 1666 alpha_c = 0.2e1_dp*A_3*t50*t64 1667 t66 = 2**(0.1e1_dp/0.3e1_dp) 1668 t69 = 1/(2*t66 - 2) 1669 t70 = 0.1e1_dp + chi 1670 t71 = t70**(0.1e1_dp/0.3e1_dp) 1671 t72 = t71*t70 1672 t73 = 0.1e1_dp - chi 1673 t74 = t73**(0.1e1_dp/0.3e1_dp) 1674 t75 = t74*t73 1675 f = (t72 + t75 - 0.2e1_dp)*t69 1676 t77 = alpha_c*f 1677 t78 = 0.9e1_dp/0.8e1_dp/t69 1678 t79 = chi**2 1679 t80 = t79**2 1680 t82 = t78*(0.1e1_dp - t80) 1681 t84 = -0.2e1_dp*A_2*t32*t46 - e_c_u_0 1682 t85 = t84*f 1683 epsilon_c_unif = e_c_u_0 + t77*t82 + t85*t80 1684 t87 = t71**2 1685 t88 = t74**2 1686 phi = t87/0.2e1_dp + t88/0.2e1_dp 1687 t91 = t90*my_rho 1688 t92 = t91**(0.1e1_dp/0.3e1_dp) 1689 t93 = t3*t92*t7 1690 t94 = SQRT(t93) 1691 k_s = 0.2e1_dp*t94 1692 t95 = 0.1e1_dp/phi 1693 t96 = my_norm_drho*t95 1694 t97 = 0.1e1_dp/k_s 1695 t98 = t97*t2 1696 t = t96*t98/0.2e1_dp 1697 t100 = 0.1e1_dp/gamma_var 1698 t101 = beta*t100 1699 t102 = epsilon_c_unif*t100 1700 t103 = phi**2 1701 t104 = t103*phi 1702 t105 = 0.1e1_dp/t104 1703 t107 = EXP(-t102*t105) 1704 t108 = t107 - 0.1e1_dp 1705 A = t101/t108 1706 t110 = gamma_var*t104 1707 t111 = t**2 1708 t112 = A*t111 1709 t113 = 0.1e1_dp + t112 1710 t115 = A**2 1711 t116 = t111**2 1712 t118 = 0.1e1_dp + t112 + t115*t116 1713 t119 = 0.1e1_dp/t118 1714 t122 = 0.1e1_dp + t101*t111*t113*t119 1715 t123 = LOG(t122) 1716 epsilon_cGGA = epsilon_c_unif + t110*t123 1717 t124 = t3*t66 1718 t125 = t90*my_rhoa 1719 t126 = t125**(0.1e1_dp/0.3e1_dp) 1720 kf_a = t124*t126 1721 ex_unif_a = -0.3e1_dp/0.4e1_dp*t7*kf_a 1722 t129 = 0.1e1_dp/kf_a 1723 t130 = my_norm_drhoa*t129 1724 t131 = 0.1e1_dp/my_rhoa 1725 s_a = t130*t131/0.2e1_dp 1726 t133 = s_a**2 1727 t135 = 0.1e1_dp/kappa 1728 t137 = 0.1e1_dp + mu*t133*t135 1729 Fx_a = 0.1e1_dp + kappa - kappa/t137 1730 t140 = my_rhoa*ex_unif_a 1731 t142 = t90*my_rhob 1732 t143 = t142**(0.1e1_dp/0.3e1_dp) 1733 kf_b = t124*t143 1734 ex_unif_b = -0.3e1_dp/0.4e1_dp*t7*kf_b 1735 t146 = 0.1e1_dp/kf_b 1736 t147 = my_norm_drhob*t146 1737 t148 = 0.1e1_dp/my_rhob 1738 s_b = t147*t148/0.2e1_dp 1739 t150 = s_b**2 1740 t153 = 0.1e1_dp + mu*t150*t135 1741 Fx_b = 0.1e1_dp + kappa - kappa/t153 1742 t156 = my_rhob*ex_unif_b 1743 1744 IF (grad_deriv >= 0) THEN 1745 e_0(ii) = e_0(ii) + & 1746 scale_ex*(0.2e1_dp*t140*Fx_a + 0.2e1_dp*t156*Fx_b) & 1747 /0.2e1_dp + scale_ec*my_rho*epsilon_cGGA 1748 END IF 1749 1750 t162 = my_rho**2 1751 t163 = 0.1e1_dp/t162 1752 t164 = t1*t163 1753 chirhoa = t2 - t164 1754 t165 = t9**2 1755 t167 = 0.1e1_dp/t165*t7 1756 rsrhoa = -t6*t167*t163/0.12e2_dp 1757 t171 = A_1*alpha_1_1 1758 t175 = t23**2 1759 t176 = 0.1e1_dp/t175 1760 t177 = t12*t176 1761 t178 = 0.1e1_dp/t15 1762 t179 = beta_1_1*t178 1763 t183 = beta_3_1*t15 1764 t187 = 0.1e1_dp/rs 1765 t190 = t179*rsrhoa/0.2e1_dp + beta_2_1*rsrhoa + 0.3e1_dp/ & 1766 0.2e1_dp*t183*rsrhoa + t22*t20*rsrhoa*t187 1767 t191 = 0.1e1_dp/t27 1768 t192 = t190*t191 1769 e_c_u_0rhoa = -0.2e1_dp*t171*rsrhoa*t28 + t177*t192 1770 t194 = A_2*alpha_1_2 1771 t198 = t41**2 1772 t199 = 0.1e1_dp/t198 1773 t200 = t32*t199 1774 t201 = beta_1_2*t178 1775 t205 = beta_3_2*t15 1776 t211 = t201*rsrhoa/0.2e1_dp + beta_2_2*rsrhoa + 0.3e1_dp/ & 1777 0.2e1_dp*t205*rsrhoa + t40*t38*rsrhoa*t187 1778 t212 = 0.1e1_dp/t45 1779 t213 = t211*t212 1780 e_c_u_1rhoa = -0.2e1_dp*t194*rsrhoa*t46 + t200*t213 1781 t215 = A_3*alpha_1_3 1782 t219 = t59**2 1783 t220 = 0.1e1_dp/t219 1784 t221 = t50*t220 1785 t222 = beta_1_3*t178 1786 t226 = beta_3_3*t15 1787 t232 = t222*rsrhoa/0.2e1_dp + beta_2_3*rsrhoa + 0.3e1_dp/ & 1788 0.2e1_dp*t226*rsrhoa + t58*t56*rsrhoa*t187 1789 t233 = 0.1e1_dp/t63 1790 t234 = t232*t233 1791 alpha_crhoa = 0.2e1_dp*t215*rsrhoa*t64 - t221*t234 1792 frhoa = (0.4e1_dp/0.3e1_dp*t71*chirhoa - 0.4e1_dp/0.3e1_dp & 1793 *t74*chirhoa)*t69 1794 t240 = alpha_crhoa*f 1795 t242 = alpha_c*frhoa 1796 t244 = t79*chi 1797 t245 = t78*t244 1798 t246 = t245*chirhoa 1799 t248 = 0.4e1_dp*t77*t246 1800 t249 = e_c_u_1rhoa - e_c_u_0rhoa 1801 t250 = t249*f 1802 t252 = t84*frhoa 1803 t254 = t244*chirhoa 1804 t256 = 0.4e1_dp*t85*t254 1805 epsilon_c_unifrhoa = e_c_u_0rhoa + t240*t82 + t242*t82 - t248 & 1806 + t250*t80 + t252*t80 + t256 1807 t257 = 0.1e1_dp/t71 1808 t259 = 0.1e1_dp/t74 1809 phirhoa = t257*chirhoa/0.3e1_dp - t259*chirhoa/0.3e1_dp 1810 t262 = t92**2 1811 k_frhoa = t3/t262*t90/0.3e1_dp 1812 t266 = 0.1e1_dp/t94 1813 k_srhoa = t266*k_frhoa*t7 1814 t268 = 0.1e1_dp/t103 1815 t269 = my_norm_drho*t268 1816 t272 = k_s**2 1817 t273 = 0.1e1_dp/t272 1818 t274 = t273*t2 1819 t277 = t97*t163 1820 t278 = t96*t277 1821 trhoa = -t269*t98*phirhoa/0.2e1_dp - t96*t274*k_srhoa/ & 1822 0.2e1_dp - t278/0.2e1_dp 1823 t280 = t108**2 1824 t281 = 0.1e1_dp/t280 1825 t282 = epsilon_c_unifrhoa*t100 1826 t284 = t103**2 1827 t285 = 0.1e1_dp/t284 1828 t286 = t285*phirhoa 1829 t289 = -t282*t105 + 0.3e1_dp*t102*t286 1830 Arhoa = -t101*t281*t289*t107 1831 t293 = gamma_var*t103 1832 t294 = t123*phirhoa 1833 t297 = t101*t 1834 t298 = t113*t119 1835 t299 = t298*trhoa 1836 t302 = Arhoa*t111 1837 t303 = A*t 1838 t305 = 0.2e1_dp*t303*trhoa 1839 t306 = t302 + t305 1840 t310 = t101*t111 1841 t311 = t118**2 1842 t312 = 0.1e1_dp/t311 1843 t313 = t113*t312 1844 t314 = A*t116 1845 t317 = t111*t 1846 t318 = t115*t317 1847 t321 = t302 + t305 + 0.2e1_dp*t314*Arhoa + 0.4e1_dp*t318*trhoa 1848 t324 = 0.2e1_dp*t297*t299 + t101*t111*t306*t119 - t310* & 1849 t313*t321 1850 t325 = 0.1e1_dp/t122 1851 t326 = t324*t325 1852 epsilon_cGGArhoa = epsilon_c_unifrhoa + 0.3e1_dp*t293*t294 + & 1853 t110*t326 1854 t329 = t126**2 1855 kf_arhoa = t124/t329*t90/0.3e1_dp 1856 ex_unif_arhoa = -0.3e1_dp/0.4e1_dp*t7*kf_arhoa 1857 t335 = kf_a**2 1858 t336 = 0.1e1_dp/t335 1859 t337 = my_norm_drhoa*t336 1860 t340 = my_rhoa**2 1861 t341 = 0.1e1_dp/t340 1862 s_arhoa = -t337*t131*kf_arhoa/0.2e1_dp - t130*t341/0.2e1_dp 1863 t344 = t137**2 1864 t346 = 0.1e1_dp/t344*mu 1865 Fx_arhoa = 0.2e1_dp*t346*s_a*s_arhoa 1866 t350 = my_rhoa*ex_unif_arhoa 1867 1868 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 1869 e_ra(ii) = e_ra(ii) + & 1870 scale_ex*(0.2e1_dp*ex_unif_a*Fx_a + 0.2e1_dp* & 1871 t350*Fx_a + 0.2e1_dp*t140*Fx_arhoa)/0.2e1_dp + scale_ec*( & 1872 epsilon_cGGA + my_rho*epsilon_cGGArhoa) 1873 END IF 1874 1875 chirhob = -t2 - t164 1876 rsrhob = rsrhoa 1877 t368 = t179*rsrhob/0.2e1_dp + beta_2_1*rsrhob + 0.3e1_dp/ & 1878 0.2e1_dp*t183*rsrhob + t22*t20*rsrhob*t187 1879 e_c_u_0rhob = -0.2e1_dp*t171*rsrhob*t28 + t177*t368*t191 1880 t382 = t201*rsrhob/0.2e1_dp + beta_2_2*rsrhob + 0.3e1_dp/ & 1881 0.2e1_dp*t205*rsrhob + t40*t38*rsrhob*t187 1882 e_c_u_1rhob = -0.2e1_dp*t194*rsrhob*t46 + t200*t382*t212 1883 t396 = t222*rsrhob/0.2e1_dp + beta_2_3*rsrhob + 0.3e1_dp/ & 1884 0.2e1_dp*t226*rsrhob + t58*t56*rsrhob*t187 1885 alpha_crhob = 0.2e1_dp*t215*rsrhob*t64 - t221*t396*t233 1886 frhob = (0.4e1_dp/0.3e1_dp*t71*chirhob - 0.4e1_dp/0.3e1_dp & 1887 *t74*chirhob)*t69 1888 t403 = alpha_crhob*f 1889 t405 = alpha_c*frhob 1890 t407 = t245*chirhob 1891 t409 = 0.4e1_dp*t77*t407 1892 t410 = e_c_u_1rhob - e_c_u_0rhob 1893 t411 = t410*f 1894 t413 = t84*frhob 1895 t415 = t244*chirhob 1896 t417 = 0.4e1_dp*t85*t415 1897 epsilon_c_unifrhob = e_c_u_0rhob + t403*t82 + t405*t82 - t409 & 1898 + t411*t80 + t413*t80 + t417 1899 phirhob = t257*chirhob/0.3e1_dp - t259*chirhob/0.3e1_dp 1900 k_frhob = k_frhoa 1901 k_srhob = t266*k_frhob*t7 1902 trhob = -t269*t98*phirhob/0.2e1_dp - t96*t274*k_srhob/ & 1903 0.2e1_dp - t278/0.2e1_dp 1904 t427 = epsilon_c_unifrhob*t100 1905 t429 = t285*phirhob 1906 t432 = -t427*t105 + 0.3e1_dp*t102*t429 1907 Arhob = -t101*t281*t432*t107 1908 t436 = t123*phirhob 1909 t439 = t298*trhob 1910 t442 = Arhob*t111 1911 t444 = 0.2e1_dp*t303*trhob 1912 t445 = t442 + t444 1913 t453 = t442 + t444 + 0.2e1_dp*t314*Arhob + 0.4e1_dp*t318*trhob 1914 t456 = 0.2e1_dp*t297*t439 + t101*t111*t445*t119 - t310* & 1915 t313*t453 1916 t457 = t456*t325 1917 epsilon_cGGArhob = epsilon_c_unifrhob + 0.3e1_dp*t293*t436 + & 1918 t110*t457 1919 t460 = t143**2 1920 kf_brhob = t124/t460*t90/0.3e1_dp 1921 ex_unif_brhob = -0.3e1_dp/0.4e1_dp*t7*kf_brhob 1922 t466 = kf_b**2 1923 t467 = 0.1e1_dp/t466 1924 t468 = my_norm_drhob*t467 1925 t471 = my_rhob**2 1926 t472 = 0.1e1_dp/t471 1927 s_brhob = -t468*t148*kf_brhob/0.2e1_dp - t147*t472/0.2e1_dp 1928 t475 = t153**2 1929 t477 = 0.1e1_dp/t475*mu 1930 Fx_brhob = 0.2e1_dp*t477*s_b*s_brhob 1931 t481 = my_rhob*ex_unif_brhob 1932 1933 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 1934 e_rb(ii) = e_rb(ii) + & 1935 scale_ex*(0.2e1_dp*ex_unif_b*Fx_b + 0.2e1_dp* & 1936 t481*Fx_b + 0.2e1_dp*t156*Fx_brhob)/0.2e1_dp + scale_ec*( & 1937 epsilon_cGGA + my_rho*epsilon_cGGArhob) 1938 END IF 1939 1940 t488 = t95*t97 1941 tnorm_drho = t488*t2/0.2e1_dp 1942 t493 = t101*t317 1943 t494 = A*tnorm_drho 1944 t502 = 0.2e1_dp*t303*tnorm_drho + 0.4e1_dp*t318*tnorm_drho 1945 t505 = 0.2e1_dp*t297*t298*tnorm_drho + 0.2e1_dp*t493* & 1946 t494*t119 - t310*t313*t502 1947 t506 = t505*t325 1948 Hnorm_drho = t110*t506 1949 1950 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 1951 e_ndr(ii) = e_ndr(ii) + & 1952 scale_ec*my_rho*Hnorm_drho 1953 END IF 1954 1955 s_anorm_drhoa = t129*t131/0.2e1_dp 1956 Fx_anorm_drhoa = 0.2e1_dp*t346*s_a*s_anorm_drhoa 1957 1958 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 1959 e_ndra(ii) = e_ndra(ii) + & 1960 scale_ex*t140*Fx_anorm_drhoa 1961 END IF 1962 1963 s_bnorm_drhob = t146*t148/0.2e1_dp 1964 Fx_bnorm_drhob = 0.2e1_dp*t477*s_b*s_bnorm_drhob 1965 1966 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 1967 e_ndrb(ii) = e_ndrb(ii) + & 1968 scale_ex*t156*Fx_bnorm_drhob 1969 END IF 1970 1971 IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN 1972 t518 = 0.1e1_dp/t162/my_rho 1973 t519 = t1*t518 1974 chirhoarhoa = -0.2e1_dp*t163 + 0.2e1_dp*t519 1975 t523 = 0.1e1_dp/t90 1976 t525 = t162**2 1977 rsrhoarhoa = -t6/t165/t8*t523/t525/0.18e2_dp + & 1978 t6*t167*t518/0.6e1_dp 1979 t536 = alpha_1_1*rsrhoa 1980 t538 = t176*t190*t191 1981 t543 = t12/t175/t23 1982 t544 = t190**2 1983 t548 = 0.1e1_dp/t18 1984 t549 = beta_1_1*t548 1985 t550 = rsrhoa**2 1986 t556 = beta_3_1*t178 1987 t561 = t20**2 1988 t563 = rs**2 1989 t564 = 0.1e1_dp/t563 1990 t576 = t175**2 1991 t578 = t12/t576 1992 t579 = t27**2 1993 t580 = 0.1e1_dp/t579 1994 e_c_u_0rhoarhoa = -0.2e1_dp*t171*rsrhoarhoa*t28 + 0.2e1_dp* & 1995 t536*t538 - 0.2e1_dp*t543*t544*t191 + t177*(-t549*t550 & 1996 /0.4e1_dp + t179*rsrhoarhoa/0.2e1_dp + beta_2_1*rsrhoarhoa + & 1997 0.3e1_dp/0.4e1_dp*t556*t550 + 0.3e1_dp/0.2e1_dp*t183* & 1998 rsrhoarhoa + t22*t561*t550*t564 + t22*t20*rsrhoarhoa* & 1999 t187 - t22*t20*t550*t564)*t191 + t578*t544*t580*t14/ & 2000 0.2e1_dp 2001 e_c_u_01rhoa = e_c_u_0rhoa 2002 t588 = alpha_1_2*rsrhoa 2003 t590 = t199*t211*t212 2004 t595 = t32/t198/t41 2005 t596 = t211**2 2006 t600 = beta_1_2*t548 2007 t606 = beta_3_2*t178 2008 t611 = t38**2 2009 t624 = t198**2 2010 t626 = t32/t624 2011 t627 = t45**2 2012 t628 = 0.1e1_dp/t627 2013 t636 = alpha_1_3*rsrhoa 2014 t638 = t220*t232*t233 2015 t643 = t50/t219/t59 2016 t644 = t232**2 2017 t648 = beta_1_3*t548 2018 t654 = beta_3_3*t178 2019 t659 = t56**2 2020 t672 = t219**2 2021 t674 = t50/t672 2022 t675 = t63**2 2023 t676 = 0.1e1_dp/t675 2024 alpha_c1rhoa = alpha_crhoa 2025 t681 = 0.1e1_dp/t87 2026 t682 = chirhoa**2 2027 t687 = 0.1e1_dp/t88 2028 frhoarhoa = (0.4e1_dp/0.9e1_dp*t681*t682 + 0.4e1_dp/ & 2029 0.3e1_dp*t71*chirhoarhoa + 0.4e1_dp/0.9e1_dp*t687*t682 - & 2030 0.4e1_dp/0.3e1_dp*t74*chirhoarhoa)*t69 2031 f1rhoa = frhoa 2032 t705 = alpha_c1rhoa*f 2033 t708 = alpha_c*f1rhoa 2034 t711 = t78*t79 2035 t726 = e_c_u_1rhoa - e_c_u_01rhoa 2036 t733 = t726*f 2037 t736 = t84*f1rhoa 2038 t745 = -0.4e1_dp*t77*t245*chirhoarhoa + (-0.2e1_dp*t194* & 2039 rsrhoarhoa*t46 + 0.2e1_dp*t588*t590 - 0.2e1_dp*t595*t596* & 2040 t212 + t200*(-t600*t550/0.4e1_dp + t201*rsrhoarhoa/ & 2041 0.2e1_dp + beta_2_2*rsrhoarhoa + 0.3e1_dp/0.4e1_dp*t606*t550 & 2042 + 0.3e1_dp/0.2e1_dp*t205*rsrhoarhoa + t40*t611*t550* & 2043 t564 + t40*t38*rsrhoarhoa*t187 - t40*t38*t550*t564)* & 2044 t212 + t626*t596*t628*t34/0.2e1_dp - e_c_u_0rhoarhoa)*f* & 2045 t80 + t249*f1rhoa*t80 + 0.4e1_dp*t250*t254 + t726*frhoa* & 2046 t80 + t84*frhoarhoa*t80 + 0.4e1_dp*t252*t254 + 0.4e1_dp* & 2047 t733*t254 + 0.4e1_dp*t736*t254 + 0.12e2_dp*t85*t79*t682 & 2048 + 0.4e1_dp*t85*t244*chirhoarhoa 2049 epsilon_c_unifrhoarhoa = e_c_u_0rhoarhoa + (0.2e1_dp*t215* & 2050 rsrhoarhoa*t64 - 0.2e1_dp*t636*t638 + 0.2e1_dp*t643*t644* & 2051 t233 - t221*(-t648*t550/0.4e1_dp + t222*rsrhoarhoa/ & 2052 0.2e1_dp + beta_2_3*rsrhoarhoa + 0.3e1_dp/0.4e1_dp*t654*t550 & 2053 + 0.3e1_dp/0.2e1_dp*t226*rsrhoarhoa + t58*t659*t550* & 2054 t564 + t58*t56*rsrhoarhoa*t187 - t58*t56*t550*t564)* & 2055 t233 - t674*t644*t676*t52/0.2e1_dp)*f*t82 + alpha_crhoa & 2056 *f1rhoa*t82 - 0.4e1_dp*t240*t246 + alpha_c1rhoa*frhoa*t82 & 2057 + alpha_c*frhoarhoa*t82 - 0.4e1_dp*t242*t246 - 0.4e1_dp* & 2058 t705*t246 - 0.4e1_dp*t708*t246 - 0.12e2_dp*t77*t711*t682 & 2059 + t745 2060 epsilon_c_unif1rhoa = e_c_u_01rhoa + t705*t82 + t708*t82 - & 2061 t248 + t733*t80 + t736*t80 + t256 2062 t750 = 0.1e1_dp/t72 2063 t755 = 0.1e1_dp/t75 2064 phirhoarhoa = -t750*t682/0.9e1_dp + t257*chirhoarhoa/ & 2065 0.3e1_dp - t755*t682/0.9e1_dp - t259*chirhoarhoa/0.3e1_dp 2066 phi1rhoa = phirhoa 2067 t763 = t90**2 2068 k_frhoarhoa = -0.2e1_dp/0.9e1_dp*t3/t262/t91*t763 2069 t767 = 0.1e1_dp/t94/t93 2070 t768 = k_frhoa**2 2071 k_s1rhoa = k_srhoa 2072 t775 = my_norm_drho*t105*t97 2073 t776 = t2*phirhoa 2074 t779 = t269*t273 2075 t785 = t269*t277*phirhoa/0.2e1_dp 2076 t789 = t2*k_srhoa 2077 t795 = t96/t272/k_s 2078 t798 = t273*t163 2079 t801 = t96*t798*k_srhoa/0.2e1_dp 2080 t812 = t96*t97*t518 2081 trhoarhoa = t775*t776*phi1rhoa + t779*t776*k_s1rhoa/ & 2082 0.2e1_dp + t785 - t269*t98*phirhoarhoa/0.2e1_dp + t779*t789 & 2083 *phi1rhoa/0.2e1_dp + t795*t789*k_s1rhoa + t801 - t96*t274* & 2084 (-t767*t768*t523/0.2e1_dp + t266*k_frhoarhoa*t7)/ & 2085 0.2e1_dp + t269*t277*phi1rhoa/0.2e1_dp + t96*t798*k_s1rhoa & 2086 /0.2e1_dp + t812 2087 t1rhoa = -t269*t98*phi1rhoa/0.2e1_dp - t96*t274*k_s1rhoa & 2088 /0.2e1_dp - t278/0.2e1_dp 2089 t820 = t101/t280/t108 2090 t821 = t107**2 2091 t822 = t289*t821 2092 t823 = epsilon_c_unif1rhoa*t100 2093 t825 = t285*phi1rhoa 2094 t828 = -t823*t105 + 0.3e1_dp*t102*t825 2095 t839 = 0.1e1_dp/t284/phi 2096 t840 = t839*phirhoa 2097 t851 = t101*t281 2098 Arhoarhoa = 0.2e1_dp*t820*t822*t828 - t101*t281*( & 2099 -epsilon_c_unifrhoarhoa*t100*t105 + 0.3e1_dp*t282*t825 + & 2100 0.3e1_dp*t823*t286 - 0.12e2_dp*t102*t840*phi1rhoa + & 2101 0.3e1_dp*t102*t285*phirhoarhoa)*t107 - t851*t289*t828* & 2102 t107 2103 A1rhoa = -t101*t281*t828*t107 2104 t858 = gamma_var*phi 2105 t865 = A1rhoa*t111 2106 t867 = 0.2e1_dp*t303*t1rhoa 2107 t868 = t865 + t867 2108 t876 = t865 + t867 + 0.2e1_dp*t314*A1rhoa + 0.4e1_dp*t318*t1rhoa 2109 t879 = 0.2e1_dp*t297*t298*t1rhoa + t101*t111*t868*t119 & 2110 - t310*t313*t876 2111 t880 = t879*t325 2112 t904 = t306*t119 2113 t908 = Arhoarhoa*t111 2114 t909 = Arhoa*t 2115 t911 = 0.2e1_dp*t909*t1rhoa 2116 t914 = 0.2e1_dp*A1rhoa*t*trhoa 2117 t917 = 0.2e1_dp*A*t1rhoa*trhoa 2118 t919 = 0.2e1_dp*t303*trhoarhoa 2119 t924 = t306*t312 2120 t936 = t113/t311/t118 2121 t944 = A*t317 2122 t953 = t115*t111 2123 t959 = t908 + t911 + t914 + t917 + t919 + 0.2e1_dp*A1rhoa*t116 & 2124 *Arhoa + 0.8e1_dp*t944*Arhoa*t1rhoa + 0.2e1_dp*t314* & 2125 Arhoarhoa + 0.8e1_dp*t944*trhoa*A1rhoa + 0.12e2_dp*t953* & 2126 trhoa*t1rhoa + 0.4e1_dp*t318*trhoarhoa 2127 t962 = 0.2e1_dp*t101*t1rhoa*t299 + 0.2e1_dp*t297*t868* & 2128 t119*trhoa - 0.2e1_dp*t297*t313*trhoa*t876 + 0.2e1_dp* & 2129 t297*t298*trhoarhoa + 0.2e1_dp*t297*t904*t1rhoa + t101* & 2130 t111*(t908 + t911 + t914 + t917 + t919)*t119 - t310*t924* & 2131 t876 - 0.2e1_dp*t297*t313*t321*t1rhoa - t310*t868*t312* & 2132 t321 + 0.2e1_dp*t310*t936*t321*t876 - t310*t313*t959 2133 t965 = t122**2 2134 t966 = 0.1e1_dp/t965 2135 t967 = t324*t966 2136 kf_arhoarhoa = -0.2e1_dp/0.9e1_dp*t124/t329/t125*t763 2137 ex_unif_a1rhoa = ex_unif_arhoa 2138 t985 = kf_arhoa**2 2139 s_a1rhoa = s_arhoa 2140 t998 = mu**2 2141 t999 = 0.1e1_dp/t344/t137*t998 2142 t1000 = t999*t133 2143 t1001 = s_arhoa*t135 2144 Fx_a1rhoa = 0.2e1_dp*t346*s_a*s_a1rhoa 2145 2146 e_ra_ra(ii) = e_ra_ra(ii) + & 2147 scale_ex*(0.2e1_dp*ex_unif_a1rhoa*Fx_a + & 2148 0.2e1_dp*ex_unif_a*Fx_a1rhoa + 0.2e1_dp*ex_unif_arhoa*Fx_a - & 2149 0.3e1_dp/0.2e1_dp*my_rhoa*t7*kf_arhoarhoa*Fx_a + 0.2e1_dp* & 2150 t350*Fx_a1rhoa + 0.2e1_dp*ex_unif_a*Fx_arhoa + 0.2e1_dp*my_rhoa & 2151 *ex_unif_a1rhoa*Fx_arhoa + 0.2e1_dp*t140*(-0.8e1_dp*t1000 & 2152 *t1001*s_a1rhoa + 0.2e1_dp*t346*s_a1rhoa*s_arhoa + 0.2e1_dp & 2153 *t346*s_a*(my_norm_drhoa/t335/kf_a*t131*t985 + t337* & 2154 t341*kf_arhoa - t337*t131*kf_arhoarhoa/0.2e1_dp + t130/ & 2155 t340/my_rhoa)))/0.2e1_dp + scale_ec*(epsilon_c_unif1rhoa + & 2156 0.3e1_dp*t293*t123*phi1rhoa + t110*t880 + epsilon_cGGArhoa + & 2157 my_rho*(epsilon_c_unifrhoarhoa + 0.6e1_dp*t858*t294*phi1rhoa + & 2158 0.3e1_dp*t293*t880*phirhoa + 0.3e1_dp*t293*t123* & 2159 phirhoarhoa + 0.3e1_dp*t293*t326*phi1rhoa + t110*t962*t325 & 2160 - t110*t967*t879)) 2161 2162 chirhoarhob = 0.2e1_dp*t519 2163 rsrhoarhob = rsrhoarhoa 2164 t1031 = t176*t368*t191 2165 t1033 = alpha_1_1*rsrhob 2166 t1038 = rsrhoa*rsrhob 2167 t1050 = rsrhob*t564*rsrhoa 2168 e_c_u_0rhoarhob = -0.2e1_dp*t171*rsrhoarhob*t28 + t536* & 2169 t1031 + t1033*t538 - 0.2e1_dp*t543*t192*t368 + t177*(-t549 & 2170 *t1038/0.4e1_dp + t179*rsrhoarhob/0.2e1_dp + beta_2_1* & 2171 rsrhoarhob + 0.3e1_dp/0.4e1_dp*t556*t1038 + 0.3e1_dp/ & 2172 0.2e1_dp*t183*rsrhoarhob + t22*t561*t1050 + t22*t20* & 2173 rsrhoarhob*t187 - t22*t20*t1050)*t191 + t578*t190*t580* & 2174 t14*t368/0.2e1_dp 2175 t1069 = t199*t382*t212 2176 t1071 = alpha_1_2*rsrhob 2177 t1104 = t220*t396*t233 2178 t1106 = alpha_1_3*rsrhob 2179 frhoarhob = (0.4e1_dp/0.9e1_dp*t681*chirhoa*chirhob + & 2180 0.4e1_dp/0.3e1_dp*t71*chirhoarhob + 0.4e1_dp/0.9e1_dp*t687 & 2181 *chirhoa*chirhob - 0.4e1_dp/0.3e1_dp*t74*chirhoarhob)* & 2182 t69 2183 t1164 = t79*chirhoa*chirhob 2184 t1193 = -0.4e1_dp*t77*t245*chirhoarhob + (-0.2e1_dp*t194* & 2185 rsrhoarhob*t46 + t588*t1069 + t1071*t590 - 0.2e1_dp*t595* & 2186 t213*t382 + t200*(-t600*t1038/0.4e1_dp + t201*rsrhoarhob/ & 2187 0.2e1_dp + beta_2_2*rsrhoarhob + 0.3e1_dp/0.4e1_dp*t606* & 2188 t1038 + 0.3e1_dp/0.2e1_dp*t205*rsrhoarhob + t40*t611*t1050 & 2189 + t40*t38*rsrhoarhob*t187 - t40*t38*t1050)*t212 + t626 & 2190 *t211*t628*t34*t382/0.2e1_dp - e_c_u_0rhoarhob)*f*t80 + & 2191 t249*frhob*t80 + 0.4e1_dp*t250*t415 + t410*frhoa*t80 + & 2192 t84*frhoarhob*t80 + 0.4e1_dp*t252*t415 + 0.4e1_dp*t411* & 2193 t254 + 0.4e1_dp*t413*t254 + 0.12e2_dp*t85*t1164 + 0.4e1_dp* & 2194 t85*t244*chirhoarhob 2195 epsilon_c_unifrhoarhob = e_c_u_0rhoarhob + (0.2e1_dp*t215* & 2196 rsrhoarhob*t64 - t636*t1104 - t1106*t638 + 0.2e1_dp*t643* & 2197 t234*t396 - t221*(-t648*t1038/0.4e1_dp + t222*rsrhoarhob/ & 2198 0.2e1_dp + beta_2_3*rsrhoarhob + 0.3e1_dp/0.4e1_dp*t654* & 2199 t1038 + 0.3e1_dp/0.2e1_dp*t226*rsrhoarhob + t58*t659*t1050 & 2200 + t58*t56*rsrhoarhob*t187 - t58*t56*t1050)*t233 - t674 & 2201 *t232*t676*t52*t396/0.2e1_dp)*f*t82 + alpha_crhoa* & 2202 frhob*t82 - 0.4e1_dp*t240*t407 + alpha_crhob*frhoa*t82 + & 2203 alpha_c*frhoarhob*t82 - 0.4e1_dp*t242*t407 - 0.4e1_dp*t403 & 2204 *t246 - 0.4e1_dp*t405*t246 - 0.12e2_dp*t77*t78*t1164 + & 2205 t1193 2206 phirhoarhob = -t750*chirhoa*chirhob/0.9e1_dp + t257* & 2207 chirhoarhob/0.3e1_dp - t755*chirhoa*chirhob/0.9e1_dp - t259 & 2208 *chirhoarhob/0.3e1_dp 2209 k_frhoarhob = k_frhoarhoa 2210 t1228 = t269*t277*phirhob/0.2e1_dp 2211 t1231 = t96*t798*k_srhob/0.2e1_dp 2212 trhoarhob = t775*t776*phirhob + t779*t776*k_srhob/ & 2213 0.2e1_dp + t785 - t269*t98*phirhoarhob/0.2e1_dp + t779*t789 & 2214 *phirhob/0.2e1_dp + t795*t789*k_srhob + t801 - t96*t274*( & 2215 -t767*k_frhoa*t523*k_frhob/0.2e1_dp + t266*k_frhoarhob* & 2216 t7)/0.2e1_dp + t1228 + t1231 + t812 2217 Arhoarhob = 0.2e1_dp*t820*t822*t432 - t101*t281*( & 2218 -epsilon_c_unifrhoarhob*t100*t105 + 0.3e1_dp*t282*t429 + & 2219 0.3e1_dp*t427*t286 - 0.12e2_dp*t102*t840*phirhob + & 2220 0.3e1_dp*t102*t285*phirhoarhob)*t107 - t851*t289*t432* & 2221 t107 2222 t1269 = t445*t119 2223 t1283 = Arhoarhob*t111 2224 t1285 = 0.2e1_dp*t909*trhob 2225 t1286 = Arhob*t 2226 t1288 = 0.2e1_dp*t1286*trhoa 2227 t1291 = 0.2e1_dp*A*trhob*trhoa 2228 t1293 = 0.2e1_dp*t303*trhoarhob 2229 t1304 = t445*t312 2230 t1327 = t1283 + t1285 + t1288 + t1291 + t1293 + 0.2e1_dp*Arhob* & 2231 t116*Arhoa + 0.8e1_dp*t944*Arhoa*trhob + 0.2e1_dp*t314* & 2232 Arhoarhob + 0.8e1_dp*t944*trhoa*Arhob + 0.12e2_dp*t953* & 2233 trhoa*trhob + 0.4e1_dp*t318*trhoarhob 2234 t1330 = 0.2e1_dp*t101*trhob*t299 + 0.2e1_dp*t297*t1269* & 2235 trhoa - 0.2e1_dp*t297*t313*trhoa*t453 + 0.2e1_dp*t297* & 2236 t298*trhoarhob + 0.2e1_dp*t297*t904*trhob + t101*t111*( & 2237 t1283 + t1285 + t1288 + t1291 + t1293)*t119 - t310*t924*t453 - & 2238 0.2e1_dp*t297*t313*t321*trhob - t310*t1304*t321 + & 2239 0.2e1_dp*t310*t936*t321*t453 - t310*t313*t1327 2240 2241 e_ra_rb(ii) = e_ra_rb(ii) + & 2242 scale_ec*(epsilon_cGGArhob + epsilon_cGGArhoa + & 2243 my_rho*(epsilon_c_unifrhoarhob + 0.6e1_dp*t858*t294*phirhob + & 2244 0.3e1_dp*t293*t457*phirhoa + 0.3e1_dp*t293*t123* & 2245 phirhoarhob + 0.3e1_dp*t293*t326*phirhob + t110*t1330*t325 & 2246 - t110*t967*t456)) 2247 2248 chirhobrhob = 0.2e1_dp*t163 + 0.2e1_dp*t519 2249 rsrhobrhob = rsrhoarhob 2250 t1342 = t368**2 2251 t1346 = rsrhob**2 2252 e_c_u_0rhobrhob = -0.2e1_dp*t171*rsrhobrhob*t28 + 0.2e1_dp* & 2253 t1033*t1031 - 0.2e1_dp*t543*t1342*t191 + t177*(-t549* & 2254 t1346/0.4e1_dp + t179*rsrhobrhob/0.2e1_dp + beta_2_1* & 2255 rsrhobrhob + 0.3e1_dp/0.4e1_dp*t556*t1346 + 0.3e1_dp/ & 2256 0.2e1_dp*t183*rsrhobrhob + t22*t561*t1346*t564 + t22*t20 & 2257 *rsrhobrhob*t187 - t22*t20*t1346*t564)*t191 + t578* & 2258 t1342*t580*t14/0.2e1_dp 2259 e_c_u_01rhob = e_c_u_0rhob 2260 t1377 = t382**2 2261 t1411 = t396**2 2262 alpha_c1rhob = alpha_crhob 2263 t1440 = chirhob**2 2264 frhobrhob = (0.4e1_dp/0.9e1_dp*t681*t1440 + 0.4e1_dp/ & 2265 0.3e1_dp*t71*chirhobrhob + 0.4e1_dp/0.9e1_dp*t687*t1440 - & 2266 0.4e1_dp/0.3e1_dp*t74*chirhobrhob)*t69 2267 f1rhob = frhob 2268 t1462 = alpha_c1rhob*f 2269 t1465 = alpha_c*f1rhob 2270 t1482 = e_c_u_1rhob - e_c_u_01rhob 2271 t1489 = t1482*f 2272 t1492 = t84*f1rhob 2273 t1501 = -0.4e1_dp*t77*t245*chirhobrhob + (-0.2e1_dp*t194* & 2274 rsrhobrhob*t46 + 0.2e1_dp*t1071*t1069 - 0.2e1_dp*t595* & 2275 t1377*t212 + t200*(-t600*t1346/0.4e1_dp + t201*rsrhobrhob & 2276 /0.2e1_dp + beta_2_2*rsrhobrhob + 0.3e1_dp/0.4e1_dp*t606* & 2277 t1346 + 0.3e1_dp/0.2e1_dp*t205*rsrhobrhob + t40*t611*t1346 & 2278 *t564 + t40*t38*rsrhobrhob*t187 - t40*t38*t1346*t564) & 2279 *t212 + t626*t1377*t628*t34/0.2e1_dp - e_c_u_0rhobrhob)*f & 2280 *t80 + t410*f1rhob*t80 + 0.4e1_dp*t411*t415 + t1482* & 2281 frhob*t80 + t84*frhobrhob*t80 + 0.4e1_dp*t413*t415 + & 2282 0.4e1_dp*t1489*t415 + 0.4e1_dp*t1492*t415 + 0.12e2_dp*t85 & 2283 *t79*t1440 + 0.4e1_dp*t85*t244*chirhobrhob 2284 epsilon_c_unifrhobrhob = e_c_u_0rhobrhob + (0.2e1_dp*t215* & 2285 rsrhobrhob*t64 - 0.2e1_dp*t1106*t1104 + 0.2e1_dp*t643* & 2286 t1411*t233 - t221*(-t648*t1346/0.4e1_dp + t222*rsrhobrhob & 2287 /0.2e1_dp + beta_2_3*rsrhobrhob + 0.3e1_dp/0.4e1_dp*t654* & 2288 t1346 + 0.3e1_dp/0.2e1_dp*t226*rsrhobrhob + t58*t659*t1346 & 2289 *t564 + t58*t56*rsrhobrhob*t187 - t58*t56*t1346*t564) & 2290 *t233 - t674*t1411*t676*t52/0.2e1_dp)*f*t82 + & 2291 alpha_crhob*f1rhob*t82 - 0.4e1_dp*t403*t407 + alpha_c1rhob* & 2292 frhob*t82 + alpha_c*frhobrhob*t82 - 0.4e1_dp*t405*t407 - & 2293 0.4e1_dp*t1462*t407 - 0.4e1_dp*t1465*t407 - 0.12e2_dp*t77 & 2294 *t711*t1440 + t1501 2295 epsilon_c_unif1rhob = e_c_u_01rhob + t1462*t82 + t1465*t82 - & 2296 t409 + t1489*t80 + t1492*t80 + t417 2297 phirhobrhob = -t750*t1440/0.9e1_dp + t257*chirhobrhob/ & 2298 0.3e1_dp - t755*t1440/0.9e1_dp - t259*chirhobrhob/0.3e1_dp 2299 phi1rhob = phirhob 2300 t1514 = k_frhob**2 2301 k_s1rhob = k_srhob 2302 t1520 = t2*phirhob 2303 t1529 = t2*k_srhob 2304 trhobrhob = t775*t1520*phi1rhob + t779*t1520*k_s1rhob/ & 2305 0.2e1_dp + t1228 - t269*t98*phirhobrhob/0.2e1_dp + t779* & 2306 t1529*phi1rhob/0.2e1_dp + t795*t1529*k_s1rhob + t1231 - t96 & 2307 *t274*(-t767*t1514*t523/0.2e1_dp + t266*k_frhoarhob*t7) & 2308 /0.2e1_dp + t269*t277*phi1rhob/0.2e1_dp + t96*t798* & 2309 k_s1rhob/0.2e1_dp + t812 2310 t1rhob = -t269*t98*phi1rhob/0.2e1_dp - t96*t274*k_s1rhob & 2311 /0.2e1_dp - t278/0.2e1_dp 2312 t1550 = epsilon_c_unif1rhob*t100 2313 t1552 = t285*phi1rhob 2314 t1555 = -t1550*t105 + 0.3e1_dp*t102*t1552 2315 Arhobrhob = 0.2e1_dp*t820*t432*t821*t1555 - t101*t281* & 2316 (-epsilon_c_unifrhobrhob*t100*t105 + 0.3e1_dp*t427*t1552 + & 2317 0.3e1_dp*t1550*t429 - 0.12e2_dp*t102*t839*phirhob* & 2318 phi1rhob + 0.3e1_dp*t102*t285*phirhobrhob)*t107 - t851* & 2319 t432*t1555*t107 2320 A1rhob = -t101*t281*t1555*t107 2321 t1588 = A1rhob*t111 2322 t1590 = 0.2e1_dp*t303*t1rhob 2323 t1591 = t1588 + t1590 2324 t1599 = t1588 + t1590 + 0.2e1_dp*t314*A1rhob + 0.4e1_dp*t318 & 2325 *t1rhob 2326 t1602 = 0.2e1_dp*t297*t298*t1rhob + t101*t111*t1591* & 2327 t119 - t310*t313*t1599 2328 t1603 = t1602*t325 2329 t1630 = Arhobrhob*t111 2330 t1632 = 0.2e1_dp*t1286*t1rhob 2331 t1635 = 0.2e1_dp*A1rhob*t*trhob 2332 t1638 = 0.2e1_dp*A*t1rhob*trhob 2333 t1640 = 0.2e1_dp*t303*trhobrhob 2334 t1674 = t1630 + t1632 + t1635 + t1638 + t1640 + 0.2e1_dp*A1rhob & 2335 *t116*Arhob + 0.8e1_dp*t944*Arhob*t1rhob + 0.2e1_dp*t314 & 2336 *Arhobrhob + 0.8e1_dp*t944*trhob*A1rhob + 0.12e2_dp*t953* & 2337 trhob*t1rhob + 0.4e1_dp*t318*trhobrhob 2338 t1677 = 0.2e1_dp*t101*t1rhob*t439 + 0.2e1_dp*t297*t1591 & 2339 *t119*trhob - 0.2e1_dp*t297*t313*trhob*t1599 + 0.2e1_dp* & 2340 t297*t298*trhobrhob + 0.2e1_dp*t297*t1269*t1rhob + t101* & 2341 t111*(t1630 + t1632 + t1635 + t1638 + t1640)*t119 - t310* & 2342 t1304*t1599 - 0.2e1_dp*t297*t313*t453*t1rhob - t310* & 2343 t1591*t312*t453 + 0.2e1_dp*t310*t936*t453*t1599 - t310* & 2344 t313*t1674 2345 t1680 = t456*t966 2346 kf_brhobrhob = -0.2e1_dp/0.9e1_dp*t124/t460/t142*t763 2347 ex_unif_b1rhob = ex_unif_brhob 2348 t1698 = kf_brhob**2 2349 s_b1rhob = s_brhob 2350 t1711 = 0.1e1_dp/t475/t153*t998 2351 t1712 = t1711*t150 2352 t1713 = s_brhob*t135 2353 Fx_b1rhob = 0.2e1_dp*t477*s_b*s_b1rhob 2354 2355 e_rb_rb(ii) = e_rb_rb(ii) + & 2356 scale_ex*(0.2e1_dp*ex_unif_b1rhob*Fx_b + & 2357 0.2e1_dp*ex_unif_b*Fx_b1rhob + 0.2e1_dp*ex_unif_brhob*Fx_b - & 2358 0.3e1_dp/0.2e1_dp*my_rhob*t7*kf_brhobrhob*Fx_b + 0.2e1_dp* & 2359 t481*Fx_b1rhob + 0.2e1_dp*ex_unif_b*Fx_brhob + 0.2e1_dp*my_rhob & 2360 *ex_unif_b1rhob*Fx_brhob + 0.2e1_dp*t156*(-0.8e1_dp*t1712 & 2361 *t1713*s_b1rhob + 0.2e1_dp*t477*s_b1rhob*s_brhob + 0.2e1_dp & 2362 *t477*s_b*(my_norm_drhob/t466/kf_b*t148*t1698 + t468* & 2363 t472*kf_brhob - t468*t148*kf_brhobrhob/0.2e1_dp + t147/ & 2364 t471/my_rhob)))/0.2e1_dp + scale_ec*(epsilon_c_unif1rhob + & 2365 0.3e1_dp*t293*t123*phi1rhob + t110*t1603 + epsilon_cGGArhob & 2366 + my_rho*(epsilon_c_unifrhobrhob + 0.6e1_dp*t858*t436*phi1rhob & 2367 + 0.3e1_dp*t293*t1603*phirhob + 0.3e1_dp*t293*t123* & 2368 phirhobrhob + 0.3e1_dp*t293*t457*phi1rhob + t110*t1677* & 2369 t325 - t110*t1680*t1602)) 2370 t1739 = t268*t97 2371 t1741 = t95*t273 2372 t1743 = t488*t163 2373 trhoanorm_drho = -t1739*t776/0.2e1_dp - t1741*t789/ & 2374 0.2e1_dp - t1743/0.2e1_dp 2375 t1748 = t101*tnorm_drho 2376 t1765 = t909*tnorm_drho 2377 t1766 = t494*trhoa 2378 t1767 = t303*trhoanorm_drho 2379 t1801 = 0.2e1_dp*t1748*t299 + 0.4e1_dp*t310*t494*t119* & 2380 trhoa - 0.2e1_dp*t297*t313*trhoa*t502 + 0.2e1_dp*t297* & 2381 t298*trhoanorm_drho + 0.2e1_dp*t297*t904*tnorm_drho + t101* & 2382 t111*(0.2e1_dp*t1765 + 0.2e1_dp*t1766 + 0.2e1_dp*t1767)* & 2383 t119 - t310*t924*t502 - 0.2e1_dp*t297*t313*t321* & 2384 tnorm_drho - 0.2e1_dp*t493*t494*t312*t321 + 0.2e1_dp*t310 & 2385 *t936*t321*t502 - t310*t313*(0.2e1_dp*t1765 + 0.2e1_dp* & 2386 t1766 + 0.2e1_dp*t1767 + 0.8e1_dp*t944*Arhoa*tnorm_drho + & 2387 0.12e2_dp*t953*trhoa*tnorm_drho + 0.4e1_dp*t318* & 2388 trhoanorm_drho) 2389 2390 e_ra_ndr(ii) = e_ra_ndr(ii) + & 2391 scale_ec*(Hnorm_drho + my_rho*(0.3e1_dp* & 2392 t293*t506*phirhoa + t110*t1801*t325 - t110*t967*t505)) 2393 2394 trhobnorm_drho = -t1739*t1520/0.2e1_dp - t1741*t1529/ & 2395 0.2e1_dp - t1743/0.2e1_dp 2396 t1829 = t1286*tnorm_drho 2397 t1830 = t494*trhob 2398 t1831 = t303*trhobnorm_drho 2399 t1865 = 0.2e1_dp*t1748*t439 + 0.4e1_dp*t310*t494*t119* & 2400 trhob - 0.2e1_dp*t297*t313*trhob*t502 + 0.2e1_dp*t297* & 2401 t298*trhobnorm_drho + 0.2e1_dp*t297*t1269*tnorm_drho + t101 & 2402 *t111*(0.2e1_dp*t1829 + 0.2e1_dp*t1830 + 0.2e1_dp*t1831)* & 2403 t119 - t310*t1304*t502 - 0.2e1_dp*t297*t313*t453* & 2404 tnorm_drho - 0.2e1_dp*t493*t494*t312*t453 + 0.2e1_dp*t310 & 2405 *t936*t453*t502 - t310*t313*(0.2e1_dp*t1829 + 0.2e1_dp* & 2406 t1830 + 0.2e1_dp*t1831 + 0.8e1_dp*t944*Arhob*tnorm_drho + & 2407 0.12e2_dp*t953*trhob*tnorm_drho + 0.4e1_dp*t318* & 2408 trhobnorm_drho) 2409 2410 e_rb_ndr(ii) = e_rb_ndr(ii) + & 2411 scale_ec*(Hnorm_drho + my_rho*(0.3e1_dp* & 2412 t293*t506*phirhob + t110*t1865*t325 - t110*t1680*t505)) 2413 2414 t1871 = tnorm_drho**2 2415 t1876 = A*t1871 2416 t1888 = t502**2 2417 t1901 = t505**2 2418 2419 e_ndr_ndr(ii) = e_ndr_ndr(ii) + & 2420 scale_ec*my_rho*(t110*(0.2e1_dp* & 2421 t101*t1871*t113*t119 + 0.10e2_dp*t310*t1876*t119 - & 2422 0.4e1_dp*t297*t313*tnorm_drho*t502 - 0.4e1_dp*t493*t494 & 2423 *t312*t502 + 0.2e1_dp*t310*t936*t1888 - t310*t313*( & 2424 0.2e1_dp*t1876 + 0.12e2_dp*t953*t1871))*t325 - t110*t1901 & 2425 *t966) 2426 2427 e_ra_ndra(ii) = e_ra_ndra(ii) + & 2428 scale_ex*(0.2e1_dp*ex_unif_a* & 2429 Fx_anorm_drhoa + 0.2e1_dp*t350*Fx_anorm_drhoa + 0.2e1_dp*t140 & 2430 *(-0.8e1_dp*t1000*t1001*s_anorm_drhoa + 0.2e1_dp*t346* & 2431 s_anorm_drhoa*s_arhoa + 0.2e1_dp*t346*s_a*(-t336*t131* & 2432 kf_arhoa/0.2e1_dp - t129*t341/0.2e1_dp)))/0.2e1_dp 2433 2434 t1922 = s_anorm_drhoa**2 2435 2436 e_ndra_ndra(ii) = e_ndra_ndra(ii) + & 2437 scale_ex*t140*(-0.8e1_dp*t999* & 2438 t133*t1922*t135 + 0.2e1_dp*t346*t1922) 2439 2440 e_rb_ndrb(ii) = e_rb_ndrb(ii) + & 2441 scale_ex*(0.2e1_dp*ex_unif_b* & 2442 Fx_bnorm_drhob + 0.2e1_dp*t481*Fx_bnorm_drhob + 0.2e1_dp*t156 & 2443 *(-0.8e1_dp*t1712*t1713*s_bnorm_drhob + 0.2e1_dp*t477* & 2444 s_bnorm_drhob*s_brhob + 0.2e1_dp*t477*s_b*(-t467*t148* & 2445 kf_brhob/0.2e1_dp - t146*t472/0.2e1_dp)))/0.2e1_dp 2446 2447 t1949 = s_bnorm_drhob**2 2448 e_ndrb_ndrb(ii) = e_ndrb_ndrb(ii) + & 2449 scale_ex*t156*(-0.8e1_dp*t1711* & 2450 t150*t1949*t135 + 0.2e1_dp*t477*t1949) 2451 END IF 2452 END IF 2453 END DO 2454!$OMP END DO 2455 END SELECT 2456 END SUBROUTINE pbe_lsd_calc 2457 2458END MODULE xc_pbe 2459