From a5047e9b1a6d8341c0b37d3394a13ab1286a2832 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 6 Oct 2019 22:35:36 +0200 Subject: [PATCH] T-matrix self-energy --- input/methods | 2 +- input/options | 2 +- src/QuAcK/QuAcK.f90 | 2 +- src/QuAcK/excitation_density_Tmatrix.f90 | 77 +++++++++++++++++ src/QuAcK/linear_response_ph.f90 | 88 ++++++++++++++++++++ src/QuAcK/linear_response_pp.f90 | 63 +++++++------- src/QuAcK/ppRPA.f90 | 40 ++++++--- src/QuAcK/print_excitation.f90 | 4 +- src/QuAcK/renormalization_factor_Tmatrix.f90 | 65 +++++++++++++++ src/QuAcK/self_energy_Tmatrix.f90 | 71 ++++++++++++++++ src/QuAcK/self_energy_Tmatrix_diag.f90 | 67 +++++++++++++++ 11 files changed, 429 insertions(+), 52 deletions(-) create mode 100644 src/QuAcK/excitation_density_Tmatrix.f90 create mode 100644 src/QuAcK/linear_response_ph.f90 create mode 100644 src/QuAcK/renormalization_factor_Tmatrix.f90 create mode 100644 src/QuAcK/self_energy_Tmatrix.f90 create mode 100644 src/QuAcK/self_energy_Tmatrix_diag.f90 diff --git a/input/methods b/input/methods index c03cfd5..94f75a0 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF MOM T F F # MP2 MP3 MP2-F12 - F F F + T F F # CCD CCSD CCSD(T) F F F # CIS TDHF ppRPA ADC diff --git a/input/options b/input/options index d1a5689..da3080c 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.00001 F 1 # CIS/TDHF/BSE: singlet triplet - T T + T F # GF: maxSCF thresh DIIS n_diis renormalization 64 0.00001 T 10 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index d81486f..1a5a238 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -408,7 +408,7 @@ program QuAcK if(doppRPA) then call cpu_time(start_ppRPA) - call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF) + call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF) call cpu_time(end_ppRPA) t_ppRPA = end_ppRPA - start_ppRPA diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 new file mode 100644 index 0000000..ff7dba6 --- /dev/null +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -0,0 +1,77 @@ +subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,X2,Y2,rho1,rho2) + +! Compute excitation densities for T-matrix self-energy + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nR,nOO,nVV + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(out) :: X1(nVV,nVV) + double precision,intent(out) :: Y1(nOO,nVV) + double precision,intent(out) :: X2(nVV,nOO) + double precision,intent(out) :: Y2(nOO,nOO) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + integer :: p + integer :: ab,cd,ij,kl + +! Output variables + + double precision,intent(out) :: rho1(nBas,nBas,nVV) + double precision,intent(out) :: rho2(nBas,nBas,nOO) + + rho1(:,:,:) = 0d0 + rho2(:,:,:) = 0d0 + + do p=nC+1,nBas-nR + do i=nC+1,nO + do ab=1,nVV + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) + enddo + enddo + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) + enddo + enddo + + enddo + enddo + + do a=nO+1,nBas-nR + do ij=1,nOO + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho2(p,a,ij) = rho2(p,a,ij) + (ERI(p,a,c,d) - ERI(p,a,d,c))*X2(cd,ij) + enddo + enddo + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho1(p,a,ij) = rho1(p,a,ij) + (ERI(p,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij) + enddo + enddo + + enddo + enddo + enddo + +end subroutine excitation_density_Tmatrix diff --git a/src/QuAcK/linear_response_ph.f90 b/src/QuAcK/linear_response_ph.f90 new file mode 100644 index 0000000..3183188 --- /dev/null +++ b/src/QuAcK/linear_response_ph.f90 @@ -0,0 +1,88 @@ +subroutine linear_response_ph(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,Ec_phRPA,Omega,XpY) + +! Compute the p-h channel of the linear response + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dRPA,TDA,BSE + integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas),rho(nBas,nBas,nS) + +! Local variables + + double precision,external :: trace_matrix + double precision,allocatable :: A(:,:),B(:,:),M(:,:),w(:) + +! Output variables + + double precision,intent(out) :: Ec_phRPA + double precision,intent(out) :: Omega(nS),XpY(nS,nS) + + +! Memory allocation + + allocate(A(nS,nS),B(nS,nS),M(2*nS,2*nS),w(2*nS)) + +! Build A and B matrices + + call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A) + if(BSE) call Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,A) + +! Tamm-Dancoff approximation + + B(:,:) = 0d0 + if(.not. TDA) then + + call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B) + if(BSE) call Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,B) + + endif + +!------------------------------------------------------------------------ +! Solve the p-h eigenproblem +!------------------------------------------------------------------------ +! +! | +A +B | | X Y | | w 0 | | X Y | +! | | | | = | | | | +! | -B -A | | Y X | | 0 -w | | Y X | +! + +! Diagonal blocks + + M(1:nS,1:nS) = +A(1:nS,1:nS) + M(nS+1:2*nS,nS+1:2*nS) = -A(1:nS,1:nS) + + +! Off-diagonal blocks + + M(1:nS,nS+1:2*nS) = -B(1:nS,1:nS) + M(nS+1:2*nS,1:nS) = +B(1:nS,1:nS) + +! Diagonalize the p-h matrix + + call diagonalize_matrix(2*nS,M(:,:),w(:)) + + Omega(1:nS) = w(nS+1:2*nS) + +! Build X+Y + + XpY(1:nS,1:nS) = M(nS+1:2*nS,1:nS) + M(nS+1:2*nS,nS+1:2*nS) + + call DA(nS,1d0/sqrt(Omega),XpY) + +! print*,'X+Y' +! call matout(nS,nS,XpY) + +! print*,'RPA excitations' + call matout(2*nS,1,w(:)) + +! Compute the RPA correlation energy + + Ec_phRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) + + print*,'Ec_phRPA = ',Ec_phRPA + +end subroutine linear_response_ph diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index 6b7339e..044c8dc 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -1,4 +1,5 @@ -subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,Ec_ppRPA) +subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & + Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA) ! Compute the p-p channel of the linear response: see Scueria et al. JCP 139, 104113 (2013) @@ -10,42 +11,35 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,E logical,intent(in) :: dRPA logical,intent(in) :: TDA logical,intent(in) :: BSE - integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS + integer,intent(in) :: ispin,nBas,nC,nO,nV,nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables - integer :: nOO - integer :: nVV double precision :: trace_matrix double precision,allocatable :: B(:,:) double precision,allocatable :: C(:,:) double precision,allocatable :: D(:,:) double precision,allocatable :: M(:,:) double precision,allocatable :: Z(:,:) - double precision,allocatable :: w(:) - double precision,allocatable :: w1(:) - double precision,allocatable :: w2(:) - double precision,allocatable :: X1(:,:) - double precision,allocatable :: Y1(:,:) - double precision,allocatable :: X2(:,:) - double precision,allocatable :: Y2(:,:) + double precision,allocatable :: Omega(:) ! Output variables + double precision,intent(out) :: Omega1(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) :: X2(nVV,nOO) + double precision,intent(out) :: Y2(nOO,nOO) double precision,intent(out) :: Ec_ppRPA -! Useful quantities - - nOO = nO*(nO-1)/2 - nVV = nV*(nV-1)/2 - ! Memory allocation - allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),w(nOO+nVV), & - w1(nVV),w2(nOO),X1(nVV,nVV),Y1(nVV,nOO),X2(nOO,nVV),Y2(nOO,nOO)) + allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Omega(nOO+nVV)) ! Build B, C and D matrices for the pp channel @@ -78,32 +72,31 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,E ! Diagonalize the p-h matrix Z(:,:) = M(:,:) - call diagonalize_matrix(nOO+nVV,Z(:,:),w(:)) + call diagonalize_matrix(nOO+nVV,Z(:,:),Omega(:)) - write(*,*) 'pp-RPA excitation energies' - call matout(nOO+nVV,1,w(:)) - write(*,*) +! write(*,*) 'pp-RPA excitation energies' +! call matout(nOO+nVV,1,Omega(:)) +! write(*,*) ! Split the various quantities in p-p and h-h parts - w1(:) = w(nOO+1:nOO+nVV) - w2(:) = w(1:nOO) + Omega1(:) = Omega(nOO+1:nOO+nVV) + Omega2(:) = Omega(1:nOO) - X1(:,:) = Z(nOO+1:nOO+nVV, 1:nVV ) - Y1(:,:) = Z(nOO+1:nOO+nVV,nVV+1:nOO+nVV) - X2(:,:) = Z( 1:nOO , 1:nVV ) - Y2(:,:) = Z( 1:nOO ,nVV+1:nOO+nVV) + X1(:,:) = Z(nOO+1:nOO+nVV,nOO+1:nOO+nVV) + Y1(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV) + X2(:,:) = Z(nOO+1:nOO+nVV, 1:nOO ) + Y2(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV) - if(minval(w1(:)) < 0d0) call print_warning('You may have instabilities in pp-RPA!!') - if(maxval(w2(:)) > 0d0) call print_warning('You may have instabilities in pp-RPA!!') + if(minval(Omega1(:)) < 0d0) call print_warning('You may have instabilities in pp-RPA!!') + if(maxval(Omega2(:)) > 0d0) call print_warning('You may have instabilities in pp-RPA!!') ! Compute the RPA correlation energy - Ec_ppRPA = 0.5d0*( sum(w1(:)) - sum(w2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) + Ec_ppRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) print*,'Ec(pp-RPA) = ',Ec_ppRPA - print*,'Ec(pp-RPA) = ',0.5d0*( sum(abs(w(:))) - trace_matrix(nVV*nOO,M(:,:))) - print*,'Ec(pp-RPA) = ',+sum(w1(:)) - trace_matrix(nVV,C(:,:)) - print*,'Ec(pp-RPA) = ',-sum(w2(:)) - trace_matrix(nOO,D(:,:)) + print*,'Ec(pp-RPA) = ',+sum(Omega1(:)) - trace_matrix(nVV,C(:,:)) + print*,'Ec(pp-RPA) = ',-sum(Omega2(:)) - trace_matrix(nOO,D(:,:)) end subroutine linear_response_pp diff --git a/src/QuAcK/ppRPA.f90 b/src/QuAcK/ppRPA.f90 index 5b82091..c5447df 100644 --- a/src/QuAcK/ppRPA.f90 +++ b/src/QuAcK/ppRPA.f90 @@ -1,4 +1,4 @@ -subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) +subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Perform pp-RPA calculation @@ -14,7 +14,6 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF double precision,intent(in) :: e(nBas) @@ -26,10 +25,15 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF logical :: TDA logical :: BSE integer :: ispin - double precision,allocatable :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) + integer :: nOO + integer :: nVV + double precision,allocatable :: Omega1(:,:) + double precision,allocatable :: X1(:,:,:) + double precision,allocatable :: Y1(:,:,:) + double precision,allocatable :: Omega2(:,:) + double precision,allocatable :: X2(:,:,:) + double precision,allocatable :: Y2(:,:,:) - double precision :: rho double precision :: Ec_ppRPA(nspin) ! Hello world @@ -40,6 +44,11 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF write(*,*)'****************************************' write(*,*) +! Useful quantities + + nOO = nO*(nO-1)/2 + nVV = nV*(nV-1)/2 + ! Initialization Ec_ppRPA(:) = 0d0 @@ -58,7 +67,8 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF ! Memory allocation - allocate(Omega(nS,nspin),XpY(nS,nS,nspin)) + allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & + Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin)) ! Singlet manifold @@ -66,9 +76,12 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF ispin = 1 - call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, & - Ec_ppRPA) -! call print_excitation('pp-RPA ',ispin,nS,Omega(:,ispin)) + call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & + Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & + Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & + Ec_ppRPA(ispin)) + call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) + call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin)) endif @@ -78,9 +91,12 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF ispin = 2 - call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, & - Ec_ppRPA) -! call print_excitation('pp-RPA ',ispin,nS,Omega(:,ispin)) + call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & + Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & + Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & + Ec_ppRPA(ispin)) + call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) + call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin)) endif diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90 index 889ed6e..d8e8a60 100644 --- a/src/QuAcK/print_excitation.f90 +++ b/src/QuAcK/print_excitation.f90 @@ -7,7 +7,7 @@ subroutine print_excitation(method,ispin,nS,Omega) ! Input variables - character*5,intent(in) :: method + character*12,intent(in) :: method integer,intent(in) :: ispin,nS double precision,intent(in) :: Omega(nS) @@ -22,7 +22,7 @@ subroutine print_excitation(method,ispin,nS,Omega) write(*,*) write(*,*)'-------------------------------------------------------------' - write(*,'(1X,A1,1X,A4,A14,A7,A9,A25)')'|',method,' calculation: ',spin_manifold,' manifold',' |' + write(*,'(1X,A1,1X,A14,A14,A7,A9,A15)')'|',method,' calculation: ',spin_manifold,' manifold',' |' write(*,*)'-------------------------------------------------------------' write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/QuAcK/renormalization_factor_Tmatrix.f90 new file mode 100644 index 0000000..8b76586 --- /dev/null +++ b/src/QuAcK/renormalization_factor_Tmatrix.f90 @@ -0,0 +1,65 @@ +subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1,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,nOO,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 + + ! Occupied part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + eps = e(p) + e(i) - Omega1(cd) + Z(p) = Z(p) - 2d0*rho1(p,i,cd)**2*eps/(eps**2 + eta**2) + enddo + enddo + enddo + enddo + + ! Virtual part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + eps = e(p) + e(a) - Omega2(kl) + Z(p) = Z(p) - 2d0*rho2(p,a,kl)**2*(eps/(eps**2 + eta**2))**2 + enddo + enddo + enddo + enddo + +! Compute renormalization factor from derivative of SigT + + Z(:) = 1d0/(1d0 - Z(:)) + +end subroutine renormalization_factor_Tmatrix diff --git a/src/QuAcK/self_energy_Tmatrix.f90 b/src/QuAcK/self_energy_Tmatrix.f90 new file mode 100644 index 0000000..25cc8e6 --- /dev/null +++ b/src/QuAcK/self_energy_Tmatrix.f90 @@ -0,0 +1,71 @@ +subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1,rho2,SigT) + +! Compute 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,q,cd,kl + double precision :: eps + +! Output variables + + double precision,intent(out) :: SigT(nBas,nBas) + +! Initialize + + SigT = 0d0 + + ! Occupied part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do i=nC+1,nO + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + eps = e(p) + e(i) - Omega1(cd) + SigT(p,q) = SigT(p,q) + 2d0*rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + enddo + enddo + + ! Virtual part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do a=nO+1,nBas-nR + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + eps = e(p) + e(a) - Omega2(kl) + SigT(p,q) = SigT(p,q) + 2d0*rho2(p,a,kl)*rho2(q,a,kl)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + enddo + enddo + +end subroutine self_energy_Tmatrix diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/QuAcK/self_energy_Tmatrix_diag.f90 new file mode 100644 index 0000000..24cd19d --- /dev/null +++ b/src/QuAcK/self_energy_Tmatrix_diag.f90 @@ -0,0 +1,67 @@ +subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1,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 + + ! Occupied part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + eps = e(p) + e(i) - Omega1(cd) + SigT(p) = SigT(p) + 2d0*rho1(p,i,cd)**2*eps/(eps**2 + eta**2) + enddo + enddo + enddo + enddo + + ! Virtual part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + eps = e(p) + e(a) - Omega2(kl) + SigT(p) = SigT(p) + 2d0*rho2(p,a,kl)**2*eps/(eps**2 + eta**2) + enddo + enddo + enddo + enddo + +end subroutine self_energy_Tmatrix_diag