10
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00

remove obselete routines

This commit is contained in:
Pierre-Francois Loos 2022-11-29 09:35:50 +01:00
parent 99aa500483
commit 16fe319ef3
28 changed files with 0 additions and 5509 deletions

View File

@ -1,515 +0,0 @@
subroutine lsdsr(rs,z,mu,excsr,vxcsrup,vxcsrdown)
ccc Hartree atomic units used
ccc for given density parameter 'rs', spin polarization 'z'
ccc and cutoff parameter 'mu'
ccc gives the complementary short-range exchange-correlation
ccc energy (i.e., xc energy of jellium minus xc energy of long-range
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
double precision vclrup,vclrdown,vxlrup,vxlrdown
double precision vxup,vxdown,vcup,vcdown
double precision pi,alpha,cf
pi=dacos(-1.d0)
alpha=(4.d0/9.d0/pi)**(1.d0/3.d0)
cf=1.d0/alpha
ex=-3.d0*cf/rs/8.d0/pi*((1.d0+z)**(4.d0/3.d0)+
$ (1.d0-z)**(4.d0/3.d0))
ex = 0d0
vxup=-(1.d0+z)**(1.d0/3.d0)*(3.d0/2.d0/pi)**(2.d0/3.d0)/rs
vxdown=-(1.d0-z)**(1.d0/3.d0)*(3.d0/2.d0/pi)**(2.d0/3.d0)/rs
vxup = 0d0
vxdown = 0d0
call ecPW(rs,z,ec,ecd,ecz)
vcup=ec-rs/3.d0*ecd-(z-1.d0)*ecz
vcdown=ec-rs/3.d0*ecd-(z+1.d0)*ecz
call exchangelr(rs,z,mu,exlr)
exlr = 0d0
call vexchangelr(rs,z,mu,vxlrup,vxlrdown)
vxlrup = 0d0
vxlrdown = 0d0
call ecorrlr(rs,z,mu,eclr)
call vcorrlr(rs,z,mu,vclrup,vclrdown)
excsr=ex+ec-(exlr+eclr)
vxcsrup=vxup+vcup-(vxlrup+vclrup)
vxcsrdown=vxdown+vcdown-(vxlrdown+vclrdown)
return
end
subroutine ecorrlr(rs,z,mu,eclr)
ccc Hartree atomic units used
ccc for given density parameter rs, spin polarization z
ccc and cutoff parameter mu
ccc gives the correlation energy of the LR gas
ccc => eclr
implicit none
double precision rs,z,mu,eclr,ec,ecd,ecz
double precision pi,alpha,cf,phi
double precision g0,dpol,d2anti,d3anti,Qrpa
double precision coe2,coe3,coe4,coe5
double precision a1,a2,a3,a4,b0
double precision q1a,q2a,q3a,t1a,t2a,t3a,adib
pi=dacos(-1.d0)
alpha=(4.d0/9.d0/pi)**(1.d0/3.d0)
cf=1.d0/alpha
phi=((1.d0+z)**(2.d0/3.d0)+(1.d0-z)**(2.d0/3.d0))/2.d0
cc parameters from the fit
adib = 0.784949d0
q1a = -0.388d0
q2a = 0.676d0
q3a = 0.547d0
t1a = -4.95d0
t2a = 1.d0
t3a = 0.31d0
b0=adib*rs
d2anti=(q1a*rs+q2a*rs**2)*exp(-abs(q3a)*rs)/rs**2
d3anti=(t1a*rs+t2a*rs**2)*exp(-abs(t3a)*rs)/rs**3
coe2=-3.d0/8.d0/rs**3*(1.d0-z**2)*(g0(rs)-0.5d0)
coe3=-(1.d0-z**2)*g0(rs)/(sqrt(2.d0*pi)*rs**3)
if(abs(z).eq.1.d0) then
coe4=-9.d0/64.d0/rs**3*(dpol(rs)
$ -cf**2*2**(5.d0/3.d0)/5.d0/rs**2)
coe5=-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*dpol(rs)
else
coe4=-9.d0/64.d0/rs**3*(((1.d0+z)/2.d0)**2*
$ dpol(rs*(2/(1.d0+z))**(1.d0/3.d0))+((1.d0-z)/2.d0)**2
$ *dpol(rs*(2.d0/(1.d0-z))**(1.d0/3.d0))+
$ (1.-z**2)*d2anti-cf**2/10.d0*((1.d0+z)**(8.d0/3.d0)
$ +(1.-z)**(8.d0/3.d0))/rs**2)
coe5=-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*(((1.d0+z)/2.d0)**2
$ *dpol(rs*(2.d0/(1.d0+z))**(1.d0/3.d0))+((1.d0-z)/2.d0)**2
$ *dpol(rs*(2.d0/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)*
$ d3anti)
endif
call ecPW(rs,z,ec,ecd,ecz)
a1=4.d0*b0**6*coe3+b0**8*coe5
a2=4.d0*b0**6*coe2+b0**8*coe4+6.d0*b0**4*ec
a3=b0**8*coe3
a4=b0**6*(b0**2*coe2+4.d0*ec)
eclr=(phi**3*Qrpa(mu*sqrt(rs)/phi)+a1*mu**3+a2*mu**4+a3*mu**5+
$ a4*mu**6+b0**8*mu**8*ec)/((1.d0+b0**2*mu**2)**4)
return
end
subroutine vcorrlr(rs,z,mu,vclrup,vclrdown)
ccc Hartree atomic units used
ccc for given density parameter rs, spin polarization z
ccc and cutoff mu it gives the correlation LSD potential for LR interaction
ccc => vclrup (spin-up electrons), vclrdown (spin-down electrons)
implicit none
double precision rs,z,mu,eclr,eclrrs,eclrz,vclrup,vclrdown
double precision ec,ecd,ecz
double precision pi,alpha,cf,phi
double precision g0,dpol,d2anti,d3anti,Qrpa
double precision g0d,dpold,d2antid,d3antid,Qrpad,x
double precision coe2,coe3,coe4,coe5
double precision coe2rs,coe3rs,coe4rs,coe5rs
double precision coe2z,coe3z,coe4z,coe5z
double precision a1,a2,a3,a4,a5,b0,a1rs,a2rs,a3rs,a4rs,a5rs,
$ b0rs,a1z,a2z,a3z,a4z,a5z,b0z
double precision q1a,q2a,q3a,t1a,t2a,t3a,adib
pi=dacos(-1.d0)
alpha=(4.d0/9.d0/pi)**(1.d0/3.d0)
cf=1.d0/alpha
phi=((1.d0+z)**(2.d0/3.d0)+(1.d0-z)**(2.d0/3.d0))/2.d0
cc parameters from the fit
adib = 0.784949d0
q1a = -0.388d0
q2a = 0.676d0
q3a = 0.547d0
t1a = -4.95d0
t2a = 1.d0
t3a = 0.31d0
b0=adib*rs
d2anti=(q1a+q2a*rs)*exp(-q3a*rs)/rs
d3anti=(t1a+t2a*rs)*exp(-t3a*rs)/rs**2
d2antid=-((q1a + q1a*q3a*rs + q2a*q3a*rs**2)/
- rs**2)*exp(-q3a*rs)
d3antid=-((rs*t2a*(1 + rs*t3a) + t1a*(2 + rs*t3a))/
- rs**3)*exp(-rs*t3a)
coe2=-3.d0/8.d0/rs**3*(1.d0-z**2)*(g0(rs)-0.5d0)
coe2rs=-3.d0/8.d0/rs**3*(1.d0-z**2)*g0d(rs)+
$ 9.d0/8.d0/rs**4*(1.d0-z**2)*(g0(rs)-0.5d0)
coe2z=-3.d0/8.d0/rs**3*(-2.d0*z)*(g0(rs)-0.5d0)
coe3=-(1.d0-z**2)*g0(rs)/(sqrt(2.d0*pi)*rs**3)
coe3rs=-(1.d0-z**2)*g0d(rs)/(sqrt(2.d0*pi)*rs**3)+
$ 3.d0*(1.d0-z**2)*g0(rs)/(sqrt(2.d0*pi)*rs**4)
coe3z=2.d0*z*g0(rs)/(sqrt(2.d0*pi)*rs**3)
if(abs(z).eq.1.d0) then
coe4=-9.d0/64.d0/rs**3*(dpol(rs)
$ -cf**2*2**(5.d0/3.d0)/5.d0/rs**2)
coe4rs=-3.d0/rs*coe4-9.d0/64.d0/rs**3*(dpold(rs)
$ +2.d0*cf**2*2**(5.d0/3.d0)/5.d0/rs**3)
coe4z=-9.d0/64.d0/rs**3*(dpol(rs)-rs/6.d0*dpold(rs)-2.d0*d2anti
$ -4.d0/15.d0/rs**2*cf**2*2.d0**(5.d0/3.d0))*z
coe5=-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*dpol(rs)
coe5rs=-3.d0/rs*coe5-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*dpold(rs)
coe5z=-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*(dpol(rs)-rs/6.d0*
$ dpold(rs)-2.d0*d3anti)*z
else
coe4=-9.d0/64.d0/rs**3*(((1.d0+z)/2.d0)**2*
$ dpol(rs*(2/(1.d0+z))**(1.d0/3.d0))+((1.d0-z)/2.d0)**2
$ *dpol(rs*(2.d0/(1.d0-z))**(1.d0/3.d0))+
$ (1.-z**2)*d2anti-cf**2/10.d0*((1.d0+z)**(8.d0/3.d0)
$ +(1.-z)**(8.d0/3.d0))/rs**2)
coe4rs=-3.d0/rs*coe4-9.d0/64.d0/rs**3*(
$ ((1.d0+z)/2.d0)**(5.d0/3.d0)*dpold(rs*(2/(1.d0+z))**
$ (1.d0/3.d0))+((1.d0-z)/2.d0)**(5.d0/3.d0)*
$ dpold(rs*(2/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)*
$ d2antid+cf**2/5.d0*((1.d0+z)**(8.d0/3.d0)
$ +(1.d0-z)**(8.d0/3.d0))/rs**3)
coe4z=-9.d0/64.d0/rs**3*(1.d0/2.d0*(1.d0+z)*
$ dpol(rs*(2/(1.d0+z))**(1.d0/3.d0))-1.d0/2.d0*(1.d0-z)*
$ dpol(rs*(2/(1.d0-z))**(1.d0/3.d0))-rs/6.d0*
$ ((1.d0+z)/2.d0)**(2.d0/3.d0)*dpold(rs*(2/(1.d0+z))
$ **(1.d0/3.d0))+rs/6.d0*((1.d0-z)/2.d0)**(2.d0/3.d0)
$ *dpold(rs*(2/(1.d0-z))**(1.d0/3.d0))-2.d0*z*d2anti-
$ 4.d0/15.d0/rs**2*cf**2*((1.d0+z)**(5.d0/3.d0)-
$ (1.d0-z)**(5.d0/3.d0)))
coe5=-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*(((1.d0+z)/2.d0)**2
$ *dpol(rs*(2.d0/(1.d0+z))**(1.d0/3.d0))+((1.d0-z)/2.d0)**2
$ *dpol(rs*(2.d0/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)*
$ d3anti)
coe5rs=-3.d0/rs*coe5-9.d0/(40.d0*sqrt(2.d0*pi)*rs**3)*(
$ ((1.d0+z)/2.d0)**(5.d0/3.d0)*dpold(rs*(2/(1.d0+z))**
$ (1.d0/3.d0))+((1.d0-z)/2.d0)**(5.d0/3.d0)*
$ dpold(rs*(2/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)*
$ d3antid)
coe5z=-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*(1.d0/2.d0*(1.d0+z)*
$ dpol(rs*(2/(1.d0+z))**(1.d0/3.d0))-1.d0/2.d0*(1.d0-z)*
$ dpol(rs*(2/(1.d0-z))**(1.d0/3.d0))-rs/6.d0*
$ ((1.d0+z)/2.d0)**(2.d0/3.d0)*dpold(rs*(2/(1.d0+z))
$ **(1.d0/3.d0))+rs/6.d0*((1.d0-z)/2.d0)**(2.d0/3.d0)
$ *dpold(rs*(2/(1.d0-z))**(1.d0/3.d0))-2.d0*z*d3anti)
endif
call ecPW(rs,z,ec,ecd,ecz)
a1=4.d0*b0**6*coe3+b0**8*coe5
a1rs=24.d0*adib*b0**5*coe3+4.d0*b0**6*coe3rs+8.d0*adib*b0**7*
$ coe5+b0**8*coe5rs
a1z=4.d0*b0**6*coe3z+b0**8*coe5z
a2=4.d0*b0**6*coe2+b0**8*coe4+6.d0*b0**4*ec
a2rs=24.d0*adib*b0**5*coe2+4.d0*b0**6*coe2rs+8.d0*adib*b0**7*
$ coe4+b0**8*coe4rs+24.d0*adib*b0**3*ec+6.d0*b0**4*ecd
a2z=4.d0*b0**6*coe2z+b0**8*coe4z+6.d0*b0**4*ecz
a3=b0**8*coe3
a3rs=8.d0*adib*b0**7*coe3+b0**8*coe3rs
a3z=b0**8*coe3z
a4=b0**6*(b0**2*coe2+4.d0*ec)
a4rs=8.d0*adib*b0**7*coe2+b0**8*coe2rs+24.d0*adib*b0**5*ec+
$ 4.d0*b0**6*ecd
a4z=b0**6*(b0**2*coe2z+4.d0*ecz)
a5=b0**8*ec
a5rs=8.d0*adib*b0**7*ec+b0**8*ecd
a5z=b0**8*ecz
x=mu*sqrt(rs)/phi
eclr=(phi**3*Qrpa(x)+a1*mu**3+a2*mu**4+a3*mu**5+
$ a4*mu**6+a5*mu**8)/((1.d0+b0**2*mu**2)**4)
eclrrs=-4.d0/(1.d0+b0**2*mu**2)*2.d0*adib*b0*mu**2*eclr+
$ 1.d0/((1.d0+b0**2*mu**2)**4)*(phi**2*mu/(2.d0*sqrt(rs))
$ *Qrpad(x)+
$ a1rs*mu**3+a2rs*mu**4+a3rs*mu**5+a4rs*mu**6+a5rs*mu**8)
if(z.eq.1.d0) then
vclrup=eclr-rs/3.d0*eclrrs
vclrdown=0.d0
elseif(z.eq.-1.d0) then
vclrup=0.d0
vclrdown=eclr-rs/3.d0*eclrrs
else
eclrz=(phi**2*((1.d0+z)**(-1.d0/3.d0)-(1.d0-z)**(-1.d0/3.d0))
$ *Qrpa(x)-phi*Qrpad(x)*mu*sqrt(rs)*((1.d0+z)**(-1.d0/3.d0)
$ -(1.d0-z)**(-1.d0/3.d0))/3.d0+
$ a1z*mu**3+a2z*mu**4+a3z*mu**5+
$ a4z*mu**6+a5z*mu**8)/((1.d0+b0**2*mu**2)**4)
vclrup=eclr-rs/3.d0*eclrrs-(z-1.d0)*eclrz
vclrdown=eclr-rs/3.d0*eclrrs-(z+1.d0)*eclrz
endif
return
end
double precision function g0(x)
ccc on-top pair-distribution function
ccc Gori-Giorgi and Perdew, PRB 64, 155102 (2001)
ccc x -> rs
implicit none
double precision C0f,D0f,E0f,F0f,x
C0f = 0.0819306d0
D0f = 0.752411d0
E0f = -0.0127713d0
F0f = 0.00185898d0
g0=(1.d0-(0.7317d0-D0f)*x+C0f*x**2+E0f*x**3+
$ F0f*x**4)*exp(-abs(D0f)*x)/2.d0
return
end
double precision function g0d(rs)
ccc derivative of on-top pair-distribution function
ccc Gori-Giorgi and Perdew, PRB 64, 155102 (2001)
implicit none
double precision Bg0,Cg0,Dg0,Eg0,Fg0,rs
Cg0 = 0.0819306d0
Fg0 = 0.752411d0
Dg0 = -0.0127713d0
Eg0 = 0.00185898d0
Bg0 =0.7317d0-Fg0
g0d=(-Bg0+2*Cg0*rs+3*Dg0*rs**2+4*Eg0*rs**3)/2.d0*exp(-Fg0*rs)
- - (Fg0*(1 - Bg0*rs + Cg0*rs**2 + Dg0*rs**3 + Eg0*rs**4))/
- 2.d0*exp(-Fg0*rs)
return
end
double precision function dpol(rs)
implicit none
double precision cf,pi,rs,p2p,p3p
pi=dacos(-1.d0)
cf=(9.d0*pi/4.d0)**(1.d0/3.d0)
p2p = 0.04d0
p3p = 0.4319d0
dpol=2.d0**(5.d0/3.d0)/5.d0*cf**2/rs**2*(1.d0+(p3p-0.454555d0)*rs)
$ /(1.d0+p3p*rs+p2p*rs**2)
return
end
double precision function dpold(rs)
implicit none
double precision cf,pi,rs,p2p,p3p
pi=dacos(-1.d0)
cf=(9.d0*pi/4.d0)**(1.d0/3.d0)
p2p = 0.04d0
p3p = 0.4319d0
dpold=2.d0**(5.d0/3.d0)/5.d0*cf**2*
- (-2. + (0.454555 - 4.*p3p)*rs +
- (-4.*p2p +
- (0.90911 - 2.*p3p)*p3p)*rs**2
- + p2p*(1.363665 - 3.*p3p)*
- rs**3)/
- (rs**3*(1. + p3p*rs + p2p*rs**2)**2)
return
end
double precision function Qrpa(x)
implicit none
double precision pi,a2,b2,c2,d2,x,Acoul
pi=dacos(-1.d0)
Acoul=2.d0*(log(2.d0)-1.d0)/pi**2
a2 = 5.84605d0
c2 = 3.91744d0
d2 = 3.44851d0
b2=d2-3.d0/(2.d0*pi*Acoul)*(4.d0/(9.d0*pi))**(1.d0/3.d0)
Qrpa=Acoul*log((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2))
return
end
double precision function Qrpad(x)
implicit none
double precision pi,a2,b2,c2,d2,x,Acoul
pi=dacos(-1.d0)
Acoul=2.d0*(log(2.d0)-1.d0)/pi**2
a2 = 5.84605d0
c2 = 3.91744d0
d2 = 3.44851d0
b2=d2-3.d0/(2.d0*pi*Acoul)*(4.d0/(9.d0*pi))**(1.d0/3.d0)
Qrpad=Acoul*((x*(b2*(2.d0 + a2*x) +
- c2*x*(3.d0 + 2.d0*a2*x) +
- d2*(-2.d0 - a2*x + c2*x**3)))/
- ((1.d0 + a2*x + d2*x**2)*
- (1.d0 + a2*x + b2*x**2 + c2*x**3)))
return
end
subroutine exchangelr(rs,z,mu,exlr)
ccc Hartree atomic units used
ccc for given density parameter rs, spin polarization z
ccc and cutoff mu it gives the exchange energy of the LR gas
ccc => exlr
implicit none
double precision rs,z,mu,exlr
double precision pi,alpha,fx,y
pi=dacos(-1.d0)
alpha=(4.d0/9.d0/pi)**(1.d0/3.d0)
if(abs(z).eq.1.d0) then
y=mu*alpha*rs/2.d0/2.d0**(1.d0/3.d0)
fx=-((-3*y + 4*y**3 +(2*y - 4*y**3)*exp(-1./(4.*y**2)) +
$ sqrt(pi)*erf(1/(2.*y)))/pi)
exlr=mu*fx
else
y=mu*alpha*rs/2.d0/(1.+z)**(1.d0/3.d0)
fx=-((-3*y + 4*y**3 +(2*y - 4*y**3)*exp(-1./(4.*y**2)) +
$ sqrt(pi)*erf(1/(2.*y)))/pi)
exlr=(1.d0+z)*mu*fx/2.d0
y=mu*alpha*rs/2.d0/(1.-z)**(1.d0/3.d0)
fx=-((-3*y + 4*y**3 +(2*y - 4*y**3)*exp(-1./(4.*y**2)) +
$ sqrt(pi)*erf(1/(2.*y)))/pi)
exlr=exlr+(1.d0-z)*mu*fx/2.d0
endif
return
end
subroutine vexchangelr(rs,z,mu,vxlrup,vxlrdown)
ccc Hartree atomic units used
ccc for given density parameter rs, spin polarization z
ccc and cutoff mu it gives the exchange LSD potential for LR interaction
ccc => vxlrup (spin-up electrons), vxlrdown (spin-down electrons)
implicit none
double precision rs,z,mu,vxlrup,vxlrdown
double precision pi,alpha,fx,fx1,y,exlr,derrs,derz
pi=dacos(-1.d0)
alpha=(4.d0/9.d0/pi)**(1.d0/3.d0)
if(z.eq.1.d0) then
vxlrup=(rs*alpha*mu**2)/
- (2**(1.d0/3.d0)*pi) - (rs*alpha*mu**2)/(2**(1.d0/3.d0)*pi)*
- exp(-2**(2.d0/3.d0)/(rs**2*alpha**2*mu**2)) -
- (mu*erf(2**(1.d0/3.d0)/(rs*alpha*mu)))/sqrt(Pi)
vxlrdown=0.d0
elseif(z.eq.-1.d0) then
vxlrdown=(rs*alpha*mu**2)/
- (2**(1.d0/3.d0)*pi) - (rs*alpha*mu**2)/(2**(1.d0/3.d0)*pi)*
- exp(-2**(2.d0/3.d0)/(rs**2*alpha**2*mu**2)) -
- (mu*erf(2**(1.d0/3.d0)/(rs*alpha*mu)))/sqrt(Pi)
vxlrup=0.d0
else
y=mu*alpha*rs/2.d0/(1.+z)**(1.d0/3.d0)
fx=-((-3*y + 4*y**3 +(2*y - 4*y**3)*exp(-1./(4.*y**2)) +
$ sqrt(pi)*erf(1/(2.*y)))/pi)
fx1=(3.d0*(1 + (-4.d0 + 4.d0*exp(-1.d0/(4.d0*y**2)))*y**2))/pi
derrs=1.d0/4.d0*(1.d0+z)**(2.d0/3.d0)*mu**2*alpha*fx1
derz=1.d0/2.d0*mu*fx-1.d0/6.d0*fx1*mu*y
vxlrup=rs/3.d0*derrs+(z-1.d0)*derz
vxlrdown=rs/3.d0*derrs+(z+1.d0)*derz
y=mu*alpha*rs/2.d0/(1.-z)**(1.d0/3.d0)
fx=-((-3*y + 4*y**3 +(2*y - 4*y**3)*exp(-1./(4.*y**2)) +
$ sqrt(pi)*erf(1/(2.*y)))/pi)
fx1=(3.d0*(1 + (-4.d0 + 4.d0*exp(-1.d0/(4.d0*y**2)))*y**2))/pi
derrs=1.d0/4.d0*(1.d0-z)**(2.d0/3.d0)*mu**2*alpha*fx1
derz=-1.d0/2.d0*mu*fx+1.d0/6.d0*fx1*mu*y
vxlrup=vxlrup+rs/3.d0*derrs+(z-1.d0)*derz
vxlrdown=vxlrdown+rs/3.d0*derrs+(z+1.d0)*derz
call exchangelr(rs,z,mu,exlr)
vxlrup=exlr-vxlrup
vxlrdown=exlr-vxlrdown
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c correlation energy and its derivative w.r.t. rs and z at mu=infinity
c Perdew & Wang PRB 45, 13244 (1992)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine ecPW(x,y,ec,ecd,ecz)
c in Hartree; ec=ec(rs,zeta)
c x -> rs; y -> zeta
ccc ecd is d/drs ec
ccc ecz is d/dz ec
implicit none
double precision pi,f02,ff,x,y,ec,ecd,ec0,ec0d,ec1,ec1d,
$ aaa,G,Gd,alfac,alfacd,ecz
pi=dacos(-1.d0)
f02=4.d0/(9.d0*(2.d0**(1.d0/3.d0)-1.d0))
ff=((1.d0+y)**(4.d0/3.d0)+(1.d0-y)**(4.d0/3.d0)-
$ 2.d0)/(2.d0**(4.d0/3.d0)-2.d0)
aaa=(1.d0-log(2.d0))/pi**2
call GPW(x,aaa,0.21370d0,7.5957d0,3.5876d0,
$ 1.6382d0,0.49294d0,G,Gd)
ec0=G
ec0d=Gd
aaa=aaa/2.d0
call GPW(x,aaa,0.20548d0,14.1189d0,6.1977d0,
$ 3.3662d0,0.62517d0,G,Gd)
ec1=G
ec1d=Gd
call GPW(x,0.016887d0,0.11125d0,10.357d0,3.6231d0,
$ 0.88026d0,0.49671d0,G,Gd)
alfac=-G
alfacd=-Gd
ec=ec0+alfac*ff/f02*(1.d0-y**4)+(ec1-ec0)*ff*y**4
ecd=ec0d+alfacd*ff/f02*(1.d0-y**4)+(ec1d-ec0d)*
$ ff*y**4
ecz=alfac*(-4.d0*y**3)*ff/f02+alfac*(1.d0-y**4)/f02*
$ 4.d0/3.d0*((1.d0+y)**(1.d0/3.d0)-(1.d0-y)**(1.d0/3.d0))/
$ (2.d0**(4.d0/3.d0)-2.d0)+(ec1-ec0)*(4.d0*y**3*ff+
$ 4.d0/3.d0*((1.d0+y)**(1.d0/3.d0)-(1.d0-y)**(1.d0/3.d0))/
$ (2.d0**(4.d0/3.d0)-2.d0)*y**4)
return
end
subroutine GPW(x,Ac,alfa1,beta1,beta2,beta3,beta4,G,Gd)
ccc Gd is d/drs G
implicit none
double precision G,Gd,Ac,alfa1,beta1,beta2,beta3,beta4,x
G=-2.d0*Ac*(1.d0+alfa1*x)*dlog(1.d0+1.d0/(2.d0*
$ Ac*(beta1*x**0.5d0+
$ beta2*x+beta3*x**1.5d0+beta4*x**2)))
Gd=(1.d0+alfa1*x)*(beta2+beta1/(2.d0*sqrt(x))+3.d0*beta3*
$ sqrt(x)/2.d0+2.d0*beta4*x)/((beta1*sqrt(x)+beta2*x+
$ beta3*x**(3.d0/2.d0)+beta4*x**2)**2*(1.d0+1.d0/
$ (2.d0*Ac*(beta1*sqrt(x)+beta2*x+beta3*x**(3.d0/2.d0)+
$ beta4*x**2))))-2.d0*Ac*alfa1*dlog(1.d0+1.d0/(2.d0*Ac*
$ (beta1*sqrt(x)+beta2*x+beta3*x**(3.d0/2.d0)+
$ beta4*x**2)))
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

View File

@ -1,10 +0,0 @@
default:
ninja
make -C ..
clean:
ninja -t clean
debug:
ninja -t clean
make -C .. debug

View File

@ -1,90 +0,0 @@
subroutine basis_correction(nBas,nO,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
ERI,e,c,P,eG0W0)
! Compute the basis set incompleteness error
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: nShell
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: CenterShell(maxShell,ncart),DShell(maxShell,maxK),ExpShell(maxShell,maxK)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: eG0W0(nBas)
! Local variables
integer :: SGn
integer :: nRad
integer :: nAng
integer :: nGrid
double precision,allocatable :: root(:,:)
double precision,allocatable :: weight(:)
double precision,allocatable :: AO(:,:)
double precision,allocatable :: dAO(:,:,:)
double precision,allocatable :: MO(:,:)
double precision,allocatable :: dMO(:,:,:)
double precision,allocatable :: rho(:)
double precision,allocatable :: f(:)
double precision,allocatable :: mu(:)
! Output variables
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Basis set incompleteness correction |'
write(*,*)'************************************************'
write(*,*)
! Set quadrature grid
SGn = 1
call read_grid(SGn,nRad,nAng,nGrid)
! Memory allocation
allocate(root(ncart,nGrid),weight(nGrid))
allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid),MO(nBas,nGrid),dMO(ncart,nBas,nGrid))
allocate(rho(nGrid),f(nGrid),mu(nGrid))
call quadrature_grid(nRad,nAng,nGrid,root,weight)
! Calculate AO values at grid points
call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO)
! Calculate MO values at grid points
call MO_values_grid(nBas,nGrid,c,AO,dAO,MO,dMO)
! Compute one-electron density at grid points
call density(nGrid,nBas,P,AO,rho)
! Compute range-sepration function
call f_grid(nBas,nO,nGrid,MO,ERI,f)
call mu_grid(nGrid,rho,f,mu)
! Compute energy correction
call ec_srlda(nGrid,weight,rho,mu)
! Compute orbital corrections
call fc_srlda(nBas,nGrid,weight,MO,rho,mu,eG0W0)
end subroutine basis_correction

File diff suppressed because it is too large Load Diff

View File

@ -1,52 +0,0 @@
subroutine ec_srlda(nGrid,weight,rho,mu)
! Compute sr-lda ec
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: mu(nGrid)
! Local variables
integer :: iG
double precision :: r,ra,rb
double precision :: rs
double precision :: ec_lda,ecmd_lda
double precision :: ec,ecmd,vcup,vcdw
! Initialization
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)
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(*,*)
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

View File

@ -1,200 +0,0 @@
!****************************************************************************
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

View File

@ -1,44 +0,0 @@
subroutine f_grid(nBas,nO,nGrid,MO,ERI,f)
! Compute f
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: nGrid
double precision,intent(in) :: MO(nBas,nGrid)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: p,q
integer :: i,j
integer :: iG
! Output variables
double precision,intent(out) :: f(nGrid)
! Initialization
f(:) = 0d0
do p=1,nBas
do q=1,nBas
do i=1,nO
do j=1,nO
do iG=1,ngrid
f(iG) = f(iG) + MO(p,iG)*MO(q,iG)*ERI(p,q,i,j)*MO(i,iG)*MO(j,iG)
end do
end do
end do
end do
end do
end subroutine f_grid

View File

@ -1,88 +0,0 @@
subroutine fc_srlda(nBas,nGrid,weight,MO,rho,mu,eG0W0)
! Compute sr-lda ec
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: MO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: mu(nGrid)
double precision,intent(in) :: eG0W0(nBas)
! Local variables
integer :: iG,p
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
allocate(de(nBas))
! Initialization
de(:) = 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)
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
de(p)= de(p) + weight(iG)*vcup*MO(p,iG)**2
end do
end if
end do
print*, 'Eigenvalues correction from srDFT (in eV)'
call matout(nBas,1,de(:)*HaToeV)
print*, 'Corrected G0W0 eigenvalues (in eV)'
call matout(nBas,1,(eG0W0(:) + de(:))*HaToeV)
end subroutine fc_srlda

