qsGF2

input/methods

@@ -9,11 +9,11 @@
 # CIS* CIS(D) CID CISD
   F    F      F   F    
 # RPA* RPAx* ppRPA 
-  T   F    F    
+  F   F    F    
-# G0F2 evGF2 G0F3 evGF3
+# G0F2 evGF2 qsGF2 G0F3 evGF3
-  F    F     F    F
+  F    F     T     F    F
 # G0W0* evGW* qsGW*
-  F    F    T
+  F    F    F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2

mol/benzene.xyz

@@ -1,5 +1,5 @@
 12
-Benzene,^1A_{1g},CC3,aug-cc-pVTZ
+
 C    0.00000000            1.39250319            0.00000000
 C   -1.20594314            0.69625160            0.00000000
 C   -1.20594314           -0.69625160            0.00000000

src/MBPT/G0F2.f90

@@ -1,5 +1,4 @@
-subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & 
+subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
-                linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform a one-shot second-order Green function calculation
 
@@ -13,8 +12,8 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
   logical,intent(in)            :: dBSE
   logical,intent(in)            :: dTDA
   logical,intent(in)            :: evDyn
-  logical,intent(in)            :: singlet_manifold
+  logical,intent(in)            :: singlet
-  logical,intent(in)            :: triplet_manifold
+  logical,intent(in)            :: triplet
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas
@@ -116,7 +115,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
 
   if(BSE) then
 
-    call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE)
+    call BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE)
 
   end if
 

src/MBPT/evGF2.f90

@@ -1,4 +1,4 @@
-subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold, &
+subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, &
                  linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform eigenvalue self-consistent second-order Green function calculation
@@ -16,8 +16,8 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold
   integer,intent(in)            :: maxSCF
   double precision,intent(in)   :: thresh
   integer,intent(in)            :: max_diis
-  logical,intent(in)            :: singlet_manifold
+  logical,intent(in)            :: singlet
-  logical,intent(in)            :: triplet_manifold
+  logical,intent(in)            :: triplet
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas
@@ -169,7 +169,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold
 
   if(BSE) then
 
-    call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE)
+    call BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE)
 
   end if
 

src/MBPT/print_qsGF2.f90

