diff --git a/include/parameters.h b/include/parameters.h index d504995..d0e1202 100644 --- a/include/parameters.h +++ b/include/parameters.h @@ -1,7 +1,7 @@ integer,parameter :: ncart = 3 integer,parameter :: nspin = 2 integer,parameter :: nsp = 3 - integer,parameter :: maxEns = 10 + integer,parameter :: maxEns = 3 integer,parameter :: maxShell = 512 integer,parameter :: maxL = 7 integer,parameter :: n1eInt = 3 @@ -17,7 +17,6 @@ double precision,parameter :: BoToAn = 0.529177249d0 double precision,parameter :: auToD = 2.5415802529d0 - double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0) - double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) - double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) + double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0) + double precision,parameter :: CxLSDA = - (3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) diff --git a/input/dft b/input/dft index 5c2cad4..78b4137 100644 --- a/input/dft +++ b/input/dft @@ -6,18 +6,18 @@ # GGA = 2: B88,G96,PBE # MGGA = 3: # Hybrid = 4: HF,B3,PBE - 4 HF + 4 BHHLYP # correlation rung: # Hartree = 0: H # LDA = 1: VWN5,eVWN5 # GGA = 2: LYP,PBE # MGGA = 3: # Hybrid = 4: HF,LYP,PBE - 4 HF + 4 BHHLYP # quadrature grid SG-n 1 # Number of states in ensemble (nEns) - 3 + 1 # occupation numbers 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/input/methods b/input/methods index ceea951..dbe88da 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - F F T F + T F F F # MP2* MP3 MP2-F12 F F F # CCD DCD CCSD CCSD(T) @@ -9,11 +9,11 @@ # CIS* CIS(D) CID CISD F F F F # RPA* RPAx* ppRPA - F T F + F F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW* - F F F + T T F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 60308cb..1d5bf49 100644 --- a/input/options +++ b/input/options @@ -9,7 +9,7 @@ # 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 - 256 0.0000001 T 5 T 0.00367493 F F F F F + 256 0.0000001 T 5 T 0.000 F F F F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn diff --git a/src/AOtoMO/exchange_matrix_AO_basis.f90 b/src/AOtoMO/exchange_matrix_AO_basis.f90 index 2f0f1d2..8a38c0e 100644 --- a/src/AOtoMO/exchange_matrix_AO_basis.f90 +++ b/src/AOtoMO/exchange_matrix_AO_basis.f90 @@ -1,4 +1,4 @@ -subroutine exchange_matrix_AO_basis(nBas,P,G,K) +subroutine exchange_matrix_AO_basis(nBas,P,ERI,K) ! Compute exchange matrix in the AO basis @@ -9,7 +9,7 @@ subroutine exchange_matrix_AO_basis(nBas,P,G,K) integer,intent(in) :: nBas double precision,intent(in) :: P(nBas,nBas) - double precision,intent(in) :: G(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables @@ -24,7 +24,7 @@ subroutine exchange_matrix_AO_basis(nBas,P,G,K) do si=1,nBas do la=1,nBas do mu=1,nBas - K(mu,nu) = K(mu,nu) - P(la,si)*G(mu,la,si,nu) + K(mu,nu) = K(mu,nu) - P(la,si)*ERI(mu,la,si,nu) enddo enddo enddo diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index a907154..9fd3d22 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,4 +1,4 @@ -subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P) +subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P,Vx) ! Perform restricted Hartree-Fock calculation @@ -55,6 +55,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T double precision,intent(out) :: e(nBas) double precision,intent(out) :: c(nBas,nBas) double precision,intent(out) :: P(nBas,nBas) + double precision,intent(out) :: Vx(nBas) ! Hello world @@ -196,7 +197,13 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) ERHF = ET + EV + EJ + EK +! Compute dipole moments + call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) call print_RHF(nBas,nO,e,C,ENuc,ET,EV,EJ,EK,ERHF,dipole) +! Compute Vx for post-HF calculations + + call exchange_potential(nBas,c,K,Vx) + end subroutine RHF diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 4123dd5..a55e8bb 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -1,4 +1,4 @@ -subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P) +subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P,Vx) ! Perform unrestricted Hartree-Fock calculation @@ -59,6 +59,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO double precision,intent(out) :: e(nBas,nspin) double precision,intent(out) :: c(nBas,nBas,nspin) double precision,intent(out) :: P(nBas,nBas,nspin) + double precision,intent(out) :: Vx(nBas,nspin) ! Hello world @@ -249,4 +250,10 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) +! Compute Vx for post-HF calculations + + do ispin=1,nspin + call exchange_potential(nBas,c(:,:,ispin),K(:,:,ispin),Vx(:,ispin)) + end do + end subroutine UHF diff --git a/src/HF/exchange_potential.f90 b/src/HF/exchange_potential.f90 new file mode 100644 index 0000000..9f72ce7 --- /dev/null +++ b/src/HF/exchange_potential.f90 @@ -0,0 +1,34 @@ +subroutine exchange_potential(nBas,c,Fx,Vx) + +! Compute the exchange potential in the MO basis + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: Fx(nBas,nBas) + +! Local variables + + integer :: mu,nu + integer :: p + +! Output variables + + double precision,intent(out) :: Vx(nBas) + +! Compute Vx + + Vx(:) = 0d0 + do p=1,nBas + do mu=1,nBas + do nu=1,nBas + Vx(p) = Vx(p) + c(mu,p)*Fx(mu,nu)*c(nu,p) + end do + end do + end do + +end subroutine exchange_potential diff --git a/src/MBPT/G0W0.f90 b/src/MBPT/G0W0.f90 index f629708..d679ce8 100644 --- a/src/MBPT/G0W0.f90 +++ b/src/MBPT/G0W0.f90 @@ -1,6 +1,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & dBSE,dTDA,evDyn,singlet,triplet,linearize,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,eGW) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGW) ! Perform G0W0 calculation @@ -29,9 +29,13 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: Vxc(nBas) double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: PHF(nBas,nBas) ! Local variables @@ -41,6 +45,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: EcGM + double precision,allocatable :: SigX(:) double precision,allocatable :: SigC(:) double precision,allocatable :: Z(:) double precision,allocatable :: OmRPA(:) @@ -100,14 +105,14 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & ! Memory allocation - allocate(SigC(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGWlin(nBas)) + allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGWlin(nBas)) !-------------------! ! Compute screening ! !-------------------! call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, & - eHF,ERI,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + eHF,ERI_MO,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) @@ -115,12 +120,13 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & ! Compute spectral weights ! !--------------------------! - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA) !------------------------! ! Compute GW self-energy ! !------------------------! + call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) !--------------------------------! @@ -133,7 +139,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & ! Solve the quasi-particle equation ! !-----------------------------------! - eGWlin(:) = eHF(:) + Z(:)*SigC(:) + eGWlin(:) = eHF(:) + Z(:)*(SigX(:) + SigC(:) - Vxc(:)) ! Linearized or graphical solution? @@ -159,7 +165,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & ! Compute the RPA correlation energy - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & + call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,OmRPA, & rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) !--------------! @@ -180,7 +186,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) if(exchange_kernel) then @@ -214,7 +220,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcAC) + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/MBPT/UG0W0.f90 b/src/MBPT/UG0W0.f90 index 5f8b8b9..5c876cc 100644 --- a/src/MBPT/UG0W0.f90 +++ b/src/MBPT/UG0W0.f90 @@ -1,6 +1,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & - linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & - dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eGW) + linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & + dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eGW) ! Perform unrestricted G0W0 calculation @@ -34,6 +34,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) @@ -42,6 +43,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: PHF(nBas,nBas,nspin) + double precision,intent(in) :: Vxc(nBas,nspin) ! Local variables @@ -52,6 +54,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev double precision :: EcGM(nspin) double precision :: EcBSE(nspin) double precision :: EcAC(nspin) + double precision,allocatable :: SigX(:,:) double precision,allocatable :: SigC(:,:) double precision,allocatable :: Z(:,:) integer :: nS_aa,nS_bb,nS_sc @@ -106,8 +109,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev nS_bb = nS(2) nS_sc = nS_aa + nS_bb - allocate(SigC(nBas,nspin),Z(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc), & - rho_RPA(nBas,nBas,nS_sc,nspin),eGWlin(nBas,nspin)) + allocate(SigX(nBas,nspin),SigC(nBas,nspin),Z(nBas,nspin),eGWlin(nBas,nspin), & + OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) !-------------------! ! Compute screening ! @@ -132,6 +135,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Compute self-energy ! !---------------------! + do is=1,nspin + call self_energy_exchange_diag(nBas,cHF(:,:,is),PHF(:,:,is),ERI,SigX(:,is)) + end do + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) !--------------------------------! @@ -144,7 +151,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Solve the quasi-particle equation ! !-----------------------------------! - eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:) + eGWlin(:,:) = eHF(:,:) + Z(:,:)*(SigX(:,:) + SigC(:,:) - Vxc(:,:)) if(linearize) then @@ -157,10 +164,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Find graphical solution of the QP equation - do is=1,nspin - call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),OmRPA, & - rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is)) - end do + do is=1,nspin + call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),OmRPA, & + rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is)) + end do end if diff --git a/src/MBPT/evGW.f90 b/src/MBPT/evGW.f90 index d471d70..882663d 100644 --- a/src/MBPT/evGW.f90 +++ b/src/MBPT/evGW.f90 @@ -1,6 +1,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI, & - dipole_int,eHF,eG0W0) + G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & + ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -31,9 +31,13 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: Vxc(nBas,nspin) double precision,intent(in) :: eG0W0(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -54,6 +58,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE double precision,allocatable :: eGW(:) double precision,allocatable :: eOld(:) double precision,allocatable :: Z(:) + double precision,allocatable :: SigX(:) double precision,allocatable :: SigC(:) double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) @@ -117,9 +122,13 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! Memory allocation - allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & + allocate(eGW(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & rho_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) +! Compute the exchange part of the self-energy + + call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) + ! Initialization nSCF = 0 @@ -142,14 +151,14 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE if(.not. GW0 .or. nSCF == 0) then - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & + call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,OmRPA, & rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) endif ! Compute spectral weights - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA) ! Compute correlation part of the self-energy @@ -233,7 +242,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE) + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE) if(exchange_kernel) then @@ -267,7 +276,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC) + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/MBPT/evUGW.f90 b/src/MBPT/evUGW.f90 index 5379763..7af60b7 100644 --- a/src/MBPT/evUGW.f90 +++ b/src/MBPT/evUGW.f90 @@ -1,6 +1,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, & - EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eG0W0) + EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -37,10 +37,13 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE integer,intent(in) :: nR(nspin) integer,intent(in) :: nS(nspin) + double precision,intent(in) :: PHF(nBas,nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin) + double precision,intent(in) :: Vxc(nBas,nspin) double precision,intent(in) :: eG0W0(nBas,nspin) double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) @@ -67,6 +70,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE double precision,allocatable :: eOld(:,:) double precision,allocatable :: Z(:,:) integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: SigX(:,:) double precision,allocatable :: SigC(:,:) double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) @@ -127,8 +131,15 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE nS_bb = nS(2) nS_sc = nS_aa + nS_bb - allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc), & - XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin),error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin)) + allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigX(nBas,nspin),SigC(nBas,nspin), & + OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin), & + error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin)) + +! Compute the exchange part of the self-energy + + do is=1,nspin + call self_energy_exchange_diag(nBas,cHF(:,:,is),PHF(:,:,is),ERI_AO,SigX(:,is)) + end do ! Initialization @@ -182,7 +193,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE ! Solve the quasi-particle equation ! !-----------------------------------! - eGW(:,:) = eHF(:,:) + SigC(:,:) + eGW(:,:) = eHF(:,:) + SigX(:,:) + SigC(:,:) - Vxc(:,:) ! Convergence criteria diff --git a/src/MBPT/self_energy_exchange.f90 b/src/MBPT/self_energy_exchange.f90 index 26db034..09653f6 100644 --- a/src/MBPT/self_energy_exchange.f90 +++ b/src/MBPT/self_energy_exchange.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_exchange(nBas,c,P,G,SigmaX) +subroutine self_energy_exchange(nBas,c,P,ERI,SigX) ! Compute exchange part of the self-energy @@ -8,18 +8,18 @@ subroutine self_energy_exchange(nBas,c,P,G,SigmaX) ! Input variables integer,intent(in) :: nBas - double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas),G(nBas,nBas,nBas,nBas) + double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas),ERI(nBas,nBas,nBas,nBas) ! Output variables - double precision,intent(out) :: SigmaX(nBas,nBas) + double precision,intent(out) :: SigX(nBas,nBas) ! Compute exchange part of the self-energy in the AO basis - call exchange_matrix_AO_basis(nBas,P,G,SigmaX) + call exchange_matrix_AO_basis(nBas,P,ERI,SigX) ! Compute exchange part of the self-energy in the MO basis - SigmaX = matmul(transpose(c),matmul(SigmaX,c)) + SigX = matmul(transpose(c),matmul(SigX,c)) end subroutine self_energy_exchange diff --git a/src/MBPT/self_energy_exchange_diag.f90 b/src/MBPT/self_energy_exchange_diag.f90 new file mode 100644 index 0000000..0e305b7 --- /dev/null +++ b/src/MBPT/self_energy_exchange_diag.f90 @@ -0,0 +1,42 @@ +subroutine self_energy_exchange_diag(nBas,c,P,ERI,SigX) + +! Compute the diagonal elements of the exchange part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: q,mu,nu + double precision,allocatable :: Fx(:,:) + +! Output variables + + double precision,intent(out) :: SigX(nBas) + +! Compute exchange part of the self-energy in the AO basis + + allocate(Fx(nBas,nBas)) + call exchange_matrix_AO_basis(nBas,P,ERI,Fx) + +! Compute exchange part of the self-energy in the MO basis + + SigX(:) = 0d0 + do q=1,nBas + do mu=1,nBas + do nu=1,nBas + SigX(q) = SigX(q) + c(mu,q)*Fx(mu,nu)*c(nu,q) + end do + end do + end do + + deallocate(Fx) + +end subroutine self_energy_exchange_diag diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index c1154f5..ec48069 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -32,6 +32,7 @@ program QuAcK double precision,allocatable :: ZNuc(:),rNuc(:,:) double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:) + double precision,allocatable :: Vxc(:,:,:) double precision,allocatable :: eG0W0(:,:) double precision,allocatable :: eG0T0(:,:) @@ -236,7 +237,7 @@ program QuAcK allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), & S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas), & - dipole_int_AO(nBas,nBas,ncart),dipole_int_MO(nBas,nBas,ncart)) + dipole_int_AO(nBas,nBas,ncart),dipole_int_MO(nBas,nBas,ncart),Vxc(nBas,nBas,nspin)) ! Read integrals @@ -280,7 +281,7 @@ program QuAcK call cpu_time(start_HF) call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF) + nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc) call cpu_time(end_HF) t_HF = end_HF - start_HF @@ -300,7 +301,7 @@ program QuAcK call cpu_time(start_HF) call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF) + nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF,Vxc) call cpu_time(end_HF) t_HF = end_HF - start_HF @@ -315,10 +316,13 @@ program QuAcK if(doKS) then + ! Switch on the unrestricted flag + unrestricted = .true. + call cpu_time(start_KS) call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, & nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & - max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EUHF,eHF,cHF,PHF) + max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EUHF,eHF,cHF,PHF,Vxc) call cpu_time(end_KS) @@ -851,12 +855,12 @@ program QuAcK if(unrestricted) then call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & - linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & - dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0) + linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & + dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0) else call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & - linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,eG0W0) + linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) end if @@ -879,14 +883,14 @@ program QuAcK call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc, & - EUHF,S,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, & - cHF,eHF,eG0W0) + EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, & + PHF,cHF,eHF,Vxc,eG0W0) else call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) end if call cpu_time(end_evGW) @@ -1005,10 +1009,10 @@ program QuAcK if(doMCMP2) then call cpu_time(start_MCMP2) - call MCMP2(doDrift,nBas,nC,nO,nV,cHF,eHF,EcMP2, & - nMC,nEq,nWalk,dt,nPrint, & - nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - Norm,EcMCMP2,Err_EcMCMP2,Var_EcMCMP2) +! call MCMP2(doDrift,nBas,nC,nO,nV,cHF,eHF,EcMP2, & +! nMC,nEq,nWalk,dt,nPrint, & +! nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & +! Norm,EcMCMP2,Err_EcMCMP2,Var_EcMCMP2) call cpu_time(end_MCMP2) t_MCMP2 = end_MCMP2 - start_MCMP2 @@ -1040,59 +1044,59 @@ program QuAcK ! Range-separeted GT/GW !------------------------------------------------------------------------ - if(doGTGW) then +! if(doGTGW) then - ! Read and transform long-range two-electron integrals +! ! Read and transform long-range two-electron integrals - allocate(ERI_ERF_AO(nBas,nBas,nBas,nBas),ERI_ERF_MO(nBas,nBas,nBas,nBas)) - call read_LR(nBas,ERI_ERF_AO) +! allocate(ERI_ERF_AO(nBas,nBas,nBas,nBas),ERI_ERF_MO(nBas,nBas,nBas,nBas)) +! call read_LR(nBas,ERI_ERF_AO) - call cpu_time(start_AOtoMO) +! call cpu_time(start_AOtoMO) - write(*,*) - write(*,*) 'AO to MO transformation for long-range ERIs... Please be patient' - write(*,*) +! write(*,*) +! write(*,*) 'AO to MO transformation for long-range ERIs... Please be patient' +! write(*,*) - call AOtoMO_integral_transform(nBas,cHF,ERI_ERF_AO,ERI_ERF_MO) +! call AOtoMO_integral_transform(nBas,cHF,ERI_ERF_AO,ERI_ERF_MO) - call cpu_time(end_AOtoMO) +! call cpu_time(end_AOtoMO) - deallocate(ERI_ERF_AO) +! deallocate(ERI_ERF_AO) - t_AOtoMO = end_AOtoMO - start_AOtoMO +! t_AOtoMO = end_AOtoMO - start_AOtoMO - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' - write(*,*) +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' +! write(*,*) - ! Long-range G0W0 calculation +! ! Long-range G0W0 calculation - call cpu_time(start_G0W0) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0W0) - call cpu_time(end_G0W0) - - t_G0W0 = end_G0W0 - start_G0W0 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_G0W0,' seconds' - write(*,*) +! call cpu_time(start_G0W0) +! call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & +! dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & +! nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0W0) +! call cpu_time(end_G0W0) +! +! t_G0W0 = end_G0W0 - start_G0W0 +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_G0W0,' seconds' +! write(*,*) - ! Short-range G0T0 calculation +! ! Short-range G0T0 calculation - ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:) +! ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:) - call cpu_time(start_G0T0) - call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, & - singlet,triplet,linGW,eta_GW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0T0) - call cpu_time(end_G0T0) - - t_G0T0 = end_G0T0 - start_G0T0 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds' - write(*,*) +! call cpu_time(start_G0T0) +! call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, & +! singlet,triplet,linGW,eta_GW, & +! nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0T0) +! call cpu_time(end_G0T0) +! +! t_G0T0 = end_G0T0 - start_G0T0 +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds' +! write(*,*) ! call matout(nBas,1,(eG0W0+eG0T0-eHF(:,1))*HaToeV) - end if +! end if !------------------------------------------------------------------------ ! Basis set correction diff --git a/src/eDFT/RB88_gga_exchange_energy.f90 b/src/eDFT/RB88_gga_exchange_energy.f90 index 04c87f6..596710d 100644 --- a/src/eDFT/RB88_gga_exchange_energy.f90 +++ b/src/eDFT/RB88_gga_exchange_energy.f90 @@ -16,7 +16,6 @@ subroutine RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Local variables integer :: iG - double precision :: alpha double precision :: beta double precision :: r,g,x @@ -26,7 +25,6 @@ subroutine RB88_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 ! Compute GGA exchange energy @@ -41,7 +39,7 @@ subroutine RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) g = 0.25d0*(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) & + Ex = Ex + weight(iG)*CxLDA*r**(4d0/3d0) & - weight(iG)*beta*x**2*r**(4d0/3d0)/(1d0 + 6d0*beta*x*asinh(x)) end if diff --git a/src/eDFT/RB88_gga_exchange_individual_energy.f90 b/src/eDFT/RB88_gga_exchange_individual_energy.f90 index 23c05bb..56c64bf 100644 --- a/src/eDFT/RB88_gga_exchange_individual_energy.f90 +++ b/src/eDFT/RB88_gga_exchange_individual_energy.f90 @@ -17,7 +17,6 @@ subroutine RB88_gga_exchange_individual_energy(nGrid,weight,rhow,drhow,rho,drho, ! Local variables integer :: iG - double precision :: alpha double precision :: beta double precision :: r,rI,g,x double precision :: ex_p,dexdr_p @@ -28,7 +27,6 @@ subroutine RB88_gga_exchange_individual_energy(nGrid,weight,rhow,drhow,rho,drho, ! Coefficients for B88 GGA exchange functional - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) beta = 0.0042d0 ! Compute GGA exchange matrix in the AO basis @@ -45,11 +43,11 @@ subroutine RB88_gga_exchange_individual_energy(nGrid,weight,rhow,drhow,rho,drho, g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2) x = sqrt(g)/r**(4d0/3d0) - dexdr_p = 4d0/3d0*r**(1d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) & + dexdr_p = 4d0/3d0*r**(1d0/3d0)*(CxLDA - beta*g**(3d0/4d0)/r**2) & + 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0) & - 2d0*3d0/4d0*beta*g**(-1d0/4d0)/r**(2d0/3d0) - ex_p = alpha*r**(4d0/3d0) & + ex_p = CxLDA*r**(4d0/3d0) & - weight(iG)*beta*x**2*r**(4d0/3d0)/(1d0 + 6d0*beta*x*asinh(x)) Ex = Ex + weight(iG)*(ex_p*rI + dexdr_p*r*rI - dexdr_p*r*r) diff --git a/src/eDFT/RB88_gga_exchange_potential.f90 b/src/eDFT/RB88_gga_exchange_potential.f90 index 1dca23b..6829644 100644 --- a/src/eDFT/RB88_gga_exchange_potential.f90 +++ b/src/eDFT/RB88_gga_exchange_potential.f90 @@ -18,7 +18,6 @@ subroutine RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Local variables integer :: mu,nu,iG - double precision :: alpha double precision :: beta double precision :: r,g,vAO,gAO @@ -28,7 +27,6 @@ subroutine RB88_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 ! Compute GGA exchange matrix in the AO basis @@ -46,7 +44,7 @@ subroutine RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) g = 0.25d0*(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) & + + vAO*(4d0/3d0*r**(1d0/3d0)*(CxLDA - beta*g**(3d0/4d0)/r**2) & + 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0)) gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) & diff --git a/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 index 235407a..14d1083 100644 --- a/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 @@ -84,12 +84,16 @@ subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)) & * (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0))) + + case default + dCxdw1 = 0.d0 + dCxdw2 = 0.d0 + end select dCxdw1 = CxLDA*dCxdw1 dCxdw2 = CxLDA*dCxdw2 - dExdw(:) = 0d0 do iG=1,nGrid diff --git a/src/eDFT/RCC_lda_exchange_energy.f90 b/src/eDFT/RCC_lda_exchange_energy.f90 index 8fa5c6c..5f55d95 100644 --- a/src/eDFT/RCC_lda_exchange_energy.f90 +++ b/src/eDFT/RCC_lda_exchange_energy.f90 @@ -69,12 +69,19 @@ subroutine RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,C Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2) select case (Cx_choice) - case(1) - Cx = CxLDA*Fx1 - case(2) - Cx = CxLDA*Fx2 - case(3) - Cx = CxLDA*Fx2*Fx1 + + case(1) + Cx = CxLDA*Fx1 + + case(2) + Cx = CxLDA*Fx2 + + case(3) + Cx = CxLDA*Fx2*Fx1 + + case default + Cx = CxLDA + end select ! Compute GIC-LDA exchange energy diff --git a/src/eDFT/RCC_lda_exchange_individual_energy.f90 b/src/eDFT/RCC_lda_exchange_individual_energy.f90 index 02f23f8..ade4f56 100644 --- a/src/eDFT/RCC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RCC_lda_exchange_individual_energy.f90 @@ -75,12 +75,19 @@ subroutine RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2) select case (Cx_choice) - case(1) - Cx = CxLDA*Fx1 - case(2) - Cx = CxLDA*Fx2 - case(3) - Cx = CxLDA*Fx2*Fx1 + + case(1) + Cx = CxLDA*Fx1 + + case(2) + Cx = CxLDA*Fx2 + + case(3) + Cx = CxLDA*Fx2*Fx1 + + case default + Cx = CxLDA + end select ! Compute LDA exchange matrix in the AO basis diff --git a/src/eDFT/RCC_lda_exchange_potential.f90 b/src/eDFT/RCC_lda_exchange_potential.f90 index d33017b..519dfa9 100644 --- a/src/eDFT/RCC_lda_exchange_potential.f90 +++ b/src/eDFT/RCC_lda_exchange_potential.f90 @@ -74,12 +74,19 @@ subroutine RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2) select case (Cx_choice) - case(1) - Cx = CxLDA*Fx1 - case(2) - Cx = CxLDA*Fx2 - case(3) - Cx = CxLDA*Fx2*Fx1 + + case(1) + Cx = CxLDA*Fx1 + + case(2) + Cx = CxLDA*Fx2 + + case(3) + Cx = CxLDA*Fx2*Fx1 + + case default + Cx = CxLDA + end select diff --git a/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 index 15e1c3b..d3620c8 100644 --- a/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 @@ -21,6 +21,9 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r double precision :: dExdw(nEns) double precision,external :: Kronecker_delta + double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) + double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) + ! Output variables double precision,intent(out) :: ExDD(nEns) diff --git a/src/eDFT/RMFL20_lda_exchange_energy.f90 b/src/eDFT/RMFL20_lda_exchange_energy.f90 index c9ca09f..3d50004 100644 --- a/src/eDFT/RMFL20_lda_exchange_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_energy.f90 @@ -15,6 +15,9 @@ subroutine RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rho(nGrid) + double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) + double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) + ! Local variables integer :: iG diff --git a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 index 75176b0..49a1299 100644 --- a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 @@ -22,6 +22,9 @@ subroutine RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,we double precision :: r,rI double precision :: e_p,dedr + double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) + double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) + ! Output variables double precision,intent(out) :: Ex diff --git a/src/eDFT/RMFL20_lda_exchange_potential.f90 b/src/eDFT/RMFL20_lda_exchange_potential.f90 index ce6cfca..4d0ca82 100644 --- a/src/eDFT/RMFL20_lda_exchange_potential.f90 +++ b/src/eDFT/RMFL20_lda_exchange_potential.f90 @@ -22,6 +22,9 @@ subroutine RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBa double precision :: Cxw double precision :: r,vAO + double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) + double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) + ! Output variables double precision,intent(out) :: Fx(nBas,nBas) diff --git a/src/eDFT/UB88_gga_exchange_energy.f90 b/src/eDFT/UB88_gga_exchange_energy.f90 index e049b78..2ecf64d 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,b + double precision :: b double precision :: r,g,x ! Output variables @@ -25,7 +25,6 @@ subroutine UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Coefficients for B88 GGA exchange functional - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) b = 0.0042d0 ! Compute GGA exchange energy @@ -40,7 +39,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)*r**(4d0/3d0)*(alpha - b*x**2/(1d0 + 6d0*b*x*asinh(x))) + Ex = Ex + weight(iG)*r**(4d0/3d0)*(CxLSDA - 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 b8f38f6..fcec63c 100644 --- a/src/eDFT/UB88_gga_exchange_potential.f90 +++ b/src/eDFT/UB88_gga_exchange_potential.f90 @@ -18,7 +18,7 @@ subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Local variables integer :: mu,nu,iG - double precision :: alpha,b + double precision :: b double precision :: vAO,gAO double precision :: r,g,x,dxdr,dxdg,f @@ -28,7 +28,6 @@ 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) b = 0.0042d0 ! Compute GGA exchange matrix in the AO basis @@ -52,9 +51,9 @@ subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) 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 & + Fx(mu,nu) = Fx(mu,nu) + vAO*( & + 4d0/3d0*r**(1d0/3d0)*(CxLSDA - 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)) & diff --git a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 index 660e98b..dbe7cca 100644 --- a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 @@ -22,7 +22,7 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr integer :: iEns,jEns integer :: iG - double precision :: r,alpha + double precision :: r double precision,allocatable :: dExdw(:) double precision,external :: Kronecker_delta @@ -30,9 +30,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr double precision :: a2,b2,c2,w2 double precision :: dCxdw1,dCxdw2 - double precision :: nEli,nElw - - ! Output variables double precision,intent(out) :: ExDD(nEns) @@ -60,13 +57,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr b2 = aCC_w2(2) c2 = aCC_w2(3) - nElw = electron_number(nGrid,weight,rhow) - - -! Cx coefficient for unrestricted Slater LDA exchange - - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - w1 = wEns(2) w2 = wEns(3) @@ -86,11 +76,15 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)) & * (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0))) + + case default + dCxdw1 = 0d0 + dCxdw2 = 0d0 + end select - - dCxdw1 = alpha*dCxdw1 - dCxdw2 = alpha*dCxdw2 + dCxdw1 = CxLSDA*dCxdw1 + dCxdw2 = CxLSDA*dCxdw2 dExdw(:) = 0d0 diff --git a/src/eDFT/UCC_lda_exchange_energy.f90 b/src/eDFT/UCC_lda_exchange_energy.f90 index b3fcbaf..0964c15 100644 --- a/src/eDFT/UCC_lda_exchange_energy.f90 +++ b/src/eDFT/UCC_lda_exchange_energy.f90 @@ -19,7 +19,7 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,C ! Local variables integer :: iG - double precision :: r,alpha + double precision :: r double precision :: a1,b1,c1,w1 double precision :: a2,b2,c2,w2 @@ -65,10 +65,6 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,C b2 = aCC_w2(2) c2 = aCC_w2(3) -! Cx coefficient for unrestricted Slater LDA exchange - - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - ! Fx1 for states N and N-1 ! Fx2 for states N and N+1 @@ -79,21 +75,20 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,C Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2) select case (Cx_choice) - case(1) - Cx = alpha*Fx1 - case(2) - Cx = alpha*Fx2 - case(3) - Cx = alpha*Fx2*Fx1 + + case(1) + Cx = CxLSDA*Fx1 + + case(2) + Cx = CxLSDA*Fx2 + + case(3) + Cx = CxLSDA*Fx2*Fx1 + + case default + Cx = CxLSDA + end select -! for two-weights ensemble -! Cx = alpha*Fx2*Fx1 - -! for left ensemble -! Cx = alpha*Fx1 - -! for right ensemble -! Cx = alpha*Fx2 ! Compute GIC-LDA exchange energy diff --git a/src/eDFT/UCC_lda_exchange_individual_energy.f90 b/src/eDFT/UCC_lda_exchange_individual_energy.f90 index 98072ca..a0a0df3 100644 --- a/src/eDFT/UCC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/UCC_lda_exchange_individual_energy.f90 @@ -22,9 +22,8 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig ! Local variables integer :: iG - double precision :: r,rI,alpha + double precision :: r,rI double precision :: e_p,dedr - double precision :: nEli,nElw double precision :: a1,b1,c1,w1 double precision :: a2,b2,c2,w2 @@ -51,10 +50,6 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig b2 = aCC_w2(2) c2 = aCC_w2(3) -! Cx coefficient for unrestricted Slater LDA exchange - - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - w1 = wEns(2) Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2) @@ -64,21 +59,19 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig select case (Cx_choice) case(1) - Cx = alpha*Fx1 + Cx = CxLSDA*Fx1 case(2) - Cx = alpha*Fx2 + Cx = CxLSDA*Fx2 case(3) - Cx = alpha*Fx2*Fx1 + Cx = CxLSDA*Fx2*Fx1 + + case default + Cx = CxLSDA end select - nEli = electron_number(nGrid,weight,rho) - - nElw = electron_number(nGrid,weight,rhow) - - ! Compute LDA exchange matrix in the AO basis Ex = 0d0 @@ -92,19 +85,12 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig e_p = Cx*r**(1d0/3d0) dedr = 1d0/3d0*Cx*r**(-2d0/3d0) - if (doNcentered) then - Ex = Ex - weight(iG)*dedr*r*r*(nEli/nElw) - else - Ex = Ex - weight(iG)*dedr*r*r - end if + Ex = Ex - weight(iG)*dedr*r*r if(rI > threshold) then - if (doNcentered) then - Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI) - else - Ex = Ex + weight(iG)*((nEli/nElw)*e_p*rI + dedr*r*rI) - end if + Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI) + endif endif diff --git a/src/eDFT/UCC_lda_exchange_potential.f90 b/src/eDFT/UCC_lda_exchange_potential.f90 index dc14909..ecb027f 100644 --- a/src/eDFT/UCC_lda_exchange_potential.f90 +++ b/src/eDFT/UCC_lda_exchange_potential.f90 @@ -21,7 +21,7 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, ! Local variables integer :: mu,nu,iG - double precision :: r,vAO,alpha + double precision :: r,vAO double precision :: a1,b1,c1,w1 double precision :: a2,b2,c2,w2 @@ -67,10 +67,6 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, b2 = aCC_w2(2) c2 = aCC_w2(3) -! Cx coefficient for unrestricted Slater LDA exchange - - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - ! Fx1 for states N and N-1 ! Fx2 for states N and N+1 @@ -81,23 +77,21 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2) select case (Cx_choice) - case(1) - Cx = alpha*Fx1 - case(2) - Cx = alpha*Fx2 - case(3) - Cx = alpha*Fx2*Fx1 + + case(1) + Cx = CxLSDA*Fx1 + + case(2) + Cx = CxLSDA*Fx2 + + case(3) + Cx = CxLSDA*Fx2*Fx1 + + case default + Cx = CxLSDA + end select -! for two-weight ensembles -! Cx = alpha*Fx2*Fx1 - -! for left ensemble -! Cx = alpha*Fx1 - -! for right ensemble -! Cx = alpha*Fx2 - ! Compute LDA exchange matrix in the AO basis Fx(:,:) = 0d0 diff --git a/src/eDFT/UG96_gga_exchange_energy.f90 b/src/eDFT/UG96_gga_exchange_energy.f90 index 1876cc5..f09d2a7 100644 --- a/src/eDFT/UG96_gga_exchange_energy.f90 +++ b/src/eDFT/UG96_gga_exchange_energy.f90 @@ -16,7 +16,7 @@ subroutine UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Local variables integer :: iG - double precision :: alpha,beta + double precision :: beta double precision :: r,g ! Output variables @@ -25,7 +25,6 @@ subroutine UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Coefficients for G96 GGA exchange functional - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) beta = 1d0/137d0 ! Compute GGA exchange energy @@ -40,7 +39,7 @@ subroutine UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 - Ex = Ex + weight(iG)*r**(4d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) + Ex = Ex + weight(iG)*r**(4d0/3d0)*(CxLSDA - beta*g**(3d0/4d0)/r**2) end if diff --git a/src/eDFT/UG96_gga_exchange_potential.f90 b/src/eDFT/UG96_gga_exchange_potential.f90 index c8d9b40..4b28a3c 100644 --- a/src/eDFT/UG96_gga_exchange_potential.f90 +++ b/src/eDFT/UG96_gga_exchange_potential.f90 @@ -18,7 +18,7 @@ subroutine UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Local variables integer :: mu,nu,iG - double precision :: alpha,beta + double precision :: beta double precision :: r,g,vAO,gAO ! Output variables @@ -27,7 +27,6 @@ subroutine UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Coefficients for G96 GGA exchange functional - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) beta = +1d0/137d0 ! Compute GGA exchange matrix in the AO basis @@ -45,7 +44,7 @@ subroutine UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) 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) & + + vAO*(4d0/3d0*r**(1d0/3d0)*(CxLSDA - beta*g**(3d0/4d0)/r**2) & + 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0)) gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) & diff --git a/src/eDFT/UPBE_gga_exchange_energy.f90 b/src/eDFT/UPBE_gga_exchange_energy.f90 index d06f537..0cc35bd 100644 --- a/src/eDFT/UPBE_gga_exchange_energy.f90 +++ b/src/eDFT/UPBE_gga_exchange_energy.f90 @@ -16,7 +16,7 @@ subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Local variables integer :: iG - double precision :: alpha,mupbe,kappa + double precision :: mupbe,kappa double precision :: r,g,s2 ! Output variables @@ -25,7 +25,6 @@ subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Coefficients for PBE exchange functional - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) mupbe = ((1d0/2d0**(1d0/3d0))/(2d0*(3d0*pi**2)**(1d0/3d0)))**2*0.21951d0 kappa = 0.804d0 @@ -41,7 +40,7 @@ subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 s2 = g/r**(8d0/3d0) - Ex = Ex + weight(iG)*alpha*r**(4d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa)) + Ex = Ex + weight(iG)*CxLSDA*r**(4d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa)) end if diff --git a/src/eDFT/UPBE_gga_exchange_potential.f90 b/src/eDFT/UPBE_gga_exchange_potential.f90 index 245eb66..866b90a 100644 --- a/src/eDFT/UPBE_gga_exchange_potential.f90 +++ b/src/eDFT/UPBE_gga_exchange_potential.f90 @@ -18,7 +18,7 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Local variables integer :: mu,nu,iG - double precision :: alpha,mupbe,kappa + double precision :: mupbe,kappa double precision :: r,g,s2,vAO,gAO ! Output variables @@ -27,7 +27,6 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Coefficients for PBE exchange functional - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) mupbe = ((1d0/2d0**(1d0/3d0))/(2d0*(3d0*pi**2)**(1d0/3d0)))**2*0.21951d0 kappa = 0.804d0 @@ -49,15 +48,15 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) Fx(mu,nu) = Fx(mu,nu) & - + vAO*4d0/3d0*alpha*r**(1d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa)) & - - vAO*8d0/3d0*alpha*r**(1d0/3d0)*mupbe*s2/(1d0 + mupbe*s2/kappa)**2 + + vAO*4d0/3d0*CxLSDA*r**(1d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa)) & + - vAO*8d0/3d0*CxLSDA*r**(1d0/3d0)*mupbe*s2/(1d0 + mupbe*s2/kappa)**2 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 + Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*CxLSDA*r**(-4d0/3d0)*mupbe/(1d0 + mupbe*s2/kappa)**2 end if diff --git a/src/eDFT/US51_lda_exchange_energy.f90 b/src/eDFT/US51_lda_exchange_energy.f90 index c0b8702..fbaaad9 100644 --- a/src/eDFT/US51_lda_exchange_energy.f90 +++ b/src/eDFT/US51_lda_exchange_energy.f90 @@ -16,16 +16,12 @@ subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex) ! Local variables integer :: iG - double precision :: alpha,r,alphaw,a2,b2,c2,a1,b1,c1 + double precision :: r ! Output variables double precision :: Ex -! Cx coefficient for Slater LDA exchange - - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - ! Compute LDA exchange energy Ex = 0d0 @@ -35,7 +31,7 @@ subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex) if(r > threshold) then - Ex = Ex + weight(iG)*alpha*r**(4d0/3d0) + Ex = Ex + weight(iG)*CxLSDA*r**(4d0/3d0) endif diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/US51_lda_exchange_individual_energy.f90 index 36a8acf..f1541a5 100644 --- a/src/eDFT/US51_lda_exchange_individual_energy.f90 +++ b/src/eDFT/US51_lda_exchange_individual_energy.f90 @@ -17,7 +17,7 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered ! Local variables integer :: iG - double precision :: r,rI,alpha + double precision :: r,rI double precision :: e,dedr double precision :: Exrr,ExrI,ExrrI @@ -27,8 +27,6 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered ! Compute LDA exchange matrix in the AO basis - alpha = - (3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - Exrr = 0d0 ExrI = 0d0 ExrrI = 0d0 @@ -40,8 +38,8 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered if(r > threshold) then - e = alpha*r**(1d0/3d0) - dedr = 1d0/3d0*alpha*r**(-2d0/3d0) + e = CxLSDA*r**(1d0/3d0) + dedr = 1d0/3d0*CxLSDA*r**(-2d0/3d0) Exrr = Exrr - weight(iG)*dedr*r*r diff --git a/src/eDFT/US51_lda_exchange_potential.f90 b/src/eDFT/US51_lda_exchange_potential.f90 index 4066379..9c2bc0d 100644 --- a/src/eDFT/US51_lda_exchange_potential.f90 +++ b/src/eDFT/US51_lda_exchange_potential.f90 @@ -16,17 +16,12 @@ subroutine US51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) ! Local variables integer :: mu,nu,iG - double precision :: alpha double precision :: r,vAO ! Output variables double precision,intent(out) :: Fx(nBas,nBas) -! Cx coefficient for Slater LDA exchange - - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - ! Compute LDA exchange matrix in the AO basis Fx(:,:) = 0d0 @@ -39,7 +34,7 @@ subroutine US51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) if(r > threshold) then vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) - Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*alpha*r**(1d0/3d0) + Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxLSDA*r**(1d0/3d0) endif diff --git a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 index beb8031..d856c41 100644 --- a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 @@ -88,7 +88,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) decdra_f = drsdra*dxdrs*decdx_f - Ec(1) = Ec(1) + weight(iG)*(ec_z + decdra_f*r)*rI + Ec(1) = Ec(1) + weight(iG)*(ec_f + decdra_f*r)*rI end if @@ -191,7 +191,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) decdra_f = drsdra*dxdrs*decdx_f - Ec(3) = Ec(3) + weight(iG)*(ec_z + decdra_f*r)*rI + Ec(3) = Ec(3) + weight(iG)*(ec_f + decdra_f*r)*rI end if diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 2f9ec5a..bbfbed6 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -1,6 +1,6 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, & nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & - max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI,dipole_int,Ew,eKS,cKS,PKS) + max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI,dipole_int,Ew,eKS,cKS,PKS,Vxc) ! exchange-correlation density-functional theory calculations @@ -50,8 +50,6 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n ! Local variables - double precision,allocatable :: c(:,:) - character(len=8) :: method integer :: x_rung,c_rung character(len=12) :: x_DFA ,c_DFA @@ -88,6 +86,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n double precision,intent(out) :: eKS(nBas,nspin) double precision,intent(out) :: cKS(nBas,nBas,nspin) double precision,intent(out) :: PKS(nBas,nBas,nspin) + double precision,intent(out) :: Vxc(nBas,nspin) ! Hello World @@ -111,7 +110,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n ! Allocate ensemble weights and MO coefficients - allocate(c(nBas,nspin),wEns(maxEns),occnum(nBas,nspin,maxEns)) + allocate(wEns(maxEns),occnum(nBas,nspin,maxEns)) call read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, & doNcentered,occnum,Cx_choice) @@ -140,73 +139,73 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n ! Compute GOK-RKS energy !------------------------------------------------------------------------ - if(method == 'GOK-RKS') then +! if(method == 'GOK-RKS') then - call cpu_time(start_KS) - call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, & - maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO(1),nV(1), & - S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum,Cx_choice) - call cpu_time(end_KS) +! call cpu_time(start_KS) +! call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, & +! maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO(1),nV(1), & +! S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum,Cx_choice) +! call cpu_time(end_KS) - t_KS = end_KS - start_KS - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GOK-RKS = ',t_KS,' seconds' - write(*,*) +! t_KS = end_KS - start_KS +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GOK-RKS = ',t_KS,' seconds' +! write(*,*) - end if +!end if !------------------------------------------------------------------------ ! Compute LIM excitation energies !------------------------------------------------------------------------ - if(method == 'LIM-RKS') then +! if(method == 'LIM-RKS') then - call cpu_time(start_KS) - call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), & - aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & - S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered) - call cpu_time(end_KS) +! call cpu_time(start_KS) +! call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), & +! aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & +! S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered) +! call cpu_time(end_KS) - t_KS = end_KS - start_KS - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for LIM-RKS = ',t_KS,' seconds' - write(*,*) +! t_KS = end_KS - start_KS +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for LIM-RKS = ',t_KS,' seconds' +! write(*,*) - end if +! end if !------------------------------------------------------------------------ ! Compute MOM excitation energies !------------------------------------------------------------------------ - if(method == 'MOM-RKS') then +! if(method == 'MOM-RKS') then - call cpu_time(start_KS) - call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), & - aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & - S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered) - call cpu_time(end_KS) +! call cpu_time(start_KS) +! call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), & +! aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & +! S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered) +! call cpu_time(end_KS) - t_KS = end_KS - start_KS - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM-RKS = ',t_KS,' seconds' - write(*,*) +! t_KS = end_KS - start_KS +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM-RKS = ',t_KS,' seconds' +! write(*,*) - end if +! end if !------------------------------------------------------------------------ ! Compute GOK-UKS energy (BROKEN) !------------------------------------------------------------------------ - if(method == 'GOK-UKS') then +! if(method == 'GOK-UKS') then - call cpu_time(start_KS) - call GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, & - nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum, & - Cx_choice,doNcentered) - call cpu_time(end_KS) +! call cpu_time(start_KS) +! call GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, & +! nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum, & +! Cx_choice,doNcentered) +! call cpu_time(end_KS) - t_KS = end_KS - start_KS - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UKS = ',t_KS,' seconds' - write(*,*) +! t_KS = end_KS - start_KS +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UKS = ',t_KS,' seconds' +! write(*,*) - end if +! end if !------------------------------------------------------------------------ ! Compute N-centered UKS energy @@ -216,7 +215,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n call cpu_time(start_KS) call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type,mix, & - nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS) + nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc) call cpu_time(end_KS) t_KS = end_KS - start_KS diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index 373e1ea..d0bf4e0 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -1,5 +1,5 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, & - nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eps,c,Pw) + nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eps,c,Pw,Vxc) ! Perform unrestricted Kohn-Sham calculation for ensembles @@ -80,7 +80,8 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig double precision,intent(out) :: Ew double precision,intent(out) :: eps(nBas,nspin) double precision,intent(out) :: Pw(nBas,nBas,nspin) - double precision,intent(out) :: c(nBAs,nBas,nspin) + double precision,intent(out) :: c(nBas,nBas,nspin) + double precision,intent(out) :: Vxc(nBas,nspin) ! Hello world @@ -386,12 +387,17 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig call print_UKS(nBas,nEns,occnum,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) +! Compute Vxc for post-HF calculations + + call xc_potential(nBas,c,Fx,Fc,Vxc) + !------------------------------------------------------------------------ ! Compute individual energies from ensemble energy !------------------------------------------------------------------------ - call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, & - AO,dAO,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum, & - Cx_choice,doNcentered) + if(nEns > 1) & + call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, & + AO,dAO,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum, & + Cx_choice,doNcentered) end subroutine eDFT_UKS diff --git a/src/eDFT/read_options_dft.f90 b/src/eDFT/read_options_dft.f90 index 77964a9..8c0fbc0 100644 --- a/src/eDFT/read_options_dft.f90 +++ b/src/eDFT/read_options_dft.f90 @@ -96,7 +96,7 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns, occnum(:,:,:) = 0d0 - do iEns=1,nEns + do iEns=1,maxEns read(1,*) read(1,*) (occnum(iBas,1,iEns),iBas=1,nBas) read(1,*) (occnum(iBas,2,iEns),iBas=1,nBas) @@ -117,9 +117,9 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns, ! Read ensemble weights for real physical (fractional number of electrons) ensemble (w1,w2) - allocate(nEl(nEns)) + allocate(nEl(maxEns)) nEl(:) = 0d0 - do iEns=1,nEns + do iEns=1,maxEns do iBas=1,nBas nEl(iEns) = nEl(iEns) + occnum(iBas,1,iEns) + occnum(iBas,2,iEns) end do diff --git a/src/eDFT/unrestricted_hybrid_correlation_energy.f90 b/src/eDFT/unrestricted_hybrid_correlation_energy.f90 index 0e030d0..4a3413e 100644 --- a/src/eDFT/unrestricted_hybrid_correlation_energy.f90 +++ b/src/eDFT/unrestricted_hybrid_correlation_energy.f90 @@ -31,7 +31,7 @@ subroutine unrestricted_hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho Ec(:) = 0d0 - case('LYP') + case('B3LYP') aC = 0.81d0 @@ -40,11 +40,13 @@ subroutine unrestricted_hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho Ec(:) = EcLDA(:) + aC*(EcGGA(:) - EcLDA(:)) + case('BHHLYP') + + call unrestricted_gga_correlation_energy('LYP ',nEns,wEns,nGrid,weight,rho,drho,Ec) + case('PBE') - call unrestricted_gga_correlation_energy('PBE ',nEns,wEns,nGrid,weight,rho,drho,EcGGA) - - Ec(:) = EcGGA(:) + call unrestricted_gga_correlation_energy('PBE ',nEns,wEns,nGrid,weight,rho,drho,Ec) case default diff --git a/src/eDFT/unrestricted_hybrid_correlation_potential.f90 b/src/eDFT/unrestricted_hybrid_correlation_potential.f90 index 41afc8f..96a9935 100644 --- a/src/eDFT/unrestricted_hybrid_correlation_potential.f90 +++ b/src/eDFT/unrestricted_hybrid_correlation_potential.f90 @@ -36,7 +36,7 @@ subroutine unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight, Fc(:,:,:) = 0d0 - case('LYP') + case('B3LYP') allocate(FcLDA(nBas,nBas,nspin),FcGGA(nBas,nBas,nspin)) @@ -47,13 +47,17 @@ subroutine unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight, Fc(:,:,:) = FcLDA(:,:,:) + aC*(FcGGA(:,:,:) - FcLDA(:,:,:)) + case('BHHLYP') + + allocate(FcGGA(nBas,nBas,nspin)) + + call unrestricted_gga_correlation_potential('LYP ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + case('PBE') allocate(FcGGA(nBas,nBas,nspin)) - call unrestricted_gga_correlation_potential('PBE ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FcGGA) - - Fc(:,:,:) = FcGGA(:,:,:) + call unrestricted_gga_correlation_potential('PBE ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) case default diff --git a/src/eDFT/unrestricted_hybrid_exchange_energy.f90 b/src/eDFT/unrestricted_hybrid_exchange_energy.f90 index d0365b4..8fee449 100644 --- a/src/eDFT/unrestricted_hybrid_exchange_energy.f90 +++ b/src/eDFT/unrestricted_hybrid_exchange_energy.f90 @@ -38,7 +38,7 @@ subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1 call unrestricted_fock_exchange_energy(nBas,P,FxHF,Ex) - case ('B3') + case ('B3LYP') a0 = 0.20d0 aX = 0.72d0 @@ -52,6 +52,13 @@ subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1 + a0*(ExHF - ExLDA) & + aX*(ExGGA - ExLDA) + case ('BHHLYP') + + call unrestricted_gga_exchange_energy('B88 ',nEns,wEns,nGrid,weight,rho,drho,ExGGA) + call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF) + + Ex = 0.5d0*ExHF + 0.5d0*ExGGA + case ('PBE') call unrestricted_gga_exchange_energy('PBE ',nEns,wEns,nGrid,weight,rho,drho,ExGGA) diff --git a/src/eDFT/unrestricted_hybrid_exchange_potential.f90 b/src/eDFT/unrestricted_hybrid_exchange_potential.f90 index f7c5d37..9875832 100644 --- a/src/eDFT/unrestricted_hybrid_exchange_potential.f90 +++ b/src/eDFT/unrestricted_hybrid_exchange_potential.f90 @@ -43,7 +43,7 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) Fx(:,:) = FxHF(:,:) - case('B3') + case('B3LYP') allocate(FxLDA(nBas,nBas),FxGGA(nBas,nBas)) @@ -59,6 +59,15 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC + a0*(FxHF(:,:) - FxLDA(:,:)) & + aX*(FxGGA(:,:) - FxLDA(:,:)) + case('BHHLYP') + + allocate(FxGGA(nBas,nBas)) + + call unrestricted_gga_exchange_potential('B88 ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) + call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) + + Fx(:,:) = 0.5d0*FxHF(:,:) + 0.5d0*FxGGA(:,:) + case('PBE') allocate(FxGGA(nBas,nBas)) diff --git a/src/eDFT/xc_potential.f90 b/src/eDFT/xc_potential.f90 new file mode 100644 index 0000000..2037ab2 --- /dev/null +++ b/src/eDFT/xc_potential.f90 @@ -0,0 +1,40 @@ +subroutine xc_potential(nBas,c,Fx,Fc,Vxc) + +! Compute the exchange-correlation potential in the MO basis + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + double precision,intent(in) :: c(nBas,nBas,nspin) + double precision,intent(in) :: Fx(nBas,nBas,nspin) + double precision,intent(in) :: Fc(nBas,nBas,nspin) + +! Local variables + + integer :: mu,nu + integer :: p + integer :: ispin + +! Output variables + + double precision,intent(out) :: Vxc(nBas,nspin) + +! Compute Vxc + + Vxc(:,:) = 0d0 + do p=1,nBas + do ispin=1,nspin + do mu=1,nBas + do nu=1,nBas + Vxc(p,ispin) = Vxc(p,ispin) & + + c(mu,p,ispin)*(Fx(mu,nu,ispin) + Fc(mu,nu,ispin))*c(nu,p,ispin) + + end do + end do + end do + end do + +end subroutine xc_potential