View File

@ -1,46 +0,0 @@
subroutine mu_grid(nGrid,rho,f,mu)
! Compute mu
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: f(nGrid)
! Local variables
integer :: iG
double precision,parameter :: thres = 1d-10
double precision :: n2
! Output variables
double precision,intent(out) :: mu(nGrid)
! Initialization
mu(:) = 0d0
do iG=1,ngrid
n2 = 0.25d0*rho(iG)**2
if(abs(n2) > thres) then
mu(iG) = f(iG)/n2
else
mu(iG) = 1d0/thres
end if
end do
mu(:) = 0.5d0*sqrt(pi)*mu(:)
end subroutine mu_grid

View File

@ -1,77 +0,0 @@
subroutine quadrature_grid(nRad,nAng,nGrid,root,weight)
! Build roots and weights of quadrature grid
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nRad,nAng,nGrid
! Local variables
integer :: i,j,k
double precision :: scale
double precision,allocatable :: Radius(:)
double precision,allocatable :: RadWeight(:)
double precision,allocatable :: XYZ(:,:)
double precision,allocatable :: XYZWeight(:)
! Output variables
double precision,intent(out) :: root(3,nGrid)
double precision,intent(out) :: weight(nGrid)
! Memory allocation
allocate(Radius(nRad),RadWeight(nRad),XYZ(3,nAng),XYZWeight(nAng))
! Find the radial grid
scale = 1d0
call EulMac(Radius,RadWeight,nRad,scale)
write(*,20)
write(*,30)
write(*,20)
do i = 1,nRad
write(*,40) i,Radius(i),RadWeight(i)
end do
write(*,20)
write(*,*)
! Find the angular grid
call Lebdev(XYZ,XYZWeight,nAng)
write(*,20)
write(*,50)
write(*,20)
do j = 1,nAng
write(*,60) j,(XYZ(k,j),k=1,3), XYZWeight(j)
end do
write(*,20)
! Form the roots and weights
k = 0
do i=1,nRad
do j=1,nAng
k = k + 1
root(:,k) = Radius(i)*XYZ(:,j)
weight(k) = RadWeight(i)*XYZWeight(j)
enddo
enddo
! Compute values of the basis functions (and the its gradient if required) at each grid point
20 format(T2,58('-'))
30 format(T20,'Radial Quadrature',/, &
T6,'I',T26,'Radius',T50,'Weight')
40 format(T3,I4,T18,F17.10,T35,F25.10)
50 format(T20,'Angular Quadrature',/, &
T6,'I',T19,'X',T29,'Y',T39,'Z',T54,'Weight')
60 format(T3,I4,T13,3F10.5,T50,F10.5)
end subroutine quadrature_grid