@@ -0,0 +1,122 @@
+subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigC,Z,EqsGF2,dipole)
+
+! Print one-electron energies and other stuff for qsGF2
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)                 :: nBas
+  integer,intent(in)                 :: nO
+  integer,intent(in)                 :: nSCF
+  double precision,intent(in)        :: ENuc
+  double precision,intent(in)        :: Conv
+  double precision,intent(in)        :: thresh
+  double precision,intent(in)        :: eHF(nBas)
+  double precision,intent(in)        :: eGF2(nBas)
+  double precision,intent(in)        :: c(nBas)
+  double precision,intent(in)        :: P(nBas,nBas) 
+  double precision,intent(in)        :: T(nBas,nBas),V(nBas,nBas)
+  double precision,intent(in)        :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas)
+  double precision,intent(in)        :: Z(nBas),SigC(nBas,nBas)
+  double precision,intent(in)        :: dipole(ncart)
+
+! Local variables
+
+  integer                            :: q,ixyz,HOMO,LUMO
+  double precision                   :: Gap,ET,EV,EJ,Ex,Ec
+  double precision,external          :: trace_matrix
+
+! Output variables
+
+  double precision,intent(out)       :: EqsGF2
+
+! HOMO and LUMO
+
+  HOMO = nO
+  LUMO = HOMO + 1
+  Gap = eGF2(LUMO)-eGF2(HOMO)
+
+! Compute energies
+
+  ET = trace_matrix(nBas,matmul(P,T))
+  EV = trace_matrix(nBas,matmul(P,V))
+  EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
+  Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
+  Ec = 0.50d0*trace_matrix(nBas,matmul(P,SigC))
+  EqsGF2 = ET + EV + EJ + Ex + Ec
+
+! Dump results
+
+  write(*,*)'-------------------------------------------------------------------------------'
+  if(nSCF < 10) then
+    write(*,'(1X,A21,I1,A1,I1,A12)')'  Self-consistent qsG',nSCF,'W',nSCF,' calculation'
+  else
+    write(*,'(1X,A21,I2,A1,I2,A12)')'  Self-consistent qsG',nSCF,'W',nSCF,' calculation'
+  endif
+  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)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|'
+  write(*,*)'-------------------------------------------------------------------------------'
+
+  do q=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)') &
+    '|',q,'|',eHF(q)*HaToeV,'|',SigC(q,q)*HaToeV,'|',Z(q),'|',eGF2(q)*HaToeV,'|'
+  enddo
+
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A10,I3)')   'Iteration ',nSCF
+  write(*,'(2X,A19,F15.5)')'max(|FPS - SPF|) = ',Conv
+  write(*,*)'-------------------------------------------'
+  write(*,'(2X,A30,F15.6,A3)') 'qsGF2 HOMO      energy:',eGF2(HOMO)*HaToeV,' eV'
+  write(*,'(2X,A30,F15.6,A3)') 'qsGF2 LUMO      energy:',eGF2(LUMO)*HaToeV,' eV'
+  write(*,'(2X,A30,F15.6,A3)') 'qsGF2 HOMO-LUMO gap   :',Gap*HaToeV,' eV'
+  write(*,*)'-------------------------------------------'
+  write(*,'(2X,A30,F15.6,A3)') '    qsGF2 total       energy:',EqsGF2 + ENuc,' au'
+  write(*,'(2X,A30,F15.6,A3)') '    qsGF2 exchange    energy:',Ex,' au'
+  write(*,*)'-------------------------------------------'
+  write(*,*)
+
+! Dump results for final iteration
+
+  if(Conv < thresh) then
+
+    write(*,*)
+    write(*,'(A50)')           '---------------------------------------'
+    write(*,'(A32)')           ' Summary              '
+    write(*,'(A50)')           '---------------------------------------'
+    write(*,'(A32,1X,F16.10,A3)') ' One-electron energy: ',ET + EV,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Kinetic      energy: ',ET,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Potential    energy: ',EV,' au'
+    write(*,'(A50)')           '---------------------------------------'
+    write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + Ex,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Hartree      energy: ',EJ,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Exchange     energy: ',Ex,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Correlation  energy: ',Ec,' au'
+    write(*,'(A50)')           '---------------------------------------'
+    write(*,'(A32,1X,F16.10,A3)') ' Electronic   energy: ',EqsGF2,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Nuclear   repulsion: ',ENuc,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' qsGF2        energy: ',ENuc + EqsGF2,' au'
+    write(*,'(A50)')           '---------------------------------------'
+    write(*,'(A35)')           ' Dipole moment (Debye)    '
+    write(*,'(10X,4A10)')      'X','Y','Z','Tot.'
+    write(*,'(10X,4F10.6)')    (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
+    write(*,'(A50)')           '-----------------------------------------'
+    write(*,*)
+ 
+    write(*,'(A50)') '---------------------------------------'
+    write(*,'(A32)') ' qsGF2 MO coefficients'
+    write(*,'(A50)') '---------------------------------------'
+    call matout(nBas,nBas,c)
+    write(*,*)
+    write(*,'(A50)') '---------------------------------------'
+    write(*,'(A32)') ' qsGF2 MO energies'
+    write(*,'(A50)') '---------------------------------------'
+    call matout(nBas,1,eGF2)
+    write(*,*)
+
+  endif
+
+
+end subroutine print_qsGF2

src/MBPT/qsGF2.f90

@@ -0,0 +1,222 @@
+subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, & 
+                eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
+                S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
+
+! Perform a quasiparticle self-consistent GF2 calculation
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: maxSCF
+  integer,intent(in)            :: max_diis
+  double precision,intent(in)   :: thresh
+  logical,intent(in)            :: BSE
+  logical,intent(in)            :: TDA
+  logical,intent(in)            :: dBSE
+  logical,intent(in)            :: dTDA
+  logical,intent(in)            :: evDyn
+  logical,intent(in)            :: singlet
+  logical,intent(in)            :: triplet
+  double precision,intent(in)   :: eta
+
+  integer,intent(in)            :: nNuc
+  double precision,intent(in)   :: ZNuc(nNuc)
+  double precision,intent(in)   :: rNuc(nNuc,ncart)
+  double precision,intent(in)   :: ENuc
+
+  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  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)   :: S(nBas,nBas)
+  double precision,intent(in)   :: T(nBas,nBas)
+  double precision,intent(in)   :: V(nBas,nBas)
+  double precision,intent(in)   :: Hc(nBas,nBas)
+  double precision,intent(in)   :: X(nBas,nBas)
+  double precision,intent(in)   :: ERI_AO(nBas,nBas,nBas,nBas)
+  double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: dipole_int_AO(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int_MO(nBas,nBas,ncart)
+
+! Local variables
+
+  integer                       :: nSCF
+  integer                       :: nBasSq
+  integer                       :: ispin
+  integer                       :: n_diis
+  double precision              :: EqsGF2
+  double precision              :: Conv
+  double precision              :: rcond
+  double precision,external     :: trace_matrix
+  double precision              :: dipole(ncart)
+  double precision              :: EcBSE(nspin)
+
+  double precision,allocatable  :: error_diis(:,:)
+  double precision,allocatable  :: F_diis(:,:)
+  double precision,allocatable  :: c(:,:)
+  double precision,allocatable  :: cp(:,:)
+  double precision,allocatable  :: eGF2(:)
+  double precision,allocatable  :: P(:,:)
+  double precision,allocatable  :: F(:,:)
+  double precision,allocatable  :: Fp(:,:)
+  double precision,allocatable  :: J(:,:)
+  double precision,allocatable  :: K(:,:)
+  double precision,allocatable  :: SigC(:,:)
+  double precision,allocatable  :: SigCp(:,:)
+  double precision,allocatable  :: SigCm(:,:)
+  double precision,allocatable  :: Z(:)
+  double precision,allocatable  :: error(:,:)
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'|       Self-consistent qsGF2 calculation      |'
+  write(*,*)'************************************************'
+  write(*,*)
+
+! Warning 
+
+  write(*,*) '!! ERIs in MO basis will be overwritten in qsGF2 !!'
+  write(*,*)
+
+! Stuff 
+
+  nBasSq = nBas*nBas
+
+! TDA 
+
+  if(TDA) then 
+    write(*,*) 'Tamm-Dancoff approximation activated!'
+    write(*,*)
+  end if
+
+! Memory allocation
+
+  allocate(eGF2(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas),        &
+           J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & 
+           error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
+
+! Initialization
+  
+  nSCF            = -1
+  n_diis          = 0
+  ispin           = 1
+  Conv            = 1d0
+  P(:,:)          = PHF(:,:)
+  eGF2(:)         = eHF(:)
+  c(:,:)          = cHF(:,:)
+  F_diis(:,:)     = 0d0
+  error_diis(:,:) = 0d0
+
+!------------------------------------------------------------------------
+! Main loop
+!------------------------------------------------------------------------
+
+  do while(Conv > thresh .and. nSCF <= maxSCF)
+
+    ! Increment
+
+    nSCF = nSCF + 1
+
+    ! Buid Coulomb matrix
+
+    call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J)
+
+    ! Compute exchange part of the self-energy 
+
+    call exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
+
+    ! AO to MO transformation of two-electron integrals
+
+    call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO)
+
+    ! Compute self-energy and renormalization factor
+
+    call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z)
+
+    ! Make correlation self-energy Hermitian and transform it back to AO basis
+   
+    SigCp = 0.5d0*(SigC + transpose(SigC))
+    SigCm = 0.5d0*(SigC - transpose(SigC))
+
+    call MOtoAO_transform(nBas,S,c,SigCp)
+ 
+    ! Solve the quasi-particle equation
+
+    F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:)
+
+    ! Compute commutator and convergence criteria
+
+    error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
+    Conv = maxval(abs(error))
+
+    ! DIIS extrapolation 
+
+    n_diis = min(n_diis+1,max_diis)
+    call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
+
+!   Reset DIIS if required
+
+    if(abs(rcond) < 1d-15) n_diis = 0
+
+    ! Diagonalize Hamiltonian in AO basis
+
+    Fp = matmul(transpose(X),matmul(F,X))
+    cp(:,:) = Fp(:,:)
+    call diagonalize_matrix(nBas,cp,eGF2)
+    c = matmul(X,cp)
+
+    ! Compute new density matrix in the AO basis
+
+    P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
+
+    ! Print results
+
+    call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
+    call print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigCp,Z,EqsGF2,dipole)
+
+  enddo
+!------------------------------------------------------------------------
+! End main loop
+!------------------------------------------------------------------------
+
+! Did it actually converge?
+
+  if(nSCF == maxSCF+1) then
+
+    write(*,*)
+    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+    write(*,*)'                 Convergence failed                 '
+    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+    write(*,*)
+
+    stop
+
+  endif
+
+! Deallocate memory
+
+  deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,error,error_diis,F_diis)
+
+! Perform BSE calculation
+
+  if(BSE) then
+
+    call BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGF2,EcBSE)
+
+    write(*,*)
+    write(*,*)'-------------------------------------------------------------------------------'
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy (singlet) =',EcBSE(1)
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy (triplet) =',EcBSE(2)
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy           =',EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 total energy                 =',ENuc + EqsGF2 + EcBSE(1) + EcBSE(2)
+    write(*,*)'-------------------------------------------------------------------------------'
+    write(*,*)
+
+  end if
+
+end subroutine qsGF2

