diff --git a/input/methods b/input/methods index 07cc57d..20a60d5 100644 --- a/input/methods +++ b/input/methods @@ -15,7 +15,7 @@ # G0W0* evGW* qsGW* ufG0W0 ufGW F F F F F # G0T0 evGT qsGT - T F F + T F F # MCMP2 F # * unrestricted version available diff --git a/input/options b/input/options index 90a6cbd..0b58cf0 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 + T T T 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/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..33217db 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)) @@ -114,16 +115,16 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & !---------------------------------------------- ispin = 1 - iblock = 1 + iblock = 3 ! iblock = 1 ! 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(:)) @@ -134,16 +135,15 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & !---------------------------------------------- ispin = 2 - iblock = 2 + iblock = 4 ! 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) @@ -155,16 +155,15 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & !---------------------------------------------- ispin = 2 - iblock = 3 + iblock = 7 ! 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 + + iblock = 3 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 = 4 + 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 = 7 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(:,:) @@ -257,20 +250,20 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & !alpha-beta block ispin = 1 - iblock = 1 + 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, & + 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 + iblock = 4 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, & + nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & Omega2aa,X2aa,Y2aa,EcRPA(ispin)) @@ -278,10 +271,10 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & !beta-beta block - iblock = 3 + iblock = 7 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, & + nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & Omega2bb,X2bb,Y2bb,EcRPA(ispin)) @@ -290,7 +283,7 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) -! call print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eG0T0,EcGM,EcRPA) + call print_UG0T0(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eG0T0,EcGM,EcRPA) ! Free memory diff --git a/src/GT/print_UG0T0.f90 b/src/GT/print_UG0T0.f90 new file mode 100644 index 0000000..b777d9e --- /dev/null +++ b/src/GT/print_UG0T0.f90 @@ -0,0 +1,69 @@ +subroutine print_UG0T0(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA) + +! Print one-electron energies and other stuff for UG0T0 + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) + double precision,intent(in) :: ENuc + double precision,intent(in) :: EUHF + double precision,intent(in) :: EcGM + double precision,intent(in) :: EcRPA(nspin) + double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: SigT(nBas,nspin) + double precision,intent(in) :: Z(nBas,nspin) + double precision,intent(in) :: eGT(nBas,nspin) + + integer :: p + integer :: ispin + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) + double precision :: Gap(nspin) + +! HOMO and LUMO + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = eGT(nO(ispin),ispin) + LUMO(ispin) = eGT(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = eGT(1,ispin) + Gap(ispin) = 0d0 + end if + end do + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + write(*,*)' Unrestricted one-shot G0T0 calculation (T-matrix self-energy) ' + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & + '|','#','|','e_HF (eV)','|','Sigma_T (eV)','|','Z','|','e_QP (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + do p=1,nBas + write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & + '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigT(p,1)*HaToeV,SigT(p,2)*HaToeV,'|', & + Z(p,1),Z(p,2),'|',eGT(p,1)*HaToeV,eGT(p,2)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F15.6,A3)') 'UG0T0 HOMO energy (eV) =',maxval(HOMO(:))*HaToeV,' eV' + write(*,'(2X,A50,F15.6,A3)') 'UG0T0 LUMO energy (eV) =',minval(LUMO(:))*HaToeV,' eV' + write(*,'(2X,A50,F15.6,A3)') 'UG0T0 HOMO-LUMO gap (eV) =',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV' + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 correlation energy (singlet) =',EcRPA(1),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 correlation energy (triplet) =',EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 correlation energy =',EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' GM@UG0T0 correlation energy =',EcGM,' au' + write(*,'(2X,A50,F20.10,A3)') ' GM@UG0T0 total energy =',ENuc + EUHF + EcGM,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +end subroutine print_UG0T0 + + 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/print_excitation.f90 b/src/LR/print_excitation.f90 index 1720300..d746cd6 100644 --- a/src/LR/print_excitation.f90 +++ b/src/LR/print_excitation.f90 @@ -20,9 +20,10 @@ subroutine print_excitation(method,ispin,nS,Omega) if(ispin == 1) spin_manifold = 'singlet' if(ispin == 2) spin_manifold = 'triplet' if(ispin == 3) spin_manifold = 'alpha-beta' - if(ispin == 4) spin_manifold = 'alpha-alpha' + if(ispin == 4) spin_manifold = 'alpha-alpha' if(ispin == 5) spin_manifold = 'spin-conserved' if(ispin == 6) spin_manifold = 'spin-flip' + if(ispin == 7) spin_manifold = 'beta-beta' write(*,*) write(*,*)'-------------------------------------------------------------' diff --git a/src/LR/unrestricted_linear_response_B_pp.f90 b/src/LR/unrestricted_linear_response_B_pp.f90 index 6f5db5b..e1c57bc 100644 --- a/src/LR/unrestricted_linear_response_B_pp.f90 +++ b/src/LR/unrestricted_linear_response_B_pp.f90 @@ -41,11 +41,7 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP eF = 0d0 -!----------------------------------------------- -! Build B matrix for spin-conserved transitions -!----------------------------------------------- - - if(ispin == 1) then + if(ispin == 3) then ! abab block @@ -67,11 +63,8 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP end if -!----------------------------------------------- -! Build B matrix for spin-flip transitions -!----------------------------------------------- - if(ispin == 2) then + if(ispin == 4) then ! aaaa block @@ -92,7 +85,7 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP end do end if - if (ispin == 3) then + if (ispin == 7) then ! bbbb block diff --git a/src/LR/unrestricted_linear_response_C_pp.f90 b/src/LR/unrestricted_linear_response_C_pp.f90 index 7e37f78..dd14c72 100644 --- a/src/LR/unrestricted_linear_response_C_pp.f90 +++ b/src/LR/unrestricted_linear_response_C_pp.f90 @@ -37,11 +37,8 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP eF = 0d0 -!----------------------------------------------- -! Build C matrix for spin-conserved transitions -!----------------------------------------------- - if(ispin == 1) then + if(ispin == 3) then ! abab block @@ -62,11 +59,7 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP end if -!----------------------------------------------- -! Build C matrix for spin-flip transitions -!----------------------------------------------- - - if(ispin == 2) then + if(ispin == 4) then ! aaaa block @@ -90,7 +83,7 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP end if - if (ispin == 3) then + if (ispin == 7) then ! bbbb block diff --git a/src/LR/unrestricted_linear_response_D_pp.f90 b/src/LR/unrestricted_linear_response_D_pp.f90 index 3d7c10e..28c4c78 100644 --- a/src/LR/unrestricted_linear_response_D_pp.f90 +++ b/src/LR/unrestricted_linear_response_D_pp.f90 @@ -38,11 +38,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH eF = 0d0 -!----------------------------------------------- -! Build D matrix for spin-conserved transitions -!----------------------------------------------- - - if(ispin == 1) then + if(ispin == 3) then ! abab block @@ -63,12 +59,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH end if - -!----------------------------------------------- -! Build D matrix for spin-flip transitions -!----------------------------------------------- - - if(ispin == 2) then + if(ispin == 4) then ! aaaa block @@ -91,7 +82,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH end do end if - if (ispin == 3) then + if (ispin == 7) then ! bbbb block 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 diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90 index 75b0dd4..752e98d 100644 --- a/src/RPA/ppURPA.f90 +++ b/src/RPA/ppURPA.f90 @@ -51,82 +51,75 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH Ec_ppURPA(:) = 0d0 EcAC(:) = 0d0 -! Spin-conserved manifold +!alpha-beta block - if(spin_conserved) then + ispin = 1 + iblock = 3 - ispin = 1 - iblock = 1 + nPab = nV(1)*nV(2) + nHab = nO(1)*nO(2) -!Spin-conserved quantities - - nPab = nV(1)*nV(2) - nHab = nO(1)*nO(2) - - nP_sc = nPab - nH_sc = nHab + nP_sc = nPab + nH_sc = nHab ! Memory allocation - allocate(Omega1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), & - Omega2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc)) + allocate(Omega1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), & + Omega2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc)) - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,e,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1sc,X1sc,Y1sc, & - Omega2sc,X2sc,Y2sc,Ec_ppURPA(ispin)) + call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,e,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Omega1sc,X1sc,Y1sc, & + Omega2sc,X2sc,Y2sc,Ec_ppURPA(ispin)) - call print_excitation('pp-RPA (N+2)',5,nP_sc,Omega1sc) - call print_excitation('pp-RPA (N-2)',5,nH_sc,Omega2sc) + call print_excitation('pp-RPA (N+2)',iblock,nP_sc,Omega1sc) + call print_excitation('pp-RPA (N-2)',iblock,nH_sc,Omega2sc) - endif +!alpha-alpha block -! Spin-flip manifold + ispin = 2 + iblock = 4 - if(spin_flip) then + nPaa = nV(1)*(nV(1)-1)/2 + nPbb = nV(2)*(nV(2)-1)/2 - ispin = 2 - iblock = 2 + nP_sf = nPaa -!Spin-flip quantities + nHaa = nO(1)*(nO(1)-1)/2 + nHbb = nO(2)*(nO(2)-1)/2 - nPaa = nV(1)*(nV(1)-1)/2 - nPbb = nV(2)*(nV(2)-1)/2 + nH_sf = nHaa - nP_sf = nPaa + allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & + Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) - nHaa = nO(1)*(nO(1)-1)/2 - nHbb = nO(2)*(nO(2)-1)/2 + call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf, & + Omega2sf,X2sf,Y2sf,Ec_ppURPA(ispin)) - nH_sf = nHaa + call print_excitation('pp-RPA (N+2)',iblock,nP_sf,Omega1sf) + call print_excitation('pp-RPA (N-2)',iblock,nH_sf,Omega2sf) - allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & - Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) + deallocate(Omega1sf,X1sf,Y1sf,Omega2sf,X2sf,Y2sf) - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf, & - Omega2sf,X2sf,Y2sf,Ec_ppURPA(ispin)) +!beta-beta block - deallocate(Omega1sf,X1sf,Y1sf,Omega2sf,X2sf,Y2sf) + iblock = 7 - iblock = 3 + nP_sf = nPbb + nH_sf = nHbb - nP_sf = nPbb - nH_sf = nHbb + allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & + Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) - allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & - Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) + call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,& + nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,& + ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf,& + Omega2sf,X2sf,Y2sf,Ec_ppURPA(ispin)) - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,& - nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,& - ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf,& - Omega2sf,X2sf,Y2sf,Ec_ppURPA(ispin)) - - call print_excitation('pp-RPA (N+2)',6,nP_sf,Omega1sf) - call print_excitation('pp-RPA (N-2)',6,nH_sf,Omega2sf) - - endif + call print_excitation('pp-RPA (N+2)',iblock,nP_sf,Omega1sf) + call print_excitation('pp-RPA (N-2)',iblock,nH_sf,Omega2sf) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -159,5 +152,4 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ! end if - end subroutine ppURPA