From 91954dba8d202587613cb0ed6274d25ca28da92d Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 17 Jul 2023 22:31:28 +0200 Subject: [PATCH] remove so GT routines --- src/GT/Bethe_Salpeter_Tmatrix_so.f90 | 78 ------------ 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 | 123 ------------------- src/LR/ppLR.f90 | 47 +++---- 6 files changed, 16 insertions(+), 438 deletions(-) delete mode 100644 src/GT/Bethe_Salpeter_Tmatrix_so.f90 delete mode 100644 src/GT/excitation_density_Tmatrix_so.f90 delete mode 100644 src/GT/renormalization_factor_Tmatrix_so.f90 delete mode 100644 src/GT/self_energy_Tmatrix_diag_so.f90 delete mode 100644 src/GT/soG0T0.f90 diff --git a/src/GT/Bethe_Salpeter_Tmatrix_so.f90 b/src/GT/Bethe_Salpeter_Tmatrix_so.f90 deleted file mode 100644 index 4eb6991..0000000 --- a/src/GT/Bethe_Salpeter_Tmatrix_so.f90 +++ /dev/null @@ -1,78 +0,0 @@ -subroutine Bethe_Salpeter_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,Om1,X1,Y1,Om2,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) :: Om1(nVV) - double precision,intent(in) :: X1(nVV,nVV) - double precision,intent(in) :: Y1(nOO,nVV) - double precision,intent(in) :: Om2(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 ppLR(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eT,ERI,Om1,X1,Y1,Om2,X2,Y2,EcRPA) - - call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,Om1,rho1,Om2,rho2,TA) - call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,Om1,rho1,Om2,rho2,TB) - -!------------------! -! Singlet manifold ! -!------------------! - - ispin = 3 - EcBSE = 0d0 - - ! Compute BSE singlet excitation energies - - call linear_response_BSE(ispin,.false.,.false.,.true.,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 diff --git a/src/GT/excitation_density_Tmatrix_so.f90 b/src/GT/excitation_density_Tmatrix_so.f90 deleted file mode 100644 index f0749f6..0000000 --- a/src/GT/excitation_density_Tmatrix_so.f90 +++ /dev/null @@ -1,80 +0,0 @@ -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 deleted file mode 100644 index ca98217..0000000 --- a/src/GT/renormalization_factor_Tmatrix_so.f90 +++ /dev/null @@ -1,63 +0,0 @@ -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 deleted file mode 100644 index 9c7834b..0000000 --- a/src/GT/self_energy_Tmatrix_diag_so.f90 +++ /dev/null @@ -1,63 +0,0 @@ -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 deleted file mode 100644 index c749389..0000000 --- a/src/GT/soG0T0.f90 +++ /dev/null @@ -1,123 +0,0 @@ -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 ppLR(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/LR/ppLR.f90 b/src/LR/ppLR.f90 index d6c189d..64abe6c 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -1,4 +1,4 @@ -subroutine ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcRPA) +subroutine ppLR(TDA,nBas,nOO,nVV,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013) @@ -7,18 +7,12 @@ subroutine ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Omega1,X1,Y1,Ome ! Input variables - integer,intent(in) :: ispin logical,intent(in) :: TDA - 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) :: lambda - double precision,intent(in) :: e(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: B(nVV,nOO) + double precision,intent(in) :: C(nVV,nVV) + double precision,intent(in) :: D(nOO,nOO) ! Local variables @@ -27,28 +21,24 @@ subroutine ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Omega1,X1,Y1,Ome double precision :: trace_matrix double precision :: EcRPA1 double precision :: EcRPA2 - double precision,allocatable :: B(:,:) - double precision,allocatable :: C(:,:) - double precision,allocatable :: D(:,:) double precision,allocatable :: M(:,:) double precision,allocatable :: Z(:,:) - double precision,allocatable :: Omega(:) + double precision,allocatable :: Om(:) ! Output variables - double precision,intent(out) :: Omega1(nVV) + double precision,intent(out) :: Om1(nVV) double precision,intent(out) :: X1(nVV,nVV) double precision,intent(out) :: Y1(nOO,nVV) - double precision,intent(out) :: Omega2(nOO) + double precision,intent(out) :: Om2(nOO) double precision,intent(out) :: X2(nVV,nOO) double precision,intent(out) :: Y2(nOO,nOO) double precision,intent(out) :: EcRPA ! Memory allocation - allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Omega(nOO+nVV)) -!write(*,*) 'nOO', nOO -!write(*,*) 'nVV', nVV + allocate(M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Om(nOO+nVV)) + !-------------------------------------------------! ! Solve the p-p eigenproblem ! !-------------------------------------------------! @@ -59,20 +49,15 @@ subroutine ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Omega1,X1,Y1,Ome ! ! !-------------------------------------------------! -! Build B, C and D matrices for the pp channel - - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D) - if(TDA) then X1(:,:) = +C(:,:) Y1(:,:) = 0d0 - if(nVV > 0) call diagonalize_matrix(nVV,X1,Omega1) + if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1) X2(:,:) = 0d0 Y2(:,:) = -D(:,:) - if(nOO > 0) call diagonalize_matrix(nOO,Y2,Omega2) + if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2) else @@ -90,19 +75,19 @@ subroutine ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Omega1,X1,Y1,Ome ! Diagonalize the p-p matrix - if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z) + if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Om,Z) ! Split the various quantities in p-p and h-h parts - call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:)) + call sort_ppRPA(nOO,nVV,Om(:),Z(:,:),Om1(:),X1(:,:),Y1(:,:),Om2(:),X2(:,:),Y2(:,:)) end if ! Compute the RPA correlation energy - EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) - EcRPA1 = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:)) - EcRPA2 = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:)) + EcRPA = 0.5d0*( sum(Om1(:)) - sum(Om2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) + EcRPA1 = +sum(Om1(:)) - trace_matrix(nVV,C(:,:)) + EcRPA2 = -sum(Om2(:)) - trace_matrix(nOO,D(:,:)) if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'