From 6e3abd89271bd9740849feff03f010213457cc56 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 16 Oct 2019 18:14:47 +0200 Subject: [PATCH] starting G0T0 --- examples/molecule.LiH | 2 +- input/methods | 6 +- src/QuAcK/G0T0.f90 | 153 +++++++++++++++++++ src/QuAcK/QuAcK.f90 | 35 +++-- src/QuAcK/excitation_density_Tmatrix.f90 | 2 +- src/QuAcK/print_G0T0.f90 | 49 ++++++ src/QuAcK/read_methods.f90 | 14 ++ src/QuAcK/renormalization_factor_Tmatrix.f90 | 2 +- src/QuAcK/self_energy_Tmatrix_diag.f90 | 2 +- 9 files changed, 247 insertions(+), 18 deletions(-) create mode 100644 src/QuAcK/G0T0.f90 create mode 100644 src/QuAcK/print_G0T0.f90 diff --git a/examples/molecule.LiH b/examples/molecule.LiH index 6878c9e..771beea 100644 --- a/examples/molecule.LiH +++ b/examples/molecule.LiH @@ -2,4 +2,4 @@ 2 2 2 0 0 # Znuc x y z Li 0. 0. 0. - H 0. 0. 3.015 + H 0. 0. 3.061356 diff --git a/input/methods b/input/methods index 7ef4e42..9ac4f33 100644 --- a/input/methods +++ b/input/methods @@ -5,10 +5,12 @@ # CCD CCSD CCSD(T) F F F # CIS TDHF ppRPA ADC - F F T F + F F F F # GF2 GF3 F F -# G0W0 evGW qsGW +# G0W0 evGW qsGW F F F +# G0T0 evGT qsGT + T F F # MCMP2 F diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 new file mode 100644 index 0000000..6f8b47e --- /dev/null +++ b/src/QuAcK/G0T0.f90 @@ -0,0 +1,153 @@ +subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nOO,nVV,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0T0) + +! Perform G0W0 calculation with a T-matrix self-energy (G0T0) + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: BSE + logical,intent(in) :: singlet_manifold + logical,intent(in) :: triplet_manifold + double precision,intent(in) :: eta + + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: PHF(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: H(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + logical :: dRPA + integer :: ispin,jspin + double precision :: EcRPA(nspin) + double precision :: EcBSE(nspin) + double precision :: EcGM + double precision,allocatable :: SigT(:) + double precision,allocatable :: Z(:) + 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(:,:,:,:) + +! Output variables + + double precision :: eG0T0(nBas) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| One-shot G0T0 calculation |' + write(*,*)'************************************************' + write(*,*) + +! Spin manifold + + ispin = 1 + +! Memory allocation + + allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & + Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin), & + rho1(nBas,nBas,nVV,nspin),rho2(nBas,nBas,nOO,nspin), & + SigT(nBas),Z(nBas)) + +! Compute linear response + + call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,eHF,ERI, & + Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & + Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & + EcRPA(ispin)) + +! Compute excitation densities for the T-matrix + + call excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI, & + X1(:,:,ispin),Y1(:,:,ispin),rho1(:,:,:,ispin), & + X2(:,:,ispin),Y2(:,:,ispin),rho2(:,:,:,ispin)) + +! Compute T-matrix version of the self-energy + + call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF, & + Omega1(:,ispin),rho1(:,:,:,ispin), & + Omega2(:,ispin),rho2(:,:,:,ispin), & + SigT) + +! Compute renormalization factor for T-matrix self-energy + + call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF, & + Omega1(:,ispin),rho1(:,:,:,ispin), & + Omega2(:,ispin),rho2(:,:,:,ispin), & + Z(:)) + +! Solve the quasi-particle equation + + eG0T0(:) = eHF(:) + Z(:)*SigT(:) + +! Dump results + + call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) + call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin)) + + call print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eG0T0,EcRPA(ispin)) + +! Perform BSE calculation + +! if(BSE) then + +! ! Singlet manifold + +! if(singlet_manifold) then + +! ispin = 1 +! EcBSE(ispin) = 0d0 + +! call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, & +! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) +! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + +! call linear_response(ispin,.false.,.false.,BSE,nBas,nC,nO,nV,nR,nS,eG0T0,ERI, & +! rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) +! call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) + +! endif + +! ! Triplet manifold + +! if(triplet_manifold) then + +! ispin = 2 +! EcBSE(ispin) = 0d0 + +! call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, & +! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) +! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + +! call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0T0,ERI, & +! rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) +! call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) + +! endif + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy (singlet) =',EcBSE(1) +! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy (triplet) =',EcBSE(2) +! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy =',EcBSE(1) + EcBSE(2) +! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! endif + +end subroutine G0T0 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 33cce74..af170fd 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -10,18 +10,20 @@ program QuAcK logical :: doCIS,doTDHF,doppRPA,doADC logical :: doGF2,doGF3 logical :: doG0W0,doevGW,doqsGW + logical :: doG0T0,doevGT,doqsGT logical :: doMCMP2,doMinMCMP2 - logical :: doeNcusp logical :: doBas integer :: nNuc,nBas,nBasCABS - integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin),nS(nspin) + integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin) + integer :: nS(nspin),nOO(nspin),nVV(nspin) double precision :: ENuc,ERHF,EUHF,Norm double precision :: EcMP2(3),EcMP3,EcMP2F12(3),EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3) double precision,allocatable :: ZNuc(:),rNuc(:,:) double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:) double precision,allocatable :: eG0W0(:) + double precision,allocatable :: eG0T0(:) integer :: nShell integer,allocatable :: TotAngMomShell(:),KShell(:) @@ -50,7 +52,9 @@ program QuAcK double precision :: start_G0W0 ,end_G0W0 ,t_G0W0 double precision :: start_evGW ,end_evGW ,t_evGW double precision :: start_qsGW ,end_qsGW ,t_qsGW - double precision :: start_eNcusp ,end_eNcusp ,t_eNcusp + double precision :: start_G0T0 ,end_G0T0 ,t_G0T0 + double precision :: start_evGT ,end_evGT ,t_evGT + double precision :: start_qsGT ,end_qsGT ,t_qsGT double precision :: start_MP2 ,end_MP2 ,t_MP2 double precision :: start_MP3 ,end_MP3 ,t_MP3 double precision :: start_MP2F12 ,end_MP2F12 ,t_MP2F12 @@ -108,6 +112,7 @@ program QuAcK doCIS,doTDHF,doppRPA,doADC, & doGF2,doGF3, & doG0W0,doevGW,doqsGW, & + doG0T0,doevGT,doqsGT, & doMCMP2) ! Read options for methods @@ -121,7 +126,6 @@ program QuAcK ! Weird stuff - doeNCusp = .false. doMinMCMP2 = .false. !------------------------------------------------------------------------ @@ -175,7 +179,7 @@ program QuAcK ! Memory allocation for one- and two-electron integrals - allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas),PHF(nBas,nBas,nspin), & + allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas),eG0T0(nBas),PHF(nBas,nBas,nspin), & S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas), & ERI_AO_basis(nBas,nBas,nBas,nBas),ERI_MO_basis(nBas,nBas,nBas,nBas)) @@ -522,16 +526,23 @@ program QuAcK end if !------------------------------------------------------------------------ -! Compute e-N cusp dressing +! Perform G0T0 calculatiom !------------------------------------------------------------------------ - if(doeNcusp) then - call cpu_time(start_eNcusp) -! call eNcusp() - call cpu_time(end_eNcusp) + eG0T0(:) = eHF(:,1) - t_eNcusp = end_eNcusp - start_eNcusp - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for e-N cusp dressing = ',t_eNcusp,' seconds' + if(doG0T0) then + + nOO(:) = nO(:)*(nO(:)+1)/2 + nVV(:) = nV(:)*(nV(:)+1)/2 + + call cpu_time(start_G0T0) + call G0T0(BSE,singlet_manifold,triplet_manifold,eta, & + nBas,nC(1),nO(1),nV(1),nR(1),nOO(1),nVV(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0T0) + call cpu_time(end_G0T0) + + t_G0T0 = end_G0T0 - start_G0T0 + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds' write(*,*) end if diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 index ff7dba6..dbabd5c 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -1,4 +1,4 @@ -subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,X2,Y2,rho1,rho2) +subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) ! Compute excitation densities for T-matrix self-energy diff --git a/src/QuAcK/print_G0T0.f90 b/src/QuAcK/print_G0T0.f90 new file mode 100644 index 0000000..53272b7 --- /dev/null +++ b/src/QuAcK/print_G0T0.f90 @@ -0,0 +1,49 @@ +subroutine print_G0T0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) + +! Print one-electron energies and other stuff for G0T0 + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas,nO + double precision,intent(in) :: ENuc + double precision,intent(in) :: EHF + double precision,intent(in) :: EcRPA + double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) + + integer :: x,HOMO,LUMO + double precision :: Gap + +! HOMO and LUMO + + HOMO = nO + LUMO = HOMO + 1 + Gap = eGW(LUMO)-eGW(HOMO) + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + write(*,*)' 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_c (eV)','|','Z','|','e_QP (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + do x=1,nBas + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'G0T0 HOMO energy (eV):',eGW(HOMO)*HaToeV + write(*,'(2X,A30,F15.6)') 'G0T0 LUMO energy (eV):',eGW(LUMO)*HaToeV + write(*,'(2X,A30,F15.6)') 'G0T0 HOMO-LUMO gap (eV):',Gap*HaToeV + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'RPA@G0T0 total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@G0T0 correlation energy =',EcRPA + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +end subroutine print_G0T0 + + diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 5866a97..a8ea50e 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -4,6 +4,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & doCIS,doTDHF,doppRPA,doADC, & doGF2,doGF3, & doG0W0,doevGW,doqsGW, & + doG0T0,doevGT,doqsGT, & doMCMP2) ! Read desired methods @@ -18,6 +19,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & logical,intent(out) :: doCIS,doTDHF,doppRPA,doADC logical,intent(out) :: doGF2,doGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW + logical,intent(out) :: doG0T0,doevGT,doqsGT logical,intent(out) :: doMCMP2 ! Local variables @@ -51,6 +53,10 @@ subroutine read_methods(doRHF,doUHF,doMOM, & doGF3 = .false. doG0W0 = .false. + doevGT = .false. + doqsGT = .false. + + doG0T0 = .false. doevGW = .false. doqsGW = .false. @@ -104,6 +110,14 @@ subroutine read_methods(doRHF,doUHF,doMOM, & if(answer2 == 'T') doevGW = .true. if(answer3 == 'T') doqsGW = .true. +! Read GT methods + + read(1,*) + read(1,*) answer1,answer2,answer3 + if(answer1 == 'T') doG0T0 = .true. + if(answer2 == 'T') doevGT = .true. + if(answer3 == 'T') doqsGT = .true. + ! Read stochastic methods read(1,*) diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/QuAcK/renormalization_factor_Tmatrix.f90 index 8b76586..d95f8c5 100644 --- a/src/QuAcK/renormalization_factor_Tmatrix.f90 +++ b/src/QuAcK/renormalization_factor_Tmatrix.f90 @@ -1,4 +1,4 @@ -subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1,rho2,Z) +subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,Z) ! Compute renormalization factor of the T-matrix self-energy diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/QuAcK/self_energy_Tmatrix_diag.f90 index 24cd19d..0c9e8db 100644 --- a/src/QuAcK/self_energy_Tmatrix_diag.f90 +++ b/src/QuAcK/self_energy_Tmatrix_diag.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1,rho2,SigT) +subroutine self_energy_Tmatrix_diag(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