diff --git a/input/methods b/input/methods index 61c5be6..07cc57d 100644 --- a/input/methods +++ b/input/methods @@ -13,9 +13,9 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - T F F F F + 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 fb295a4..e50dfd1 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.00001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T T T T + F T F T T # 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 diff --git a/src/GT/Bethe_Salpeter_Tmatrix.f90 b/src/GT/Bethe_Salpeter_Tmatrix.f90 index 1b415c2..bab1596 100644 --- a/src/GT/Bethe_Salpeter_Tmatrix.f90 +++ b/src/GT/Bethe_Salpeter_Tmatrix.f90 @@ -85,8 +85,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TAs) - if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TBs) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,Omega1s,rho1s,Omega2s,rho2s,TAs) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,Omega1s,rho1s,Omega2s,rho2s,TBs) ! print*,'ab block of TA' ! call matout(nS,nS,TAs) @@ -105,8 +105,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TAt) - if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TBt) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,Omega1t,rho1t,Omega2t,rho2t,TAt) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,Omega1t,rho1t,Omega2t,rho2t,TBt) ! print*,'aa block of TA' ! call matout(nS,nS,TAt) diff --git a/src/GT/dynamic_Tmatrix_A.f90 b/src/GT/dynamic_Tmatrix_A.f90 index 91a3075..78a0460 100644 --- a/src/GT/dynamic_Tmatrix_A.f90 +++ b/src/GT/dynamic_Tmatrix_A.f90 @@ -62,12 +62,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*eps/(eps**2 + eta**2) end do -! do kl=1,nOO -! eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) -! chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*eps/(eps**2 + eta**2) -! end do + do kl=1,nOO + eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) + chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*eps/(eps**2 + eta**2) + end do - TA(ia,jb) = TA(ia,jb) + 1d0*lambda*chi + TA(ia,jb) = TA(ia,jb) + lambda*chi chi = 0d0 @@ -76,12 +76,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do -! do kl=1,nOO -! eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) -! chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 -! end do + do kl=1,nOO + eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) + chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do - ZA(ia,jb) = ZA(ia,jb) - 1d0*lambda*chi + ZA(ia,jb) = ZA(ia,jb) - lambda*chi end do end do diff --git a/src/GT/dynamic_Tmatrix_B.f90 b/src/GT/dynamic_Tmatrix_B.f90 new file mode 100644 index 0000000..7e0aa90 --- /dev/null +++ b/src/GT/dynamic_Tmatrix_B.f90 @@ -0,0 +1,91 @@ +subroutine dynamic_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,TB,ZB) + +! Compute the off-diagonal dynamic part of the Bethe-Salpeter equation matrices for GT + + 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) :: lambda + double precision,intent(in) :: eGT(nBas) + double precision,intent(in) :: OmBSE + + + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,cd,kl + +! Output variables + + double precision,intent(out) :: TB(nS,nS) + double precision,intent(out) :: ZB(nS,nS) + +! Initialization + + TB(:,:) = 0d0 + ZB(:,:) = 0d0 + +! Build dynamic A matrix + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + chi = 0d0 + + do cd=1,nVV + eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) + chi = chi + rho1(i,j,cd)*rho1(b,a,cd)*eps/(eps**2 + eta**2) + end do + + do kl=1,nOO + eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) + chi = chi + rho2(i,j,kl)*rho2(b,a,kl)*eps/(eps**2 + eta**2) + end do + + TB(ia,jb) = TB(ia,jb) + lambda*chi + + chi = 0d0 + + do cd=1,nVV + eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) + chi = chi + rho1(i,j,cd)*rho1(b,a,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + + do kl=1,nOO + eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) + chi = chi + rho2(i,a,kl)*rho2(b,a,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + + ZB(ia,jb) = ZB(ia,jb) - lambda*chi + + end do + end do + end do + end do + +end subroutine dynamic_Tmatrix_B diff --git a/src/GT/evGT.f90 b/src/GT/evGT.f90 index c600825..1c5e518 100644 --- a/src/GT/evGT.f90 +++ b/src/GT/evGT.f90 @@ -46,7 +46,6 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ! Local variables - logical :: linear_mixing integer :: nSCF integer :: n_diis double precision :: rcond diff --git a/src/GT/qsGT.f90 b/src/GT/qsGT.f90 index e2a1cf0..1e4c349 100644 --- a/src/GT/qsGT.f90 +++ b/src/GT/qsGT.f90 @@ -70,7 +70,6 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T integer :: nOOs,nOOt integer :: nVVs,nVVt - logical :: print_W = .false. double precision,allocatable :: error_diis(:,:) double precision,allocatable :: F_diis(:,:) double precision,allocatable :: c(:,:) @@ -212,22 +211,47 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO, & X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & - Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) + if(regularize) then - call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & - Omega1s,rho1s,Omega2s,rho2s,Z) + call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & + Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) + + call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & + Omega1s,rho1s,Omega2s,rho2s,Z) + + else + + call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & + Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) + + call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & + Omega1s,rho1s,Omega2s,rho2s,Z) + + end if iblock = 4 call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO, & X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & - Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) + if(regularize) then - call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & - Omega1t,rho1t,Omega2t,rho2t,Z) + call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & + Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) + + call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & + Omega1t,rho1t,Omega2t,rho2t,Z) + + else + + + call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & + Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) + + call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & + Omega1t,rho1t,Omega2t,rho2t,Z) + + end if Z(:) = 1d0/(1d0 - Z(:)) diff --git a/src/GT/static_Tmatrix_A.f90 b/src/GT/static_Tmatrix_A.f90 index 2e1586a..0ac807e 100644 --- a/src/GT/static_Tmatrix_A.f90 +++ b/src/GT/static_Tmatrix_A.f90 @@ -1,4 +1,4 @@ -subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TA) +subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,Omega2,rho2,TA) ! Compute the OOVV block of the static T-matrix for the resonant block @@ -7,7 +7,6 @@ subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome ! Input variables - integer,intent(in) :: ispin double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -18,7 +17,6 @@ subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome integer,intent(in) :: nOO integer,intent(in) :: nVV double precision,intent(in) :: lambda - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: Omega1(nVV) double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: Omega2(nOO) diff --git a/src/GT/static_Tmatrix_B.f90 b/src/GT/static_Tmatrix_B.f90 index 01ef3ab..eac09fd 100644 --- a/src/GT/static_Tmatrix_B.f90 +++ b/src/GT/static_Tmatrix_B.f90 @@ -1,4 +1,4 @@ -subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TB) +subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,Omega2,rho2,TB) ! Compute the OVVO block of the static T-matrix for the coupling block @@ -7,7 +7,6 @@ subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome ! Input variables - integer,intent(in) :: ispin double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -18,7 +17,6 @@ subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome integer,intent(in) :: nOO integer,intent(in) :: nVV double precision,intent(in) :: lambda - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: Omega1(nVV) double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: Omega2(nOO)