View File

@ -1,169 +0,0 @@
subroutine read_F12_integrals(nBas,S,C,F,Y,FC)
! Read one- and two-electron integrals from files
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: S(nBas,nBas)
! Local variables
logical :: debug
integer :: mu,nu,la,si,ka,ta
double precision :: ERI,F12,Yuk,F13C12,ExpS
! Output variables
double precision,intent(out) :: C(nBas,nBas,nBas,nBas)
double precision,intent(out) :: F(nBas,nBas,nBas,nBas)
double precision,intent(out) :: Y(nBas,nBas,nBas,nBas)
double precision,intent(out) :: FC(nBas,nBas,nBas,nBas,nBas,nBas)
debug = .false.
! Open file with integrals
open(unit=21,file='int/ERI.dat')
open(unit=22,file='int/F12.dat')
open(unit=23,file='int/Yuk.dat')
open(unit=31,file='int/3eInt_Type1.dat')
! Read 1/r12 integrals
C = 0d0
do
read(21,*,end=21) mu,nu,la,si,ERI
! <12|34>
C(mu,nu,la,si) = ERI
! <32|14>
C(la,nu,mu,si) = ERI
! <14|32>
C(mu,si,la,nu) = ERI
! <34|12>
C(la,si,mu,nu) = ERI
! <41|23>
C(si,mu,nu,la) = ERI
! <23|41>
C(nu,la,si,mu) = ERI
! <21|43>
C(nu,mu,si,la) = ERI
! <43|21>
C(si,la,nu,mu) = ERI
enddo
21 close(unit=21)
! Read f12 integrals
F = 0d0
do
read(22,*,end=22) mu,nu,la,si,F12
! <12|34>
F(mu,nu,la,si) = F12
! <32|14>
F(la,nu,mu,si) = F12
! <14|32>
F(mu,si,la,nu) = F12
! <34|12>
F(la,si,mu,nu) = F12
! <41|23>
F(si,mu,nu,la) = F12
! <23|41>
F(nu,la,si,mu) = F12
! <21|43>
F(nu,mu,si,la) = F12
! <43|21>
F(si,la,nu,mu) = F12
enddo
22 close(unit=22)
! Read f12/r12 integrals
Y = 0d0
do
read(23,*,end=23) mu,nu,la,si,Yuk
! <12|34>
Y(mu,nu,la,si) = Yuk
! <32|14>
Y(la,nu,mu,si) = Yuk
! <14|32>
Y(mu,si,la,nu) = Yuk
! <34|12>
Y(la,si,mu,nu) = Yuk
! <41|23>
Y(si,mu,nu,la) = Yuk
! <23|41>
Y(nu,la,si,mu) = Yuk
! <21|43>
Y(nu,mu,si,la) = Yuk
! <43|21>
Y(si,la,nu,mu) = Yuk
enddo
23 close(unit=23)
! Read f13/r12 integrals
FC = 0d0
do
read(31,*,end=31) mu,nu,la,si,ka,ta,F13C12
FC(mu,nu,la,si,ka,ta) = F13C12
enddo
31 close(unit=31)
! Print results
if(debug) then
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Electron repulsion integrals'
write(*,'(A28)') '----------------------'
do la=1,nBas
do si=1,nBas
call matout(nBas,nBas,C(1,1,la,si))
enddo
enddo
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'F12 integrals'
write(*,'(A28)') '----------------------'
do la=1,nBas
do si=1,nBas
call matout(nBas,nBas,F(1,1,la,si))
enddo
enddo
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Yukawa integrals'
write(*,'(A28)') '----------------------'
do la=1,nBas
do si=1,nBas
call matout(nBas,nBas,Y(1,1,la,si))
enddo
enddo
write(*,*)
endif
! Read exponent of Slater geminal
open(unit=4,file='input/geminal')
read(4,*) ExpS
close(unit=4)
! Transform two-electron integrals
! do mu=1,nBas
! do nu=1,nBas
! do la=1,nBas
! do si=1,nBas
! F(mu,nu,la,si) = (S(mu,la)*S(nu,si) - F(mu,nu,la,si))/ExpS
! Y(mu,nu,la,si) = (C(mu,nu,la,si) - Y(mu,nu,la,si))/ExpS
! enddo
! enddo
! enddo
! enddo
end subroutine read_F12_integrals

