From 44dfbef76634cb3da027005d0e86b4aa2fe936e7 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 12 Feb 2021 16:47:40 +0100 Subject: [PATCH] OK with GGA exchange functionals --- input/dft | 8 +- input/methods | 2 +- input/options | 2 +- src/eDFT/UB88_gga_exchange_energy.f90 | 7 +- src/eDFT/UB88_gga_exchange_potential.f90 | 27 ++++-- src/eDFT/ULYP_gga_correlation_energy.f90 | 73 ++++++++++++++++ src/eDFT/ULYP_gga_correlation_potential.f90 | 85 +++++++++++++++++++ src/eDFT/UPBE_gga_exchange_potential.f90 | 1 - .../unrestricted_gga_correlation_energy.f90 | 18 ++-- 9 files changed, 193 insertions(+), 30 deletions(-) create mode 100644 src/eDFT/ULYP_gga_correlation_energy.f90 create mode 100644 src/eDFT/ULYP_gga_correlation_potential.f90 diff --git a/input/dft b/input/dft index 42c6dab..1b97ce1 100644 --- a/input/dft +++ b/input/dft @@ -4,13 +4,13 @@ # Hartree = 0 # LDA = 1: S51,CC-S51 # GGA = 2: B88,G96,PBE -# Hybrid = 4 +# Hybrid = 4: B3LYP,PBE0 # Hartree-Fock = 666 - 2 PBE + 2 B88 # correlation rung: -# Hartree = 0 +# Hartree = 0: H # LDA = 1: VWN5,eVWN5 -# GGA = 2: +# GGA = 2: LYP,PBE # Hybrid = 4: # Hartree-Fock = 666 0 H diff --git a/input/methods b/input/methods index 94ad703..1da5665 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - F F T F + T F T F # MP2* MP3 MP2-F12 F F F # CCD DCD CCSD CCSD(T) diff --git a/input/options b/input/options index 9d7b8f9..a8c1b79 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.00001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - T T T T T + F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 diff --git a/src/eDFT/UB88_gga_exchange_energy.f90 b/src/eDFT/UB88_gga_exchange_energy.f90 index 2f6b14d..e049b78 100644 --- a/src/eDFT/UB88_gga_exchange_energy.f90 +++ b/src/eDFT/UB88_gga_exchange_energy.f90 @@ -16,7 +16,7 @@ subroutine UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Local variables integer :: iG - double precision :: alpha,beta + double precision :: alpha,b double precision :: r,g,x ! Output variables @@ -26,7 +26,7 @@ subroutine UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Coefficients for B88 GGA exchange functional alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - beta = 0.0042d0 + b = 0.0042d0 ! Compute GGA exchange energy @@ -40,8 +40,7 @@ subroutine UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 x = sqrt(g)/r**(4d0/3d0) - Ex = Ex + weight(iG)*alpha*r**(4d0/3d0) & - - weight(iG)*beta*x**2*r**(4d0/3d0)/(1d0 + 6d0*beta*x*asinh(x)) + Ex = Ex + weight(iG)*r**(4d0/3d0)*(alpha - b*x**2/(1d0 + 6d0*b*x*asinh(x))) end if diff --git a/src/eDFT/UB88_gga_exchange_potential.f90 b/src/eDFT/UB88_gga_exchange_potential.f90 index 14244ff..b8f38f6 100644 --- a/src/eDFT/UB88_gga_exchange_potential.f90 +++ b/src/eDFT/UB88_gga_exchange_potential.f90 @@ -18,8 +18,9 @@ subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Local variables integer :: mu,nu,iG - double precision :: alpha,beta - double precision :: r,g,vAO,gAO + double precision :: alpha,b + double precision :: vAO,gAO + double precision :: r,g,x,dxdr,dxdg,f ! Output variables @@ -28,7 +29,7 @@ subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Coefficients for B88 GGA exchange functional alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - beta = 0.0042d0 + b = 0.0042d0 ! Compute GGA exchange matrix in the AO basis @@ -42,19 +43,27 @@ subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) if(r > threshold) then - g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) - Fx(mu,nu) = Fx(mu,nu) & - + vAO*(4d0/3d0*r**(1d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) & - + 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0)) + + g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 + x = sqrt(g)/r**(4d0/3d0) + dxdr = - 4d0*sqrt(g)/(3d0*r**(7d0/3d0))/x + dxdg = + 1d0/(2d0*sqrt(g)*r**(4d0/3d0))/x + + f = b*x**2/(1d0 + 6d0*b*x*asinh(x)) + + Fx(mu,nu) = Fx(mu,nu) + vAO*( & + 4d0/3d0*r**(1d0/3d0)*(alpha - f) & + - 2d0*r**(4d0/3d0)*dxdr*f & + + r**(4d0/3d0)*dxdr*(6d0*b*x*asinh(x) + 6d0*b*x**2/sqrt(1d0+x**2))*f/(1d0 + 6d0*b*x*asinh(x)) ) gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) & + drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) & + drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG)) - gAO = weight(iG)*gAO - Fx(mu,nu) = Fx(mu,nu) - 2d0*gAO*3d0/4d0*beta*g**(-1d0/4d0)/r**(2d0/3d0) + Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*r**(4d0/3d0)*dxdg*( & + - 2d0*f + (6d0*b*x*asinh(x) + 6d0*b*x**2/sqrt(1d0+x**2))*f/(1d0 + 6d0*b*x*asinh(x)) ) end if diff --git a/src/eDFT/ULYP_gga_correlation_energy.f90 b/src/eDFT/ULYP_gga_correlation_energy.f90 new file mode 100644 index 0000000..974d860 --- /dev/null +++ b/src/eDFT/ULYP_gga_correlation_energy.f90 @@ -0,0 +1,73 @@ +subroutine ULYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec) + +! Compute unrestricted LYP GGA correlation energy + + implicit none + + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid,nspin) + double precision,intent(in) :: drho(ncart,nGrid,nspin) + +! Local variables + + integer :: iG + double precision :: ra,rb,r + double precision :: ga,gab,gb,g + + double precision :: a,b,c,d + double precision :: Cf,omega,delta + +! Output variables + + double precision :: Ec(nsp) + +! Parameters of the functional + + a = 0.04918d0 + b = 0.132d0 + c = 0.2533d0 + d = 0.349d0 + + Cf = 3d0/10d0*(3d0*pi**2)**(2d0/3d0) + +! Initialization + + Ec(:) = 0d0 + + do iG=1,nGrid + + ra = max(0d0,rho(iG,1)) + rb = max(0d0,rho(iG,2)) + r = ra + rb + + if(r > threshold) then + + ga = drho(1,iG,1)**2 + drho(2,iG,1)**2 + drho(3,iG,1)**2 + gb = drho(1,iG,2)**2 + drho(2,iG,2)**2 + drho(3,iG,2)**2 + gab = drho(1,iG,1)*drho(1,iG,2) + drho(2,iG,1)*drho(2,iG,2) + drho(3,iG,1)*drho(3,iG,2) + g = ga + gab + gb + + omega = exp(-c*r**(-1d0/3d0))/(1d0 + d*r**(-1d0/3d0))*r**(-11d0/3d0) + delta = c*r**(-1d0/3d0) + d*r**(-1d0/3d0)/(1d0 + d*r**(-1d0/3d0)) + + Ec(2) = Ec(2) - weight(iG)*4d0*a/(1d0 + d*r**(-1d0/3d0))*ra*rb/r & + - weight(iG)*a*b*omega*ra*rb*( & + 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) & + + (47d0/18d0 - 7d0*delta/18d0)*g & + - (5d0/2d0 - delta/18d0)*(ga + gb) & + - (delta - 11d0)/9d0*(ra/r*ga + rb/r*gb) ) & + - weight(iG)*a*b*omega*( & + - 2d0*r**2/3d0*g & + + (2d0*r**2/3d0 - ra**2)*gb & + + (2d0*r**2/3d0 - rb**2)*ga ) + + end if + + end do + +end subroutine ULYP_gga_correlation_energy diff --git a/src/eDFT/ULYP_gga_correlation_potential.f90 b/src/eDFT/ULYP_gga_correlation_potential.f90 new file mode 100644 index 0000000..18022a8 --- /dev/null +++ b/src/eDFT/ULYP_gga_correlation_potential.f90 @@ -0,0 +1,85 @@ +subroutine ULYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + +! Compute LYP correlation potential + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: dAO(ncart,nBas,nGrid) + double precision,intent(in) :: rho(nGrid,nspin) + double precision,intent(in) :: drho(ncart,nGrid,nspin) + +! Local variables + + integer :: mu,nu,iG + double precision :: vAO,gaAO,gbAO + double precision :: ra,rb,r + double precision :: ga,gab,gb,g + + double precision :: a,b,c,d + double precision :: Cf,omega,delta + +! Output variables + + double precision,intent(out) :: Fc(nBas,nBas) + +! Prameter of the functional + + a = 0.04918d0 + b = 0.132d0 + c = 0.2533d0 + d = 0.349d0 + + Cf = 3d0/10d0*(3d0*pi**2)**(2d0/3d0) + +! Compute matrix elements in the AO basis + + Fc(:,:) = 0d0 + + do mu=1,nBas + do nu=1,nBas + do iG=1,nGrid + + ra = max(0d0,rho(iG,1)) + rb = max(0d0,rho(iG,2)) + r = ra + rb + + if(r > threshold) then + + ga = drho(1,iG,1)**2 + drho(2,iG,1)**2 + drho(3,iG,1)**2 + gb = drho(1,iG,2)**2 + drho(2,iG,2)**2 + drho(3,iG,2)**2 + gab = drho(1,iG,1)*drho(1,iG,2) + drho(2,iG,1)*drho(2,iG,2) + drho(3,iG,1)*drho(3,iG,2) + g = ga + gab + gb + + omega = exp(-c*r**(-1d0/3d0))/(1d0 + d*r**(-1d0/3d0))*r**(-11d0/3d0) + delta = c*r**(-1d0/3d0) + d*r**(-1d0/3d0)/(1d0 + d*r**(-1d0/3d0)) + + vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) + + Fc(mu,nu) = Fc(mu,nu) + vAO + + gaAO = drho(1,iG,1)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) & + + drho(2,iG,1)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) & + + drho(3,iG,1)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG)) + gaAO = weight(iG)*gaAO + + gbAO = drho(1,iG,2)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) & + + drho(2,iG,2)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) & + + drho(3,iG,2)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG)) + gbAO = weight(iG)*gbAO + + Fc(mu,nu) = Fc(mu,nu) + 2d0*gaAO + gbAO + + end if + + end do + end do + end do + +end subroutine ULYP_gga_correlation_potential diff --git a/src/eDFT/UPBE_gga_exchange_potential.f90 b/src/eDFT/UPBE_gga_exchange_potential.f90 index f05e0b6..245eb66 100644 --- a/src/eDFT/UPBE_gga_exchange_potential.f90 +++ b/src/eDFT/UPBE_gga_exchange_potential.f90 @@ -55,7 +55,6 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) & + drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) & + drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG)) - gAO = weight(iG)*gAO Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*alpha*r**(-4d0/3d0)*mupbe/(1d0 + mupbe*s2/kappa)**2 diff --git a/src/eDFT/unrestricted_gga_correlation_energy.f90 b/src/eDFT/unrestricted_gga_correlation_energy.f90 index 1ba444a..34072e4 100644 --- a/src/eDFT/unrestricted_gga_correlation_energy.f90 +++ b/src/eDFT/unrestricted_gga_correlation_energy.f90 @@ -1,6 +1,6 @@ subroutine unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) -! Compute unrstricted GGA correlation energy +! Compute unrestricted GGA correlation energy implicit none include 'parameters.h' @@ -24,19 +24,17 @@ subroutine unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,dr double precision :: Ec(nsp) -! Coefficients for ??? GGA exchange functional + select case (DFA) -! Compute GGA exchange energy + case ('LYP') - Ec(:) = 0d0 + call ULYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec) - do iG=1,nGrid + case default - ra = rho(iG,1) - rb = rho(iG,2) - ga = drho(1,iG,1)**2 + drho(2,iG,1)**2 + drho(3,iG,1)**2 - gb = drho(1,iG,2)**2 + drho(2,iG,2)**2 + drho(3,iG,2)**2 + call print_warning('!!! GGA correlation energy not available !!!') + stop - enddo + end select end subroutine unrestricted_gga_correlation_energy