From df5345e570d24f0e35fb52e85f59e20404e06b87 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 26 Jan 2022 13:05:16 +0100 Subject: [PATCH] spinorbital BSE@G0T0 --- input/methods | 4 +- input/options | 6 +- mol/h2.xyz | 2 +- src/GT/Bethe_Salpeter_Tmatrix_so.f90 | 79 ++++++++++++ src/GT/G0T0.f90 | 4 + src/GT/excitation_density_Tmatrix_so.f90 | 80 ++++++++++++ src/GT/renormalization_factor_Tmatrix_so.f90 | 63 ++++++++++ src/GT/self_energy_Tmatrix_diag_so.f90 | 63 ++++++++++ src/GT/soG0T0.f90 | 124 +++++++++++++++++++ src/GW/self_energy_exchange_diag.f90 | 3 +- src/HF/RHF.f90 | 2 +- src/HF/UHF.f90 | 2 +- src/HF/mo_fock_exchange_potential.f90 | 17 ++- src/QuAcK/QuAcK.f90 | 1 + 14 files changed, 436 insertions(+), 14 deletions(-) create mode 100644 src/GT/Bethe_Salpeter_Tmatrix_so.f90 create mode 100644 src/GT/excitation_density_Tmatrix_so.f90 create mode 100644 src/GT/renormalization_factor_Tmatrix_so.f90 create mode 100644 src/GT/self_energy_Tmatrix_diag_so.f90 create mode 100644 src/GT/soG0T0.f90 diff --git a/input/methods b/input/methods index 979b54d..07cc57d 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 pCCD DCD CCSD CCSD(T) @@ -15,7 +15,7 @@ # G0W0* evGW* qsGW* ufG0W0 ufGW F F F F F # G0T0 evGT qsGT - F F F + T F F # MCMP2 F # * unrestricted version available diff --git a/input/options b/input/options index bcd3046..bac7e06 100644 --- a/input/options +++ b/input/options @@ -7,14 +7,14 @@ # spin: TDA singlet triplet spin_conserved spin_flip F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm reg - 256 0.00001 T 5 T 0.0 3 F + 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.00001 T 5 T 0.0 F F F F F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg 256 0.00001 T 5 T 0.0 F F # ACFDT: AC Kx XBS F F F # BSE: BSE dBSE dTDA evDyn - T T T F + T F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/mol/h2.xyz b/mol/h2.xyz index a4e936a..3c8e04d 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 2.000000 +H 0. 0. 1.5 diff --git a/src/GT/Bethe_Salpeter_Tmatrix_so.f90 b/src/GT/Bethe_Salpeter_Tmatrix_so.f90 new file mode 100644 index 0000000..ded548e --- /dev/null +++ b/src/GT/Bethe_Salpeter_Tmatrix_so.f90 @@ -0,0 +1,79 @@ +subroutine Bethe_Salpeter_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,Omega1,X1,Y1,Omega2,X2,Y2,rho1,rho2, & + ERI,eT,eGT,EcBSE) + +! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel + + 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 + + integer,intent(in) :: nOO + integer,intent(in) :: nVV + + double precision,intent(in) :: eT(nBas) + double precision,intent(in) :: eGT(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: ispin + + double precision :: EcRPA + double precision,allocatable :: TA(:,:),TB(:,:) + double precision,allocatable :: OmBSE(:) + double precision,allocatable :: XpY_BSE(:,:) + double precision,allocatable :: XmY_BSE(:,:) + +! Output variables + + double precision,intent(out) :: EcBSE + +! Memory allocation + + allocate(TA(nS,nS),TB(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) + +!------------------! +! Compute T-matrix ! +!------------------! + + ispin = 4 + + call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eT,ERI, & + Omega1,X1,Y1,Omega2,X2,Y2,EcRPA) + + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,Omega1,rho1,Omega2,rho2,TA) + call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,Omega1,rho1,Omega2,rho2,TB) + +!------------------! +! Singlet manifold ! +!------------------! + + ispin = 3 + EcBSE = 0d0 + + ! Compute BSE singlet excitation energies + + call linear_response_Tmatrix(ispin,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & + EcBSE,OmBSE,XpY_BSE,XmY_BSE) + + call print_excitation('BSE@GT ',ispin,nS,OmBSE) + +end subroutine Bethe_Salpeter_Tmatrix_so diff --git a/src/GT/G0T0.f90 b/src/GT/G0T0.f90 index 11f4415..993b94f 100644 --- a/src/GT/G0T0.f90 +++ b/src/GT/G0T0.f90 @@ -184,8 +184,12 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing if(linearize) then +! eG0T0(:) = eHF(:) + Z(:)*SigT(:) eG0T0(:) = eHF(:) + Z(:)*(SigX(:) + SigT(:) - Vxc(:)) + call matout(nBas,1,SigX) + call matout(nBas,1,Vxc) + else eG0T0(:) = eHF(:) + SigX(:) + SigT(:) - Vxc(:) diff --git a/src/GT/excitation_density_Tmatrix_so.f90 b/src/GT/excitation_density_Tmatrix_so.f90 new file mode 100644 index 0000000..f0749f6 --- /dev/null +++ b/src/GT/excitation_density_Tmatrix_so.f90 @@ -0,0 +1,80 @@ +subroutine excitation_density_Tmatrix_so(nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) + +! Compute excitation densities for T-matrix self-energy + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + +! Local variables + + integer :: j,k,l + integer :: b,c,d + integer :: p,q + integer :: ab,cd,ij,kl + double precision,external :: Kronecker_delta + +! Output variables + + double precision,intent(out) :: rho1(nBas,nBas,nVV) + double precision,intent(out) :: rho2(nBas,nBAs,nOO) + +! Initialization + + rho1(:,:,:) = 0d0 + rho2(:,:,:) = 0d0 + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + + do ab=1,nVV + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) + end do + end do + + end do + + do ij=1,nOO + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) + end do + end do + + end do + end do + + end do + +end subroutine excitation_density_Tmatrix_so diff --git a/src/GT/renormalization_factor_Tmatrix_so.f90 b/src/GT/renormalization_factor_Tmatrix_so.f90 new file mode 100644 index 0000000..ca98217 --- /dev/null +++ b/src/GT/renormalization_factor_Tmatrix_so.f90 @@ -0,0 +1,63 @@ +subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,Z) + +! Compute renormalization factor of the T-matrix self-energy + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: i,j,k,l,a,b,c,d,p,cd,kl + double precision :: eps + +! Output variables + + double precision,intent(out) :: Z(nBas) + +! Initialize + + Z(:) = 0d0 + +!---------------------------------------------- +! T-matrix renormalization factor in the spinorbital basis +!---------------------------------------------- + +! Occupied part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do cd=1,nVV + eps = e(p) + e(i) - Omega1(cd) + Z(p) = Z(p) + (rho1(p,i,cd)/eps)**2 + enddo + enddo + enddo + +! Virtual part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do kl=1,nOO + eps = e(p) + e(a) - Omega2(kl) + Z(p) = Z(p) + (rho2(p,a,kl)/eps)**2 + enddo + enddo + enddo + +! Compute renormalization factor from derivative of SigT + + Z(:) = 1d0/(1d0 + Z(:)) + +end subroutine renormalization_factor_Tmatrix_so diff --git a/src/GT/self_energy_Tmatrix_diag_so.f90 b/src/GT/self_energy_Tmatrix_diag_so.f90 new file mode 100644 index 0000000..9c7834b --- /dev/null +++ b/src/GT/self_energy_Tmatrix_diag_so.f90 @@ -0,0 +1,63 @@ +subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT) + +! Compute diagonal of the correlation part of the T-matrix 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) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: i,j,k,l,a,b,c,d,p,cd,kl + double precision :: eps + +! Output variables + + double precision,intent(out) :: SigT(nBas) + +! Initialize + + SigT(:) = 0d0 + +!---------------------------------------------- +! T-matrix self-energy in the spinorbital basis +!---------------------------------------------- + +! Occupied part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do cd=1,nVV + eps = e(p) + e(i) - Omega1(cd) + SigT(p) = SigT(p) + rho1(p,i,cd)**2/eps + enddo + enddo + enddo + +! Virtual part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do kl=1,nOO + eps = e(p) + e(a) - Omega2(kl) + SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps + enddo + enddo + enddo + +end subroutine self_energy_Tmatrix_diag_so diff --git a/src/GT/soG0T0.f90 b/src/GT/soG0T0.f90 new file mode 100644 index 0000000..540ed47 --- /dev/null +++ b/src/GT/soG0T0.f90 @@ -0,0 +1,124 @@ +subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF) + +! Perform G0W0 calculation with a T-matrix self-energy (G0T0) in the spinorbital basis + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: ispin + integer :: nOO + integer :: nVV + double precision :: EcRPA + double precision :: EcGM + double precision :: EcBSE + integer :: nBas2,nC2,nO2,nV2,nR2,nS2 + double precision,allocatable :: Omega1(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + double precision,allocatable :: rho1(:,:,:) + double precision,allocatable :: Omega2(:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) + double precision,allocatable :: rho2(:,:,:) + double precision,allocatable :: SigT(:) + double precision,allocatable :: Z(:) + double precision,allocatable :: eG0T0(:) + double precision,allocatable :: seHF(:) + double precision,allocatable :: sERI(:,:,:,:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| One-shot soG0T0 calculation |' + write(*,*)'************************************************' + write(*,*) + +! Define occupied and virtual spaces + + nBas2 = 2*nBas + nO2 = 2*nO + nV2 = 2*nV + nC2 = 2*nC + nR2 = 2*nR + nS2 = nO2*nV2 + + +! Spatial to spin orbitals + + allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) + + call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) + call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI) + +! Dimensions of the rr-RPA linear reponse matrices + + nOO = nO2*(nO2 - 1)/2 + nVV = nV2*(nV2 - 1)/2 + +! Memory allocation + + allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & + Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + rho1(nBas2,nBas2,nVV),rho2(nBas2,nBas2,nOO), & + eG0T0(nBas2),SigT(nBas2),Z(nBas2)) + +!---------------------------------------------- +! Spinorbital basis +!---------------------------------------------- + + ispin = 4 + +! Compute linear response + + call linear_response_pp(ispin,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,1d0,seHF,sERI, & + Omega1,X1,Y1,Omega2,X2,Y2,EcRPA) + + call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1) + call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2) + +! Compute excitation densities for the T-matrix + + call excitation_density_Tmatrix_so(nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI,X1,Y1,rho1,X2,Y2,rho2) + +!---------------------------------------------- +! Compute T-matrix version of the self-energy +!---------------------------------------------- + + call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF, & + Omega1,rho1,Omega2,rho2,SigT) + +! Compute renormalization factor for T-matrix self-energy + + call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF, & + Omega1,rho1,Omega2,rho2,Z) + +!---------------------------------------------- +! Solve the quasi-particle equation +!---------------------------------------------- + + eG0T0(:) = seHF(:) + Z(:)*SigT(:) + +!---------------------------------------------- +! Dump results +!---------------------------------------------- + + call print_G0T0(nBas2,nO2,seHF,ENuc,ERHF,SigT,Z,eG0T0,EcGM,EcRPA) + + + call Bethe_Salpeter_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nS2,nOO,nVV,Omega1,X1,Y1,Omega2,X2,Y2,rho1,rho2, & + sERI,seHF,eG0T0,EcBSE) + +end subroutine soG0T0 diff --git a/src/GW/self_energy_exchange_diag.f90 b/src/GW/self_energy_exchange_diag.f90 index 0e305b7..eb0cebc 100644 --- a/src/GW/self_energy_exchange_diag.f90 +++ b/src/GW/self_energy_exchange_diag.f90 @@ -14,7 +14,8 @@ subroutine self_energy_exchange_diag(nBas,c,P,ERI,SigX) ! Local variables - integer :: q,mu,nu + integer :: mu,nu + integer :: q double precision,allocatable :: Fx(:,:) ! Output variables diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index a66f2aa..78e1dad 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -205,6 +205,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T ! Compute Vx for post-HF calculations - call mo_fock_exchange_potential(nBas,c,K,Vx) + call mo_fock_exchange_potential(nBas,c,P,ERI,Vx) end subroutine RHF diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index afa1fad..c6801c2 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -253,7 +253,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO ! Compute Vx for post-HF calculations do ispin=1,nspin - call mo_fock_exchange_potential(nBas,c(:,:,ispin),K(:,:,ispin),Vx(:,ispin)) + call mo_fock_exchange_potential(nBas,c(:,:,ispin),P(:,:,ispin),ERI,Vx(:,ispin)) end do end subroutine UHF diff --git a/src/HF/mo_fock_exchange_potential.f90 b/src/HF/mo_fock_exchange_potential.f90 index 56203c3..d07fcce 100644 --- a/src/HF/mo_fock_exchange_potential.f90 +++ b/src/HF/mo_fock_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine mo_fock_exchange_potential(nBas,c,Fx,Vx) +subroutine mo_fock_exchange_potential(nBas,c,P,ERI,Vx) ! Compute the exchange potential in the MO basis @@ -9,12 +9,14 @@ subroutine mo_fock_exchange_potential(nBas,c,Fx,Vx) integer,intent(in) :: nBas double precision,intent(in) :: c(nBas,nBas) - double precision,intent(in) :: Fx(nBas,nBas) + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables integer :: mu,nu - integer :: p + integer :: q + double precision,allocatable :: Fx(:,:) ! Output variables @@ -22,13 +24,18 @@ subroutine mo_fock_exchange_potential(nBas,c,Fx,Vx) ! Compute Vx + allocate(Fx(nBas,nBas)) + call exchange_matrix_AO_basis(nBas,P,ERI,Fx) + Vx(:) = 0d0 - do p=1,nBas + do q=1,nBas do mu=1,nBas do nu=1,nBas - Vx(p) = Vx(p) + c(mu,p)*Fx(mu,nu)*c(nu,p) + Vx(q) = Vx(q) + c(mu,q)*Fx(mu,nu)*c(nu,q) end do end do end do + deallocate(Fx) + end subroutine mo_fock_exchange_potential diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 2959cec..67ff676 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -1159,6 +1159,7 @@ program QuAcK else +! call soG0T0(eta_GT,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,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)