View File

@ -1,47 +0,0 @@
subroutine read_grid(SGn,nRad,nAng,nGrid)
! Read grid type
implicit none
! Input variables
integer,intent(in) :: SGn
! Output variables
integer,intent(out) :: nRad
integer,intent(out) :: nAng
integer,intent(out) :: nGrid
write(*,*)'----------------------------------------------------------'
write(*,'(A22,I1)')' Quadrature grid: SG-',SGn
write(*,*)'----------------------------------------------------------'
select case (SGn)
case(0)
nRad = 23
nAng = 170
case(1)
nRad = 50
nAng = 194
case(2)
nRad = 75
nAng = 302
case(3)
nRad = 99
nAng = 590
case default
call print_warning('!!! Quadrature grid not available !!!')
stop
end select
nGrid = nRad*nAng
end subroutine read_grid

Binary file not shown.

View File

@ -1,108 +0,0 @@
subroutine AO_values(doDrift,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,r,AO,dAO)
! Compute values of the AOs and their derivatives (if required)
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: doDrift
integer,intent(in) :: nBas,nShell,nWalk
double precision,intent(in) :: CenterShell(maxShell,3)
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
double precision,intent(in) :: r(nWalk,2,3)
! Local variables
integer :: atot,nShellFunction,a(3)
integer,allocatable :: ShellFunction(:,:)
double precision :: rASq,xA,yA,zA,norm_coeff,prim
integer :: iSh,iShF,iK,iW,iEl,iBas,ixyz
! Output variables
double precision,intent(out) :: AO(nWalk,2,nBas),dAO(nWalk,2,3,nBas)
! Initialization
AO = 0d0
if(doDrift) dAO = 0d0
iBas = 0
!------------------------------------------------------------------------
! Loops over shells
!------------------------------------------------------------------------
do iSh=1,nShell
atot = TotAngMomShell(iSh)
nShellFunction = (atot*atot + 3*atot + 2)/2
allocate(ShellFunction(1:nShellFunction,1:3))
call generate_shell(atot,nShellFunction,ShellFunction)
do iShF=1,nShellFunction
iBas = iBas + 1
a(1) = ShellFunction(iShF,1)
a(2) = ShellFunction(iShF,2)
a(3) = ShellFunction(iShF,3)
do iW=1,nWalk
do iEl=1,2
xA = r(iW,iEl,1) - CenterShell(iSh,1)
yA = r(iW,iEl,2) - CenterShell(iSh,2)
zA = r(iW,iEl,3) - CenterShell(iSh,3)
! Calculate distance for exponential
rASq = xA**2 + yA**2 + zA**2
!------------------------------------------------------------------------
! Loops over contraction degrees
!-------------------------------------------------------------------------
do iK=1,KShell(iSh)
! Calculate the exponential part
prim = DShell(iSh,iK)*norm_coeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq)
AO(iW,iEl,iBas) = AO(iW,iEl,iBas) + prim
if(doDrift) then
prim = -2d0*ExpShell(iSh,iK)*prim
do ixyz=1,3
dAO(iW,iEl,ixyz,iBas) = dAO(iW,iEl,ixyz,iBas) + prim
enddo
endif
enddo
if(doDrift) then
dAO(iW,iEl,1,iBas) = xA**(a(1)+1)*yA**a(2)*zA**a(3)*dAO(iW,iEl,1,iBas)
if(a(1) > 0) dAO(iW,iEl,1,iBas) = dAO(iW,iEl,1,iBas) + dble(a(1))*xA**(a(1)-1)*yA**a(2)*zA**a(3)*AO(iW,iEl,iBas)
dAO(iW,iEl,2,iBas) = xA**a(1)*yA**(a(2)+1)*zA**a(3)*dAO(iW,iEl,2,iBas)
if(a(2) > 0) dAO(iW,iEl,2,iBas) = dAO(iW,iEl,2,iBas) + dble(a(2))*xA**a(1)*yA**(a(2)-1)*zA**a(3)*AO(iW,iEl,iBas)
dAO(iW,iEl,3,iBas) = xA**a(1)*yA**a(2)*zA**(a(3)+1)*dAO(iW,iEl,3,iBas)
if(a(3) > 0) dAO(iW,iEl,3,iBas) = dAO(iW,iEl,3,iBas) + dble(a(3))*xA**a(1)*yA**a(2)*zA**(a(3)-1)*AO(iW,iEl,iBas)
endif
! Calculate polynmial part
AO(iW,iEl,iBas) = xA**a(1)*yA**a(2)*zA**a(3)*AO(iW,iEl,iBas)
enddo
enddo
enddo
deallocate(ShellFunction)
enddo
!------------------------------------------------------------------------
! End loops over shells
!------------------------------------------------------------------------
end subroutine AO_values

