From feb1a1bb40c940dce0fa57c5ccecbd3256d89838 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 15 Oct 2021 13:51:04 +0200 Subject: [PATCH] working on BSE@GT --- src/MBPT/Bethe_Salpeter.f90 | 7 +- src/MBPT/Bethe_Salpeter_Tmatrix.f90 | 156 ++++++++++++++++++++++++++++ src/MBPT/static_Tmatrix_TA.f90 | 64 ++++++++++++ src/MBPT/static_Tmatrix_TB.f90 | 64 ++++++++++++ src/MBPT/static_screening_WB.f90 | 7 +- 5 files changed, 296 insertions(+), 2 deletions(-) create mode 100644 src/MBPT/Bethe_Salpeter_Tmatrix.f90 create mode 100644 src/MBPT/static_Tmatrix_TA.f90 create mode 100644 src/MBPT/static_Tmatrix_TB.f90 diff --git a/src/MBPT/Bethe_Salpeter.f90 b/src/MBPT/Bethe_Salpeter.f90 index 4d8621d..a71949a 100644 --- a/src/MBPT/Bethe_Salpeter.f90 +++ b/src/MBPT/Bethe_Salpeter.f90 @@ -16,7 +16,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, logical,intent(in) :: triplet double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + 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) :: eW(nBas) double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) diff --git a/src/MBPT/Bethe_Salpeter_Tmatrix.f90 b/src/MBPT/Bethe_Salpeter_Tmatrix.f90 new file mode 100644 index 0000000..d62cbb8 --- /dev/null +++ b/src/MBPT/Bethe_Salpeter_Tmatrix.f90 @@ -0,0 +1,156 @@ +subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eT,eGT,EcBSE) + +! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA_T + logical,intent(in) :: TDA + logical,intent(in) :: singlet + logical,intent(in) :: triplet + + 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) :: eT(nBas) + double precision,intent(in) :: eGT(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: nOOs + integer :: nOOt + integer :: nVVs + integer :: nVVt + integer :: ispin + integer :: iblock + integer :: dERI + integer :: xERI + + double precision :: EcRPA(nspin) + double precision,allocatable :: Omega1s(:),Omega1t(:) + double precision,allocatable :: X1s(:,:),X1t(:,:) + double precision,allocatable :: Y1s(:,:),Y1t(:,:) + double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:) + double precision,allocatable :: Omega2s(:),Omega2t(:) + double precision,allocatable :: X2s(:,:),X2t(:,:) + double precision,allocatable :: Y2s(:,:),Y2t(:,:) + double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) + + 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(nspin) + +! Memory allocation + + allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & + rho1s(nBas,nO,nVVs),rho2s(nBas,nV,nOOs), & + Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & + rho1t(nBas,nO,nVVt),rho2t(nBas,nV,nOOt)) + + allocate(TA(nS,nS),TB(nS,nS),OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) + +! Dimensions of the pp-RPA linear reponse matrices + + nOOs = nO*nO + nVVs = nV*nV + + nOOt = nO*(nO - 1)/2 + nVVt = nV*(nV - 1)/2 + +! Initialize T matrix + + TA(:,:) = 0d0 + TB(:,:) = 0d0 + +!---------------------------------------------- +! Compute T-matrix for alpha-beta block +!---------------------------------------------- + + ispin = 1 + iblock = 3 + dERI = +1d0 + xERI = +0d0 + + call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eT(:),ERI(:,:,:,:), & + Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin)) + + call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & + X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:)) + + call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) + if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) + +!---------------------------------------------- +! Compute T-matrix for alpha-alpha block +!---------------------------------------------- + + ispin = 2 + iblock = 4 + dERI = +1d0 + xERI = -1d0 + + call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eT(:),ERI(:,:,:,:), & + Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin)) + + call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & + X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) + + call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) + if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) + +!------------------- +! Singlet manifold +!------------------- + + if(singlet) then + + ispin = 1 + EcBSE(ispin) = 0d0 + + ! Compute BSE singlet excitation energies + + call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,OmRPA, & + rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + + end if + +!------------------- +! Triplet manifold +!------------------- + + if(triplet) then + + ispin = 2 + EcBSE(ispin) = 0d0 + + ! Compute BSE triplet excitation energies + + call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,OmRPA, & + rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + + end if + +end subroutine Bethe_Salpeter_Tmatrix diff --git a/src/MBPT/static_Tmatrix_TA.f90 b/src/MBPT/static_Tmatrix_TA.f90 new file mode 100644 index 0000000..37f88ed --- /dev/null +++ b/src/MBPT/static_Tmatrix_TA.f90 @@ -0,0 +1,64 @@ +subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TA) + +! Compute the OOVV block of the static T-matrix for the resonant block + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: eta + 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,nO,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho2(nBas,nV,nOO) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kl,cd + +! Output variables + + double precision,intent(out) :: TA(nS,nS) + + 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 = Omega1(cd)**2 + eta**2 + chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/eps + enddo + + do kl=1,nOO + eps = Omega2(kl)**2 + eta**2 + chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/eps + enddo + + TA(ia,jb) = TA(ia,jb) + lambda*ERI(i,b,j,a) - 2d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine static_Tmatrix_TA diff --git a/src/MBPT/static_Tmatrix_TB.f90 b/src/MBPT/static_Tmatrix_TB.f90 new file mode 100644 index 0000000..b804f7e --- /dev/null +++ b/src/MBPT/static_Tmatrix_TB.f90 @@ -0,0 +1,64 @@ +subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TB) + +! Compute the OVVO block of the static T-matrix for the coupling block + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: eta + 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,nO,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho2(nBas,nV,nOO) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kl,cd + +! Output variables + + double precision,intent(out) :: TB(nS,nS) + + 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 = Omega1(cd)**2 + eta**2 + chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/eps + enddo + + do kl=1,nOO + eps = Omega2(kl)**2 + eta**2 + chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/eps + enddo + + TB(ia,jb) = TB(ia,jb) + lambda*ERI(i,j,b,a) - 2d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine static_Tmatrix_TB diff --git a/src/MBPT/static_screening_WB.f90 b/src/MBPT/static_screening_WB.f90 index 4419660..ee7ac48 100644 --- a/src/MBPT/static_screening_WB.f90 +++ b/src/MBPT/static_screening_WB.f90 @@ -7,7 +7,12 @@ subroutine static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WB) ! Input variables - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + 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) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)