From 76a3b9cfffa19eb1bf5ec84bef9a09afb46c32ef Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 29 Nov 2022 12:11:09 +0100 Subject: [PATCH] SRG qsGW --- input/methods | 8 +- input/options | 4 +- src/GW/SRG_qsGW.f90 | 323 +++++++++++++++++++++++++ src/GW/renormalization_factor_SRG.f90 | 38 +++ src/GW/self_energy_correlation_SRG.f90 | 80 ++++++ src/QuAcK/QuAcK.f90 | 32 ++- src/QuAcK/read_methods.f90 | 30 +-- 7 files changed, 493 insertions(+), 22 deletions(-) create mode 100644 src/GW/SRG_qsGW.f90 create mode 100644 src/GW/renormalization_factor_SRG.f90 create mode 100644 src/GW/self_energy_correlation_SRG.f90 diff --git a/input/methods b/input/methods index 36fcf41..c4629c5 100644 --- a/input/methods +++ b/input/methods @@ -7,13 +7,13 @@ # drCCD rCCD crCCD lCCD F F F F # CIS* CIS(D) CID CISD FCI - T F F F F + F F F F F # RPA* RPAx* crRPA ppRPA - F T F F + F F F F # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F -# G0W0* evGW* qsGW* ufG0W0 ufGW - T F F F F +# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW + F F F T F F # G0T0 evGT qsGT F F F # * unrestricted version available diff --git a/input/options b/input/options index 9683845..0f1eb47 100644 --- a/input/options +++ b/input/options @@ -9,10 +9,10 @@ # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 3 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg - 256 0.00001 T 5 T 0.0 F F F F F F + 256 0.0000001 F 5 T 0.0 F F F F F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg 10 0.00001 T 5 T 0.0 F F # ACFDT: AC Kx XBS F T T # BSE: BSE dBSE dTDA evDyn ppBSE BSE2 - T F T F F T + F F T F F T diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 new file mode 100644 index 0000000..1607859 --- /dev/null +++ b/src/GW/SRG_qsGW.f90 @@ -0,0 +1,323 @@ +subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_W,TDA, & + dBSE,dTDA,evDyn,singlet,triplet,eta,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 + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: BSE + logical,intent(in) :: BSE2 + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: evDyn + logical,intent(in) :: singlet + logical,intent(in) :: triplet + double precision,intent(in) :: eta + + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: PHF(nBas,nBas) + 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_AO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart) + +! Local variables + + integer :: nSCF + integer :: nBasSq + integer :: ispin + integer :: n_diis + double precision :: ET + double precision :: EV + double precision :: EJ + double precision :: Ex + double precision :: EqsGW + double precision :: EcRPA + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision :: EcGM + double precision :: Conv + double precision :: rcond + double precision,external :: trace_matrix + double precision :: dipole(ncart) + + logical :: print_W = .false. + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: F_diis(:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + double precision,allocatable :: c(:,:) + double precision,allocatable :: cp(:,:) + double precision,allocatable :: eGW(:) + double precision,allocatable :: eOld(:) + double precision,allocatable :: P(:,:) + double precision,allocatable :: F(:,:) + double precision,allocatable :: Fp(:,:) + double precision,allocatable :: J(:,:) + double precision,allocatable :: K(:,:) + double precision,allocatable :: SigC(:,:) + double precision,allocatable :: Z(:) + double precision,allocatable :: error(:,:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Self-consistent SRG-qsGW calculation |' + write(*,*)'************************************************' + write(*,*) + +! Warning + + write(*,*) '!! ERIs in MO basis will be overwritten in qsGW !!' + write(*,*) + +! Stuff + + nBasSq = nBas*nBas + +! TDA for W + + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Memory allocation + + allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & + J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & + rho_RPA(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + +! Initialization + + nSCF = -1 + n_diis = 0 + ispin = 1 + Conv = 1d0 + P(:,:) = PHF(:,:) + eGW(:) = eHF(:) + eOld(:) = eHF(:) + c(:,:) = cHF(:,:) + F_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + rcond = 0d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF <= maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Buid Coulomb matrix + + call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J) + + ! Compute exchange part of the self-energy + + call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) + + ! AO to MO transformation of two-electron integrals + + call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO) + + ! Compute linear response + + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) + if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,OmRPA) + + ! Compute correlation part of the self-energy + + call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA) + call self_energy_correlation_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + + ! Make correlation self-energy Hermitian and transform it back to AO basis + + call MOtoAO_transform(nBas,S,c,SigC) + + ! Solve the quasi-particle equation + + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigC(:,:) + + ! Compute commutator and convergence criteria + + error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) + + ! 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 + + ! Diagonalize Hamiltonian in AO basis + + Fp = matmul(transpose(X),matmul(F,X)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nBas,cp,eGW) + c = matmul(X,cp) + SigC = matmul(transpose(c),matmul(SigC,c)) + + ! Compute new density matrix in the AO basis + + P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + + ! Save quasiparticles energy for next cycle + + Conv = maxval(abs(eGW - eOld)) + eOld(:) = eGW(:) + + !------------------------------------------------------------------------ + ! Compute total energy + !------------------------------------------------------------------------ + + ! Kinetic energy + + ET = trace_matrix(nBas,matmul(P,T)) + + ! Potential energy + + EV = trace_matrix(nBas,matmul(P,V)) + + ! Coulomb energy + + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + + ! Exchange energy + + Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) + + ! Total energy + + EqsGW = ET + EV + EJ + Ex + + ! Print results + + call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole) + + enddo +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF+1) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +! Deallocate memory + + deallocate(c,cp,P,F,Fp,J,K,SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis) + +! Perform BSE calculation + + if(BSE) then + + call Bethe_Salpeter(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & + eGW,eGW,EcBSE) + + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(2) + + end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + + if(doACFDT) then + + write(*,*) '------------------------------------------------------' + write(*,*) 'Adiabatic connection version of BSE correlation energy' + write(*,*) '------------------------------------------------------' + write(*,*) + + if(doXBS) then + + write(*,*) '*** scaled screening version (XBS) ***' + write(*,*) + + end if + + 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(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcAC(1) + EcAC(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + + end if + +end subroutine SRG_qsGW diff --git a/src/GW/renormalization_factor_SRG.f90 b/src/GW/renormalization_factor_SRG.f90 new file mode 100644 index 0000000..f292b18 --- /dev/null +++ b/src/GW/renormalization_factor_SRG.f90 @@ -0,0 +1,38 @@ +subroutine renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) + +! Compute renormalization factor for GW + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: p,i,a,jb + double precision :: eps + +! Output variables + + double precision,intent(out) :: Z(nBas) + +! Initialize + + Z(:) = 0d0 + +! Compute renormalization factor from derivative of SigC + + Z(:) = 1d0/(1d0 - Z(:)) + +end subroutine renormalization_factor_SRG diff --git a/src/GW/self_energy_correlation_SRG.f90 b/src/GW/self_energy_correlation_SRG.f90 new file mode 100644 index 0000000..3235ac5 --- /dev/null +++ b/src/GW/self_energy_correlation_SRG.f90 @@ -0,0 +1,80 @@ +subroutine self_energy_correlation_SRG(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) + +! Compute correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,j,a,b + integer :: p,q,r + integer :: m + double precision :: Dpim,Dqim,Dpam,Dqam + +! Output variables + + double precision,intent(out) :: EcGM + double precision,intent(out) :: SigC(nBas,nBas) + +! Initialize + + SigC(:,:) = 0d0 + +!--------------------! +! SRG-GW self-energy ! +!--------------------! + +! Occupied part of the correlation self-energy + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do i=nC+1,nO + do m=1,nS + Dpim = e(p) - e(i) + Omega(m) + Dqim = e(q) - e(i) + Omega(m) + SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,m)*rho(q,i,m)*(Dpim + Dqim)/(Dpim**2 + Dqim**2) + end do + end do + end do + end do + +! Virtual part of the correlation self-energy + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do a=nO+1,nBas-nR + do m=1,nS + Dpam = e(p) - e(a) - Omega(m) + Dqam = e(q) - e(a) - Omega(m) + SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,m)*rho(q,a,m)*(Dpam + Dqam)/(Dpam**2 + Dqam**2) + end do + end do + end do + end do + +! Galitskii-Migdal correlation energy + + EcGM = 0d0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! do m=1,nS +! EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) +! end do +! end do +! end do + +end subroutine self_energy_correlation_SRG diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 9ddfc32..6700edf 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -15,7 +15,7 @@ program QuAcK logical :: doRPA,doRPAx,docrRPA,doppRPA logical :: doADC logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 - logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW + logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW logical :: doG0T0,doevGT,doqsGT integer :: nNuc,nBas,nBasCABS @@ -167,7 +167,7 @@ program QuAcK doRPA,doRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2, & doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW, & + doG0W0,doevGW,doqsGW,doSRGqsGW, & doufG0W0,doufGW, & doG0T0,doevGT,doqsGT) @@ -1080,6 +1080,34 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Perform SRG-qsGW calculation +!------------------------------------------------------------------------ + + if(doSRGqsGW) then + + call cpu_time(start_qsGW) + + if(unrestricted) then + + print*,'Unrestricted version of SRG-qsGW NYI' + + else + + call SRG_qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn, & + singlet,triplet,eta_GW,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) + + end if + + call cpu_time(end_qsGW) + + t_qsGW = end_qsGW - start_qsGW + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_qsGW,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Perform ufG0W0 calculatiom !------------------------------------------------------------------------ diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index d025073..c46cb19 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -6,7 +6,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doRPA,doRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2, & doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW, & + doG0W0,doevGW,doqsGW,doSRGqsGW, & doufG0W0,doufGW, & doG0T0,doevGT,doqsGT) @@ -23,12 +23,12 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD,doFCI logical,intent(out) :: doRPA,doRPAx,docrRPA,doppRPA logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 - logical,intent(out) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW + logical,intent(out) :: doG0W0,doevGW,doqsGW,doSRGqsGW,doufG0W0,doufGW logical,intent(out) :: doG0T0,doevGT,doqsGT ! Local variables - character(len=1) :: answer1,answer2,answer3,answer4,answer5 + character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6 ! Open file with method specification @@ -73,11 +73,12 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doG0F3 = .false. doevGF3 = .false. - doG0W0 = .false. - doevGW = .false. - doqsGW = .false. - doufG0W0 = .false. - doufGW = .false. + doG0W0 = .false. + doevGW = .false. + doqsGW = .false. + doSRGqsGW = .false. + doufG0W0 = .false. + doufGW = .false. doG0T0 = .false. doevGT = .false. @@ -149,12 +150,13 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Read GW methods read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5 - if(answer1 == 'T') doG0W0 = .true. - if(answer2 == 'T') doevGW = .true. - if(answer3 == 'T') doqsGW = .true. - if(answer4 == 'T') doufG0W0 = .true. - if(answer5 == 'T') doufGW = .true. + read(1,*) answer1,answer2,answer3,answer4,answer5,answer6 + if(answer1 == 'T') doG0W0 = .true. + if(answer2 == 'T') doevGW = .true. + if(answer3 == 'T') doqsGW = .true. + if(answer4 == 'T') doSRGqsGW = .true. + if(answer5 == 'T') doufG0W0 = .true. + if(answer6 == 'T') doufGW = .true. ! Read GT methods