View File

@ -1,101 +0,0 @@
subroutine AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
nGrid,root,AO,dAO)
! Compute values of the AOs and their derivatives with respect to the cartesian coordinates
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nShell
double precision,intent(in) :: CenterShell(maxShell,3)
integer,intent(in) :: TotAngMomShell(maxShell)
integer,intent(in) :: KShell(maxShell)
double precision,intent(in) :: DShell(maxShell,maxK)
double precision,intent(in) :: ExpShell(maxShell,maxK)
double precision,intent(in) :: root(3,nGrid)
integer,intent(in) :: nGrid
! Local variables
integer :: atot,nShellFunction,a(3)
integer,allocatable :: ShellFunction(:,:)
double precision :: rASq,xA,yA,zA,norm_coeff,prim
integer :: iSh,iShF,iK,iG,iBas
! Output variables
double precision,intent(out) :: AO(nBas,nGrid)
double precision,intent(out) :: dAO(3,nBas,nGrid)
! Initialization
iBas = 0
AO(:,:) = 0d0
dAO(:,:,:) = 0d0
!------------------------------------------------------------------------
! Loops over shells
!------------------------------------------------------------------------
do iSh=1,nShell
atot = TotAngMomShell(iSh)
nShellFunction = (atot*atot + 3*atot + 2)/2
allocate(ShellFunction(1:nShellFunction,1:3))
call generate_shell(atot,nShellFunction,ShellFunction)
do iShF=1,nShellFunction
iBas = iBas + 1
a(:) = ShellFunction(iShF,:)
do iG=1,nGrid
xA = root(1,iG) - CenterShell(iSh,1)
yA = root(2,iG) - CenterShell(iSh,2)
zA = root(3,iG) - CenterShell(iSh,3)
! Calculate distance for exponential
rASq = xA**2 + yA**2 + zA**2
!------------------------------------------------------------------------
! Loops over contraction degrees
!-------------------------------------------------------------------------
do iK=1,KShell(iSh)
! Calculate the exponential part
prim = DShell(iSh,iK)*norm_coeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq)
AO(iBas,iG) = AO(iBas,iG) + prim
prim = -2d0*ExpShell(iSh,iK)*prim
dAO(:,iBas,iG) = dAO(:,iBas,iG) + prim
enddo
dAO(1,iBas,iG) = xA**(a(1)+1)*yA**a(2)*zA**a(3)*dAO(1,iBas,iG)
if(a(1) > 0) dAO(1,iBas,iG) = dAO(1,iBas,iG) + dble(a(1))*xA**(a(1)-1)*yA**a(2)*zA**a(3)*AO(iBas,iG)
dAO(2,iBas,iG) = xA**a(1)*yA**(a(2)+1)*zA**a(3)*dAO(2,iBas,iG)
if(a(2) > 0) dAO(2,iBas,iG) = dAO(2,iBas,iG) + dble(a(2))*xA**a(1)*yA**(a(2)-1)*zA**a(3)*AO(iBas,iG)
dAO(3,iBas,iG) = xA**a(1)*yA**a(2)*zA**(a(3)+1)*dAO(3,iBas,iG)
if(a(3) > 0) dAO(3,iBas,iG) = dAO(3,iBas,iG) + dble(a(3))*xA**a(1)*yA**a(2)*zA**(a(3)-1)*AO(iBas,iG)
! Calculate polynmial part
AO(iBas,iG) = xA**a(1)*yA**a(2)*zA**a(3)*AO(iBas,iG)
enddo
enddo
deallocate(ShellFunction)
enddo
!------------------------------------------------------------------------
! End loops over shells
!------------------------------------------------------------------------
end subroutine AO_values_grid

View File

