From 4d6e969a54bc411b3ae1f779986099b69c1978a6 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 15 Jul 2019 12:07:54 +0200 Subject: [PATCH] ECMD --- input/basis | 34 ++++-- input/molecule | 4 +- input/weight | 34 ++++-- src/QuAcK/LSD_SR.f | 23 +--- src/QuAcK/basis_correction.f90 | 2 +- src/QuAcK/ec_srlda.f90 | 21 ++-- src/QuAcK/ecmd_lda.f90 | 200 +++++++++++++++++++++++++++++++++ src/QuAcK/fc_srlda.f90 | 36 +++++- 8 files changed, 303 insertions(+), 51 deletions(-) create mode 100644 src/QuAcK/ecmd_lda.f90 diff --git a/input/basis b/input/basis index b9ca7b5..eb58543 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,29 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 +1 6 +S 8 1.00 + 17880.0000000 0.0007380 + 2683.0000000 0.0056770 + 611.5000000 0.0288830 + 173.5000000 0.1085400 + 56.6400000 0.2909070 + 20.4200000 0.4483240 + 7.8100000 0.2580260 + 1.6530000 0.0150630 +S 8 1.00 + 17880.0000000 -0.0001720 + 2683.0000000 -0.0013570 + 611.5000000 -0.0067370 + 173.5000000 -0.0276630 + 56.6400000 -0.0762080 + 20.4200000 -0.1752270 + 7.8100000 -0.1070380 + 1.6530000 0.5670500 S 1 1.00 - 0.2976000 1.0000000 + 0.4869000 1.0000000 +P 3 1.00 + 28.3900000 0.0460870 + 6.2700000 0.2401810 + 1.6950000 0.5087440 P 1 1.00 - 1.2750000 1.0000000 + 0.4317000 1.0000000 +D 1 1.00 + 2.2020000 1.0000000 diff --git a/input/molecule b/input/molecule index c78e87e..edeba31 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 1 5 5 0 0 # Znuc x y z - He 0.0 0.0 0.0 + Ne 0.0 0.0 0.0 diff --git a/input/weight b/input/weight index b9ca7b5..eb58543 100644 --- a/input/weight +++ b/input/weight @@ -1,9 +1,29 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 +1 6 +S 8 1.00 + 17880.0000000 0.0007380 + 2683.0000000 0.0056770 + 611.5000000 0.0288830 + 173.5000000 0.1085400 + 56.6400000 0.2909070 + 20.4200000 0.4483240 + 7.8100000 0.2580260 + 1.6530000 0.0150630 +S 8 1.00 + 17880.0000000 -0.0001720 + 2683.0000000 -0.0013570 + 611.5000000 -0.0067370 + 173.5000000 -0.0276630 + 56.6400000 -0.0762080 + 20.4200000 -0.1752270 + 7.8100000 -0.1070380 + 1.6530000 0.5670500 S 1 1.00 - 0.2976000 1.0000000 + 0.4869000 1.0000000 +P 3 1.00 + 28.3900000 0.0460870 + 6.2700000 0.2401810 + 1.6950000 0.5087440 P 1 1.00 - 1.2750000 1.0000000 + 0.4317000 1.0000000 +D 1 1.00 + 2.2020000 1.0000000 diff --git a/src/QuAcK/LSD_SR.f b/src/QuAcK/LSD_SR.f index e9f627d..72ffe09 100644 --- a/src/QuAcK/LSD_SR.f +++ b/src/QuAcK/LSD_SR.f @@ -1,25 +1,3 @@ -ccc example: use the subroutine lsdsr to compute the complementary -ccc short-range exchange-correlation energy 'excsr' and -ccc the corresponding up and down potentials 'vxcsrup','vxcsrdown' -ccc at polarization z=0.7, cutoff mu=0.5, and for 0.2 < rs < 20, -ccc and write them on a file -c program testex -c implicit none -c double precision z,rs,mu -c double precision excsr,vxcsrup,vxcsrdown -c integer i -c open(9,file='testex',status='unknown') -c z=0.7d0 -c mu=0.5d0 -c do i=1,100 -c rs=0.2*i -c call lsdsr(rs,z,mu,excsr,vxcsrup,vxcsrdown) -c write(9,*) rs,excsr,vxcsrup,vxcsrdown -c enddo -c stop -c end - - subroutine lsdsr(rs,z,mu,excsr,vxcsrup,vxcsrdown) ccc Hartree atomic units used ccc for given density parameter 'rs', spin polarization 'z' @@ -30,6 +8,7 @@ ccc interacting electron gas) => 'excsr' ccc and the corresponding exchange-correlation potentials for ccc spin-up and spin-down electrons => 'vxcsrup','vxcsrdown' ccc from Paziani, Moroni, Gori-Giorgi, and Bachelet, cond-mat/0601353 + implicit none double precision rs,z,mu,excsr,vxcsrup,vxcsrdown double precision eclr,exlr,ec,ecd,ecz,ex diff --git a/src/QuAcK/basis_correction.f90 b/src/QuAcK/basis_correction.f90 index ce4351c..44b7ec7 100644 --- a/src/QuAcK/basis_correction.f90 +++ b/src/QuAcK/basis_correction.f90 @@ -50,7 +50,7 @@ subroutine basis_correction(nBas,nO,nShell,CenterShell,TotAngMomShell,KShell,DSh ! Set quadrature grid - SGn = 0 + SGn = 1 call read_grid(SGn,nRad,nAng,nGrid) diff --git a/src/QuAcK/ec_srlda.f90 b/src/QuAcK/ec_srlda.f90 index 3dccbb7..96ffd03 100644 --- a/src/QuAcK/ec_srlda.f90 +++ b/src/QuAcK/ec_srlda.f90 @@ -15,31 +15,38 @@ subroutine ec_srlda(nGrid,weight,rho,mu) ! Local variables integer :: iG - double precision :: r + double precision :: r,ra,rb double precision :: rs - double precision :: ecsr - double precision :: ec,vcup,vcdw + double precision :: ec_lda,ecmd_lda + double precision :: ec,ecmd,vcup,vcdw ! Initialization - ecsr = 0d0 + ec_lda = 0d0 + ecmd_lda = 0d0 do iG=1,ngrid r = max(0d0,rho(iG)) + ra = 0.5d0*r + rb = 0.5d0*r if(r > threshold) then rs = (4d0*pi*r/3d0)**(-1d0/3d0) call lsdsr(rs,0d0,mu(iG),ec,vcup,vcdw) - - ecsr = ecsr + weight(iG)*ec*r + call ESRC_MD_LDAERF(mu(iG),ra,rb,.true.,ecmd) + ec_lda = ec_lda + weight(iG)*ec*r + ecmd_lda = ecmd_lda + weight(iG)*ecmd*r end if end do - write(*,'(A32,1X,F16.10)') 'ecsr = ',ecsr + write(*,*) + write(*,'(A32,1X,F16.10)') 'Ec(sr-LDA) = ',ec_lda + write(*,'(A32,1X,F16.10)') 'Ecmd(sr-LDA) = ',ec_lda + ecmd_lda + write(*,*) end subroutine ec_srlda diff --git a/src/QuAcK/ecmd_lda.f90 b/src/QuAcK/ecmd_lda.f90 new file mode 100644 index 0000000..a7fa705 --- /dev/null +++ b/src/QuAcK/ecmd_lda.f90 @@ -0,0 +1,200 @@ +!**************************************************************************** + subroutine ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,e) +!***************************************************************************** +! Short-range spin-dependent LDA correlation functional with multideterminant reference +! for OEP calculations from Section V of +! Paziani, Moroni, Gori-Giorgi and Bachelet, PRB 73, 155111 (2006) +! +! Input: rhot : total density +! rhos : spin density +! mu : Interation parameter +! dospin : use spin density +! +! Ouput: e : energy +! +! Created: 26-08-11, J. Toulouse +!***************************************************************************** + implicit none + + double precision, intent(in) :: rho_a,rho_b,mu + logical, intent(in) :: dospin + double precision, intent(out):: e + + double precision :: e1 + double precision :: rhoa,rhob + double precision :: rhot, rhos + rhoa=max(rho_a,1.0d-15) + rhob=max(rho_b,1.0d-15) + rhot = rhoa + rhob + rhos = rhoa - rhob + + call ec_only_lda_sr(mu,rho_a,rho_b,e1) + if(isnan(e1))then + print*,'e1 is NaN' + print*,mu,rho_a,rho_b + stop + endif + e1 = 0d0 + call DELTA_LRSR_LDAERF (rhot,rhos,mu,dospin,e) + if(isnan(e))then + print*,'e is NaN' + print*,mu,rhot,rhos + stop + endif + e = e1 + e + + end + +!**************************************************************************** + subroutine DELTA_LRSR_LDAERF (rhot,rhos,mu,dospin,e) +!***************************************************************************** +! LDA approximation to term Delta_(LR-SR) from Eq. 42 of +! Paziani, Moroni, Gori-Giorgi and Bachelet, PRB 73, 155111 (2006) +! +! Input: rhot : total density +! rhos : spin density +! mu : Interation parameter +! dospin : use spin density +! +! Ouput: e : energy +! +! Warning: not tested for z != 0 +! +! Created: 26-08-11, J. Toulouse +!***************************************************************************** + implicit none + + double precision rhot, rhos, mu + logical dospin + double precision e + + double precision f13, f83, pi, rsfac, alpha2 + double precision rs, rs2, rs3 + + double precision rhoa, rhob, z, z2, onepz, onemz, zp, zm, phi8 + double precision g0s,g0f + double precision bd2, bd3 + double precision c45, c4, c5 + double precision bc2, bc4, bc3t, bc5t, d0 + double precision delta2,delta3,delta4,delta5,delta6 + double precision delta + + parameter(f13 = 0.333333333333333d0) + parameter(f83 = 2.6666666666666665d0) + parameter(pi = 3.141592653589793d0) + parameter(rsfac = 0.620350490899400d0) + parameter(alpha2 = 0.2715053589826032d0) + + rs = rsfac/(rhot**f13) + rs2 = rs*rs + rs3 = rs2*rs + +! Spin-unpolarized case + if (.not.dospin) then + z = 0.d0 + +! Spin-polarized case + else + rhoa=max((rhot+rhos)*.5d0,1.0d-15) + rhob=max((rhot-rhos)*.5d0,1.0d-15) + z=min((rhoa-rhob)/(rhoa+rhob),0.9999999999d0) + endif + + z2=z*z + + bd2=dexp(-0.547d0*rs)*(-0.388d0*rs+0.676*rs2)/rs2 + bd3=dexp(-0.31d0*rs)*(-4.95d0*rs+rs2)/rs3 + + onepz=1.d0+z + onemz=1.d0-z + phi8=0.5d0*(onepz**f83+onemz**f83) + + zp=onepz/2.d0 + zm=onemz/2.d0 + c45=(zp**2)*g0s(rs*zp**(-f13))+(zm**2)*g0s(rs*zm**(-f13)) + c4=c45+(1.d0-z2)*bd2-phi8/(5.d0*alpha2*rs2) + c5=c45+(1.d0-z2)*bd3 + + bc2=-3.d0*(1-z2)*(g0f(rs)-0.5d0)/(8.d0*rs3) + bc4=-9.d0*c4/(64.d0*rs3) + bc3t=-(1-z2)*g0f(rs)*(2.d0*dsqrt(2.d0)-1)/(2.d0*dsqrt(pi)*rs3) + bc5t = -3.d0*c5*(3.d0-dsqrt(2.d0))/(20.d0*dsqrt(2.d0*pi)*rs3) + + d0=(0.70605d0+0.12927d0*z2)*rs + delta2=0.073867d0*(rs**(1.5d0)) + delta3=4*(d0**6)*bc3t+(d0**8)*bc5t; + delta4=4*(d0**6)*bc2+(d0**8)*bc4; + delta5=(d0**8)*bc3t; + delta6=(d0**8)*bc2; + delta=(delta2*(mu**2)+delta3*(mu**3)+delta4*(mu**4)+delta5*(mu**5)+delta6*(mu**6))/((1+(d0**2)*(mu**2))**4) + + +! multiply by rhot to get energy density +! e=delta*rhot + e=delta + + end + +!***************************************************************************** + double precision function g0s(rs) +!***************************************************************************** +! g"(0,rs,z=1) from Eq. 32 of +! Paziani, Moroni, Gori-Giorgi and Bachelet, PRB 73, 155111 (2006) +! +! Created: 26-08-11, J. Toulouse +!***************************************************************************** + implicit none + double precision rs + double precision rs2, f53, alpha2 + parameter(f53 = 1.6666666666666667d0) + parameter(alpha2 = 0.2715053589826032d0) + rs2=rs*rs + g0s=(2.d0**f53)*(1.d0-0.02267d0*rs)/((5.d0*alpha2*rs2)*(1.d0+0.4319d0*rs+0.04d0*rs2)) + end + + double precision function g0f(x) +!cc on-top pair-distribution function +!cc Gori-Giorgi and Perdew, PRB 64, 155102 (2001) +!cc x -> rs + implicit none + double precision C0f,D0f,E0f,F0f,x + C0f = 0.0819306d0 + D0f = 0.752411d0 + E0f = -0.0127713d0 + F0f = 0.00185898d0 + g0f=(1.d0-(0.7317d0-D0f)*x+C0f*x**2+E0f*x**3+ & + F0f*x**4)*exp(-abs(D0f)*x)/2.d0 + return + end + +subroutine ec_only_lda_sr(mu,rho_a,rho_b,ec) + implicit none + include 'parameters.h' + + double precision, intent(out) :: ec + double precision, intent(in) :: mu,rho_a,rho_b + +! Double precision numbers + + double precision :: rsfac,rho,rs,rhoa,rhob,z + double precision :: eccoul, ecd, ecz, ecdd, eczd + double precision :: eclr + rsfac = (3.0d0/(4.0d0*pi))**(1d0/3d0) + + ec = 0.d0 +! Test on density + rho = rho_a + rho_b + if (dabs(rho).ge.1.d-12) then + + rs=rsfac/(rho**(1d0/3d0)) + rhoa=max(rho_a,1.0d-15) + rhob=max(rho_b,1.0d-15) + z=(rhoa-rhob)/(rhoa+rhob) + + call ecPW(rs,z,eccoul,ecd,ecz,ecdd,eczd) + call ecorrlr(rs,z,mu,eclr) + ec=(eccoul-eclr)*rho + + endif + +end diff --git a/src/QuAcK/fc_srlda.f90 b/src/QuAcK/fc_srlda.f90 index 139642b..4af54ae 100644 --- a/src/QuAcK/fc_srlda.f90 +++ b/src/QuAcK/fc_srlda.f90 @@ -18,9 +18,10 @@ subroutine fc_srlda(nBas,nGrid,weight,MO,rho,mu,eG0W0) ! Local variables integer :: iG,p - double precision :: r - double precision :: rs - double precision :: ec,vcup,vcdw + double precision :: r,ra,rb,rap,ram + double precision :: rs,rsp,rsm + double precision :: ec,ecp,ecm,vcup,vcdw + double precision,parameter :: delta = 1d-6 double precision,allocatable :: de(:) ! Memory allocation @@ -33,13 +34,38 @@ subroutine fc_srlda(nBas,nGrid,weight,MO,rho,mu,eG0W0) do iG=1,ngrid - r = max(0d0,rho(iG)) + r = max(0d0,rho(iG)) + ra = 0.5d0*r + rb = 0.5d0*r if(r > threshold) then rs = (4d0*pi*r/3d0)**(-1d0/3d0) - call lsdsr(rs,0d0,mu(iG),ec,vcup,vcdw) +! call lsdsr(rs,0d0,mu(iG),ec,vcup,vcdw) + if(abs(ra) > delta) then + + rap = ra + delta + ram = ra - delta + + rsp = (4d0*pi*rap/3d0)**(-1d0/3d0) + rsm = (4d0*pi*ram/3d0)**(-1d0/3d0) + +! call lsdsr(rsp,0d0,mu(iG),ecp,vcup,vcdw) +! call lsdsr(rsm,0d0,mu(iG),ecm,vcup,vcdw) + call lsdsr(rs,0d0,mu(iG),ec,vcup,vcdw) + + +! call ESRC_MD_LDAERF(mu(iG),rap,rb,.true.,ecp) +! call ESRC_MD_LDAERF(mu(iG),ram,rb,.true.,ecm) + +! vcup = (ecp - ecm)/(2d0*delta) + + else + + vcup = 0d0 + + end if do p=1,nBas