diff --git a/mol/h2.xyz b/mol/h2.xyz index 9e53116..fc52308 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 0.741 +H 0. 0. 1.0 diff --git a/src/GT/UG0T0.f90 b/src/GT/UG0T0.f90 index 0c5f7a5..5e13e1e 100644 --- a/src/GT/UG0T0.f90 +++ b/src/GT/UG0T0.f90 @@ -106,6 +106,7 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & Omega2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), & Omega1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & + Omega2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & SigX(nBas,nspin),SigT(nBas,nspin),Z(nBas,nspin)) @@ -120,10 +121,10 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & ! Compute linear response call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,eHF,ERI_aaaa, & + nPab,nHaa,nHab,nHbb,nHab,1d0,eHF,ERI_aaaa, & ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & Omega2ab,X2ab,Y2ab,EcRPA(ispin)) - + ! EcRPA(ispin) = 1d0*EcRPA(ispin) call print_excitation('pp-RPA (N+2)',iblock,nPab,Omega1ab(:)) @@ -139,11 +140,10 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & ! Compute linear response call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nP_sf,nHaa,nHab,nHbb,nH_sc,1d0,eHF,ERI_aaaa, & + nPaa,nHaa,nHab,nHbb,nHaa,1d0,eHF,ERI_aaaa, & ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & Omega2aa,X2aa,Y2aa,EcRPA(ispin)) - - + ! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin) @@ -160,11 +160,10 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & ! Compute linear response call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nP_sc,nHaa,nHab,nHbb,nH_sf,1d0,eHF,ERI_aaaa, & + nPbb,nHaa,nHab,nHbb,nHbb,1d0,eHF,ERI_aaaa, & ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & Omega2bb,X2bb,Y2bb,EcRPA(ispin)) - ! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin) @@ -180,44 +179,38 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & Z(:,:) = 0d0 !alpha-beta block - ispin = 1 + iblock = 1 call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHab,nPab, & ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & rho1ab,X2ab,Y2ab,rho2ab) - - call unrestricted_self_energy_Tmatrix_diag(ispin,eta,nBas,nC,nO,nV,nR,nHab,nPab,eHF, & - Omega1ab,rho1ab,Omega2ab,rho2ab,EcGM,SigT) - - call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHab,nPab,eHF, & - Omega1ab,rho1ab,Omega2ab,rho2ab,Z) !alpha-alpha block - ispin = 2 - iblock = 2 + iblock = 2 + call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa, & ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & rho1aa,X2aa,Y2aa,rho2aa) - call unrestricted_self_energy_Tmatrix_diag(ispin,eta,nBas,nC,nO,nV,nR,nHaa,nPaa,eHF, & - Omega1aa,rho1aa,Omega2aa,rho2aa,EcGM,SigT) - - call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nPaa,eHF, & - Omega1aa,rho1aa,Omega2aa,rho2aa,Z) -!beta-beta block - ispin = 2 - iblock = 3 +!beta-beta block + + iblock = 3 call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb, & ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & rho1bb,X2bb,Y2bb,rho2bb) - call unrestricted_self_energy_Tmatrix_diag(ispin,eta,nBas,nC,nO,nV,nR,nHbb,nPbb,eHF, & - Omega1bb,rho1bb,Omega2bb,rho2bb,EcGM,SigT) + call unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& + nPab,nPbb,eHF,Omega1aa,Omega1ab,Omega1bb,& + rho1aa,rho1ab,rho1bb,Omega2aa,Omega2ab,& + Omega2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) - call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHbb,nPbb,eHF, & - Omega1bb,rho1bb,Omega2bb,rho2bb,Z) + call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,& + nPaa,nPab,nPbb,eHF,Omega1aa,Omega1ab,& + Omega1bb,rho1aa,rho1ab,rho1bb, & + Omega2aa,Omega2ab,Omega2bb,rho2aa, & + rho2ab,rho2bb,Z) Z(:,:) = 1d0/(1d0 - Z(:,:)) @@ -241,7 +234,7 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & ! call matout(nBas,1,SigX) ! call matout(nBas,1,Vxc) - + call matout(nBas,1,eG0T0(:,1)) else eG0T0(:,:) = eHF(:,:) + SigX(:,:) + SigT(:,:) - Vxc(:,:) @@ -259,20 +252,20 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & ispin = 1 iblock = 1 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & - Omega2ab,X2ab,Y2ab,EcRPA(ispin)) + !call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + ! nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, & + ! ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & + ! Omega2ab,X2ab,Y2ab,EcRPA(ispin)) !alpha-alpha block ispin = 2 iblock = 2 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & - Omega2aa,X2aa,Y2aa,EcRPA(ispin)) + !call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + ! nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & + ! ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & + ! Omega2aa,X2aa,Y2aa,EcRPA(ispin)) Ecaa = EcRPA(2) @@ -280,10 +273,10 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & iblock = 3 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & - Omega2bb,X2bb,Y2bb,EcRPA(ispin)) + !call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + ! nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & + ! ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & + ! Omega2bb,X2bb,Y2bb,EcRPA(ispin)) Ecbb = EcRPA(2) EcRPA(2) = Ecaa + Ecbb diff --git a/src/GT/unrestricted_excitation_density_Tmatrix.f90 b/src/GT/unrestricted_excitation_density_Tmatrix.f90 index 178b17d..73a14fa 100644 --- a/src/GT/unrestricted_excitation_density_Tmatrix.f90 +++ b/src/GT/unrestricted_excitation_density_Tmatrix.f90 @@ -52,7 +52,7 @@ subroutine unrestricted_excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nH,nP, do ab=1,nP cd = 0 do c=nO(1)+1,nBas-nR(1) - do d=nO(2),nBas-nR(1) + do d=nO(2)+1,nBas-nR(1) cd = cd + 1 rho1(p,q,ab) = rho1(p,q,ab) & diff --git a/src/GT/unrestricted_renormalization_factor_Tmatrix.f90 b/src/GT/unrestricted_renormalization_factor_Tmatrix.f90 index 94635f7..918f22d 100644 --- a/src/GT/unrestricted_renormalization_factor_Tmatrix.f90 +++ b/src/GT/unrestricted_renormalization_factor_Tmatrix.f90 @@ -1,4 +1,8 @@ -subroutine unrestricted_renormalization_factor_Tmatrix(ispin,eta,nBas,nC,nO,nV,nR,nH,nP,e,Omega1,rho1,Omega2,rho2,Z) +subroutine unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, & + nPaa,nPab,nPbb,e,Omega1aa,Omega1ab, & + Omega1bb,rho1aa,rho1ab,rho1bb, & + Omega2aa,Omega2ab,Omega2bb,rho2aa, & + rho2ab,rho2bb,Z) ! Compute renormalization factor of the T-matrix self-energy @@ -9,13 +13,15 @@ subroutine unrestricted_renormalization_factor_Tmatrix(ispin,eta,nBas,nC,nO,nV,n double precision,intent(in) :: eta integer,intent(in) :: nBas,nC(nspin),nO(nspin),nV(nspin),nR(nspin) - integer,intent(in) :: nH,ispin - integer,intent(in) :: nP + integer,intent(in) :: nHaa,nHab,nHbb + integer,intent(in) :: nPaa,nPab,nPbb double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: Omega1(nP) - double precision,intent(in) :: rho1(nBas,nBas,nP,nspin) - double precision,intent(in) :: Omega2(nH) - double precision,intent(in) :: rho2(nBas,nBas,nH,nspin) + double precision,intent(in) :: Omega1aa(nPaa),Omega1ab(nPab),Omega1bb(nPbb) + double precision,intent(in) :: rho1aa(nBas,nBas,nPaa),rho1ab(nBas,nBas,nPab) + double precision,intent(in) :: rho1bb(nBas,nBas,nPbb) + double precision,intent(in) :: Omega2aa(nHaa),Omega2ab(nHab),Omega2bb(nHbb) + double precision,intent(in) :: rho2aa(nBas,nBas,nHaa),rho2ab(nBas,nBas,nHab) + double precision,intent(in) :: rho2bb(nBas,nBas,nHbb) ! Local variables @@ -25,57 +31,81 @@ subroutine unrestricted_renormalization_factor_Tmatrix(ispin,eta,nBas,nC,nO,nV,n ! Output variables double precision,intent(out) :: Z(nBas,nspin) - -! Occupied part of the T-matrix self-energy - if(ispin==1) then - - do p=nC(1)+1,nBas-nR(1) - do i=nC(1)+1,nO(1) - do cd=1,nP - eps = e(p,1) + e(i,1) - Omega1(cd) - Z(p,1) = Z(p,1) - rho1(p,i,cd,1)**2*(eps/(eps**2 + eta**2))**2 - enddo - enddo - enddo - -! Virtual part of the T-matrix self-energy - - do p=nC(1)+1,nBas-nR(1) - do a=nO(1)+1,nBas-nR(1) - do kl=1,nH - eps = e(p,1) + e(a,1) - Omega2(kl) - Z(p,1) = Z(p,1) - rho2(p,a,kl,1)**2*(eps/(eps**2 + eta**2))**2 - enddo - enddo - enddo - - end if +!spin up part ! Occupied part of the T-matrix self-energy - - if(ispin==2) then - do p=nC(2)+1,nBas-nR(2) - do i=nC(2)+1,nO(2) - do cd=1,nP - eps = e(p,2) + e(i,2) - Omega1(cd) - Z(p,2) = Z(p,2) - rho1(p,i,cd,2)**2*(eps/(eps**2 + eta**2))**2 - enddo + do p=nC(1)+1,nBas-nR(1) + do i=nC(1)+1,nO(1) + do cd=1,nPaa + eps = e(p,1) + e(i,1) - Omega1aa(cd) + Z(p,1) = Z(p,1) - rho1aa(p,i,cd)**2*(eps/(eps**2 + eta**2))**2 enddo enddo + + do i=nC(1)+1,nO(1) + do cd=1,nPab + eps = e(p,1) + e(i,1) - Omega1ab(cd) + Z(p,1) = Z(p,1) - rho1ab(p,i,cd)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + enddo ! Virtual part of the T-matrix self-energy - do p=nC(2)+1,nBas-nR(2) - do a=nO(2)+1,nBas-nR(2) - do kl=1,nH - eps = e(p,2) + e(a,2) - Omega2(kl) - Z(p,2) = Z(p,2) - rho2(p,a,kl,2)**2*(eps/(eps**2 + eta**2))**2 - enddo + do p=nC(1)+1,nBas-nR(1) + do a=nO(1)+1,nBas-nR(1) + do kl=1,nHaa + eps = e(p,1) + e(a,1) - Omega2aa(kl) + Z(p,1) = Z(p,1) - rho2aa(p,a,kl)**2*(eps/(eps**2 + eta**2))**2 enddo enddo - end if + do a=nO(2)+1,nBas-nR(2) + do kl=1,nHab + eps = e(p,1) + e(a,1) - Omega2ab(kl) + Z(p,1) = Z(p,1) - rho2ab(p,a,kl)**2*(eps/(eps**2 + eta**2))**2 + enddo + enddo + enddo + +!spin down part + +! Occupied part of the T-matrix self-energy + + do p=nC(2)+1,nBas-nR(2) + do i=nC(2)+1,nO(2) + do cd=1,nPbb + eps = e(p,2) + e(i,2) - Omega1bb(cd) + Z(p,2) = Z(p,2) - rho1bb(p,i,cd)**2*(eps/(eps**2 + eta**2))**2 + enddo + enddo + + do i=nC(1)+1,nO(1) + do cd=1,nPab + eps = e(p,2) + e(i,2) - Omega1ab(cd) + Z(p,2) = Z(p,2) - rho1ab(p,i,cd)**2*(eps/(eps**2 + eta**2))**2 + enddo + enddo + enddo + +! Virtual part of the T-matrix self-energy + + do p=nC(2)+1,nBas-nR(2) + do a=nO(2)+1,nBas-nR(2) + do kl=1,nHbb + eps = e(p,2) + e(a,2) - Omega2bb(kl) + Z(p,2) = Z(p,2) - rho2bb(p,a,kl)**2*(eps/(eps**2 + eta**2))**2 + enddo + enddo + + do a=nO(1)+1,nBas-nR(1) + do kl=1,nHab + eps = e(p,2) + e(a,2) - Omega2ab(kl) + Z(p,2) = Z(p,2) - rho2ab(p,a,kl)**2*(eps/(eps**2 + eta**2))**2 + enddo + enddo + enddo end subroutine unrestricted_renormalization_factor_Tmatrix diff --git a/src/GT/unrestricted_self_energy_Tmatrix_diag.f90 b/src/GT/unrestricted_self_energy_Tmatrix_diag.f90 index b5311ce..fea06a4 100644 --- a/src/GT/unrestricted_self_energy_Tmatrix_diag.f90 +++ b/src/GT/unrestricted_self_energy_Tmatrix_diag.f90 @@ -1,4 +1,7 @@ -subroutine unrestricted_self_energy_Tmatrix_diag(ispin,eta,nBas,nC,nO,nV,nR,nH,nP,e,Omega1,rho1,Omega2,rho2,EcGM,SigT) +subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& + nPab,nPbb,e,Omega1aa,Omega1ab,Omega1bb,& + rho1aa,rho1ab,rho1bb,Omega2aa,Omega2ab,& + Omega2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) ! Compute diagonal of the correlation part of the T-matrix self-energy @@ -13,13 +16,15 @@ subroutine unrestricted_self_energy_Tmatrix_diag(ispin,eta,nBas,nC,nO,nV,nR,nH,n integer,intent(in) :: nO(nspin) integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) - integer,intent(in) :: nH,ispin - integer,intent(in) :: nP + integer,intent(in) :: nHaa,nHab,nHbb + integer,intent(in) :: nPaa,nPab,nPbb double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: Omega1(nP) - double precision,intent(in) :: rho1(nBas,nBas,nP,nspin) - double precision,intent(in) :: Omega2(nH) - double precision,intent(in) :: rho2(nBas,nBas,nH,nspin) + double precision,intent(in) :: Omega1aa(nPaa),Omega1ab(nPab),Omega1bb(nPbb) + double precision,intent(in) :: rho1aa(nBas,nBas,nPaa),rho1ab(nBas,nBas,nPab) + double precision,intent(in) :: rho1bb(nBas,nBas,nPbb) + double precision,intent(in) :: Omega2aa(nHaa),Omega2ab(nHab),Omega2bb(nHbb) + double precision,intent(in) :: rho2aa(nBas,nBas,nHaa),rho2ab(nBas,nBas,nHab) + double precision,intent(in) :: rho2bb(nBas,nBas,nHbb) ! Local variables @@ -35,113 +40,162 @@ subroutine unrestricted_self_energy_Tmatrix_diag(ispin,eta,nBas,nC,nO,nV,nR,nH,n ! Occupied part of the T-matrix self-energy !---------------------------------------------- - if(ispin==1) then +!spin up part - do p=nC(1)+1,nBas-nR(1) - do i=nC(1)+1,nO(1) - do cd=1,nP - eps = e(p,1) + e(i,1) - Omega1(cd) - SigT(p,1) = SigT(p,1) + rho1(p,i,cd,1)**2*eps/(eps**2 + eta**2) - enddo + do p=nC(1)+1,nBas-nR(1) + do i=nC(1)+1,nO(1) + do cd=1,nPaa + eps = e(p,1) + e(i,1) - Omega1aa(cd) + SigT(p,1) = SigT(p,1) + rho1aa(p,i,cd)**2*eps/(eps**2 + eta**2) enddo enddo - end if + do i=nC(2)+1,nO(2) + do cd=1,nPab + eps = e(p,1) + e(i,1) - Omega1ab(cd) + SigT(p,1) = SigT(p,1) + rho1ab(p,i,cd)**2*eps/(eps**2 + eta**2) + end do + end do + enddo -!beta part +!spin down part - if(ispin==2) then - - do p=nC(2)+1,nBas-nR(2) - do i=nC(2)+1,nO(2) - do cd=1,nP - eps = e(p,2) + e(i,2) - Omega1(cd) - SigT(p,2) = SigT(p,2) + rho1(p,i,cd,2)**2*eps/(eps**2 + eta**2) - enddo + do p=nC(2)+1,nBas-nR(2) + do i=nC(2)+1,nO(2) + do cd=1,nPbb + eps = e(p,2) + e(i,2) - Omega1bb(cd) + SigT(p,2) = SigT(p,2) + rho1bb(p,i,cd)**2*eps/(eps**2 + eta**2) enddo enddo - - end if + + do i=nC(2)+1,nO(2) + do cd=1,nPab + eps = e(p,2) + e(i,2) - Omega1ab(cd) + SigT(p,2) = SigT(p,2) + rho1ab(p,i,cd)**2*eps/(eps**2 + eta**2) + end do + end do + enddo !---------------------------------------------- ! Virtual part of the T-matrix self-energy !---------------------------------------------- - !alpha part +! spin up part - if(ispin==1) then - - do p=nC(1)+1,nBas-nR(1) - do a=nO(1)+1,nBas-nR(1) - do kl=1,nH - eps = e(p,1) + e(a,1) - Omega2(kl) - SigT(p,1) = SigT(p,1) + rho2(p,a,kl,1)**2*eps/(eps**2 + eta**2) - enddo - enddo - enddo - - end if - - !alpha part - - if(ispin==2) then - - do p=nC(2)+1,nBas-nR(2) - do a=nO(2)+1,nBas-nR(2) - do kl=1,nH - eps = e(p,2) + e(a,2) - Omega2(kl) - SigT(p,2) = SigT(p,2) + rho2(p,a,kl,2)**2*eps/(eps**2 + eta**2) - enddo - enddo - enddo - - end if - -!---------------------------------------------- -! Galitskii-Migdal correlation energy -!---------------------------------------------- - if(ispin==1) then - - do i=nC(1)+1,nO(1) - do j=nC(1)+1,nO(1) - do cd=1,nP - eps = e(i,1) + e(j,1) - Omega1(cd) - EcGM(1) = EcGM(1) + rho1(i,j,cd,1)*rho1(i,j,cd,1)*eps/(eps**2 + eta**2) - enddo - enddo - enddo - + do p=nC(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1) - do b=nO(1)+1,nBas-nR(1) - do kl=1,nH - eps = e(a,1) + e(b,1) - Omega2(kl) - EcGM(1) = EcGM(1) - rho2(a,b,kl,1)*rho2(a,b,kl,1)*eps/(eps**2 + eta**2) - enddo + do kl=1,nHaa + eps = e(p,1) + e(a,1) - Omega2aa(kl) + SigT(p,1) = SigT(p,1) + rho2aa(p,a,kl)**2*eps/(eps**2 + eta**2) enddo - enddo + end do + + do a=nO(1)+1,nBas-nR(1) + do kl=1,nHab + eps = e(p,1) + e(a,1) - Omega2ab(kl) + SigT(p,1) = SigT(p,1) + rho2ab(p,a,kl)**2*eps/(eps**2 + eta**2) + end do + end do + enddo - end if +!spin down part - if(ispin==2) then - - do i=nC(2)+1,nO(2) - do j=nC(2)+1,nO(2) - do cd=1,nP - eps = e(i,2) + e(j,2) - Omega1(cd) - EcGM(2) = EcGM(2) + rho1(i,j,cd,2)*rho1(i,j,cd,2)*eps/(eps**2 + eta**2) - enddo + do p=nC(2)+1,nBas-nR(2) + do a=nO(2)+1,nBas-nR(2) + do kl=1,nHbb + eps = e(p,2) + e(a,2) - Omega2bb(kl) + SigT(p,2) = SigT(p,2) + rho2bb(p,a,kl)**2*eps/(eps**2 + eta**2) enddo enddo do a=nO(2)+1,nBas-nR(2) - do b=nO(2)+1,nBas-nR(2) - do kl=1,nH - eps = e(a,2) + e(b,2) - Omega2(kl) - EcGM(2) = EcGM(2) - rho2(a,b,kl,2)*rho2(a,b,kl,2)*eps/(eps**2 + eta**2) - enddo + do kl=1,nHab + eps = e(p,2) + e(a,2) - Omega2ab(kl) + SigT(p,2) = SigT(p,2) + rho2ab(p,a,kl)**2*eps/(eps**2 + eta**2) + end do + end do + enddo + + + +!---------------------------------------------- +! Galitskii-Migdal correlation energy +!---------------------------------------------- + +!spin up part + + do i=nC(1)+1,nO(1) + do j=nC(1)+1,nO(1) + do cd=1,nPaa + eps = e(i,1) + e(j,1) - Omega1aa(cd) + EcGM(1) = EcGM(1) + rho1aa(i,j,cd)*rho1aa(i,j,cd)*eps/(eps**2 + eta**2) enddo enddo + enddo + + do i=nC(1)+1,nO(1) + do j=nC(2)+1,nO(2) + do cd=1,nPab + eps = e(i,1) + e(j,1) - Omega1ab(cd) + EcGM(1) = EcGM(1) + rho1ab(i,j,cd)*rho1ab(i,j,cd)*eps/(eps**2 + eta**2) + end do + end do + end do - end if + do a=nO(1)+1,nBas-nR(1) + do b=nO(1)+1,nBas-nR(1) + do kl=1,nHaa + eps = e(a,1) + e(b,1) - Omega2aa(kl) + EcGM(1) = EcGM(1) - rho2aa(a,b,kl)*rho2aa(a,b,kl)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + + do a=nO(1)+1,nBas-nR(1) + do b=nO(1)+1,nBas-nR(1) + do kl=1,nHab + eps = e(a,1) + e(b,1) - Omega2ab(kl) + EcGM(1) = EcGM(1) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + +! spin down part + + do i=nC(2)+1,nO(2) + do j=nC(2)+1,nO(2) + do cd=1,nPbb + eps = e(i,2) + e(j,2) - Omega1bb(cd) + EcGM(2) = EcGM(2) + rho1bb(i,j,cd)*rho1bb(i,j,cd)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + + do i=nC(1)+1,nO(1) + do j=nC(2)+1,nO(2) + do cd=1,nPab + eps = e(i,2) + e(j,2) - Omega1ab(cd) + EcGM(2) = EcGM(2) + rho1ab(i,j,cd)*rho1ab(i,j,cd)*eps/(eps**2 + eta**2) + end do + end do + end do + + do a=nO(1)+1,nBas-nR(1) + do b=nO(2)+1,nBas-nR(2) + do kl=1,nHab + eps = e(a,2) + e(b,2) - Omega2ab(kl) + EcGM(2) = EcGM(2) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + + do a=nO(2)+1,nBas-nR(2) + do b=nO(2)+1,nBas-nR(2) + do kl=1,nHbb + eps = e(a,2) + e(b,2) - Omega2bb(kl) + EcGM(2) = EcGM(2) - rho2bb(a,b,kl)*rho2bb(a,b,kl)*eps/(eps**2 + eta**2) + enddo + enddo + enddo end subroutine unrestricted_self_energy_Tmatrix_diag diff --git a/src/LR/unrestricted_linear_response_pp.f90 b/src/LR/unrestricted_linear_response_pp.f90 index 89acb7f..72aa75b 100644 --- a/src/LR/unrestricted_linear_response_pp.f90 +++ b/src/LR/unrestricted_linear_response_pp.f90 @@ -65,8 +65,8 @@ subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab, call unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,& lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,C) - call unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,& - nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,& + call unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,& + nHaa,nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,& ERI_bbbb,B) call unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,& @@ -90,7 +90,7 @@ subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab, call sort_ppRPA(nHt,nPt,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),& Y2(:,:)) - + ! Compute the RPA correlation energy EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nPt,C(:,:)) & @@ -99,5 +99,5 @@ subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab, EcRPA2 = -sum(Omega2(:)) - trace_matrix(nHt,D(:,:)) if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' - + end subroutine unrestricted_linear_response_pp