@ -1,65 +0,0 @@
subroutine Green_function(nBas,nO,nV,nWalk,nWP,cO,cV,eO_Quad,eV_Quad,AO, &
o1MO,o2MO,v1MO,v2MO,o11,o12,o21,o22,v11,v12,v21,v22)
! Calculate the Green functions
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
integer,intent(in) :: nBas,nO,nV,nWalk,nWP
double precision,intent(in) :: AO(nWalk,2,nBas),cO(nBas,nO),cV(nBas,nV), &
eO_Quad(nQuad,nO),eV_Quad(nQuad,nV)
! Local variables
integer :: kW,lW,klW,i,a,q
double precision :: o1MO(nWalk,nO),o2MO(nWalk,nO),v1MO(nWalk,nV),v2MO(nWalk,nV)
! Output variables
double precision,intent(out) :: o11(nQuad,nWP),o12(nQuad,nWP),o21(nQuad,nWP),o22(nQuad,nWP)
double precision,intent(out) :: v11(nQuad,nWP),v12(nQuad,nWP),v21(nQuad,nWP),v22(nQuad,nWP)
! Calculate occupied and virtual MOs
o1MO = matmul(AO(:,1,:),cO)
o2MO = matmul(AO(:,2,:),cO)
v1MO = matmul(AO(:,1,:),cV)
v2MO = matmul(AO(:,2,:),cV)
! Compute occupied Green functions
o11 = 0d0
o12 = 0d0
o21 = 0d0
o22 = 0d0
v11 = 0d0
v12 = 0d0
v21 = 0d0
v22 = 0d0
do q=1,nQuad
klW = 0
do kW=1,nWalk-1
do lW=kW+1,nWalk
klW = klW + 1
do i=1,nO
o11(q,klW) = o11(q,klW) + o1MO(kW,i)*o1MO(lW,i)*eO_Quad(q,i)
o12(q,klW) = o12(q,klW) + o1MO(kW,i)*o2MO(lW,i)*eO_Quad(q,i)
o21(q,klW) = o21(q,klW) + o2MO(kW,i)*o1MO(lW,i)*eO_Quad(q,i)
o22(q,klW) = o22(q,klW) + o2MO(kW,i)*o2MO(lW,i)*eO_Quad(q,i)
enddo
do a=1,nV
v11(q,klW) = v11(q,klW) + v1MO(kW,a)*v1MO(lW,a)*eV_Quad(q,a)
v12(q,klW) = v12(q,klW) + v1MO(kW,a)*v2MO(lW,a)*eV_Quad(q,a)
v21(q,klW) = v21(q,klW) + v2MO(kW,a)*v1MO(lW,a)*eV_Quad(q,a)
v22(q,klW) = v22(q,klW) + v2MO(kW,a)*v2MO(lW,a)*eV_Quad(q,a)
enddo
enddo
enddo
enddo
end subroutine Green_function

View File

@ -1,344 +0,0 @@
subroutine MCMP2(doDrift,nBas,nC,nO,nV,c,e,EcMP2, &
nMC,nEq,nWalk,dt,nPrint, &
nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
Norm, &
EcMCMP2,Err_EcMCMP2,Var_EcMCMP2)
! Perform Monte Carlo MP2 calculation
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: doDrift
integer,intent(in) :: nBas,nC,nO,nV,nMC,nEq,nWalk,nPrint
double precision,intent(inout):: dt
double precision,intent(in) :: EcMP2(3)
double precision,intent(in) :: c(nBas,nBas),e(nBas)
integer,intent(in) :: nShell
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: CenterShell(maxShell,3),DShell(maxShell,maxK),ExpShell(maxShell,maxK)
! Local variables
logical :: AcPh,EqPh,Accept,dump
double precision :: start_Eq,end_Eq,t_Eq,start_Ac,end_Ac,t_Ac
integer :: nWP
double precision :: Norm,NormSq,nData,tau
double precision,allocatable :: chi1(:,:,:),chi2(:,:,:),eta(:)
double precision,allocatable :: cO(:,:),cV(:,:),eO(:),eV(:),P(:,:),eO_Quad(:,:),eV_Quad(:,:)
double precision,allocatable :: r(:,:,:), r12(:), gAO(:,:,:), g(:,:), w(:)
double precision,allocatable :: rp(:,:,:),r12p(:),gAOp(:,:,:), gp(:,:),wp(:)
double precision,allocatable :: o1MO(:,:),o2MO(:,:),v1MO(:,:),v2MO(:,:)
double precision,allocatable :: o11(:,:),o12(:,:),o21(:,:),o22(:,:)
double precision,allocatable :: v11(:,:),v12(:,:),v21(:,:),v22(:,:)
double precision,allocatable :: fd_Quad(:,:),fx_Quad(:,:),fd(:),fx(:),fdx(:)
double precision,allocatable :: dgAO(:,:,:,:),dg(:,:,:),dgAOp(:,:,:,:),dgp(:,:,:)
double precision,allocatable :: F(:,:,:),Fp(:,:,:),T(:),Tp(:)
double precision :: acceptance,D
double precision :: eloc_MP2(3),mean_MP2(3),variance_MP2(3)
integer :: iW,kW,lW,klW,iMC,q
! Output variables
double precision,intent(out) :: EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3)
! Number of distinct walker pairs
nWP = nWalk*(nWalk-1)/2
! Diffusion coefficient
D = 0.5d0
! Do diffusion-drift moves?
if(doDrift) then
write(*,*)
write(*,*) '*** Diffusion-drift algorithm ***'
write(*,*)
else
write(*,*)
write(*,*) '*** Diffusion-only algorithm ***'
write(*,*)
endif
! Print results
dump = .true.
if(dump) open(unit=13,file='results/data')
!------------------------------------------------------------------------
! Memory allocation
!------------------------------------------------------------------------
allocate(cO(nBas,nO),cV(nBas,nV),eO(nO),eV(nV), &
eO_Quad(nQuad,nO),eV_Quad(nQuad,nV), &
P(nBas,nBas),r(nWalk,2,3),rp(nWalk,2,3), &
chi1(nWalk,2,3),chi2(nWalk,2,3),eta(nWalk), &
r12(nWalk),r12p(nWalk),w(nWalk),wp(nWalk), &
g(nWalk,2),gp(nWalk,2),gAO(nWalk,2,nBas),gAOp(nWalk,2,nBas), &
dg(nWalk,2,3),dgp(nWalk,2,3),dgAO(nWalk,2,3,nBas),dgAOp(nWalk,2,3,nBas), &
o1MO(nWalk,nO),v1MO(nWalk,nV),o2MO(nWalk,nO),v2MO(nWalk,nV), &
o11(nQuad,nWP),v11(nQuad,nWP),o12(nQuad,nWP),v12(nQuad,nWP), &
o21(nQuad,nWP),v21(nQuad,nWP),o22(nQuad,nWP),v22(nQuad,nWP), &
fd_Quad(nQuad,nWP),fd(nWP),fx_Quad(nQuad,nWP),fx(nWP),fdx(nWP), &
T(nWalk),Tp(nWalk),F(nWalk,2,3),Fp(nWalk,2,3))
! Split MOs into occupied and virtual sets
eO(1:nO) = e(nC+1:nC+nO)
eV(1:nV) = e(nC+nO+1:nBas)
do q=1,nQuad
tau = 1d0/rQuad(q)
eO_Quad(q,1:nO) = exp(+eO(1:nO)*(tau-1d0))*sqrt(tau)
eV_Quad(q,1:nV) = exp(-eV(1:nV)*(tau-1d0))*sqrt(tau)
enddo
cO(1:nBas,1:nO) = c(1:nBas,nC+1:nC+nO)
cV(1:nBas,1:nV) = c(1:nBas,nC+nO+1:nBas)
! Compute norm of the trial wave function
call norm_trial(nBas,nO,cO,P,Norm,NormSq)
!------------------------------------------------------------------------
! Initialize MC-MP2 calculation
!------------------------------------------------------------------------
! Initialize electron coordinates
call random_number(r)
r = 2d0*r - 1d0
! Compute initial interelectronic distances
call rij(nWalk,r,r12)
! Compute initial AO values and their derivatives (if required)
call AO_values(doDrift,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,r,gAO,dgAO)
! Compute initial weight function
call density(doDrift,nBas,nWalk,P,gAO,dgAO,g,dg)
! Compute initial weights
w(1:nWalk) = g(1:nWalk,1)*g(1:nWalk,2)/r12(1:nWalk)
! Compute initial quantum force
if(doDrift) call drift(nWalk,r,r12,g,dg,F)
! Equilibration or Accumulation?
AcPh = .false.
EqPh = .true.
! Initialization
nData = 0d0
acceptance = 0d0
mean_MP2 = 0d0
variance_MP2 = 0d0
T = 1d0
Tp = 1d0
!------------------------------------------------------------------------
! Start main Monte Carlo loop
!------------------------------------------------------------------------
call cpu_time(start_Eq)
do iMC=1,nEq+nMC
! Timings
if(iMC == nEq + 1) then
AcPh = .true.
EqPh = .false.
write(*,*) 'Time step value at the end of equilibration: dt = ',dt
write(*,*)
call cpu_time(end_Eq)
t_Eq = end_Eq - start_Eq
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for equilibration = ',t_Eq,' seconds'
write(*,*)
call cpu_time(start_Ac)
endif
! Optimize time step to reach 50% acceptance
if(EqPh .and. mod(iMC,100) == 0) call optimize_timestep(nWalk,iMC,acceptance,dt)
! Move electrons
call random_number(chi1)
call random_number(chi2)
! Diffusion
rp(:,:,:) = r(:,:,:) + sqrt(2d0*D*dt)*sqrt(-2d0*log(chi1(:,:,:)))*cos(2d0*pi*chi2(:,:,:))
! Drift
if(doDrift) rp(:,:,:) = rp(:,:,:) + D*dt*F(:,:,:)
! Compute new interelectronic distances
call rij(nWalk,rp,r12p)
! Compute new AO values and their derivatives (if required)
call AO_values(doDrift,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,rp,gAOp,dgAOp)
call Density(doDrift,nBas,nWalk,P,gAOp,dgAOp,gp,dgp)
! Compute new weights
wp(1:nWalk) = gp(1:nWalk,1)*gp(1:nWalk,2)/r12p(1:nWalk)
! Compute new quantum force and transition probability
if(doDrift) then
call Drift(nWalk,rp,r12p,gp,dgp,Fp)
call transition_probability(nWalk,dt,D,r,rp,F,Fp,T,Tp)
endif
! Move for walkers
call random_number(eta)
do iW=1,nWalk
Accept = (wp(iW)*Tp(iW))/(w(iW)*T(iW)) > eta(iW)
if(Accept) then
acceptance = acceptance + 1d0
r(iW,1:2,1:3) = rp(iW,1:2,1:3)
gAO(iW,1:2,1:nBas) = gAOp(iW,1:2,1:nBas)
r12(iW) = r12p(iW)
w(iW) = wp(iW)
if(doDrift) F(iW,1:2,1:3) = Fp(iW,1:2,1:3)
endif
enddo
! Accumulation phase
if(AcPh) then
nData = nData + 1d0
! Calculate Green functions
call Green_function(nBas,nO,nV,nWalk,nWP,cO,cV,eO_Quad,eV_Quad,gAO, &
o1MO,o2MO,v1MO,v2MO,o11,o12,o21,o22,v11,v12,v21,v22)
! Compute local energy
fd_Quad = o11*o22*v11*v22 + o12*o21*v12*v21
fx_Quad = o11*o22*v12*v21 + o12*o21*v11*v22
fd = matmul(wQuad,fd_Quad)
fx = matmul(wQuad,fx_Quad)
eloc_MP2 = 0d0
klW = 0
do kW=1,nWalk-1
do lW=kW+1,nWalk
klW = klW + 1
eloc_MP2(2) = eloc_MP2(2) + fd(klW)/(r12(kW)*r12(lW)*w(kW)*w(lW))
eloc_MP2(3) = eloc_MP2(3) + fx(klW)/(r12(kW)*r12(lW)*w(kW)*w(lW))
enddo
enddo
eloc_MP2(2) = -2d0*eloc_MP2(2)/dble(2*nWP)
eloc_MP2(3) = eloc_MP2(3)/dble(2*nWP)
fdx = -2d0*fd + fx
eloc_MP2(1) = eloc_MP2(2) + eloc_MP2(3)
! Accumulate results
mean_MP2 = mean_MP2 + eloc_MP2
variance_MP2 = variance_MP2 + eloc_MP2*eloc_MP2
! Print results
if(mod(iMC,nPrint) == 0) then
call compute_error(nData,mean_MP2,variance_MP2,Err_EcMCMP2)
EcMCMP2 = mean_MP2/nData
Var_EcMCMP2 = variance_MP2/nData
EcMCMP2 = Norm*EcMCMP2
Var_EcMCMP2 = Norm*Var_EcMCMP2
Err_EcMCMP2 = Norm*Err_EcMCMP2
write(*,*)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,I15)') 'Number of data points ','|',int(nData)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10I15)') 'acceptance ','|',int(100*acceptance/dble(nWalk*iMC))
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'MP2 correlation energy Total ','|',EcMCMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',EcMCMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',EcMCMP2(3)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'Statistical error Total ','|',Err_EcMCMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',Err_EcMCMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',Err_EcMCMP2(3)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'Variance Total ','|',Var_EcMCMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',Var_EcMCMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',Var_EcMCMP2(3)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'Dev. wrt deterministic Total ','|',EcMCMP2(1) - EcMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',EcMCMP2(2) - EcMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',EcMCMP2(3) - EcMP2(3)
write(*,*)'-------------------------------------------------------'
if(dump) write(13,*) int(nData),EcMCMP2(1),Err_EcMCMP2(1)
endif
endif
!------------------------------------------------------------------------
! End main Monte Carlo loop
!------------------------------------------------------------------------
enddo
! Timing
call cpu_time(end_Ac)
t_Ac = end_Ac - start_Ac
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for accumulation = ',t_Ac,' seconds'
write(*,*)
! Close files
if(dump) close(unit=13)
end subroutine MCMP2

