From 2c7653b7bc44525078f070178e2122e704181157 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 17 Jul 2023 13:35:24 +0200 Subject: [PATCH] add HF module --- input/methods | 8 +- input/options | 2 +- src/GT/G0T0eh.f90 | 11 +- src/GT/G0T0pp.f90 | 11 +- src/GT/GTeh_excitation_density.f90 | 20 ++- src/GT/GTeh_self_energy_diag.f90 | 28 ++-- src/GT/dynamic_Tmatrix_A.f90 | 16 +- src/GT/evGTeh.f90 | 5 +- src/GT/evGTpp.f90 | 8 +- src/GW/G0W0.f90 | 9 +- src/GW/evGW.f90 | 13 +- src/GW/qsGW.f90 | 6 +- src/HF/HF.f90 | 136 +++++++++++++++ src/HF/MOM.f90 | 195 ---------------------- src/HF/RMOM.f90 | 152 ++++++++--------- src/QuAcK/QuAcK.f90 | 254 ++++++++++------------------- src/QuAcK/read_methods.f90 | 24 +-- 17 files changed, 370 insertions(+), 528 deletions(-) create mode 100644 src/HF/HF.f90 delete mode 100644 src/HF/MOM.f90 diff --git a/input/methods b/input/methods index 9b32421..8848a96 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ -# RHF UHF KS MOM - T F F F +# RHF UHF RMOM UMOM KS + T F F F F # MP2* MP3 F F # CCD pCCD DCD CCSD CCSD(T) @@ -11,9 +11,9 @@ # phRPA* phRPAx* crRPA ppRPA F F F F # G0F2* evGF2* qsGF2* G0F3 evGF3 - T F F F F + F F F F F # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW F F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh - F F F T F F + F F F F F F # * unrestricted version available diff --git a/input/options b/input/options index 3fbc216..ce371d3 100644 --- a/input/options +++ b/input/options @@ -15,4 +15,4 @@ # ACFDT: AC Kx XBS F T T # BSE: phBSE phBSE2 ppBSE dBSE dTDA evDyn - T F F T F F + F F F T F F diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90 index fb37dfe..b5b5299 100644 --- a/src/GT/G0T0eh.f90 +++ b/src/GT/G0T0eh.f90 @@ -1,6 +1,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & - ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGT) + ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) ! Perform ehG0T0 calculation @@ -58,15 +58,14 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rhoL_RPA(:,:,:) - double precision,allocatable :: rhoR_RPA(:,:,:) + double precision,allocatable :: rhoL_RPA(:,:,:,:) + double precision,allocatable :: rhoR_RPA(:,:,:,:) + double precision,allocatable :: eGT(:) double precision,allocatable :: eGTlin(:) ! Output variables - double precision :: eGT(nBas) - ! Hello world write(*,*) @@ -100,7 +99,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD ! Memory allocation allocate(Sig(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & - rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS),eGTlin(nBas)) + rhoL_RPA(nBas,nBas,nS,2),rhoR_RPA(nBas,nBas,nS,2),eGT(nBas),eGTlin(nBas)) !-------------------! ! Compute screening ! diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index c0f29fa..a57ccb0 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -1,5 +1,5 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGT) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -62,12 +62,11 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp double precision,allocatable :: SigX(:) double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) + double precision,allocatable :: eGT(:) double precision,allocatable :: eGTlin(:) ! Output variables - double precision,intent(out) :: eGT(nBas) - ! Hello world write(*,*) @@ -88,11 +87,11 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp allocate(Om1ab(nVVab),X1ab(nVVab,nVVab),Y1ab(nOOab,nVVab), & Om2ab(nOOab),X2ab(nVVab,nOOab),Y2ab(nOOab,nOOab), & - rho1ab(nBas,nBas,nVVab),rho2ab(nBas,nBas,nOOab), & + rho1ab(nBas,nBas,nVVab),rho2ab(nBas,nBas,nOOab), & Om1aa(nVVaa),X1aa(nVVaa,nVVaa),Y1aa(nOOaa,nVVaa), & Om2aa(nOOaa),X2aa(nVVaa,nOOaa),Y2aa(nOOaa,nOOaa), & - rho1aa(nBas,nBas,nVVaa),rho2aa(nBas,nBas,nOOaa), & - SigX(nBas),Sig(nBas),Z(nBas),eGTlin(nBas)) + rho1aa(nBas,nBas,nVVaa),rho2aa(nBas,nBas,nOOaa), & + SigX(nBas),Sig(nBas),Z(nBas),eGT(nBas),eGTlin(nBas)) !---------------------------------------------- ! alpha-beta block diff --git a/src/GT/GTeh_excitation_density.f90 b/src/GT/GTeh_excitation_density.f90 index 40c415e..11da30e 100644 --- a/src/GT/GTeh_excitation_density.f90 +++ b/src/GT/GTeh_excitation_density.f90 @@ -14,17 +14,17 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR) ! Local variables integer :: m,jb,p,q,j,b - double precision :: X,Y + double precision :: X,Y,Xt,Yt ! Output variables - double precision,intent(out) :: rhoL(nBas,nBas,nS) - double precision,intent(out) :: rhoR(nBas,nBas,nS) + double precision,intent(out) :: rhoL(nBas,nBas,nS,2) + double precision,intent(out) :: rhoR(nBas,nBas,nS,2) ! Initialization - rhoL(:,:,:) = 0d0 - rhoR(:,:,:) = 0d0 + rhoL(:,:,:,:) = 0d0 + rhoR(:,:,:,:) = 0d0 !$OMP PARALLEL & !$OMP SHARED(nC,nBas,nR,nO,nS,rhoL,rhoR,ERI,XpY,XmY) & @@ -42,8 +42,14 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR) X = 0.5d0*(XpY(m,jb) + XmY(m,jb)) Y = 0.5d0*(XpY(m,jb) - XmY(m,jb)) - rhoL(p,q,m) = rhoL(p,q,m) + ERI(p,j,b,q)*X + ERI(p,b,j,q)*Y - rhoR(p,q,m) = rhoR(p,q,m) + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*X + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*Y + Xt = 0.5d0*(XpY(jb,m) + XmY(jb,m)) + Yt = 0.5d0*(XpY(jb,m) - XmY(jb,m)) + + rhoL(p,q,m,1) = rhoL(p,q,m,1) + ERI(p,j,b,q)*X + ERI(p,b,j,q)*Y + rhoL(p,q,m,2) = rhoL(p,q,m,2) + ERI(p,b,j,q)*X + ERI(p,j,b,q)*Y + + rhoR(p,q,m,1) = rhoR(p,q,m,1) + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*X + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*Y + rhoR(p,q,m,2) = rhoR(p,q,m,2) + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*X + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*Y enddo enddo diff --git a/src/GT/GTeh_self_energy_diag.f90 b/src/GT/GTeh_self_energy_diag.f90 index f48c06c..41df1a1 100644 --- a/src/GT/GTeh_self_energy_diag.f90 +++ b/src/GT/GTeh_self_energy_diag.f90 @@ -16,12 +16,12 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig integer,intent(in) :: nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rhoL(nBas,nBas,nS) - double precision,intent(in) :: rhoR(nBas,nBas,nS) + double precision,intent(in) :: rhoL(nBas,nBas,nS,2) + double precision,intent(in) :: rhoR(nBas,nBas,nS,2) ! Local variables - integer :: i,a,p,q,m + integer :: i,a,p,m double precision :: num,eps ! Output variables @@ -46,7 +46,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig do m=1,nS eps = e(p) - e(i) + Om(m) - num = rhoL(i,p,m)*rhoR(i,p,m) + num = rhoL(p,i,m,2)*rhoR(i,p,m,1) Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 @@ -61,7 +61,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig do m=1,nS eps = e(p) - e(a) - Om(m) - num = rhoL(p,a,m)*rhoR(p,a,m) + num = rhoL(p,a,m,1)*rhoR(a,p,m,2) Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 @@ -72,17 +72,17 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig ! GM correlation energy EcGM = 0d0 - do i=nC+1,nO - do a=nO+1,nBas-nR - do m=1,nS +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! do m=1,nS - eps = e(a) - e(i) + Om(m) - num = rhoL(i,a,m)*rhoR(i,a,m) - EcGM = EcGM - num*eps/(eps**2 + eta**2) +! eps = e(a) - e(i) + Om(m) +! num = rhoL(i,a,m)*rhoR(i,a,m) +! EcGM = EcGM - num*eps/(eps**2 + eta**2) - end do - end do - end do +! end do +! end do +! end do ! Compute renormalization factor from derivative diff --git a/src/GT/dynamic_Tmatrix_A.f90 b/src/GT/dynamic_Tmatrix_A.f90 index 65361d7..a5c8ec9 100644 --- a/src/GT/dynamic_Tmatrix_A.f90 +++ b/src/GT/dynamic_Tmatrix_A.f90 @@ -1,4 +1,4 @@ -subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,TA,ZA) +subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Om1,Om2,rho1,rho2,OmBSE,TA,ZA) ! Compute the dynamic part of the Bethe-Salpeter equation matrices for GT @@ -23,8 +23,8 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O double precision,intent(in) :: OmBSE - double precision,intent(in) :: Omega1(nVV) - double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: Om1(nVV) + double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: rho2(nBas,nBas,nOO) @@ -47,7 +47,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O ! Build dynamic A matrix jb = 0 -!$omp parallel do default(private) shared(TA,ZA,Omega1,Omega2,OmBSE,eGT,rho1,rho2,nO,nBas,nVV,nOO,chi,eps,eta,nC,nR,lambda) +!$omp parallel do default(private) shared(TA,ZA,Om1,Om2,OmBSE,eGT,rho1,rho2,nO,nBas,nVV,nOO,chi,eps,eta,nC,nR,lambda) do j=nC+1,nO do b=nO+1,nBas-nR jb = (b-nO) + (j-1)*(nBas-nO) @@ -60,12 +60,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O chi = 0d0 do cd=1,nVV - eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) + eps = + OmBSE - Om1(cd) + (eGT(i) + eGT(j)) chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2) end do do kl=1,nOO - eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) + eps = + OmBSE + Om2(kl) - (eGT(a) + eGT(b)) chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2) end do @@ -74,12 +74,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O chi = 0d0 do cd=1,nVV - eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) + eps = + OmBSE - Om1(cd) + (eGT(i) + eGT(j)) chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do kl=1,nOO - eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) + eps = + OmBSE + Om2(kl) - (eGT(a) + eGT(b)) chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do diff --git a/src/GT/evGTeh.f90 b/src/GT/evGTeh.f90 index 07ee8b1..b11ce7c 100644 --- a/src/GT/evGTeh.f90 +++ b/src/GT/evGTeh.f90 @@ -1,6 +1,6 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, & - cHF,eHF,Vxc,eG0T0) + cHF,eHF,Vxc) ! Perform self-consistent eigenvalue-only ehGT calculation @@ -41,7 +41,6 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: Vxc(nBas) - double precision,intent(in) :: eG0T0(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) @@ -121,7 +120,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, Conv = 1d0 e_diis(:,:) = 0d0 error_diis(:,:) = 0d0 - eGT(:) = eG0T0(:) + eGT(:) = eHF(:) eOld(:) = eGT(:) Z(:) = 1d0 rcond = 0d0 diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90 index e6a6cae..e05e59e 100644 --- a/src/GT/evGTpp.f90 +++ b/src/GT/evGTpp.f90 @@ -1,6 +1,5 @@ -subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & - BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nBas, & - nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0) +subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & + singlet,triplet,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) @@ -41,7 +40,6 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & 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) :: eG0T0(nBas) ! Local variables @@ -114,7 +112,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & Conv = 1d0 e_diis(:,:) = 0d0 error_diis(:,:) = 0d0 - eGT(:) = eG0T0(:) + eGT(:) = eHF(:) eOld(:) = eGT(:) Z(:) = 1d0 rcond = 0d0 diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 951c3fa..fdff910 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -1,6 +1,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & - ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGW) + ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) ! Perform G0W0 calculation @@ -61,6 +61,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD double precision,allocatable :: XmY_RPA(:,:) double precision,allocatable :: rho_RPA(:,:,:) + double precision,allocatable :: eGW(:) double precision,allocatable :: eGWlin(:) integer :: nBas2 @@ -70,12 +71,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD integer :: nR2 integer :: nS2 - double precision,allocatable :: seHF(:),seGW(:),sERI(:,:,:,:) - ! Output variables - double precision :: eGW(nBas) - ! Hello world write(*,*) @@ -115,7 +112,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD ! Memory allocation - allocate(SigC(nBas),SigX(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),eGW(nBas),eGWlin(nBas)) !-------------------! ! Compute screening ! diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 7d914ec..1547f0f 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -1,6 +1,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, & - cHF,eHF,Vxc,eG0W0) + cHF,eHF,Vxc) ! Perform self-consistent eigenvalue-only GW calculation @@ -42,7 +42,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: Vxc(nBas) - double precision,intent(in) :: eG0W0(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) @@ -137,7 +136,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, Conv = 1d0 e_diis(:,:) = 0d0 error_diis(:,:) = 0d0 - eGW(:) = eG0W0(:) + eGW(:) = eHF(:) eOld(:) = eGW(:) Z(:) = 1d0 rcond = 0d0 @@ -312,10 +311,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcppBSE(1) - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcppBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (singlet) =',EcppBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (triplet) =',3d0*EcppBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index d395d04..98abeb3 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -1,6 +1,6 @@ -subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) +subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn, & + singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO, & + ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GW calculation diff --git a/src/HF/HF.f90 b/src/HF/HF.f90 new file mode 100644 index 0000000..baf5a21 --- /dev/null +++ b/src/HF/HF.f90 @@ -0,0 +1,136 @@ +subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,guess_type,mix,level_shift, & + nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF) + +! Hartree-Fock module + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: doRHF + logical,intent(in) :: doUHF + logical,intent(in) :: doRMOM + logical,intent(in) :: doUMOM + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + integer,intent(in) :: guess_type + double precision,intent(in) :: thresh + double precision,intent(in) :: level_shift + logical,intent(in) :: mix + + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: nSCF + integer :: n_diis + double precision :: ET + double precision :: EV + double precision :: EJ + double precision :: EK + double precision :: dipole(ncart) + + double precision :: Conv + double precision :: Gap + double precision :: rcond + double precision,external :: trace_matrix + double precision,allocatable :: error(:,:) + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: F_diis(:,:) + double precision,allocatable :: J(:,:) + double precision,allocatable :: K(:,:) + double precision,allocatable :: cp(:,:) + double precision,allocatable :: Fp(:,:) + +! Output variables + + logical,intent(out) :: unrestricted + + double precision,intent(out) :: EHF + double precision,intent(out) :: epsHF(nBas) + double precision,intent(out) :: cHF(nBas,nBas) + double precision,intent(out) :: PHF(nBas,nBas) + double precision,intent(out) :: vHF(nBas) + double precision,intent(out) :: F(nBas,nBas) + +!------------------------------------------------------------------------ +! Compute RHF energy +!------------------------------------------------------------------------ + + if(doRHF) then + + ! Check that RHF calculation is worth doing... + + if(nO(1) /= nO(2)) then + write(*,*) ' !!! The system does not appear to be closed shell !!!' + write(*,*) + stop + end if + + call RHF(maxSCF,thresh,n_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF) + + end if + +!------------------------------------------------------------------------ +! Compute UHF energy +!------------------------------------------------------------------------ + + if(doUHF) then + + ! Switch on the unrestricted flag + unrestricted = .true. + + call UHF(maxSCF,thresh,n_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF) + + end if + +!------------------------------------------------------------------------ +! Restricted maximum overlap method +!------------------------------------------------------------------------ + + if(doRMOM) then + + ! Check that RMOM calculation is worth doing... + + if(nO(1) /= nO(2)) then + write(*,*) ' !!! The system does not appear to be closed shell !!!' + write(*,*) + stop + end if + +! call RMOM(maxSCF,thresh,n_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, & +! nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF) + + end if + +!------------------------------------------------------------------------ +! Unrestricted maximum overlap method +!------------------------------------------------------------------------ + + if(doUMOM) then + + ! Switch on the unrestricted flag + unrestricted = .true. + +! call UMOM(maxSCF,thresh,n_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, & +! nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF) + + end if + +end subroutine diff --git a/src/HF/MOM.f90 b/src/HF/MOM.f90 deleted file mode 100644 index 6de6630..0000000 --- a/src/HF/MOM.f90 +++ /dev/null @@ -1,195 +0,0 @@ -subroutine MOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P) - -! Maximum overlap method - - implicit none - -! Input variables - - integer,intent(in) :: maxSCF,max_diis - double precision,intent(in) :: thresh - - integer,intent(in) :: nBas,nO - double precision,intent(in) :: ENuc - double precision,intent(in) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),X(nBas,nBas) - -! Local variables - - integer :: iBas,jBas - integer :: nSCF,nBasSq,n_diis - double precision :: ET,EV,EJ,EK,Conv,Gap - double precision :: rcond - double precision,external :: trace_matrix - double precision,allocatable :: error(:,:),error_diis(:,:),F_diis(:,:) - double precision,allocatable :: J(:,:),K(:,:),cp(:,:),F(:,:),Fp(:,:) - double precision,allocatable :: cG(:,:),ON(:) - -! Output variables - - double precision,intent(inout):: ERHF,c(nBas,nBas),e(nBas),P(nBas,nBas) - -! Hello world - - write(*,*) - write(*,*)'************************************************' - write(*,*)'| Maximum overlap method |' - write(*,*)'************************************************' - write(*,*) - -! Useful quantities - - nBasSq = nBas*nBas - -! Memory allocation - - allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), & - cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas), & - cG(nBas,nBas),ON(nBas), & - error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) - -! Set up guess orbitals - - cG(:,:) = c(:,:) - -! Set up occupation numbers - - ON(1:nO) = 1d0 - ON(nO+1:nBas) = 0d0 - -! HOMO-LUMO transition - - ON(nO) = 0d0 - ON(nO+1) = 1d0 - - write(*,*) - write(*,*) ' --- Initial MO occupations --- ' - write(*,*) - call matout(nBas,1,ON) - write(*,*) - -! Compute density matrix - - call density_matrix(nBas,ON,c,P) - -! Initialization - - n_diis = 0 - F_diis(:,:) = 0d0 - error_diis(:,:) = 0d0 - Conv = 1d0 - nSCF = 0 - -!------------------------------------------------------------------------ -! Main SCF loop -!------------------------------------------------------------------------ - write(*,*) - write(*,*)'----------------------------------------------------' - write(*,*)'| MOM calculation |' - write(*,*)'----------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & - '|','#','|','HF energy','|','Conv','|','HL Gap','|' - write(*,*)'----------------------------------------------------' - - do while(Conv > thresh .and. nSCF < maxSCF) - -! Increment - - nSCF = nSCF + 1 - -! Build Fock matrix - - call Coulomb_matrix_AO_basis(nBas,P,ERI,J) - call exchange_matrix_AO_basis(nBas,P,ERI,K) - - F(:,:) = Hc(:,:) + J(:,:) + 0.5*K(:,:) - -! Check convergence - - error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) - Conv = maxval(abs(error)) - -! DIIS extrapolation - - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) - -! Reset DIIS if required - - if(abs(rcond) < 1d-15) n_diis = 0 - -! Diagonalize Fock matrix - - Fp = matmul(transpose(X),matmul(F,X)) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,e) - c = matmul(X,cp) - -! MOM overlap - - call MOM_overlap(nBas,nO,S,cG,c,ON) - -! Density matrix - - call density_matrix(nBas,ON,c,P) - -! Compute HF energy - - ERHF = trace_matrix(nBas,matmul(P,Hc)) & - + 0.5d0*trace_matrix(nBas,matmul(P,J)) & - + 0.5d0*trace_matrix(nBas,matmul(P,K)) - -! Compute HOMO-LUMO gap - - if(nBas > nO) then - - Gap = e(nO+1) - e(nO) - - else - - Gap = 0d0 - - endif - -! Dump results - - write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',ERHF+ENuc,'|',Conv,'|',Gap,'|' - - enddo - write(*,*)'----------------------------------------------------' -!------------------------------------------------------------------------ -! End of SCF loop -!------------------------------------------------------------------------ - -! Did it actually converge? - - if(nSCF == maxSCF) then - - write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*) - - stop - - endif - - write(*,*) - write(*,*) ' --- Final MO occupations --- ' - write(*,*) - call matout(nBas,1,ON) - write(*,*) - -! Compute HF energy - - ET = trace_matrix(nBas,matmul(P,T)) - EV = trace_matrix(nBas,matmul(P,V)) - EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) - EK = 0.5d0*trace_matrix(nBas,matmul(P,K)) - ERHF = ET + EV + EJ + EK - - call print_RHF(nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF) - -end subroutine MOM diff --git a/src/HF/RMOM.f90 b/src/HF/RMOM.f90 index c086b36..0106bf9 100644 --- a/src/HF/RMOM.f90 +++ b/src/HF/RMOM.f90 @@ -1,68 +1,40 @@ -subroutine RMOM(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) +subroutine RMOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P) -! Perform restricted Hartree-Fock calculation with MOM algorithm +! Maximum overlap method implicit none - include 'parameters.h' ! Input variables - integer,intent(in) :: maxSCF,max_diis,guess_type + integer,intent(in) :: maxSCF,max_diis double precision,intent(in) :: thresh - integer,intent(in) :: nBas - integer,intent(in) :: nO - integer,intent(in) :: nNuc - double precision,intent(in) :: ZNuc(nNuc) - double precision,intent(in) :: rNuc(nNuc,ncart) + integer,intent(in) :: nBas,nO double precision,intent(in) :: ENuc - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas) - double precision,intent(in) :: V(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),X(nBas,nBas) ! Local variables - integer :: nSCF - integer :: nBasSq - integer :: n_diis - double precision :: ET - double precision :: EV - double precision :: EJ - double precision :: EK - double precision :: dipole(ncart) - - double precision :: Conv - double precision :: Gap + integer :: iBas,jBas + integer :: nSCF,nBasSq,n_diis + double precision :: ET,EV,EJ,EK,Conv,Gap double precision :: rcond double precision,external :: trace_matrix - double precision,allocatable :: error(:,:) - double precision,allocatable :: error_diis(:,:) - double precision,allocatable :: F_diis(:,:) - double precision,allocatable :: J(:,:) - double precision,allocatable :: K(:,:) - double precision,allocatable :: cp(:,:) - double precision,allocatable :: F(:,:) - double precision,allocatable :: Fp(:,:) - double precision,allocatable :: ON(:) + double precision,allocatable :: error(:,:),error_diis(:,:),F_diis(:,:) + double precision,allocatable :: J(:,:),K(:,:),cp(:,:),F(:,:),Fp(:,:) + double precision,allocatable :: cG(:,:),ON(:) ! Output variables - double precision,intent(out) :: ERHF - 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) + double precision,intent(inout):: ERHF,c(nBas,nBas),e(nBas),P(nBas,nBas) ! Hello world write(*,*) - write(*,*)'*********************************************' - write(*,*)'| Restricted Maximum Overlap Method |' - write(*,*)'*********************************************' + write(*,*)'************************************************' + write(*,*)'| Maximum overlap method |' + write(*,*)'************************************************' write(*,*) ! Useful quantities @@ -71,22 +43,35 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S, ! Memory allocation - allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), & - cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas),ON(nBas), & + allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), & + cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas), & + cG(nBas,nBas),ON(nBas), & error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) -! Guess coefficients and eigenvalues - - call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) - -! ON(:) = 0d0 -! do i=1,nO -! ON(i) = 1d0 -! ON(i) = dble(2*i-1) -! end do - -! call density_matrix(nBas,ON,c,P) +! Set up guess orbitals + cG(:,:) = c(:,:) + +! Set up occupation numbers + + ON(1:nO) = 1d0 + ON(nO+1:nBas) = 0d0 + +! HOMO-LUMO transition + + ON(nO) = 0d0 + ON(nO+1) = 1d0 + + write(*,*) + write(*,*) ' --- Initial MO occupations --- ' + write(*,*) + call matout(nBas,1,ON) + write(*,*) + +! Compute density matrix + + call density_matrix(nBas,ON,c,P) + ! Initialization n_diis = 0 @@ -100,7 +85,7 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S, !------------------------------------------------------------------------ write(*,*) write(*,*)'----------------------------------------------------' - write(*,*)'| RHF calculation |' + write(*,*)'| MOM calculation |' write(*,*)'----------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & '|','#','|','HF energy','|','Conv','|','HL Gap','|' @@ -117,7 +102,7 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S, call Coulomb_matrix_AO_basis(nBas,P,ERI,J) call exchange_matrix_AO_basis(nBas,P,ERI,K) - F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + F(:,:) = Hc(:,:) + J(:,:) + 0.5*K(:,:) ! Check convergence @@ -127,11 +112,11 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S, ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - if(abs(rcond) > 1d-7) then - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) - else - n_diis = 0 - end if + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + +! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 ! Diagonalize Fock matrix @@ -140,17 +125,19 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S, call diagonalize_matrix(nBas,cp,e) c = matmul(X,cp) +! MOM overlap + + call MOM_overlap(nBas,nO,S,cG,c,ON) + ! Density matrix - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + call density_matrix(nBas,ON,c,P) -! call density_matrix(nBas,ON,c,P) - ! Compute HF energy - ERHF = trace_matrix(nBas,matmul(P,Hc)) & - + 0.5d0*trace_matrix(nBas,matmul(P,J)) & - + 0.25d0*trace_matrix(nBas,matmul(P,K)) + ERHF = trace_matrix(nBas,matmul(P,Hc)) & + + 0.5d0*trace_matrix(nBas,matmul(P,J)) & + + 0.5d0*trace_matrix(nBas,matmul(P,K)) ! Compute HOMO-LUMO gap @@ -189,21 +176,20 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S, endif + write(*,*) + write(*,*) ' --- Final MO occupations --- ' + write(*,*) + call matout(nBas,1,ON) + write(*,*) + ! Compute HF energy - ET = trace_matrix(nBas,matmul(P,T)) - EV = trace_matrix(nBas,matmul(P,V)) - EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) - EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) + ET = trace_matrix(nBas,matmul(P,T)) + EV = trace_matrix(nBas,matmul(P,V)) + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + EK = 0.5d0*trace_matrix(nBas,matmul(P,K)) ERHF = ET + EV + EJ + EK -! Compute dipole moments + call print_RHF(nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF) - 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 RMOM +end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 8c4d86c..9112742 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -3,9 +3,8 @@ program QuAcK implicit none include 'parameters.h' - logical :: doSph logical :: unrestricted = .false. - logical :: doRHF,doUHF,doMOM + logical :: doHF,doRHF,doUHF,doRMOM,doUMOM logical :: dostab logical :: doKS logical :: doMP2,doMP3 @@ -32,9 +31,6 @@ program QuAcK double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:) double precision,allocatable :: Vxc(:,:) - double precision,allocatable :: eG0W0(:,:) - double precision,allocatable :: eG0T0(:,:) - logical :: doACFDT logical :: exchange_kernel logical :: doXBS @@ -122,14 +118,12 @@ program QuAcK write(*,*) ! Spherium calculation? - - doSph = .false. call wall_time(start_QuAcK) ! Which calculations do you want to do? - call read_methods(doRHF,doUHF,doKS,doMOM, & + call read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, & doMP2,doMP3, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_crCCD,do_lCCD, & @@ -168,7 +162,7 @@ program QuAcK ! nS = number of single excitation ! = nO*nV - call read_molecule(nNuc,nEl(:),nO(:),nC(:),nR(:)) + call read_molecule(nNuc,nEl,nO,nC,nR) allocate(ZNuc(nNuc),rNuc(nNuc,ncart)) ! Read geometry @@ -193,24 +187,16 @@ program QuAcK ! Memory allocation for one- and two-electron integrals - 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),Vxc(nBas,nspin),F_AO(nBas,nBas)) + allocate(cHF(nBas,nBas,nspin),eHF(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),Vxc(nBas,nspin),F_AO(nBas,nBas)) ! Read integrals call wall_time(start_int) - if(doSph) then - - call read_integrals_sph(nBas,S,T,V,Hc,ERI_AO) - - else - - call read_integrals(nBas,S,T,V,Hc,ERI_AO) - call read_dipole_integrals(nBas,dipole_int_AO) - - end if + call read_integrals(nBas,S,T,V,Hc,ERI_AO) + call read_dipole_integrals(nBas,dipole_int_AO) call wall_time(end_int) @@ -224,52 +210,27 @@ program QuAcK call orthogonalization_matrix(ortho_type,nBas,S,X) !------------------------------------------------------------------------ -! Compute RHF energy +! Hartree-Fock module !------------------------------------------------------------------------ - if(doRHF) then + doHF = doRHF .or. doUHF .or. doRMOM .or. doUMOM - ! Check that RHF calculation is worth doing... - - if(nO(1) /= nO(2)) then - write(*,*) ' !!! The system does not appear to be closed shell !!!' - write(*,*) - stop - end if + if(doHF) then call wall_time(start_HF) - call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,F_AO,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc) + call HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF_HF,thresh_HF,n_diis_HF, & + guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F_AO, & + ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc) call wall_time(end_HF) t_HF = end_HF - start_HF - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for HF = ',t_HF,' seconds' write(*,*) end if !------------------------------------------------------------------------ -! Compute UHF energy -!------------------------------------------------------------------------ - - if(doUHF) then - - ! Switch on the unrestricted flag - unrestricted = .true. - - call cpu_time(start_HF) - call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - 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 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UHF = ',t_HF,' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Compute KS energy +! KS module !------------------------------------------------------------------------ if(doKS) then @@ -290,33 +251,6 @@ program QuAcK end if -!------------------------------------------------------------------------ -! Maximum overlap method -!------------------------------------------------------------------------ - - if(doMOM) then - - call cpu_time(start_HF) - - if(unrestricted) then - -! call UMOM() - - else - - call MOM(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,Vxc) - - end if - - call cpu_time(end_HF) - - t_HF = end_HF - start_HF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM = ',t_HF,' seconds' - write(*,*) - - end if - !------------------------------------------------------------------------ ! AO to MO integral transform for post-HF methods !------------------------------------------------------------------------ @@ -327,81 +261,70 @@ program QuAcK write(*,*) 'AO to MO transformation... Please be patient' write(*,*) - if(doSph) then + if(unrestricted) then - allocate(ERI_MO(nBas,nBas,nBas,nBas)) - ERI_MO(:,:,:,:) = ERI_AO(:,:,:,:) - print*,'!!! MO = AO !!!' - deallocate(ERI_AO) + ! Read and transform dipole-related integrals + + allocate(dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart)) + dipole_int_aa(:,:,:) = dipole_int_AO(:,:,:) + dipole_int_bb(:,:,:) = dipole_int_AO(:,:,:) + do ixyz=1,ncart + call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz)) + call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz)) + end do + + ! Memory allocation + + allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas)) + + ! 4-index transform for (aa|aa) block + + bra1 = 1 + bra2 = 1 + ket1 = 1 + ket2 = 1 + call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aaaa) + + ! 4-index transform for (aa|bb) block + + bra1 = 1 + bra2 = 1 + ket1 = 2 + ket2 = 2 + call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aabb) + + ! 4-index transform for (bb|bb) block + + bra1 = 2 + bra2 = 2 + ket1 = 2 + ket2 = 2 + call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_bbbb) else - if(unrestricted) then + ! Memory allocation + + allocate(ERI_MO(nBas,nBas,nBas,nBas)) + allocate(F_MO(nBas,nBas)) - ! Read and transform dipole-related integrals - - allocate(dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart)) - dipole_int_aa(:,:,:) = dipole_int_AO(:,:,:) - dipole_int_bb(:,:,:) = dipole_int_AO(:,:,:) - do ixyz=1,ncart - call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz)) - call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz)) - end do + ! Read and transform dipole-related integrals + + dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:) + do ixyz=1,ncart + call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz)) + end do - ! Memory allocation - - allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas)) - - ! 4-index transform for (aa|aa) block - - bra1 = 1 - bra2 = 1 - ket1 = 1 - ket2 = 1 - call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aaaa) - - ! 4-index transform for (aa|bb) block - - bra1 = 1 - bra2 = 1 - ket1 = 2 - ket2 = 2 - call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aabb) - - ! 4-index transform for (bb|bb) block - - bra1 = 2 - bra2 = 2 - ket1 = 2 - ket2 = 2 - call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_bbbb) + ! 4-index transform + + bra1 = 1 + bra2 = 1 + ket1 = 1 + ket2 = 1 + call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO) - else - - ! Memory allocation - - allocate(ERI_MO(nBas,nBas,nBas,nBas)) - allocate(F_MO(nBas,nBas)) - - ! Read and transform dipole-related integrals - - dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:) - do ixyz=1,ncart - call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz)) - end do - - ! 4-index transform - - bra1 = 1 - bra2 = 1 - ket1 = 1 - ket2 = 1 - call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO) - - F_MO(:,:) = F_AO(:,:) - call AOtoMO_transform(nBas,cHF,F_MO) - - end if + F_MO(:,:) = F_AO(:,:) + call AOtoMO_transform(nBas,cHF,F_MO) end if @@ -917,8 +840,6 @@ program QuAcK ! Perform G0W0 calculatiom !------------------------------------------------------------------------ - eG0W0(:,:) = eHF(:,:) - if(doG0W0) then call cpu_time(start_GW) @@ -926,11 +847,11 @@ program QuAcK call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,dophBSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & linGW,eta_GW,regGW,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) + dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc) else call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet, & - linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) + linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc) end if call cpu_time(end_GW) @@ -953,13 +874,13 @@ program QuAcK call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,dophBSE,TDA_W,TDA, & dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,regGW,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) + PHF,cHF,eHF,Vxc) else call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX, & dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet,linGW,eta_GW,regGW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc) end if call cpu_time(end_GW) @@ -1042,7 +963,7 @@ program QuAcK write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_GW,' seconds' write(*,*) - if(dophBSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) +! if(dophBSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) end if @@ -1061,7 +982,7 @@ program QuAcK write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_GW,' seconds' write(*,*) - if(dophBSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) +! if(dophBSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) end if @@ -1069,8 +990,6 @@ program QuAcK ! Perform G0T0pp calculatiom !------------------------------------------------------------------------ - eG0T0(:,:) = eHF(:,:) - if(doG0T0pp) then call cpu_time(start_GT) @@ -1081,14 +1000,13 @@ program QuAcK call UG0T0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn, & spin_conserved,spin_flip,linGT,eta_GT,regGT,nBas,nC,nO,nV, & nR,nS,ENuc,EUHF,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & - dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0T0) + dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc) else ! call soG0T0(eta_GT,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet, & - linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO, & - PHF,cHF,eHF,Vxc,eG0T0) + linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc) end if @@ -1114,14 +1032,14 @@ program QuAcK dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,& eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_AO, & ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa, & - dipole_int_bb,PHF,cHF,eHF,Vxc,eG0T0) + dipole_int_bb,PHF,cHF,eHF,Vxc) else call evGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO, & - PHF,cHF,eHF,Vxc,eG0T0) + PHF,cHF,eHF,Vxc) end if @@ -1168,8 +1086,6 @@ program QuAcK ! Perform G0T0eh calculatiom !------------------------------------------------------------------------ - eG0T0(:,:) = eHF(:,:) - if(doG0T0eh) then call cpu_time(start_GT) @@ -1181,7 +1097,7 @@ program QuAcK else call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet, & - linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0T0) + linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc) end if @@ -1206,7 +1122,7 @@ program QuAcK call evGTeh(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet,linGT,eta_GT,regGT, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0T0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc) end if call cpu_time(end_GT) diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 105a02f..e36695d 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -1,4 +1,4 @@ -subroutine read_methods(doRHF,doUHF,doKS,doMOM, & +subroutine read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, & doMP2,doMP3, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_crCCD,do_lCCD, & @@ -17,7 +17,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Input variables - logical,intent(out) :: doRHF,doUHF,doKS,doMOM + logical,intent(out) :: doRHF,doUHF,doRMOM,doUMOM,doKS logical,intent(out) :: doMP2,doMP3 logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD @@ -38,10 +38,11 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Set all the booleans to false - doRHF = .false. - doUHF = .false. - doKS = .false. - doMOM = .false. + doRHF = .false. + doUHF = .false. + doRMOM = .false. + doUMOM = .false. + doKS = .false. doMP2 = .false. doMP3 = .false. @@ -91,11 +92,12 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Read mean-field methods read(1,*) - read(1,*) answer1,answer2,answer3,answer4 - if(answer1 == 'T') doRHF = .true. - if(answer2 == 'T') doUHF = .true. - if(answer3 == 'T') doKS = .true. - if(answer4 == 'T') doMOM = .true. + read(1,*) answer1,answer2,answer3,answer4,answer5 + if(answer1 == 'T') doRHF = .true. + if(answer2 == 'T') doUHF = .true. + if(answer3 == 'T') doRMOM = .true. + if(answer4 == 'T') doUMOM = .true. + if(answer5 == 'T') doKS = .true. ! Read MPn methods