From 2feb62101d5eeee6c129568d2dfdb42c930b6a30 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 9 Jul 2019 23:09:32 +0200 Subject: [PATCH] GW correction --- input/basis | 29 ++------ input/weight | 29 ++------ src/QuAcK/LSD_SR.f | 89 ++++++++++++++++++++++++- src/QuAcK/MO_values_grid.f90 | 4 +- src/QuAcK/QuAcK.f90 | 88 +++++++++--------------- src/QuAcK/density.f90 | 9 +-- src/QuAcK/ec_srlda.f90 | 12 ++-- src/QuAcK/f_grid.f90 | 24 +------ src/QuAcK/fc_srlda.f90 | 35 +++++++--- src/QuAcK/srlda.f90 | 125 ----------------------------------- 10 files changed, 169 insertions(+), 275 deletions(-) delete mode 100644 src/QuAcK/srlda.f90 diff --git a/input/basis b/input/basis index 05207b0..b9ca7b5 100644 --- a/input/basis +++ b/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 diff --git a/input/weight b/input/weight index 05207b0..b9ca7b5 100644 --- a/input/weight +++ b/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 diff --git a/src/QuAcK/LSD_SR.f b/src/QuAcK/LSD_SR.f index 795cab2..e9f627d 100644 --- a/src/QuAcK/LSD_SR.f +++ b/src/QuAcK/LSD_SR.f @@ -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 diff --git a/src/QuAcK/MO_values_grid.f90 b/src/QuAcK/MO_values_grid.f90 index a0eeb35..ee0864b 100644 --- a/src/QuAcK/MO_values_grid.f90 +++ b/src/QuAcK/MO_values_grid.f90 @@ -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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index b222703..e198060 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -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 !------------------------------------------------------------------------ diff --git a/src/QuAcK/density.f90 b/src/QuAcK/density.f90 index e6bfd22..8e6a574 100644 --- a/src/QuAcK/density.f90 +++ b/src/QuAcK/density.f90 @@ -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 - end subroutine density diff --git a/src/QuAcK/ec_srlda.f90 b/src/QuAcK/ec_srlda.f90 index 7bcac27..3dccbb7 100644 --- a/src/QuAcK/ec_srlda.f90 +++ b/src/QuAcK/ec_srlda.f90 @@ -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 diff --git a/src/QuAcK/f_grid.f90 b/src/QuAcK/f_grid.f90 index edf476f..63bd3f8 100644 --- a/src/QuAcK/f_grid.f90 +++ b/src/QuAcK/f_grid.f90 @@ -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 diff --git a/src/QuAcK/fc_srlda.f90 b/src/QuAcK/fc_srlda.f90 index 891ceb5..139642b 100644 --- a/src/QuAcK/fc_srlda.f90 +++ b/src/QuAcK/fc_srlda.f90 @@ -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 diff --git a/src/QuAcK/srlda.f90 b/src/QuAcK/srlda.f90 deleted file mode 100644 index aa8609e..0000000 --- a/src/QuAcK/srlda.f90 +++ /dev/null @@ -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