View File

@ -1,44 +0,0 @@
subroutine MO_values_grid(nBas,nGrid,c,AO,dAO,MO,dMO)
! Compute values of the MOs and their derivatives with respect to the cartesian coordinates
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nGrid
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
! Local variables
integer :: p,mu,iG
! Output variables
double precision,intent(out) :: MO(nBas,nGrid)
double precision,intent(out) :: dMO(ncart,nBas,nGrid)
! Initialization
MO(:,:) = 0d0
dMO(:,:,:) = 0d0
do p=1,nBas
do mu=1,nBas
do iG=1,ngrid
MO(p,iG) = MO(p,iG) + c(mu,p)*AO(mu,iG)
dMO(1,p,iG) = dMO(1,p,iG) + c(mu,p)*dAO(1,mu,iG)
dMO(2,p,iG) = dMO(2,p,iG) + c(mu,p)*dAO(2,mu,iG)
dMO(3,p,iG) = dMO(3,p,iG) + c(mu,p)*dAO(3,mu,iG)
end do
end do
end do
end subroutine MO_values_grid

View File

@ -1,10 +0,0 @@
default:
ninja
make -C ..
clean:
ninja -t clean
debug:
ninja -t clean
make -C .. debug

View File

@ -1,121 +0,0 @@
subroutine MinMCMP2(nBas,nEl,nC,nO,nV,c,e,EcMP2, &
nMC,nEq,nWalk,dt,nPrint, &
nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
TrialType,Norm,cTrial,gradient,hessian)
! Minimize the variance of MC-MP2
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nEl,nC,nO,nV,nMC,nEq,nWalk,nPrint
double precision,intent(in) :: EcMP2(3),dt
double precision,intent(in) :: c(nBas,nBas),e(nBas)
integer,intent(in) :: nShell
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: CenterShell(maxShell,3),DShell(maxShell,maxK),ExpShell(maxShell,maxK)
! Local variables
logical :: debug,varmin,mincvg
double precision :: thresh
double precision,allocatable :: max_gradient(:),energy_MCMP2(:),variance_MCMP2(:),error_MCMP2(:),NormTr(:)
double precision :: EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3)
integer :: it,nIt,i
! Output variables
integer,intent(in) :: TrialType
double precision,intent(inout):: Norm,cTrial(nBas),gradient(nBas),hessian(nBas,nBas)
! Debuging mode
! debug = .true.
debug = .false.
! Minimization parameters
varmin = .true.
mincvg = .false.
nIt = 10
thresh = 1d-5
allocate(max_gradient(nIt),energy_MCMP2(nIt),variance_MCMP2(nIt),error_MCMP2(nIt),normTr(nIt))
if(TrialType == 1) then
! Use HOMO as guess
cTrial(1:nBas) = c(1:nBas,nEl/2)
! Normalization factor will be computed later
endif
!------------------------------------------------------------------------
! Start MC-MP2 variance minimization
!------------------------------------------------------------------------
it = 0
do while (it < nIt .and. .not.(mincvg))
it = it + 1
write(*,*) '**********************************************************************'
write(*,*) ' Variance minimization of MC-MP2 energy '
write(*,*) '**********************************************************************'
write(*,*) ' Iteration n.',it
write(*,*) '**********************************************************************'
write(*,*)
write(*,*) ' Trial wave function coefficients at iteration n.',it
call matout(nBas,1,cTrial)
write(*,*)
call MCMP2(varmin,nBas,nEl,nC,nO,nV,c,e,EcMP2, &
nMC,nEq,nWalk,dt,nPrint, &
nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
TrialType,Norm,cTrial,gradient,hessian, &
EcMCMP2,Err_EcMCMP2,Var_EcMCMP2)
! Newton update of the coefficients
call Newton(nBas,gradient,hessian,cTrial)
! Check for convergence
max_gradient(it) = maxval(abs(gradient))
energy_MCMP2(it) = EcMCMP2(1)
variance_MCMP2(it) = Var_EcMCMP2(1)
error_MCMP2(it) = Err_EcMCMP2(1)
NormTr(it) = Norm
write(*,*)
write(*,*) 'Maximum gradient at iteration n.',it,':',max_gradient(it)
write(*,*)
if(max_gradient(it) < thresh) then
write(*,*) ' Miracle! Variance minimization of MC-MP2 has converged!'
mincvg = .true.
endif
enddo
write(*,*)
write(*,*) '********************************'
write(*,*) 'Summary of variance minimization'
write(*,*) '********************************'
write(*,*)
write(*,'(A3,A20,A20,A20,A20,A20,A20)') &
'It.','Gradient','Ec(MC-MPC2)','Variance','Error','Ec(MC-MP2)-Ec(MP2)','Norm'
write(*,'(I3,4X,F16.10,4X,F16.10,4X,F16.10,4X,F16.10,4X,F16.10,4X,F16.10)') &
(i,max_gradient(i),energy_MCMP2(i),variance_MCMP2(i),error_MCMP2(i),energy_MCMP2(i)-EcMP2(1),NormTr(i),i=1,it)
write(*,*)
!------------------------------------------------------------------------
! End MC-MP2 variance minimization
!------------------------------------------------------------------------
end subroutine MinMCMP2

View File

