|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE xc_lyp 00013 USE bibliography, ONLY: Lee1988,& 00014 cite_reference 00015 USE f77_blas 00016 USE input_section_types, ONLY: section_vals_type,& 00017 section_vals_val_get 00018 USE kinds, ONLY: dp 00019 USE mathconstants, ONLY: pi 00020 USE timings, ONLY: timeset,& 00021 timestop 00022 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00023 xc_dset_get_derivative 00024 USE xc_derivative_types, ONLY: xc_derivative_get,& 00025 xc_derivative_type 00026 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00027 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 00028 xc_rho_set_type 00029 #include "cp_common_uses.h" 00030 00031 IMPLICIT NONE 00032 PRIVATE 00033 00034 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00035 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_lyp' 00036 REAL(kind=dp), PARAMETER, PRIVATE :: a=0.04918_dp, b=0.132_dp, 00037 c=0.2533_dp,d=0.349_dp 00038 00039 PUBLIC :: lyp_lda_info, lyp_lsd_info, lyp_lda_eval, lyp_lsd_eval 00040 CONTAINS 00041 00042 ! ***************************************************************************** 00054 SUBROUTINE lyp_lda_info(reference,shortform, needs, max_deriv, error) 00055 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00056 TYPE(xc_rho_cflags_type), 00057 INTENT(inout), OPTIONAL :: needs 00058 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00059 TYPE(cp_error_type), INTENT(inout) :: error 00060 00061 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_lda_info', 00062 routineP = moduleN//':'//routineN 00063 00064 IF ( PRESENT ( reference ) ) THEN 00065 reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LDA version}" 00066 END IF 00067 IF ( PRESENT ( shortform ) ) THEN 00068 shortform = "Lee-Yang-Parr correlation energy functional (LDA)" 00069 END IF 00070 IF (PRESENT(needs)) THEN 00071 needs%rho=.TRUE. 00072 needs%rho_1_3=.TRUE. 00073 needs%norm_drho=.TRUE. 00074 END IF 00075 IF (PRESENT(max_deriv)) max_deriv=3 00076 00077 END SUBROUTINE lyp_lda_info 00078 00079 ! ***************************************************************************** 00091 SUBROUTINE lyp_lsd_info(reference,shortform, needs, max_deriv, error) 00092 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00093 TYPE(xc_rho_cflags_type), 00094 INTENT(inout), OPTIONAL :: needs 00095 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00096 TYPE(cp_error_type), INTENT(inout) :: error 00097 00098 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_lsd_info', 00099 routineP = moduleN//':'//routineN 00100 00101 IF ( PRESENT ( reference ) ) THEN 00102 reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LSD version}" 00103 END IF 00104 IF ( PRESENT ( shortform ) ) THEN 00105 shortform = "Lee-Yang-Parr correlation energy functional (LSD)" 00106 END IF 00107 IF (PRESENT(needs)) THEN 00108 needs%rho_spin=.TRUE. 00109 needs%norm_drho_spin=.TRUE. 00110 needs%norm_drho=.TRUE. 00111 END IF 00112 IF (PRESENT(max_deriv)) max_deriv=2 00113 00114 END SUBROUTINE lyp_lsd_info 00115 00116 ! ***************************************************************************** 00132 SUBROUTINE lyp_lda_eval(rho_set,deriv_set,grad_deriv,lyp_params,error) 00133 TYPE(xc_rho_set_type), POINTER :: rho_set 00134 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00135 INTEGER, INTENT(in) :: grad_deriv 00136 TYPE(section_vals_type), POINTER :: lyp_params 00137 TYPE(cp_error_type), INTENT(inout) :: error 00138 00139 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_lda_eval', 00140 routineP = moduleN//':'//routineN 00141 00142 INTEGER :: handle, npoints, stat 00143 INTEGER, DIMENSION(:, :), POINTER :: bo 00144 LOGICAL :: failure 00145 REAL(kind=dp) :: epsilon_norm_drho, 00146 epsilon_rho, sc 00147 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, 00148 e_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, 00149 e_rho_rho, e_rho_rho_rho, norm_drho, rho, rho_1_3 00150 TYPE(xc_derivative_type), POINTER :: deriv 00151 00152 CALL timeset(routineN,handle) 00153 00154 failure=.FALSE. 00155 NULLIFY(bo) 00156 00157 CALL section_vals_val_get(lyp_params,"scale_c",r_val=sc,error=error) 00158 CALL cite_reference(Lee1988) 00159 00160 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00161 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00162 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00163 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00164 IF (.NOT. failure) THEN 00165 CALL xc_rho_set_get(rho_set,rho_1_3=rho_1_3,rho=rho,& 00166 norm_drho=norm_drho,local_bounds=bo,rho_cutoff=epsilon_rho,& 00167 drho_cutoff=epsilon_norm_drho,error=error) 00168 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00169 00170 ! meaningful default for the arrays we don't need: let us make compiler 00171 ! and debugger happy... 00172 IF (cp_debug) THEN 00173 ALLOCATE(dummy(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 00174 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00175 ELSE 00176 dummy=> rho 00177 END IF 00178 00179 e_0 => dummy 00180 e_rho => dummy 00181 e_ndrho => dummy 00182 e_rho_rho => dummy 00183 e_ndrho_rho => dummy 00184 e_ndrho_ndrho => dummy 00185 e_rho_rho_rho => dummy 00186 e_ndrho_rho_rho => dummy 00187 e_ndrho_ndrho_rho => dummy 00188 00189 IF (grad_deriv>=0) THEN 00190 deriv => xc_dset_get_derivative(deriv_set,"",& 00191 allocate_deriv=.TRUE., error=error) 00192 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00193 END IF 00194 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 00195 deriv => xc_dset_get_derivative(deriv_set,"(rho)",& 00196 allocate_deriv=.TRUE.,error=error) 00197 CALL xc_derivative_get(deriv,deriv_data=e_rho,error=error) 00198 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",& 00199 allocate_deriv=.TRUE.,error=error) 00200 CALL xc_derivative_get(deriv,deriv_data=e_ndrho,error=error) 00201 END IF 00202 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 00203 deriv => xc_dset_get_derivative(deriv_set,"(rho)(rho)",& 00204 allocate_deriv=.TRUE.,error=error) 00205 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho,error=error) 00206 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(rho)",& 00207 allocate_deriv=.TRUE.,error=error) 00208 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_rho,error=error) 00209 deriv => xc_dset_get_derivative(deriv_set,& 00210 "(norm_drho)(norm_drho)", allocate_deriv=.TRUE.,error=error) 00211 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho,error=error) 00212 END IF 00213 IF (grad_deriv>=3.OR.grad_deriv==-3) THEN 00214 deriv => xc_dset_get_derivative(deriv_set,"(rho)(rho)(rho)",& 00215 allocate_deriv=.TRUE.,error=error) 00216 CALL xc_derivative_get(deriv,deriv_data=e_rho_rho_rho,error=error) 00217 deriv => xc_dset_get_derivative(deriv_set,& 00218 "(norm_drho)(rho)(rho)",allocate_deriv=.TRUE.,error=error) 00219 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_rho_rho,error=error) 00220 deriv => xc_dset_get_derivative(deriv_set,& 00221 "(norm_drho)(norm_drho)(rho)",allocate_deriv=.TRUE.,error=error) 00222 CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho_rho,error=error) 00223 !FM deriv => xc_dset_get_derivative(deriv_set,& 00224 !FM "(norm_drho)(norm_drho)(norm_drho)", allocate_deriv=.TRUE.,& 00225 !FM error=error) 00226 !FM call xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho_ndrho,error=error) 00227 END IF 00228 IF (grad_deriv>3.OR.grad_deriv<-3) THEN 00229 CALL cp_unimplemented_error(fromWhere=routineP, & 00230 message="derivatives bigger than 3 not implemented", & 00231 error=error, error_level=cp_failure_level) 00232 END IF 00233 00234 !$omp parallel default(none) & 00235 !$omp shared(rho, rho_1_3, norm_drho, e_0, e_rho, e_ndrho) & 00236 !$omp shared(e_rho_rho, e_ndrho_rho, e_ndrho_ndrho) & 00237 !$omp shared(e_rho_rho_rho, e_ndrho_rho_rho) & 00238 !$omp shared(e_ndrho_ndrho_rho, grad_deriv, npoints) & 00239 !$omp shared(epsilon_rho, epsilon_norm_drho, sc, error) 00240 00241 CALL lyp_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho,& 00242 e_0=e_0,e_rho=e_rho,e_ndrho=e_ndrho,e_rho_rho=e_rho_rho,& 00243 e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, & 00244 e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho,& 00245 e_ndrho_ndrho_rho=e_ndrho_ndrho_rho,& 00246 grad_deriv=grad_deriv,& 00247 npoints=npoints,epsilon_rho=epsilon_rho,epsilon_norm_drho=epsilon_norm_drho,sc=sc,& 00248 error=error) 00249 00250 !$omp end parallel 00251 00252 IF (cp_debug) THEN 00253 DEALLOCATE(dummy,stat=stat) 00254 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00255 ELSE 00256 NULLIFY(dummy) 00257 END IF 00258 END IF 00259 CALL timestop(handle) 00260 END SUBROUTINE lyp_lda_eval 00261 00262 ! ***************************************************************************** 00278 SUBROUTINE lyp_lda_calc(rho, rho_1_3, norm_drho,& 00279 e_0,e_rho,e_ndrho,e_rho_rho,e_ndrho_rho,& 00280 e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho,& 00281 grad_deriv,npoints,epsilon_rho, epsilon_norm_drho,sc,error) 00282 INTEGER, INTENT(in) :: npoints, grad_deriv 00283 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_rho, 00284 e_ndrho_rho_rho, e_rho_rho_rho, e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, 00285 e_ndrho, e_rho, e_0 00286 REAL(kind=dp), DIMENSION(1:npoints), 00287 INTENT(in) :: norm_drho, rho_1_3, rho 00288 REAL(kind=dp), INTENT(in) :: epsilon_rho, 00289 epsilon_norm_drho, sc 00290 TYPE(cp_error_type), INTENT(inout) :: error 00291 00292 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_lda_calc', 00293 routineP = moduleN//':'//routineN 00294 00295 INTEGER :: ii 00296 REAL(kind=dp) :: cf, epsilon_rho13, my_ndrho, my_rho, my_rho_1_3, t1, 00297 t102, t103, t105, t106, t11, t112, t123, t124, t127, t128, t13, t133, 00298 t134, t14, t161, t165, t166, t173, t176, t184, t185, t189, t192, t193, 00299 t196, t199, t2, t200, t201, t202, t203, t215, t22, t227, t228, t231, 00300 t245, t246, t248, t26, t265, t278, t279, t3, t332, t37, t38, t39, t4, 00301 t40, t41, t44, t45, t48, t5, t52, t6, t61, t62, t69, t7, t70, t77, t78, 00302 t80, t87, t88, t89, t93, t94, t95, t98, t99 00303 00304 epsilon_rho13=epsilon_rho**(1.0_dp/3.0_dp) 00305 cf=0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp) 00306 SELECT CASE(grad_deriv) 00307 CASE (1) 00308 !$omp do 00309 DO ii=1,npoints 00310 my_rho = rho(ii) 00311 IF (my_rho>epsilon_rho) THEN 00312 my_rho_1_3 = rho_1_3(ii) 00313 my_ndrho = norm_drho(ii) 00314 t1 = my_rho_1_3 ** 2 00315 t2 = t1 * my_rho 00316 t3 = 0.1e1_dp / t2 00317 t4 = a * t3 00318 t5 = my_rho ** 2 00319 t6 = t5 * my_rho 00320 t7 = my_rho_1_3 * t6 00321 t11 = 0.1e1_dp / my_rho_1_3 00322 t13 = EXP(-c * t11) 00323 t14 = b * t13 00324 t22 = my_ndrho ** 2 00325 t26 = my_rho_1_3 * t22 00326 t37 = -0.72e2_dp * t7 - 0.72e2_dp * t6 * d - 0.72e2_dp * t14 *& 00327 & t7 * cf - 0.72e2_dp * t14 * t6 * cf * d + 0.3e1_dp * t14& 00328 & * t1 * t22 + 0.10e2_dp * t14 * t26 * d + 0.7e1_dp * t14 & 00329 &* t26 * c + 0.7e1_dp * t14 * t22 * c * d 00330 t38 = my_rho_1_3 + d 00331 t39 = t38 ** 2 00332 t40 = 0.1e1_dp / t39 00333 t41 = t37 * t40 00334 00335 e_0(ii) = e_0(ii)& 00336 + ( t4 * t41 / 0.72e2_dp ) * sc 00337 t44 = 0.1e1_dp / t1 / t5 00338 t45 = a * t44 00339 t48 = my_rho_1_3 * t5 00340 t52 = b * c 00341 t61 = t13 * cf 00342 t62 = t61 * d 00343 t69 = 0.1e1_dp / t1 00344 t70 = t69 * t13 00345 t77 = 0.1e1_dp / my_rho 00346 t78 = t52 * t77 00347 t80 = t13 * t22 * d 00348 t87 = c ** 2 00349 t88 = b * t87 00350 t89 = t77 * t13 00351 t93 = my_rho_1_3 * my_rho 00352 t94 = 0.1e1_dp / t93 00353 t95 = t88 * t94 00354 t98 = -0.240e3_dp * t48 - 0.216e3_dp * t5 * d - 0.24e2_dp * t52& 00355 & * t5 * t13 * cf - 0.240e3_dp * t14 * t48 * cf -& 00356 & 0.24e2_dp * t52 * t2 * t62 - 0.216e3_dp * t14 * t5 * cf & 00357 &* d + 0.10e2_dp / 0.3e1_dp * t52 * t70 * t22 + 0.2e1_dp *& 00358 & t14 * t11 * t22 + 0.10e2_dp / 0.3e1_dp * t78 * t80 +& 00359 & 0.10e2_dp / 0.3e1_dp * t14 * t69 * t22 * d + 0.7e1_dp /& 00360 & 0.3e1_dp * t88 * t89 * t22 + 0.7e1_dp / 0.3e1_dp * t95 *& 00361 & t80 00362 t99 = t98 * t40 00363 t102 = 0.1e1_dp / t48 00364 t103 = a * t102 00365 t105 = 0.1e1_dp / t39 / t38 00366 t106 = t37 * t105 00367 00368 e_rho(ii) = e_rho(ii) & 00369 - ( 0.5e1_dp / 0.216e3_dp * t45 * t41 - t4 * t99 / 0.72e2_dp& 00370 & + t103 * t106 / 0.108e3_dp) * sc 00371 00372 t112 = my_rho_1_3 * my_ndrho 00373 t123 = 0.6e1_dp * t14 * t1 * my_ndrho + 0.20e2_dp * t14 * t112 & 00374 &* d + 0.14e2_dp * t14 * t112 * c + 0.14e2_dp * t14 *& 00375 & my_ndrho * c * d 00376 t124 = t123 * t40 00377 00378 e_ndrho(ii) = e_ndrho(ii)& 00379 + ( t4 * t124 / 0.72e2_dp ) * sc 00380 END IF 00381 END DO 00382 !$end do 00383 CASE default 00384 !$omp do 00385 DO ii=1,npoints 00386 my_rho = rho(ii) 00387 IF (my_rho>epsilon_rho) THEN 00388 my_rho_1_3 = rho_1_3(ii) 00389 my_ndrho = norm_drho(ii) 00390 t1 = my_rho_1_3 ** 2 00391 t2 = t1 * my_rho 00392 t3 = 0.1e1_dp / t2 00393 t4 = a * t3 00394 t5 = my_rho ** 2 00395 t6 = t5 * my_rho 00396 t7 = my_rho_1_3 * t6 00397 t11 = 0.1e1_dp / my_rho_1_3 00398 t13 = EXP(-c * t11) 00399 t14 = b * t13 00400 t22 = my_ndrho ** 2 00401 t26 = my_rho_1_3 * t22 00402 t37 = -0.72e2_dp * t7 - 0.72e2_dp * t6 * d - 0.72e2_dp * t14 *& 00403 & t7 * cf - 0.72e2_dp * t14 * t6 * cf * d + 0.3e1_dp * t14& 00404 & * t1 * t22 + 0.10e2_dp * t14 * t26 * d + 0.7e1_dp * t14 & 00405 &* t26 * c + 0.7e1_dp * t14 * t22 * c * d 00406 t38 = my_rho_1_3 + d 00407 t39 = t38 ** 2 00408 t40 = 0.1e1_dp / t39 00409 t41 = t37 * t40 00410 00411 IF (grad_deriv>=0) THEN 00412 e_0(ii) = e_0(ii)& 00413 + ( t4 * t41 / 0.72e2_dp ) * sc 00414 END IF 00415 00416 t44 = 0.1e1_dp / t1 / t5 00417 t45 = a * t44 00418 t48 = my_rho_1_3 * t5 00419 t52 = b * c 00420 t61 = t13 * cf 00421 t62 = t61 * d 00422 t69 = 0.1e1_dp / t1 00423 t70 = t69 * t13 00424 t77 = 0.1e1_dp / my_rho 00425 t78 = t52 * t77 00426 t80 = t13 * t22 * d 00427 t87 = c ** 2 00428 t88 = b * t87 00429 t89 = t77 * t13 00430 t93 = my_rho_1_3 * my_rho 00431 t94 = 0.1e1_dp / t93 00432 t95 = t88 * t94 00433 t98 = -0.240e3_dp * t48 - 0.216e3_dp * t5 * d - 0.24e2_dp * t52& 00434 & * t5 * t13 * cf - 0.240e3_dp * t14 * t48 * cf -& 00435 & 0.24e2_dp * t52 * t2 * t62 - 0.216e3_dp * t14 * t5 * cf & 00436 &* d + 0.10e2_dp / 0.3e1_dp * t52 * t70 * t22 + 0.2e1_dp *& 00437 & t14 * t11 * t22 + 0.10e2_dp / 0.3e1_dp * t78 * t80 +& 00438 & 0.10e2_dp / 0.3e1_dp * t14 * t69 * t22 * d + 0.7e1_dp /& 00439 & 0.3e1_dp * t88 * t89 * t22 + 0.7e1_dp / 0.3e1_dp * t95 *& 00440 & t80 00441 t99 = t98 * t40 00442 t102 = 0.1e1_dp / t48 00443 t103 = a * t102 00444 t105 = 0.1e1_dp / t39 / t38 00445 t106 = t37 * t105 00446 t112 = my_rho_1_3 * my_ndrho 00447 t123 = 0.6e1_dp * t14 * t1 * my_ndrho + 0.20e2_dp * t14 * t112 & 00448 &* d + 0.14e2_dp * t14 * t112 * c + 0.14e2_dp * t14 *& 00449 & my_ndrho * c * d 00450 t124 = t123 * t40 00451 00452 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 00453 e_rho(ii) = e_rho(ii) & 00454 - ( 0.5e1_dp / 0.216e3_dp * t45 * t41 - t4 * t99 / 0.72e2_dp& 00455 & + t103 * t106 / 0.108e3_dp ) * sc 00456 e_ndrho(ii) = e_ndrho(ii)& 00457 + ( t4 * t124 / 0.72e2_dp ) * sc 00458 END IF 00459 00460 t127 = 0.1e1_dp / t1 / t6 00461 t128 = a * t127 00462 t133 = 0.1e1_dp / t7 00463 t134 = a * t133 00464 t161 = t3 * t13 00465 t165 = 0.1e1_dp / t5 00466 t166 = t165 * t13 00467 t173 = t52 * t165 00468 t176 = t88 * t102 00469 t184 = b * t87 * c 00470 t185 = t102 * t13 00471 t189 = t184 * t44 00472 t192 = -0.560e3_dp * t93 - 0.432e3_dp * my_rho * d - 0.128e3_dp& 00473 & * t52 * my_rho * t13 * cf - 0.8e1_dp * t88 * t1 * t13 *& 00474 & cf - 0.560e3_dp * t14 * t93 * cf - 0.112e3_dp * t52 * t1& 00475 & * t62 - 0.8e1_dp * t88 * my_rho_1_3 * t62 - 0.432e3_dp *& 00476 & t14 * my_rho * cf * d - 0.14e2_dp / 0.9e1_dp * t52 *& 00477 & t161 * t22 - 0.11e2_dp / 0.9e1_dp * t88 * t166 * t22 -& 00478 & 0.2e1_dp / 0.3e1_dp * t14 * t94 * t22 - 0.20e2_dp /& 00479 & 0.9e1_dp * t173 * t80 - 0.2e1_dp * t176 * t80 -& 00480 & 0.20e2_dp / 0.9e1_dp * t14 * t3 * t22 * d + 0.7e1_dp /& 00481 & 0.9e1_dp * t184 * t185 * t22 + 0.7e1_dp / 0.9e1_dp *& 00482 & t189 * t80 00483 t193 = t192 * t40 00484 t196 = t98 * t105 00485 t199 = 0.1e1_dp / t6 00486 t200 = a * t199 00487 t201 = t39 ** 2 00488 t202 = 0.1e1_dp / t201 00489 t203 = t37 * t202 00490 t215 = t13 * my_ndrho * d 00491 t227 = 0.20e2_dp / 0.3e1_dp * t52 * t70 * my_ndrho + 0.4e1_dp *& 00492 & t14 * t11 * my_ndrho + 0.20e2_dp / 0.3e1_dp * t78 * t215& 00493 & + 0.20e2_dp / 0.3e1_dp * t14 * t69 * my_ndrho * d +& 00494 & 0.14e2_dp / 0.3e1_dp * t88 * t89 * my_ndrho + 0.14e2_dp & 00495 &/ 0.3e1_dp * t95 * t215 00496 t228 = t227 * t40 00497 t231 = t123 * t105 00498 t245 = 0.6e1_dp * t14 * t1 + 0.20e2_dp * t14 * my_rho_1_3 * d +& 00499 & 0.14e2_dp * t14 * my_rho_1_3 * c + 0.14e2_dp * t14 * c *& 00500 & d 00501 t246 = t245 * t40 00502 00503 IF (grad_deriv>=2 .OR.grad_deriv==-2) THEN 00504 e_rho_rho(ii) = e_rho_rho(ii)& 00505 + ( 0.5e1_dp / 0.81e2_dp * t128 * t41 - 0.5e1_dp / 0.108e3_dp& 00506 & * t45 * t99 + t134 * t106 / 0.27e2_dp + t4 * t193 /& 00507 & 0.72e2_dp - t103 * t196 / 0.54e2_dp + t200 * t203 /& 00508 & 0.108e3_dp ) * sc 00509 e_ndrho_rho(ii) = e_ndrho_rho(ii)& 00510 - ( 0.5e1_dp / 0.216e3_dp * t45 * t124 - t4 * t228 /& 00511 & 0.72e2_dp + t103 * t231 / 0.108e3_dp ) * sc 00512 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii)& 00513 + ( t4 * t246 / 0.72e2_dp ) * sc 00514 END IF 00515 00516 t248 = t5 ** 2 00517 t265 = 0.1e1_dp / t248 00518 t278 = t87 ** 2 00519 t279 = b * t278 00520 t332 = -0.432e3_dp * d - 0.2240e4_dp / 0.3e1_dp * my_rho_1_3 -& 00521 & 0.74e2_dp / 0.27e2_dp * t184 * t127 * t80 + 0.100e3_dp /& 00522 & 0.27e2_dp * t14 * t44 * t22 * d + 0.7e1_dp / 0.27e2_dp *& 00523 & t279 * t127 * t13 * t22 - 0.8e1_dp / 0.3e1_dp * t184 *& 00524 & t70 * cf - 0.2240e4_dp / 0.3e1_dp * t14 * my_rho_1_3 *& 00525 & cf - 0.48e2_dp * t88 * t11 * t13 * cf - 0.944e3_dp /& 00526 & 0.3e1_dp * t52 * t61 - 0.40e2_dp * t88 * t69 * t62 -& 00527 & 0.656e3_dp / 0.3e1_dp * t52 * t11 * t62 + 0.7e1_dp /& 00528 & 0.27e2_dp * t279 * t265 * t80 + 0.64e2_dp / 0.27e2_dp *& 00529 & t52 * t44 * t13 * t22 - 0.8e1_dp / 0.3e1_dp * t184 * t77& 00530 & * t62 - 0.432e3_dp * t14 * cf * d + 0.52e2_dp /& 00531 & 0.27e2_dp * t88 * t199 * t13 * t22 - 0.20e2_dp /& 00532 & 0.9e1_dp * t184 * t133 * t13 * t22 + 0.8e1_dp / 0.9e1_dp& 00533 & * t14 * t102 * t22 + 0.100e3_dp / 0.27e2_dp * t52 * t199& 00534 & * t80 + 0.106e3_dp / 0.27e2_dp * t88 * t133 * t80 00535 00536 IF (grad_deriv>=3 .OR. grad_deriv==-3) THEN 00537 e_rho_rho_rho(ii) = e_rho_rho_rho(ii)& 00538 - ( 0.55e2_dp / 0.243e3_dp * a / t1 / t248 * t41 - 0.5e1_dp & 00539 &/ 0.27e2_dp * t128 * t99 + 0.40e2_dp / 0.243e3_dp * a /& 00540 & my_rho_1_3 / t248 * t106 + 0.5e1_dp / 0.72e2_dp * t45 *& 00541 & t193 - t134 * t196 / 0.9e1_dp + 0.7e1_dp / 0.108e3_dp *& 00542 & a * t265 * t203 - t4 * t332 * t40 / 0.72e2_dp + t103 *& 00543 & t192 * t105 / 0.36e2_dp - t200 * t98 * t202 / 0.36e2_dp & 00544 &+ t128 * t37 / t201 / t38 / 0.81e2_dp ) * sc 00545 e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii)& 00546 + ( 0.5e1_dp / 0.81e2_dp * t128 * t124 - 0.5e1_dp /& 00547 & 0.108e3_dp * t45 * t228 + t134 * t231 / 0.27e2_dp + t4 *& 00548 & (-0.28e2_dp / 0.9e1_dp * t52 * t161 * my_ndrho -& 00549 & 0.22e2_dp / 0.9e1_dp * t88 * t166 * my_ndrho - 0.4e1_dp & 00550 &/ 0.3e1_dp * t14 * t94 * my_ndrho - 0.40e2_dp / 0.9e1_dp & 00551 &* t173 * t215 - 0.4e1_dp * t176 * t215 - 0.40e2_dp /& 00552 & 0.9e1_dp * t14 * t3 * my_ndrho * d + 0.14e2_dp /& 00553 & 0.9e1_dp * t184 * t185 * my_ndrho + 0.14e2_dp / 0.9e1_dp& 00554 & * t189 * t215) * t40 / 0.72e2_dp - t103 * t227 * t105 /& 00555 & 0.54e2_dp + t200 * t123 * t202 / 0.108e3_dp ) * sc 00556 e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii)& 00557 - ( 0.5e1_dp / 0.216e3_dp * t45 * t246 - t4 * (0.20e2_dp /& 00558 & 0.3e1_dp * t52 * t70 + 0.4e1_dp * t14 * t11 + 0.20e2_dp & 00559 &/ 0.3e1_dp * t52 * t89 * d + 0.20e2_dp / 0.3e1_dp * t14 *& 00560 & t69 * d + 0.14e2_dp / 0.3e1_dp * t88 * t89 + 0.14e2_dp /& 00561 & 0.3e1_dp * t88 * t94 * t13 * d) * t40 / 0.72e2_dp + t103& 00562 & * t245 * t105 / 0.108e3_dp ) * sc 00563 END IF 00564 END IF 00565 00566 !FM t1 = my_rho_1_3 ** 2 00567 !FM t2 = t1 * my_rho 00568 !FM t3 = 0.1e1_dp / t2 00569 !FM t4 = a * t3 00570 !FM t5 = my_rho ** 2 00571 !FM t6 = t5 * my_rho 00572 !FM t7 = my_rho_1_3 * t6 00573 !FM t11 = 0.1e1_dp / my_rho_1_3 00574 !FM t13 = exp(-c * t11) 00575 !FM t14 = b * t13 00576 !FM t22 = my_ndrho ** 2 00577 !FM t26 = my_rho_1_3 * t22 00578 !FM t37 = -0.72e2_dp *( t7 - t14 *& 00579 !FM & t7 * cf - t6 * d * (1.0_dp+ t14 * cf)) + t14 *(t22*(0.3e1_dp & 00580 !FM & * t1 + 0.7e1_dp * c * d)& 00581 !FM + 0.10e2_dp * t26 * d + 0.7e1_dp & 00582 !FM &* t26 * c ) 00583 !FM t38 = my_rho_1_3 + d 00584 !FM t39 = t38 ** 2 00585 !FM t40 = 0.1e1_dp / t39 00586 !FM t41 = t37 * t40 00587 !FM 00588 !FM IF (grad_deriv>=0.AND.my_rho>epsilon_rho) THEN 00589 !FM e_0(ii) = e_0(ii)& 00590 !FM + t4 * t41 / 0.72e2_dp 00591 !FM END IF 00592 !FM 00593 !FM t44 = 0.1e1_dp / (t1 * t5) 00594 !FM t45 = a * t44 00595 !FM t48 = my_rho_1_3 * t5 00596 !FM t52 = b * c 00597 !FM t61 = t13 * cf 00598 !FM t62 = t61 * d 00599 !FM t69 = 0.1e1_dp / t1 00600 !FM t70 = t69 * t13 00601 !FM t77 = 0.1e1_dp / my_rho 00602 !FM t78 = t52 * t77 00603 !FM t80 = t13 * t22 * d 00604 !FM t87 = c ** 2 00605 !FM t88 = b * t87 00606 !FM t89 = t77 * t13 00607 !FM t93 = my_rho_1_3 * my_rho 00608 !FM t94 = 0.1e1_dp / t93 00609 !FM t95 = t88 * t94 00610 !FM t98 = -0.216e3_dp * t5 *d -0.240e3_dp * t48(1.0_dp+ t14 * cf) -& 00611 !FM & 0.24e2_dp * t52 * (t2 * t62 +t13*t5*cf) & 00612 !FM - 0.216e3_dp * t14 * t5 * cf & 00613 !FM &* d + t22 *(0.10e2_dp / 0.3e1_dp * t52 * t70 + 0.2e1_dp *& 00614 !FM & t14 * t11 + 0.10e2_dp / 0.3e1_dp * t14 * t69 * d + 0.7e1_dp /& 00615 !FM & 0.3e1_dp * t88 * t89 ) +(0.10e2_dp / 0.3e1_dp * t78 +& 00616 !FM & 0.7e1_dp / 0.3e1_dp * t95 ) *& 00617 !FM & t80 00618 !FM t99 = t98 * t40 00619 !FM t102 = 0.1e1_dp / t48 00620 !FM t103 = a * t102 00621 !FM t105 = 0.1e1_dp / (t39 * t38) 00622 !FM t106 = t37 * t105 00623 !FM t112 = my_rho_1_3 * my_ndrho 00624 !FM t123 = t14 *(0.6e1_dp t1 * my_ndrho + t112 * 0.20e2_dp & 00625 !FM &* d + 0.14e2_dp * c *(t112 + t14 *& 00626 !FM & my_ndrho * d)) 00627 !FM t124 = t123 * t40 00628 !FM 00629 !FM IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 00630 !FM e_rho(ii) = e_rho(ii) & 00631 !FM -0.5e1_dp / 0.216e3_dp * t45 * t41 + t4 * t99 / 0.72e2_dp& 00632 !FM & - t103 * t106 / 0.108e3_dp 00633 !FM e_ndrho(ii) = e_ndrho(ii)& 00634 !FM +t4 * t124 / 0.72e2_dp 00635 !FM END IF 00636 !FM 00637 !FM t127 = 0.1e1_dp / (t1 * t6) 00638 !FM t128 = a * t127 00639 !FM t133 = 0.1e1_dp / t7 00640 !FM t134 = a * t133 00641 !FM t161 = t3 * t13 00642 !FM t165 = 0.1e1_dp / t5 00643 !FM t166 = t165 * t13 00644 !FM t173 = t52 * t165 00645 !FM t176 = t88 * t102 00646 !FM t184 = b * t87 * c 00647 !FM t185 = t102 * t13 00648 !FM t189 = t184 * t44 00649 !FM t192 = -0.560e3_dp * t93 - 0.432e3_dp * my_rho * d +cf*(& 00650 !FM -t13*(0.128e3_dp& 00651 !FM & * t52 * my_rho + 0.8e1_dp * t88 * t1 )& 00652 !FM & - 0.560e3_dp * t14 * t93 ) - 0.112e3_dp * t52 * t1& 00653 !FM & * t62 - 0.8e1_dp * t88 * my_rho_1_3 * t62 - 0.432e3_dp *& 00654 !FM & t14 * my_rho * cf * d - 0.14e2_dp / 0.9e1_dp * t52 *& 00655 !FM & t161 * t22 - 0.11e2_dp / 0.9e1_dp * t88 * t166 * t22 -& 00656 !FM & 0.2e1_dp / 0.3e1_dp * t14 * t94 * t22 - 0.20e2_dp /& 00657 !FM & 0.9e1_dp * t173 * t80 - 0.2e1_dp * t176 * t80 -& 00658 !FM & 0.20e2_dp / 0.9e1_dp * t14 * t3 * t22 * d + 0.7e1_dp /& 00659 !FM & 0.9e1_dp * t184 * t185 * t22 + 0.7e1_dp / 0.9e1_dp *& 00660 !FM & t189 * t80 00661 !FM t193 = t192 * t40 00662 !FM t196 = t98 * t105 00663 !FM t199 = 0.1e1_dp / t6 00664 !FM t200 = a * t199 00665 !FM t201 = t39 ** 2 00666 !FM t202 = 0.1e1_dp / t201 00667 !FM t203 = t37 * t202 00668 !FM t215 = t13 * my_ndrho * d 00669 !FM t227 = 0.20e2_dp / 0.3e1_dp * t52 * t70 * my_ndrho + 0.4e1_dp *& 00670 !FM & t14 * t11 * my_ndrho + 0.20e2_dp / 0.3e1_dp * t78 * t215& 00671 !FM & + 0.20e2_dp / 0.3e1_dp * t14 * t69 * my_ndrho * d +& 00672 !FM & 0.14e2_dp / 0.3e1_dp * t88 * t89 * my_ndrho + 0.14e2_dp & 00673 !FM &/ 0.3e1_dp * t95 * t215 00674 !FM t228 = t227 * t40 00675 !FM t231 = t123 * t105 00676 !FM t245 = 0.6e1_dp * t14 * t1 + 0.20e2_dp * t14 * my_rho_1_3 * d +& 00677 !FM & 0.14e2_dp * t14 * my_rho_1_3 * c + 0.14e2_dp * t14 * c *& 00678 !FM & d 00679 !FM t246 = t245 * t40 00680 !FM 00681 !FM IF (grad_deriv>=2 .OR.grad_deriv==-2) THEN 00682 !FM e_rho_rho(ii) = e_rho_rho(ii)& 00683 !FM +0.5e1_dp / 0.81e2_dp * t128 * t41 - 0.5e1_dp / 0.108e3_dp& 00684 !FM & * t45 * t99 + t134 * t106 / 0.27e2_dp + t4 * t193 /& 00685 !FM & 0.72e2_dp - t103 * t196 / 0.54e2_dp + t200 * t203 /& 00686 !FM & 0.108e3_dp 00687 !FM e_ndrho_rho(ii) = e_ndrho_rho(ii)& 00688 !FM -0.5e1_dp / 0.216e3_dp * t45 * t124 + t4 * t228 /& 00689 !FM & 0.72e2_dp - t103 * t231 / 0.108e3_dp 00690 !FM e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii)& 00691 !FM +t4 * t246 / 0.72e2_dp 00692 !FM END IF 00693 !FM 00694 !FM t248 = t5 ** 2 00695 !FM t265 = 0.1e1_dp / t248 00696 !FM t278 = t87 ** 2 00697 !FM t279 = b * t278 00698 !FM t332 = -0.432e3_dp * d - 0.2240e4_dp / 0.3e1_dp * my_rho_1_3 -& 00699 !FM & 0.74e2_dp / 0.27e2_dp * t184 * t127 * t80 + 0.100e3_dp /& 00700 !FM & 0.27e2_dp * t14 * t44 * t22 * d + 0.7e1_dp / 0.27e2_dp *& 00701 !FM & t279 * t127 * t13 * t22 - 0.8e1_dp / 0.3e1_dp * t184 *& 00702 !FM & t70 * cf - 0.2240e4_dp / 0.3e1_dp * t14 * my_rho_1_3 *& 00703 !FM & cf - 0.48e2_dp * t88 * t11 * t13 * cf - 0.944e3_dp /& 00704 !FM & 0.3e1_dp * t52 * t61 - 0.40e2_dp * t88 * t69 * t62 -& 00705 !FM & 0.656e3_dp / 0.3e1_dp * t52 * t11 * t62 + 0.7e1_dp /& 00706 !FM & 0.27e2_dp * t279 * t265 * t80 + 0.64e2_dp / 0.27e2_dp *& 00707 !FM & t52 * t44 * t13 * t22 - 0.8e1_dp / 0.3e1_dp * t184 * t77& 00708 !FM & * t62 - 0.432e3_dp * t14 * cf * d + 0.52e2_dp /& 00709 !FM & 0.27e2_dp * t88 * t199 * t13 * t22 - 0.20e2_dp /& 00710 !FM & 0.9e1_dp * t184 * t133 * t13 * t22 + 0.8e1_dp / 0.9e1_dp& 00711 !FM & * t14 * t102 * t22 + 0.100e3_dp / 0.27e2_dp * t52 * t199& 00712 !FM & * t80 + 0.106e3_dp / 0.27e2_dp * t88 * t133 * t80 00713 !FM 00714 !FM IF (grad_deriv>=3 .OR. grad_deriv==-3) THEN 00715 !FM e_rho_rho_rho(ii) = e_rho_rho_rho(ii)& 00716 !FM -0.55e2_dp / 0.243e3_dp * a / t1 / t248 * t41 + 0.5e1_dp & 00717 !FM &/ 0.27e2_dp * t128 * t99 - 0.40e2_dp / 0.243e3_dp * a /& 00718 !FM & my_rho_1_3 / t248 * t106 - 0.5e1_dp / 0.72e2_dp * t45 *& 00719 !FM & t193 + t134 * t196 / 0.9e1_dp - 0.7e1_dp / 0.108e3_dp *& 00720 !FM & a * t265 * t203 + t4 * t332 * t40 / 0.72e2_dp - t103 *& 00721 !FM & t192 * t105 / 0.36e2_dp + t200 * t98 * t202 / 0.36e2_dp & 00722 !FM &- t128 * t37 / t201 / t38 / 0.81e2_dp 00723 !FM e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii)& 00724 !FM +0.5e1_dp / 0.81e2_dp * t128 * t124 - 0.5e1_dp /& 00725 !FM & 0.108e3_dp * t45 * t228 + t134 * t231 / 0.27e2_dp + t4 *& 00726 !FM & (-0.28e2_dp / 0.9e1_dp * t52 * t161 * my_ndrho -& 00727 !FM & 0.22e2_dp / 0.9e1_dp * t88 * t166 * my_ndrho - 0.4e1_dp & 00728 !FM &/ 0.3e1_dp * t14 * t94 * my_ndrho - 0.40e2_dp / 0.9e1_dp & 00729 !FM &* t173 * t215 - 0.4e1_dp * t176 * t215 - 0.40e2_dp /& 00730 !FM & 0.9e1_dp * t14 * t3 * my_ndrho * d + 0.14e2_dp /& 00731 !FM & 0.9e1_dp * t184 * t185 * my_ndrho + 0.14e2_dp / 0.9e1_dp& 00732 !FM & * t189 * t215) * t40 / 0.72e2_dp - t103 * t227 * t105 /& 00733 !FM & 0.54e2_dp + t200 * t123 * t202 / 0.108e3_dp 00734 !FM e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii)& 00735 !FM -0.5e1_dp / 0.216e3_dp * t45 * t246 + t4 * (0.20e2_dp /& 00736 !FM & 0.3e1_dp * t52 * t70 + 0.4e1_dp * t14 * t11 + 0.20e2_dp & 00737 !FM &/ 0.3e1_dp * t52 * t89 * d + 0.20e2_dp / 0.3e1_dp * t14 *& 00738 !FM & t69 * d + 0.14e2_dp / 0.3e1_dp * t88 * t89 + 0.14e2_dp /& 00739 !FM & 0.3e1_dp * t88 * t94 * t13 * d) * t40 / 0.72e2_dp - t103& 00740 !FM & * t245 * t105 / 0.108e3_dp 00741 !FM END IF 00742 00743 END DO 00744 00745 !$omp end do 00746 00747 END SELECT 00748 END SUBROUTINE lyp_lda_calc 00749 00750 ! ***************************************************************************** 00760 SUBROUTINE lyp_lsd_eval(rho_set,deriv_set,grad_deriv,lyp_params,error) 00761 TYPE(xc_rho_set_type), POINTER :: rho_set 00762 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00763 INTEGER, INTENT(in) :: grad_deriv 00764 TYPE(section_vals_type), POINTER :: lyp_params 00765 TYPE(cp_error_type), INTENT(inout) :: error 00766 00767 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_lsd_eval', 00768 routineP = moduleN//':'//routineN 00769 00770 INTEGER :: handle, npoints, stat 00771 INTEGER, DIMENSION(:, :), POINTER :: bo 00772 LOGICAL :: failure 00773 REAL(kind=dp) :: epsilon_drho, epsilon_rho, sc 00774 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, 00775 e_ndr_ndr, e_ndr_ra, e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, 00776 e_ndra_rb, e_ndrb, e_ndrb_ndrb, e_ndrb_ra, e_ndrb_rb, e_ra, e_ra_ra, 00777 e_ra_rb, e_rb, e_rb_rb, norm_drho, norm_drhoa, norm_drhob, rhoa, rhob 00778 TYPE(xc_derivative_type), POINTER :: deriv 00779 00780 CALL timeset(routineN,handle) 00781 failure=.FALSE. 00782 NULLIFY(deriv, bo) 00783 00784 CALL section_vals_val_get(lyp_params,"scale_c",r_val=sc,error=error) 00785 CALL cite_reference(Lee1988) 00786 00787 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00788 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00789 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00790 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00791 IF (.NOT. failure) THEN 00792 CALL xc_rho_set_get(rho_set,& 00793 rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, & 00794 norm_drhob=norm_drhob, norm_drho=norm_drho, & 00795 rho_cutoff=epsilon_rho,& 00796 drho_cutoff=epsilon_drho, local_bounds=bo, error=error) 00797 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00798 00799 ! meaningful default for the arrays we don't need: let us make compiler 00800 ! and debugger happy... 00801 IF (cp_debug) THEN 00802 ALLOCATE(dummy(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 00803 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00804 ELSE 00805 dummy=> rhoa 00806 END IF 00807 e_0 => dummy 00808 e_ra => dummy 00809 e_rb => dummy 00810 e_ndra_ra => dummy 00811 e_ndra_rb => dummy 00812 e_ndrb_ra => dummy 00813 e_ndrb_rb => dummy 00814 e_ndr_ndr => dummy 00815 e_ndra_ndra => dummy 00816 e_ndrb_ndrb => dummy 00817 e_ndr => dummy 00818 e_ndra => dummy 00819 e_ndrb => dummy 00820 e_ra_ra => dummy 00821 e_ra_rb => dummy 00822 e_rb_rb => dummy 00823 e_ndr_ra => dummy 00824 e_ndr_rb => dummy 00825 00826 IF (grad_deriv>=0) THEN 00827 deriv => xc_dset_get_derivative(deriv_set,"",& 00828 allocate_deriv=.TRUE., error=error) 00829 CALL xc_derivative_get(deriv, deriv_data=e_0,error=error) 00830 END IF 00831 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 00832 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)",& 00833 allocate_deriv=.TRUE.,error=error) 00834 CALL xc_derivative_get(deriv,deriv_data=e_ra,error=error) 00835 deriv => xc_dset_get_derivative(deriv_set,"(rhob)",& 00836 allocate_deriv=.TRUE.,error=error) 00837 CALL xc_derivative_get(deriv,deriv_data=e_rb,error=error) 00838 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",& 00839 allocate_deriv=.TRUE.,error=error) 00840 CALL xc_derivative_get(deriv,deriv_data=e_ndr,error=error) 00841 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)",& 00842 allocate_deriv=.TRUE.,error=error) 00843 CALL xc_derivative_get(deriv,deriv_data=e_ndra,error=error) 00844 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)",& 00845 allocate_deriv=.TRUE.,error=error) 00846 CALL xc_derivative_get(deriv,deriv_data=e_ndrb,error=error) 00847 END IF 00848 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 00849 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)(rhoa)",& 00850 allocate_deriv=.TRUE.,error=error) 00851 CALL xc_derivative_get(deriv,deriv_data=e_ra_ra,error=error) 00852 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)(rhob)",& 00853 allocate_deriv=.TRUE.,error=error) 00854 CALL xc_derivative_get(deriv,deriv_data=e_ra_rb,error=error) 00855 deriv => xc_dset_get_derivative(deriv_set,"(rhob)(rhob)",& 00856 allocate_deriv=.TRUE.,error=error) 00857 CALL xc_derivative_get(deriv,deriv_data=e_rb_rb,error=error) 00858 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(rhoa)",& 00859 allocate_deriv=.TRUE.,error=error) 00860 CALL xc_derivative_get(deriv,deriv_data=e_ndr_ra,error=error) 00861 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(rhob)",& 00862 allocate_deriv=.TRUE.,error=error) 00863 CALL xc_derivative_get(deriv,deriv_data=e_ndr_rb,error=error) 00864 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(rhoa)",& 00865 allocate_deriv=.TRUE.,error=error) 00866 CALL xc_derivative_get(deriv,deriv_data=e_ndra_ra,error=error) 00867 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(rhob)",& 00868 allocate_deriv=.TRUE.,error=error) 00869 CALL xc_derivative_get(deriv,deriv_data=e_ndra_rb,error=error) 00870 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)(rhob)",& 00871 allocate_deriv=.TRUE.,error=error) 00872 CALL xc_derivative_get(deriv,deriv_data=e_ndrb_rb,error=error) 00873 deriv => xc_dset_get_derivative(deriv_set,& 00874 "(norm_drho)(norm_drho)", allocate_deriv=.TRUE.,error=error) 00875 CALL xc_derivative_get(deriv,deriv_data=e_ndr_ndr,error=error) 00876 deriv => xc_dset_get_derivative(deriv_set,& 00877 "(norm_drhoa)(norm_drhoa)", allocate_deriv=.TRUE.,error=error) 00878 CALL xc_derivative_get(deriv,deriv_data=e_ndra_ndra,error=error) 00879 deriv => xc_dset_get_derivative(deriv_set,& 00880 "(norm_drhob)(norm_drhob)", allocate_deriv=.TRUE.,error=error) 00881 CALL xc_derivative_get(deriv,deriv_data=e_ndrb_ndrb,error=error) 00882 END IF 00883 00884 !$omp parallel default(none) & 00885 !$omp shared(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob) & 00886 !$omp shared(e_0, e_ra, e_rb, e_ndra_ra, e_ndra_rb, e_ndrb_ra) & 00887 !$omp shared(e_ndrb_rb, e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb) & 00888 !$omp shared(e_ndr, e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb) & 00889 !$omp shared(e_ndr_ra, e_ndr_rb, grad_deriv) & 00890 !$omp shared(npoints, epsilon_rho, epsilon_drho, sc, error) 00891 00892 CALL lyp_lsd_calc(& 00893 rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa,& 00894 norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb,& 00895 & e_ndra_ra=e_ndra_ra, e_ndra_rb=e_ndra_rb, e_ndrb_ra& 00896 &=e_ndrb_ra, e_ndrb_rb=e_ndrb_rb, e_ndr_ndr=e_ndr_ndr,& 00897 e_ndra_ndra=e_ndra_ndra, e_ndrb_ndrb=e_ndrb_ndrb, e_ndr=e_ndr,& 00898 e_ndra=e_ndra, e_ndrb=e_ndrb, e_ra_ra=e_ra_ra, & 00899 e_ra_rb=e_ra_rb, e_rb_rb=e_rb_rb, e_ndr_ra=e_ndr_ra,& 00900 e_ndr_rb=e_ndr_rb,& 00901 grad_deriv=grad_deriv, npoints=npoints, & 00902 epsilon_rho=epsilon_rho,epsilon_drho=epsilon_drho,sc=sc,error=error) 00903 00904 !$omp end parallel 00905 00906 IF (cp_debug) THEN 00907 DEALLOCATE(dummy,stat=stat) 00908 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00909 ELSE 00910 NULLIFY(dummy) 00911 END IF 00912 END IF 00913 CALL timestop(handle) 00914 END SUBROUTINE lyp_lsd_eval 00915 00916 ! ***************************************************************************** 00931 SUBROUTINE lyp_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob,& 00932 e_0, e_ra, e_rb, e_ndra_ra, e_ndra_rb, e_ndrb_ra, e_ndrb_rb, e_ndr_ndr,& 00933 e_ndra_ndra, e_ndrb_ndrb, e_ndr,& 00934 e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, e_ndr_ra, e_ndr_rb,& 00935 grad_deriv,npoints,epsilon_rho,epsilon_drho,sc,error) 00936 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob, norm_drho, 00937 norm_drhoa, norm_drhob 00938 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_ra, e_rb, e_ndra_ra, 00939 e_ndra_rb, e_ndrb_ra, e_ndrb_rb, e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb, 00940 e_ndr, e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, e_ndr_ra, e_ndr_rb 00941 INTEGER, INTENT(in) :: grad_deriv, npoints 00942 REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_drho, sc 00943 TYPE(cp_error_type), INTENT(inout) :: error 00944 00945 CHARACTER(len=*), PARAMETER :: routineN = 'lyp_lsd_calc', 00946 routineP = moduleN//':'//routineN 00947 REAL(kind=dp), PARAMETER :: small = 1.0e-20, t_1_9 = 1.0_dp/9.0_dp, 00948 t_2_3 = 2.0_dp/3.0_dp, t_4_27 = 4.0_dp/27.0_dp, t_4_3 = 4.0_dp/3.0_dp, 00949 t_4_9 = 4.0_dp/9.0_dp, t_8_27 = 8.0_dp/27.0_dp 00950 00951 INTEGER :: ii 00952 LOGICAL :: failure 00953 REAL(kind=dp) :: cf, my_ndrho, my_ndrhoa, my_ndrhob, my_rho, my_rhoa, 00954 my_rhob, t1, t100, t101, t102, t103, t104, t107, t108, t109, t112, 00955 t115, t118, t12, t120, t124, t129, t13, t130, t132, t135, t138, t14, 00956 t140, t142, t145, t15, t153, t155, t159, t16, t164, t168, t17, t171, 00957 t176, t18, t181, t182, t185, t189, t194, t195, t198, t2, t20, t202, 00958 t205, t206, t21, t210, t214, t215, t218, t22, t222, t228, t23, t232, 00959 t236, t238, t24, t243, t249, t25, t252, t253, t255, t257, t26, t260, 00960 t265, t268, t27, t270, t29, t3, t30, t304, t31, t310, t313, t316, t319, 00961 t322, t326, t329, t332, t334, t341, t35, t354, t356, t360, t363 00962 REAL(kind=dp) :: t37, t373, t376, t381, t39, t391, t4, t40, t408, t41, 00963 t415, t419, t422, t430, t44, t445, t449, t45, t452, t46, t467, t47, 00964 t48, t483, t487, t49, t490, t5, t50, t505, t515, t519, t52, t520, t54, 00965 t56, t57, t6, t61, t62, t64, t66, t7, t72, t75, t78, t8, t82, t85, t86, 00966 t87, t88, t90, t92, t94, t95, t98, t99 00967 00968 failure=.FALSE. 00969 cf=0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp) 00970 00971 !$omp do 00972 00973 DO ii=1,npoints 00974 my_rhoa=MAX(rhoa(ii),0.0_dp) 00975 my_rhob=MAX(rhob(ii),0.0_dp) 00976 my_rho=my_rhoa+my_rhob 00977 IF (my_rho>epsilon_rho) THEN 00978 my_ndrho=norm_drho(ii) 00979 my_ndrhoa=norm_drhoa(ii) 00980 my_ndrhob=norm_drhob(ii) 00981 00982 t1 = my_rho ** (0.1e1_dp / 0.3e1_dp) 00983 t2 = 0.1e1_dp / t1 00984 t3 = d * t2 00985 t4 = 0.1e1_dp + t3 00986 t5 = 0.1e1_dp / t4 00987 t6 = a * t5 00988 t7 = my_rhoa * my_rhob 00989 t8 = 0.1e1_dp / my_rho 00990 t12 = a * b 00991 t13 = c * t2 00992 t14 = EXP(-t13) 00993 t15 = t12 * t14 00994 t16 = my_rho ** 2 00995 t17 = t16 * my_rho 00996 t18 = t1 ** 2 00997 t20 = 0.1e1_dp / t18 / t17 00998 t21 = t5 * t20 00999 t22 = 2 ** (0.1e1_dp / 0.3e1_dp) 01000 t23 = t22 ** 2 01001 t24 = t23 * cf 01002 t25 = my_rhoa ** 2 01003 t26 = my_rhoa ** (0.1e1_dp / 0.3e1_dp) 01004 t27 = t26 ** 2 01005 t29 = my_rhob ** 2 01006 t30 = my_rhob ** (0.1e1_dp / 0.3e1_dp) 01007 t31 = t30 ** 2 01008 t35 = 0.8e1_dp * t24 * (t27 * t25 + t31 * t29) 01009 t37 = t3 * t5 01010 t39 = 0.47e2_dp / 0.18e2_dp - 0.7e1_dp / 0.18e2_dp * t13 -& 01011 & 0.7e1_dp / 0.18e2_dp * t37 01012 t40 = my_ndrho ** 2 01013 t41 = t39 * t40 01014 t44 = 0.5e1_dp / 0.2e1_dp - t13 / 0.18e2_dp - t37 / 0.18e2_dp 01015 t45 = my_ndrhoa ** 2 01016 t46 = my_ndrhob ** 2 01017 t47 = t45 + t46 01018 t48 = t44 * t47 01019 t49 = t13 + t37 - 0.11e2_dp 01020 t50 = my_rhoa * t8 01021 t52 = my_rhob * t8 01022 t54 = t50 * t45 + t52 * t46 01023 t56 = t49 * t54 / 0.9e1_dp 01024 t57 = t35 + t41 - t48 - t56 01025 t61 = 0.2e1_dp / 0.3e1_dp * t16 01026 t62 = t61 - t25 01027 t64 = t61 - t29 01028 t66 = t7 * t57 - 0.2e1_dp / 0.3e1_dp * t16 * t40 + t62 * t46 +& 01029 & t64 * t45 01030 01031 IF (grad_deriv>=0 .AND. my_rho>epsilon_rho) THEN 01032 e_0(ii) = e_0(ii)& 01033 - ( 0.4e1_dp * t6 * t7 * t8 + t15 * t21 * t66 ) * sc 01034 END IF 01035 !-------- 01036 01037 t72 = t27 * my_rhoa 01038 t75 = t49 * t8 01039 t78 = 0.64e2_dp / 0.3e1_dp * t24 * t72 - t75 * t45 / 0.9e1_dp 01040 t82 = my_rhob * t57 + t7 * t78 - 0.2e1_dp * my_rhoa * t46 01041 t85 = t4 ** 2 01042 t86 = 0.1e1_dp / t85 01043 t87 = a * t86 01044 t88 = t87 * my_rhoa 01045 t90 = 0.1e1_dp / t1 / t16 01046 t92 = my_rhob * t90 * d 01047 t94 = 0.4e1_dp / 0.3e1_dp * t88 * t92 01048 t95 = 0.1e1_dp / t16 01049 t98 = 0.4e1_dp * t6 * t7 * t95 01050 t99 = t12 * c 01051 t100 = t16 ** 2 01052 t101 = t100 * my_rho 01053 t102 = 0.1e1_dp / t101 01054 t103 = t102 * t14 01055 t104 = t5 * t66 01056 t107 = t99 * t103 * t104 / 0.3e1_dp 01057 t108 = t86 * t102 01058 t109 = t66 * d 01059 t112 = t15 * t108 * t109 / 0.3e1_dp 01060 t115 = t5 / t18 / t100 01061 t118 = 0.11e2_dp / 0.3e1_dp * t15 * t115 * t66 01062 t120 = 0.1e1_dp / t1 / my_rho 01063 t124 = d ** 2 01064 t129 = c * t120 + d * t120 * t5 - t124 / t18 / my_rho * t86 01065 t130 = 0.7e1_dp / 0.54e2_dp * t129 01066 t132 = t129 / 0.54e2_dp 01067 t135 = -t129 / 0.3e1_dp 01068 t138 = my_rhoa * t95 01069 t140 = my_rhob * t95 01070 t142 = -t138 * t45 - t140 * t46 01071 t145 = t130 * t40 - t132 * t47 - t135 * t54 / 0.9e1_dp - t49 *& 01072 & t142 / 0.9e1_dp 01073 t153 = t7 * t145 - 0.4e1_dp / 0.3e1_dp * my_rho * t40 +& 01074 & 0.4e1_dp / 0.3e1_dp * my_rho * t46 + 0.4e1_dp / 0.3e1_dp& 01075 & * my_rho * t45 01076 t155 = t15 * t21 * t153 01077 01078 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 01079 e_ra(ii)=e_ra(ii)& 01080 - ( 0.4e1_dp * t6 * t52 + t15 * t21 * t82 + t94 - t98 +& 01081 & t107 + t112 - t118 + t155 ) * sc 01082 END IF 01083 01084 t159 = t31 * my_rhob 01085 t164 = 0.64e2_dp / 0.3e1_dp * t24 * t159 - t75 * t46 / 0.9e1_dp 01086 t168 = my_rhoa * t57 + t7 * t164 - 0.2e1_dp * my_rhob * t45 01087 01088 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 01089 e_rb(ii)=e_rb(ii)& 01090 - ( 0.4e1_dp * t6 * t50 + t15 * t21 * t168 + t94 - t98 +& 01091 & t107 + t112 - t118 + t155 ) * sc 01092 END IF 01093 01094 t171 = t39 * my_ndrho 01095 t176 = 0.2e1_dp * t7 * t171 - 0.4e1_dp / 0.3e1_dp * t16 * my_ndrho 01096 01097 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 01098 e_ndr(ii)=e_ndr(ii)& 01099 - ( t15 * t21 * t176 ) * sc 01100 END IF 01101 01102 t181 = t49 * my_rhoa 01103 t182 = t8 * my_ndrhoa 01104 t185 = -0.2e1_dp * t44 * my_ndrhoa - 0.2e1_dp / 0.9e1_dp * t181 * t182 01105 t189 = t7 * t185 + 0.2e1_dp * t64 * my_ndrhoa 01106 01107 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 01108 e_ndra(ii) = e_ndra(ii)& 01109 - ( t15 * t21 * t189 ) * sc 01110 END IF 01111 01112 t194 = t49 * my_rhob 01113 t195 = t8 * my_ndrhob 01114 t198 = -0.2e1_dp * t44 * my_ndrhob - 0.2e1_dp / 0.9e1_dp * t194 * t195 01115 t202 = t7 * t198 + 0.2e1_dp * t62 * my_ndrhob 01116 01117 IF (grad_deriv>=1.OR.grad_deriv==-1) THEN 01118 e_ndrb(ii) = e_ndrb(ii)& 01119 - ( t15 * t21 * t202 ) * sc 01120 END IF 01121 !------- 01122 01123 t205 = t100 * t16 01124 t206 = 0.1e1_dp / t205 01125 t210 = 0.26e2_dp / 0.9e1_dp * t99 * t206 * t14 * t104 01126 t214 = 0.2e1_dp / 0.3e1_dp * t99 * t103 * t5 * t153 01127 t215 = c ** 2 01128 t218 = 0.1e1_dp / t1 / t205 01129 t222 = t12 * t215 * t218 * t14 * t104 / 0.9e1_dp 01130 t228 = 0.2e1_dp / 0.9e1_dp * t12 * c * t218 * t14 * t86 * t109 01131 t232 = 0.26e2_dp / 0.9e1_dp * t15 * t86 * t206 * t109 01132 t236 = 0.2e1_dp / 0.3e1_dp * t15 * t108 * t153 * d 01133 t238 = 0.1e1_dp / t85 / t4 01134 t243 = 0.2e1_dp / 0.9e1_dp * t15 * t238 * t218 * t66 * t124 01135 t249 = 0.154e3_dp / 0.9e1_dp * t15 * t5 / t18 / t101 * t66 01136 t252 = 0.22e2_dp / 0.3e1_dp * t15 * t115 * t153 01137 t253 = t6 * t140 01138 t255 = t87 * t92 01139 t257 = c * t90 01140 t260 = d * t90 * t5 01141 t265 = t124 / t18 / t16 * t86 01142 t268 = 0.1e1_dp / t17 01143 t270 = t124 * d * t268 * t238 01144 t304 = t15 * t21 * (t7 * ((-0.14e2_dp / 0.81e2_dp * t257 -& 01145 & 0.14e2_dp / 0.81e2_dp * t260 + 0.7e1_dp / 0.27e2_dp *& 01146 & t265 - 0.7e1_dp / 0.81e2_dp * t270) * t40 - (-0.2e1_dp /& 01147 & 0.81e2_dp * t257 - 0.2e1_dp / 0.81e2_dp * t260 + t265 /& 01148 & 0.27e2_dp - t270 / 0.81e2_dp) * t47 - (0.4e1_dp /& 01149 & 0.9e1_dp * t257 + 0.4e1_dp / 0.9e1_dp * t260 - 0.2e1_dp & 01150 &/ 0.3e1_dp * t265 + 0.2e1_dp / 0.9e1_dp * t270) * t54 /& 01151 & 0.9e1_dp - 0.2e1_dp / 0.9e1_dp * t135 * t142 - t49 *& 01152 & (0.2e1_dp * my_rhoa * t268 * t45 + 0.2e1_dp * my_rhob *& 01153 & t268 * t46) / 0.9e1_dp) - 0.4e1_dp / 0.3e1_dp * t40 +& 01154 & 0.4e1_dp / 0.3e1_dp * t46 + 0.4e1_dp / 0.3e1_dp * t45) 01155 t310 = 0.40e2_dp / 0.9e1_dp * t88 * my_rhob / t1 / t17 * d 01156 t313 = my_rhob * t20 01157 t316 = 0.8e1_dp / 0.9e1_dp * a * t238 * my_rhoa * t313 * t124 01158 t319 = 0.8e1_dp * t6 * t7 * t268 01159 t322 = t99 * t103 * t5 * t82 01160 t326 = t15 * t108 * t82 * d 01161 t329 = t15 * t115 * t82 01162 t332 = t135 * t8 01163 t334 = t49 * t95 01164 t341 = t15 * t21 * (my_rhob * t145 + t7 * (-t332 * t45 /& 01165 & 0.9e1_dp + t334 * t45 / 0.9e1_dp)) 01166 01167 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 01168 e_ra_ra(ii) = e_ra_ra(ii)& 01169 + ( t210 - t214 - t222 - t228 + t232 - t236 - t243 - t249 +& 01170 & t252 + 0.8e1_dp * t253 - 0.8e1_dp / 0.3e1_dp * t255 -& 01171 & t304 + t310 - t316 - t319 - 0.2e1_dp / 0.3e1_dp * t322 -& 01172 & 0.2e1_dp / 0.3e1_dp * t326 + 0.22e2_dp / 0.3e1_dp * t329& 01173 & - 0.2e1_dp * t341 - t15 * t21 * (0.2e1_dp * my_rhob *& 01174 & t78 + 0.320e3_dp / 0.9e1_dp * t72 * my_rhob * t24 -& 01175 & 0.2e1_dp * t46) ) * sc 01176 END IF 01177 01178 t354 = t99 * t103 * t5 * t168 01179 t356 = t6 * t138 01180 t360 = t87 * my_rhoa * t90 * d 01181 t363 = t15 * t115 * t168 01182 t373 = t15 * t21 * (my_rhoa * t145 + t7 * (-t332 * t46 /& 01183 & 0.9e1_dp + t334 * t46 / 0.9e1_dp)) 01184 t376 = t15 * t108 * t168 * d 01185 01186 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 01187 e_rb_rb(ii)=e_rb_rb(ii)& 01188 + ( t210 - t214 - t222 - t228 + t232 - t236 - t243 - t249 +& 01189 & t252 - t304 + t310 - t316 - t319 + 0.8e1_dp * t356 -& 01190 & 0.8e1_dp / 0.3e1_dp * t360 + 0.22e2_dp / 0.3e1_dp *& 01191 & t363 - 0.2e1_dp * t373 - 0.2e1_dp / 0.3e1_dp * t354 -& 01192 & 0.2e1_dp / 0.3e1_dp * t376 - t15 * t21 * (0.2e1_dp *& 01193 & my_rhoa * t164 + 0.320e3_dp / 0.9e1_dp * my_rhoa *& 01194 & t159 * t24 - 0.2e1_dp * t45) ) * sc 01195 END IF 01196 01197 t381 = -t354 / 0.3e1_dp + 0.4e1_dp * t356 - 0.4e1_dp / 0.3e1_dp& 01198 & * t360 + 0.11e2_dp / 0.3e1_dp * t363 - t373 - t341 -& 01199 & t376 / 0.3e1_dp + t310 - 0.4e1_dp * t6 * t8 + 0.4e1_dp *& 01200 & t253 + t210 - t214 - t222 01201 t391 = -t228 + t232 - t236 - t243 - t249 + t252 - 0.4e1_dp /& 01202 & 0.3e1_dp * t255 - t304 - t316 - t319 - t322 / 0.3e1_dp -& 01203 & t326 / 0.3e1_dp + 0.11e2_dp / 0.3e1_dp * t329 - t15 *& 01204 & t21 * (t35 + t41 - t48 - t56 + my_rhob * t164 + my_rhoa & 01205 &* t78) 01206 01207 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 01208 e_ra_rb(ii)=e_ra_rb(ii)& 01209 + ( t381 + t391 ) * sc 01210 END IF 01211 01212 t408 = t12 * t14 * t5 01213 t415 = t99 * t103 * t5 * t176 / 0.3e1_dp 01214 t419 = t15 * t108 * t176 * d / 0.3e1_dp 01215 t422 = 0.11e2_dp / 0.3e1_dp * t15 * t115 * t176 01216 t430 = t15 * t21 * (0.2e1_dp * t7 * t130 * my_ndrho - 0.8e1_dp & 01217 &/ 0.3e1_dp * my_rho * my_ndrho) 01218 01219 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 01220 e_ndr_ra(ii) = e_ndr_ra(ii)& 01221 - ( 0.2e1_dp * t408 * t313 * t171 + t415 + t419 - t422 +& 01222 & t430 ) * sc 01223 e_ndr_rb(ii) = e_ndr_rb(ii)& 01224 - ( 0.2e1_dp * t408 * t20 * my_rhoa * t171 + t415 + t419 -& 01225 & t422 + t430 ) * sc 01226 END IF 01227 01228 t445 = t99 * t103 * t5 * t189 / 0.3e1_dp 01229 t449 = t15 * t108 * t189 * d / 0.3e1_dp 01230 t452 = 0.11e2_dp / 0.3e1_dp * t15 * t115 * t189 01231 t467 = t15 * t21 * (t7 * (-0.2e1_dp * t132 * my_ndrhoa -& 01232 & 0.2e1_dp / 0.9e1_dp * t135 * my_rhoa * t182 + 0.2e1_dp /& 01233 & 0.9e1_dp * t181 * t95 * my_ndrhoa) + 0.8e1_dp / 0.3e1_dp& 01234 & * my_rho * my_ndrhoa) 01235 01236 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 01237 e_ndra_ra(ii)=e_ndra_ra(ii)& 01238 - ( t15 * t21 * (my_rhob * t185 - 0.2e1_dp / 0.9e1_dp * t7& 01239 & * t75 * my_ndrhoa) + t445 + t449 - t452 + t467 ) * sc 01240 e_ndra_rb(ii)=e_ndra_rb(ii)& 01241 - ( t15 * t21 * (my_rhoa * t185 - 0.4e1_dp * my_rhob *& 01242 & my_ndrhoa) + t445 + t449 - t452 + t467 ) * sc 01243 END IF 01244 01245 t483 = t99 * t103 * t5 * t202 / 0.3e1_dp 01246 t487 = t15 * t108 * t202 * d / 0.3e1_dp 01247 t490 = 0.11e2_dp / 0.3e1_dp * t15 * t115 * t202 01248 t505 = t15 * t21 * (t7 * (-0.2e1_dp * t132 * my_ndrhob -& 01249 & 0.2e1_dp / 0.9e1_dp * t135 * my_rhob * t195 + 0.2e1_dp /& 01250 & 0.9e1_dp * t194 * t95 * my_ndrhob) + 0.8e1_dp / 0.3e1_dp& 01251 & * my_rho * my_ndrhob) 01252 01253 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 01254 e_ndrb_ra(ii) =e_ndrb_ra(ii)& 01255 - ( t15 * t21 * (my_rhob * t198 - 0.4e1_dp * my_rhoa *& 01256 & my_ndrhob) + t483 + t487 - t490 + t505 ) * sc 01257 e_ndrb_rb(ii)=e_ndrb_rb(ii)& 01258 - ( t15 * t21 * (my_rhoa * t198 - 0.2e1_dp / 0.9e1_dp * t7& 01259 & * t75 * my_ndrhob) + t483 + t487 - t490 + t505 ) * sc 01260 END IF 01261 01262 t515 = 0.4e1_dp / 0.3e1_dp * t16 01263 01264 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 01265 e_ndr_ndr(ii) = e_ndr_ndr(ii)& 01266 - ( t15 * t21 * (0.2e1_dp * t7 * t39 - t515) ) * sc 01267 END IF 01268 01269 t519 = t13 / 0.9e1_dp 01270 t520 = t37 / 0.9e1_dp 01271 01272 IF (grad_deriv>=2.OR.grad_deriv==-2) THEN 01273 e_ndra_ndra(ii) = e_ndra_ndra(ii)& 01274 - ( t15 * t21 * (t7 * (-0.5e1_dp + t519 + t520 - 0.2e1_dp & 01275 &/ 0.9e1_dp * t181 * t8) + t515 - 0.2e1_dp * t29) ) * sc 01276 01277 e_ndrb_ndrb(ii) = e_ndrb_ndrb(ii)& 01278 - ( t15 * t21 * (t7 * (-0.5e1_dp + t519 + t520 - 0.2e1_dp & 01279 &/ 0.9e1_dp * t194 * t8) + t515 - 0.2e1_dp * t25) ) * sc 01280 END IF 01281 END IF 01282 END DO 01283 !$omp end do 01284 01285 END SUBROUTINE lyp_lsd_calc 01286 01287 END MODULE xc_lyp
1.7.3