src/MBPT/self_energy_GF2.f90

@@ -0,0 +1,73 @@
+subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
+
+! Compute GF2 self-energy and its renormalization factor
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: eGF2(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer                       :: i,j,a,b
+  integer                       :: p,q
+  double precision              :: eps
+  double precision              :: num
+
+! Output variables
+
+  double precision,intent(out)  :: SigC(nBas,nBas)
+  double precision,intent(out)  :: Z(nBas)
+
+! Initialize 
+
+  SigC(:,:) = 0d0
+  Z(:)     = 0d0
+
+! Compute GF2 self-energy
+
+  do p=nC+1,nBas-nR
+    do q=nC+1,nBas-nR
+      do i=nC+1,nO
+        do j=nC+1,nO
+          do a=nO+1,nBas-nR
+
+            eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
+            num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j)
+
+            SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2)
+            if(p == q) Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  do p=nC+1,nBas-nR
+    do q=nC+1,nBas-nR
+      do i=nC+1,nO
+        do a=nO+1,nBas-nR
+          do b=nO+1,nBas-nR
+
+            eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
+            num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b)
+
+            SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2)
+            if(p == q) Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  Z(:) = 1d0/(1d0 - Z(:))
+
+end subroutine self_energy_GF2

src/QuAcK/QuAcK.f90

