mirror of
https://github.com/pfloos/quack
synced 2025-01-11 05:28:32 +01:00
GW correction
This commit is contained in:
parent
7377e5e649
commit
2feb62101d
29
input/basis
29
input/basis
@ -1,24 +1,9 @@
|
||||
1 10
|
||||
S 4 1.00
|
||||
528.5000000 0.0009400
|
||||
79.3100000 0.0072140
|
||||
18.0500000 0.0359750
|
||||
5.0850000 0.1277820
|
||||
1 3
|
||||
S 3 1.00
|
||||
38.3600000 0.0238090
|
||||
5.7700000 0.1548910
|
||||
1.2400000 0.4699870
|
||||
S 1 1.00
|
||||
1.6090000 1.0000000
|
||||
S 1 1.00
|
||||
0.5363000 1.0000000
|
||||
S 1 1.00
|
||||
0.1833000 1.0000000
|
||||
0.2976000 1.0000000
|
||||
P 1 1.00
|
||||
5.9940000 1.0000000
|
||||
P 1 1.00
|
||||
1.7450000 1.0000000
|
||||
P 1 1.00
|
||||
0.5600000 1.0000000
|
||||
D 1 1.00
|
||||
4.2990000 1.0000000
|
||||
D 1 1.00
|
||||
1.2230000 1.0000000
|
||||
F 1 1.00
|
||||
2.6800000 1.0000000
|
||||
1.2750000 1.0000000
|
||||
|
||||
29
input/weight
29
input/weight
@ -1,24 +1,9 @@
|
||||
1 10
|
||||
S 4 1.00
|
||||
528.5000000 0.0009400
|
||||
79.3100000 0.0072140
|
||||
18.0500000 0.0359750
|
||||
5.0850000 0.1277820
|
||||
1 3
|
||||
S 3 1.00
|
||||
38.3600000 0.0238090
|
||||
5.7700000 0.1548910
|
||||
1.2400000 0.4699870
|
||||
S 1 1.00
|
||||
1.6090000 1.0000000
|
||||
S 1 1.00
|
||||
0.5363000 1.0000000
|
||||
S 1 1.00
|
||||
0.1833000 1.0000000
|
||||
0.2976000 1.0000000
|
||||
P 1 1.00
|
||||
5.9940000 1.0000000
|
||||
P 1 1.00
|
||||
1.7450000 1.0000000
|
||||
P 1 1.00
|
||||
0.5600000 1.0000000
|
||||
D 1 1.00
|
||||
4.2990000 1.0000000
|
||||
D 1 1.00
|
||||
1.2230000 1.0000000
|
||||
F 1 1.00
|
||||
2.6800000 1.0000000
|
||||
1.2750000 1.0000000
|
||||
|
||||
@ -1,3 +1,23 @@
|
||||
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)
|
||||
@ -22,7 +42,7 @@ ccc from Paziani, Moroni, Gori-Giorgi, and Bachelet, cond-mat/0601353
|
||||
|
||||
ex=-3.d0*cf/rs/8.d0/pi*((1.d0+z)**(4.d0/3.d0)+
|
||||
$ (1.d0-z)**(4.d0/3.d0))
|
||||
ex=0d0
|
||||
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
|
||||
@ -38,8 +58,10 @@ ccc from Paziani, Moroni, Gori-Giorgi, and Bachelet, cond-mat/0601353
|
||||
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)
|
||||
@ -447,3 +469,68 @@ ccc => vxlrup (spin-up electrons), vxlrdown (spin-down electrons)
|
||||
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
|
||||
|
||||
@ -11,7 +11,7 @@ subroutine MO_values_grid(nBas,nGrid,c,AO,dAO,MO,dMO)
|
||||
integer,intent(in) :: nGrid
|
||||
double precision,intent(in) :: c(nBas,nBas)
|
||||
double precision,intent(in) :: AO(nBas,nGrid)
|
||||
double precision,intent(in) :: dAO(3,nBas,nGrid)
|
||||
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
|
||||
|
||||
! Local variables
|
||||
|
||||
@ -20,7 +20,7 @@ subroutine MO_values_grid(nBas,nGrid,c,AO,dAO,MO,dMO)
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: MO(nBas,nGrid)
|
||||
double precision,intent(out) :: dMO(3,nBas,nGrid)
|
||||
double precision,intent(out) :: dMO(ncart,nBas,nGrid)
|
||||
|
||||
! Initialization
|
||||
|
||||
|
||||
@ -12,6 +12,8 @@ program QuAcK
|
||||
logical :: doG0W0,doevGW,doqsGW
|
||||
logical :: doMCMP2,doMinMCMP2
|
||||
logical :: doeNcusp
|
||||
logical :: doBas
|
||||
|
||||
integer :: nNuc,nBas,nBasCABS
|
||||
integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin),nS(nspin)
|
||||
double precision :: ENuc,ERHF,EUHF,Norm
|
||||
@ -53,6 +55,7 @@ program QuAcK
|
||||
double precision :: start_MP2F12 ,end_MP2F12 ,t_MP2F12
|
||||
double precision :: start_MCMP2 ,end_MCMP2 ,t_MCMP2
|
||||
double precision :: start_MinMCMP2,end_MinMCMP2,t_MinMCMP2
|
||||
double precision :: start_Bas ,end_Bas ,t_Bas
|
||||
|
||||
integer :: maxSCF_HF,n_diis_HF
|
||||
double precision :: thresh_HF
|
||||
@ -76,20 +79,6 @@ program QuAcK
|
||||
double precision :: dt
|
||||
logical :: doDrift
|
||||
|
||||
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(:)
|
||||
|
||||
! Hello World
|
||||
|
||||
write(*,*)
|
||||
@ -226,7 +215,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute RHF energy
|
||||
@ -242,7 +231,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UHF = ',t_HF,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Maximum overlap method
|
||||
@ -259,7 +248,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM = ',t_MOM,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! AO to MO integral transform for post-HF methods
|
||||
@ -303,7 +292,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP2,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute MP3 energy
|
||||
@ -319,7 +308,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP3 = ',t_MP3,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute MP2-F12 energy
|
||||
@ -343,7 +332,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2-F12 = ',t_MP2F12,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform CCD calculation
|
||||
@ -359,7 +348,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform CCSD or CCSD(T) calculation
|
||||
@ -391,7 +380,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CIS = ',t_CIS,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute TDHF excitations
|
||||
@ -407,7 +396,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for TDHF = ',t_TDHF,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute ADC excitations
|
||||
@ -423,7 +412,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ADC = ',t_ADC,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute GF2 electronic binding energies
|
||||
@ -439,7 +428,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF2,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute GF3 electronic binding energies
|
||||
@ -455,7 +444,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF3,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform G0W0 calculatiom
|
||||
@ -474,7 +463,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_G0W0,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform evGW calculation
|
||||
@ -491,7 +480,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_evGW,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform qsGW calculation
|
||||
@ -509,7 +498,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_qsGW,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute e-N cusp dressing
|
||||
@ -524,7 +513,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for e-N cusp dressing = ',t_eNcusp,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Information for Monte Carlo calculations
|
||||
@ -558,7 +547,7 @@ program QuAcK
|
||||
TrialType = 0
|
||||
allocate(cTrial(nBas),gradient(nBas),hessian(nBas,nBas))
|
||||
|
||||
endif
|
||||
end if
|
||||
!------------------------------------------------------------------------
|
||||
! Compute MC-MP2 energy
|
||||
!------------------------------------------------------------------------
|
||||
@ -576,7 +565,7 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MC-MP2 = ',t_MCMP2,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Minimize MC-MP2 variance
|
||||
@ -595,39 +584,26 @@ program QuAcK
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MC-MP2 variance minimization = ',t_MinMCMP2,' seconds'
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Basis set correction
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Construct quadrature grid
|
||||
!------------------------------------------------------------------------
|
||||
doBas = .true.
|
||||
|
||||
SGn = 1
|
||||
if(doBas) then
|
||||
|
||||
call read_grid(SGn,nRad,nAng,nGrid)
|
||||
call cpu_time(start_Bas)
|
||||
call basis_correction(nBas,nO,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
|
||||
ERI_MO_basis,eHF,cHF,PHF,eG0W0)
|
||||
call cpu_time(end_Bas)
|
||||
|
||||
allocate(root(ncart,nGrid),weight(nGrid))
|
||||
|
||||
call quadrature_grid(nRad,nAng,nGrid,root,weight)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Calculate AO values at grid points
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid),MO(nBas,nGrid),dMO(ncart,nBas,nGrid))
|
||||
allocate(rho(nGrid),f(nGrid),mu(nGrid))
|
||||
|
||||
call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO)
|
||||
call density(nGrid,nBas,PHF(:,:,1),AO(:,:),rho(:))
|
||||
call MO_values_grid(nBas,nGrid,cHF(:,:,1),AO,dAO,MO,dMO)
|
||||
call f_grid(nBas,nO(1),nGrid,weight,MO,ERI_MO_basis,f)
|
||||
call mu_grid(nGrid,rho,f,mu)
|
||||
call ec_srlda(nGrid,weight,rho,mu)
|
||||
call fc_srlda(nEl(1),nBas,nGrid,weight,MO,rho,mu)
|
||||
t_Bas = end_Bas - start_Bas
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for basis set correction = ',t_Bas,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
!------------------------------------------------------------------------
|
||||
! End of QuAcK
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
@ -7,8 +7,6 @@ subroutine density(nGrid,nBas,P,AO,rho)
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,parameter :: thresh = 1d-15
|
||||
|
||||
integer,intent(in) :: nGrid
|
||||
integer,intent(in) :: nBas
|
||||
double precision,intent(in) :: P(nBas,nBas)
|
||||
@ -23,16 +21,15 @@ subroutine density(nGrid,nBas,P,AO,rho)
|
||||
double precision,intent(out) :: rho(nGrid)
|
||||
|
||||
rho(:) = 0d0
|
||||
|
||||
do iG=1,nGrid
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! do iG=1,nGrid
|
||||
! rho(iG) = max(rho(iG),thresh)
|
||||
! enddo
|
||||
rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine density
|
||||
|
||||
@ -19,29 +19,27 @@ subroutine ec_srlda(nGrid,weight,rho,mu)
|
||||
double precision :: rs
|
||||
double precision :: ecsr
|
||||
double precision :: ec,vcup,vcdw
|
||||
double precision,parameter :: thres = 1d-15
|
||||
|
||||
! Initialization
|
||||
|
||||
ec = 0d0
|
||||
ecsr = 0d0
|
||||
|
||||
do iG=1,ngrid
|
||||
|
||||
r = max(0d0,rho(iG))
|
||||
|
||||
if(r > thres) then
|
||||
if(r > threshold) then
|
||||
|
||||
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
|
||||
|
||||
! call srlda(rs,mu(iG),ecsr)
|
||||
call lsdsr(rs,0d0,mu(iG),ecsr,vcup,vcdw)
|
||||
call lsdsr(rs,0d0,mu(iG),ec,vcup,vcdw)
|
||||
|
||||
ec = ec + weight(iG)*ecsr*r
|
||||
ecsr = ecsr + weight(iG)*ec*r
|
||||
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
print*, 'ec = ',ec
|
||||
write(*,'(A32,1X,F16.10)') 'ecsr = ',ecsr
|
||||
|
||||
end subroutine ec_srlda
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
subroutine f_grid(nBas,nO,nGrid,weight,MO,ERI,f)
|
||||
subroutine f_grid(nBas,nO,nGrid,MO,ERI,f)
|
||||
|
||||
! Compute f
|
||||
|
||||
@ -10,7 +10,6 @@ subroutine f_grid(nBas,nO,nGrid,weight,MO,ERI,f)
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nGrid
|
||||
double precision,intent(in) :: weight(nGrid)
|
||||
double precision,intent(in) :: MO(nBas,nGrid)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
@ -19,7 +18,6 @@ subroutine f_grid(nBas,nO,nGrid,weight,MO,ERI,f)
|
||||
integer :: p,q
|
||||
integer :: i,j
|
||||
integer :: iG
|
||||
double precision :: toto
|
||||
|
||||
! Output variables
|
||||
|
||||
@ -29,26 +27,6 @@ subroutine f_grid(nBas,nO,nGrid,weight,MO,ERI,f)
|
||||
|
||||
f(:) = 0d0
|
||||
|
||||
do p=1,nBas
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do iG=1,ngrid
|
||||
|
||||
f(iG) = f(iG) + MO(i,iG)*MO(p,iG)*ERI(i,j,p,j)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
toto=0d0
|
||||
do iG=1,nGrid
|
||||
toto = toto + weight(iG)*f(iG)
|
||||
end do
|
||||
print*,'toto=',toto
|
||||
|
||||
f(:) = 0d0
|
||||
|
||||
do p=1,nBas
|
||||
do q=1,nBas
|
||||
do i=1,nO
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
subroutine fc_srlda(nEl,nBas,nGrid,weight,MO,rho,mu)
|
||||
subroutine fc_srlda(nBas,nGrid,weight,MO,rho,mu,eG0W0)
|
||||
|
||||
! Compute sr-lda ec
|
||||
|
||||
@ -7,43 +7,56 @@ subroutine fc_srlda(nEl,nBas,nGrid,weight,MO,rho,mu)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nEl
|
||||
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
|
||||
integer :: iG,p
|
||||
double precision :: r
|
||||
double precision :: rs
|
||||
double precision :: ecsr,vcup,vcdw
|
||||
double precision :: IP
|
||||
double precision,parameter :: thres = 1d-15
|
||||
double precision :: ec,vcup,vcdw
|
||||
double precision,allocatable :: de(:)
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(de(nBas))
|
||||
|
||||
! Initialization
|
||||
|
||||
IP = 0d0
|
||||
de(:) = 0d0
|
||||
|
||||
do iG=1,ngrid
|
||||
|
||||
r = max(0d0,rho(iG))
|
||||
|
||||
if(r > thres) then
|
||||
if(r > threshold) then
|
||||
|
||||
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
|
||||
|
||||
call lsdsr(rs,0d0,mu(iG),ecsr,vcup,vcdw)
|
||||
call lsdsr(rs,0d0,mu(iG),ec,vcup,vcdw)
|
||||
|
||||
IP = IP + weight(iG)*vcup*MO(nEl,iG)**2
|
||||
do p=1,nBas
|
||||
|
||||
de(p)= de(p) + weight(iG)*vcup*MO(p,iG)**2
|
||||
|
||||
end do
|
||||
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
print*, 'IP = ',IP*HaToeV
|
||||
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
|
||||
|
||||
@ -1,125 +0,0 @@
|
||||
subroutine srlda(rs,mu,ecsr)
|
||||
|
||||
! Correlation energy of a spin unpolarized uniform electron gas
|
||||
! with short-range interaction erfc(mu*r)/r
|
||||
! See Zecca et al. PRB 70, 205127 (2004)
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: rs
|
||||
double precision,intent(in) :: mu
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: ec
|
||||
double precision :: cf
|
||||
double precision :: b1
|
||||
double precision :: b2
|
||||
double precision :: b3
|
||||
double precision :: b4
|
||||
double precision :: a0
|
||||
double precision :: bb
|
||||
double precision :: m1
|
||||
|
||||
! Ouput variables
|
||||
|
||||
double precision,intent(out) :: ecsr
|
||||
|
||||
! Compute PW LDA correlation energy
|
||||
|
||||
call ecPW(rs,0d0,ec)
|
||||
|
||||
! Define various stuff
|
||||
|
||||
cf = (9d0*pi/4d0)**(1d0/3d0)
|
||||
bb = 1.27329d0
|
||||
m1 = 0.0357866d0
|
||||
a0 = ec
|
||||
b3 = bb*rs**(7d0/2d0)
|
||||
b2 = -3d0/2d0/pi/cf*rs/a0
|
||||
b1 = (b3-1d0/sqrt(3d0*pi)*rs**(3d0/2d0)/a0)/b2
|
||||
b4 = -a0*b1*rs**3/m1
|
||||
|
||||
! Compute short-range correlation energy
|
||||
|
||||
ecsr = a0*(1d0 + b1*mu)/(1d0 + b1*mu+b2*mu**2 + b3*mu**3 + b4*mu**4)
|
||||
|
||||
end subroutine srlda
|
||||
|
||||
!==========================================================================================
|
||||
|
||||
subroutine ecPW(x,y,ec)
|
||||
|
||||
! Correlation energy of the 3D electron gas of density rs and spin polarization z
|
||||
! Perdew & Wang, PRB 45, 13244 (1992)
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: x
|
||||
double precision,intent(in) :: y
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: f02
|
||||
double precision :: ff
|
||||
double precision :: aaa
|
||||
double precision :: G
|
||||
double precision :: ec0
|
||||
double precision :: ec1
|
||||
double precision :: alfac
|
||||
|
||||
! Ouput variables
|
||||
|
||||
double precision,intent(out) :: ec
|
||||
|
||||
f02 = 4d0/(9d0*(2d0**(1d0/3d0) - 1d0))
|
||||
|
||||
ff = ((1d0+y)**(4d0/3d0) + (1d0-y)**(4d0/3d0)-2d0)/(2d0**(4d0/3d0) - 2d0)
|
||||
|
||||
aaa = (1d0 - log(2d0))/pi**2
|
||||
|
||||
call GPW(x,aaa,0.21370d0,7.5957d0,3.5876d0,1.6382d0,0.49294d0,G)
|
||||
ec0 = G
|
||||
|
||||
aaa=aaa/2d0
|
||||
call GPW(x,aaa,0.20548d0,14.1189d0,6.1977d0,3.3662d0,0.62517d0,G)
|
||||
ec1 = G
|
||||
|
||||
call GPW(x,0.016887d0,0.11125d0,10.357d0,3.6231d0,0.88026d0,0.49671d0,G)
|
||||
alfac = -G
|
||||
|
||||
ec = ec0 + alfac*ff/f02*(1d0 - y**4) + (ec1 - ec0)*ff*y**4
|
||||
|
||||
end subroutine ecPW
|
||||
|
||||
!==========================================================================================
|
||||
|
||||
subroutine GPW(x,Ac,alfa1,beta1,beta2,beta3,beta4,G)
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: Ac
|
||||
double precision,intent(in) :: alfa1
|
||||
double precision,intent(in) :: beta1
|
||||
double precision,intent(in) :: beta2
|
||||
double precision,intent(in) :: beta3
|
||||
double precision,intent(in) :: beta4
|
||||
double precision,intent(in) :: x
|
||||
|
||||
! Ouput variables
|
||||
|
||||
double precision,intent(out) :: G
|
||||
|
||||
G = -2d0*Ac*(1d0 + alfa1*x)*log(1d0 &
|
||||
+ 1d0/(2d0*Ac*(beta1*sqrt(x) + beta2*x + beta3*x*sqrt(x) + beta4*x**2)))
|
||||
|
||||
end subroutine GPW
|
||||
Loading…
Reference in New Issue
Block a user