@ -1,67 +0,0 @@
subroutine NDrift(nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,P,r,w,F)
! Compute quantum force numerically
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nWalk,nBas
double precision,intent(in) :: P(nBas,nBas),r(nWalk,2,3),w(nWalk)
integer,intent(in) :: nShell
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: CenterShell(maxShell,3),DShell(maxShell,maxK),ExpShell(maxShell,maxK)
! Local variables
integer :: iW,iEl,ixyz
double precision :: delta
double precision :: wp,wm
double precision,allocatable :: rp(:,:,:),rm(:,:,:),r12p(:),r12m(:)
double precision,allocatable :: gAOp(:,:,:),dgAOp(:,:,:,:),gAOm(:,:,:),dgAOm(:,:,:,:)
double precision,allocatable :: gp(:,:),dgp(:,:,:),gm(:,:),dgm(:,:,:)
! Output variables
double precision,intent(out) :: F(nWalk,2,3)
allocate(rp(nWalk,2,3),rm(nWalk,2,3),r12p(nWalk),r12m(nWalk), &
gAOp(nWalk,2,nBas),dgAOp(nWalk,2,3,nBas),gAOm(nWalk,2,nBas),dgAOm(nWalk,2,3,nBas), &
gp(nWalk,2),dgp(nWalk,2,3),gm(nWalk,2),dgm(nWalk,2,3))
do iW=1,nWalk
do iEl=1,2
do ixyz=1,3
delta = 1d-6
rp = r
rm = r
rp(iW,iEl,ixyz) = r(iW,iEl,ixyz) + delta
rm(iW,iEl,ixyz) = r(iW,iEl,ixyz) - delta
call AO_values(.false.,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,rp,gAOp,dgAOp)
call AO_values(.false.,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,rm,gAOm,dgAOm)
call Density(.false.,nBas,nWalk,P,gAOp,dgAOp,gp,dgp)
call Density(.false.,nBas,nWalk,P,gAOm,dgAOm,gm,dgm)
call rij(nWalk,rp,r12p)
call rij(nWalk,rm,r12m)
wp = gp(iW,1)*gp(iW,2)/r12p(iW)
wm = gm(iW,1)*gm(iW,2)/r12m(iW)
F(iW,iEl,ixyz) = (wp - wm)/(2d0*delta*w(iw))
enddo
enddo
enddo
! print*,'NF',F
end subroutine NDrift

View File

@ -1,67 +0,0 @@
subroutine Newton(nWSq,gradient,hessian,cWeight)
! Calculate the Green functions
implicit none
! Input variables
integer,intent(in) :: nWSq
double precision,intent(in) :: gradient(nWSq),hessian(nWSq,nWSq)
! Local variables
integer :: info
integer,allocatable :: ipiv(:)
double precision,allocatable :: scr(:),eigval(:),eigvec(:,:)
! Output variables
double precision,intent(inout):: cWeight(nWSq)
! Memory allocation
allocate(ipiv(nWSq),scr(3*nWsq),eigval(nWSq),eigvec(nWSq,nWSq))
! Compute eigenvectors and eigenvalues
eigvec = hessian
call dsyev('V','U',nWSq,eigvec,nWSq,eigval,scr,3*nWSq,info)
if(info /= 0)then
write(*,*) ' Problem with dsyev!'
stop
endif
write(*,*)
write(*,*) 'Eigenvalues of hessian'
call matout(nWSq,1,eigval)
write(*,*)
! write(*,*) 'Eigenvectors of hessian'
! call matout(nWSq,1,eigval)
! write(*,*)
! Compute inverse of the hessian
call dgetrf(nWSq,nWSq,hessian,nWSq,ipiv,info)
if(info /= 0) then
write(*,*) ' Problem in dgetrf!'
stop
endif
call dgetri(nWSq,hessian,nWSq,ipiv,scr,nWSq,info)
if(info /= 0) then
write(*,*) ' Problem in dgetri!'
stop
endif
print*,'inverse hessian'
call matout(nWSq,nWSq,hessian)
! Compute new coefficients
cWeight = cWeight - matmul(hessian,gradient)
end subroutine Newton

View File

@ -1,50 +0,0 @@
subroutine drift(nWalk,r,r12,g,dg,F)
! Compute quantum force
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nWalk
double precision,intent(in) :: r(nWalk,2,3),r12(nWalk),g(nWalk,2),dg(nWalk,2,3)
! Local variables
logical :: smoothDrift
double precision :: rij,rijSq,w,wSq,damp
integer :: iW
! Output variables
double precision,intent(out) :: F(nWalk,2,3)
! Compute
smoothDrift = .false.
w = 0.1d0
wSq = w*w
do iW=1,nWalk
rij = r12(iW)
rijSq = rij*rij
F(iW,1,1:3) = dg(iW,1,1:3)/g(iW,1)
F(iW,2,1:3) = dg(iW,2,1:3)/g(iW,2)
if(smoothDrift) then
damp = 1d0 + 2d0*w/sqrt(pi)*rij*exp(-wSq*rijSq)/erfc(w*rij)
else
damp = 1d0
endif
F(iW,1,1:3) = F(iW,1,1:3) - damp*(r(iW,2,1:3) - r(iW,1,1:3))/rijSq
F(iW,2,1:3) = F(iW,2,1:3) - damp*(r(iW,2,1:3) - r(iW,1,1:3))/rijSq
enddo
! print*,' F',F
end subroutine drift

View File

@ -1,30 +0,0 @@
subroutine generate_shell(atot,nShellFunction,ShellFunction)
implicit none
! Input variables
integer,intent(in) :: atot,nShellFunction
! Local variables
integer :: ax,ay,az,ia
! Output variables
integer,intent(out) :: ShellFunction(nShellFunction,3)
ia = 0
do ax=atot,0,-1
do az=0,atot
ay = atot - ax - az
if(ay >= 0) then
ia = ia + 1
ShellFunction(ia,1) = ax
ShellFunction(ia,2) = ay
ShellFunction(ia,3) = az
endif
enddo
enddo
end subroutine generate_shell

View File

@ -1,25 +0,0 @@
subroutine initialize_random_generator(iSeed)
! Initialize random number generator
implicit none
! Input variables
integer,intent(in) :: iSeed
! Local variables
integer,allocatable :: Seed(:)
integer :: nSeed
call random_seed(size = nSeed)
allocate(Seed(nSeed))
call random_seed(get=Seed)
if(iSeed /= 0) then
Seed = 0
Seed(1) = iSeed
endif
call random_seed(put=Seed)
end subroutine initialize_random_generator

View File

@ -1,53 +0,0 @@
subroutine norm_trial(nBas,nO,c,P,Norm,NormSq)
! Initialize weight function
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(inout):: c(nBas,nO),P(nBas,nBas)
! Local variables
double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),G(:,:,:,:)
integer :: mu,nu,la,si
! Output variables
double precision,intent(inout):: Norm,NormSq
! Memory allocation for one- and two-electron integrals
allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),G(nBas,nBas,nBas,nBas))
! Read integrals
call read_integrals(nBas,S,T,V,Hc,G)
! Compute normalization factor
P = 2d0*matmul(c,transpose(c))
Norm = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
Norm = Norm + P(mu,nu)*P(la,si)*G(mu,la,nu,si)
enddo
enddo
enddo
enddo
Norm = Norm*Norm
NormSq = Norm*Norm
write(*,*)
write(*,*) 'Normalization of trial wave function: ',Norm
write(*,*)
end subroutine norm_trial

View File

@ -1,28 +0,0 @@
subroutine optimize_timestep(nWalk,iMC,Acc,dt)
! Optimize dt to get 50% of accepted moves
implicit none
! Input variables
integer,intent(in) :: nWalk,iMC
double precision,intent(inout):: Acc,dt
! Local variables
double precision :: TotAcc,Current_Acc,Target_Acc,delta
TotAcc = Acc/dble(nWalk)
Current_Acc = 100d0*TotAcc/dble(iMC)
Target_Acc = 50.0d0
delta = dt*abs(Target_Acc - Current_Acc)/100.d0
if(Current_Acc > Target_Acc + 0.5d0)then
dt = dt + delta
elseif(Current_Acc < Target_Acc - 0.5d0)then
dt = dt - delta
endif
end subroutine optimize_timestep

View File

@ -1,41 +0,0 @@
subroutine transition_probability(nWalk,dt,D,r,rp,F,Fp,T,Tp)
! Compute transition probability
implicit none
! Input variables
integer,intent(in) :: nWalk
double precision,intent(in) :: dt,D
double precision,intent(in) :: r(nWalk,1:2,1:3), F(nWalk,1:2,1:3)
double precision,intent(in) :: rp(nWalk,1:2,1:3),Fp(nWalk,1:2,1:3)
! Local variables
integer :: iW,iEl,ixyz
! Output variables
double precision,intent(out) :: T(nWalk),Tp(nWalk)
! Initialize
T = 0d0
Tp = 0d0
! Compute
do iW=1,nWalk
do iEl=1,2
do ixyz=1,3
T(iW) = T(iW) + (rp(iW,iEl,ixyz) - r(iW,iEl,ixyz) - D*dt*F(iW,iEl,ixyz))**2
Tp(iW) = Tp(iW) + (r(iW,iEl,ixyz) - rp(iW,iEl,ixyz) - D*dt*Fp(iW,iEl,ixyz))**2
enddo
enddo
enddo
T(:) = exp(-0.25d0*T(:)/(D*dt))
Tp(:) = exp(-0.25d0*Tp(:)/(D*dt))
end subroutine transition_probability