@@ -14,7 +14,7 @@ program QuAcK
   logical                       :: doCIS,doCIS_D,doCID,doCISD
   logical                       :: doRPA,doRPAx,doppRPA
   logical                       :: doADC
-  logical                       :: doG0F2,doevGF2,doG0F3,doevGF3
+  logical                       :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
   logical                       :: doG0W0,doevGW,doqsGW
   logical                       :: doG0T0,doevGT,doqsGT
   logical                       :: doMCMP2,doMinMCMP2
@@ -164,7 +164,8 @@ program QuAcK
                     do_drCCD,do_rCCD,do_lCCD,do_pCCD, &
                     doCIS,doCIS_D,doCID,doCISD,       & 
                     doRPA,doRPAx,doppRPA,             &
-                    doG0F2,doevGF2,doG0F3,doevGF3,    &
+                    doG0F2,doevGF2,doqsGF2,           & 
+                    doG0F3,doevGF3,                   &
                     doG0W0,doevGW,doqsGW,             &
                     doG0T0,doevGT,doqsGT,             &
                     doMCMP2)
@@ -839,6 +840,25 @@ program QuAcK
 
   end if
 
+!------------------------------------------------------------------------
+! Perform qsGF2 calculation
+!------------------------------------------------------------------------
+
+  if(doqsGF2) then 
+
+    call cpu_time(start_GF2)
+
+    call qsGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GF,nNuc,ZNuc,rNuc,ENuc, & 
+               nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
+
+    call cpu_time(end_GF2)
+
+    t_GF2 = end_GF2 - start_GF2
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF2,' seconds'
+    write(*,*)
+
+  end if
+
 !------------------------------------------------------------------------
 ! Compute G0F3 electronic binding energies
 !------------------------------------------------------------------------

src/QuAcK/read_methods.f90

@@ -4,7 +4,8 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM,           &
                         do_drCCD,do_rCCD,do_lCCD,do_pCCD, &
                         doCIS,doCIS_D,doCID,doCISD,       & 
                         doRPA,doRPAx,doppRPA,             & 
-                        doG0F2,doevGF2,doG0F3,doevGF3,    & 
+                        doG0F2,doevGF2,doqsGF2,           &
+                        doG0F3,doevGF3,                   & 
                         doG0W0,doevGW,doqsGW,             & 
                         doG0T0,doevGT,doqsGT,             & 
                         doMCMP2)
@@ -21,7 +22,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM,           &
   logical,intent(out)           :: do_drCCD,do_rCCD,do_lCCD,do_pCCD
   logical,intent(out)           :: doCIS,doCIS_D,doCID,doCISD
   logical,intent(out)           :: doRPA,doRPAx,doppRPA
-  logical,intent(out)           :: doG0F2,doevGF2,doG0F3,doevGF3  
+  logical,intent(out)           :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3  
   logical,intent(out)           :: doG0W0,doevGW,doqsGW
   logical,intent(out)           :: doG0T0,doevGT,doqsGT
   logical,intent(out)           :: doMCMP2
@@ -66,6 +67,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM,           &
 
   doG0F2  = .false.
   doevGF2 = .false.
+  doqsGF2 = .false.
   doG0F3  = .false.
   doevGF3 = .false.
 
@@ -132,11 +134,12 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM,           &
 ! Read Green function methods
 
   read(1,*) 
-  read(1,*) answer1,answer2,answer3,answer4
+  read(1,*) answer1,answer2,answer3,answer4,answer5
   if(answer1 == 'T') doG0F2  = .true.
   if(answer2 == 'T') doevGF2 = .true.
-  if(answer3 == 'T') doG0F3  = .true.
+  if(answer3 == 'T') doqsGF2 = .true.
-  if(answer4 == 'T') doevGF3 = .true.
+  if(answer4 == 'T') doG0F3  = .true.
+  if(answer5 == 'T') doevGF3 = .true.
 
 ! Read GW methods