LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 03:23:29 +01:00
Code Releases Activity

modifs in ec_lyp

Browse Source
This commit is contained in:
Emmanuel Giner LCT 2019-06-08 16:39:12 +02:00
parent 71c186e67a
commit ff4f6d1c98
1 changed files with 16 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
src/dft_utils_one_e/ec_lyp.irp.f
View File

@ -42,7 +42,7 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
double precision :: a,b,c,d,c_f,omega,delta
double precision :: rho_13,rho_inv_13,rho_83,rho_113,rho_inv_113,denom
double precision :: thr,huge_num
double precision :: thr,huge_num,rho_inv
double precision :: cst_2_113,cst_8_3,rho_2,rho_a_2,rho_b_2
double precision :: tmp1,tmp2,tmp3,tmp4
double precision :: big1,big2,big3
@ -57,24 +57,34 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
c = 0.2533d0
d = 0.349d0
thr = 1d-15
thr = 1d-10
huge_num = 1.d0/thr
if(rho.lt.0.d0)then
print*,'pb !! rho.lt.0.d0'
stop
endif
rho_13 = rho**(1d0/3d0)
rho_113 = rho**(11d0/3d0)
if(abs(rho_13) < thr) then
if(dabs(rho_13) < thr) then
rho_inv_13 = huge_num
else
rho_inv_13 = 1.d0/rho_13
endif
if (abs(rho_113) < thr) then
if (dabs(rho_113) < thr) then
rho_inv_113 = huge_num
else
rho_inv_113 = 1d0/rho_113
endif
if (dabs(rho) < thr) then
rho_inv = huge_num
else
rho_inv = 1d0/rho
endif
! Useful quantities to predefine
denom = 1d0/(1d0 + d*rho_inv_13)
@ -91,12 +101,12 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
! first term in the equation (2) of Preuss CPL, 1989
big1 = 4d0*denom*rho_a*rho_b/rho
big1 = 4d0*denom*rho_a*rho_b*rho_inv
tmp1 = cst_2_113*c_f*(rho_a**cst_8_3 + rho_b**cst_8_3)
tmp2 = (47d0/18d0 - 7d0/18d0*delta)*grad_rho_2
tmp3 = - (5d0/2d0 - 1.d0/18d0*delta)*(grad_rho_a_2 + grad_rho_b_2)
tmp4 = - (delta - 11d0)/9d0*(rho_a/rho*grad_rho_a_2 + rho_b/rho*grad_rho_b_2)
tmp4 = - (delta - 11d0)/9d0*(rho_a*rho_inv*grad_rho_a_2 + rho_b*rho_inv*grad_rho_b_2)
big2 = rho_a*rho_b*(tmp1 + tmp2 + tmp3 + tmp4)
tmp1 = -2d0/3d0*rho_2*grad_rho_2