fix bug in HF diis

input/methods

@@ -13,7 +13,7 @@
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F      F      F    F
 # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
-  F     T     F     F        F      F
+  F     F     F     F        F      F
 # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
-  F      F      F      F      F      F
+  F      F      F      F      T      F
 # * unrestricted version available

src/GT/GTeh_self_energy_diag.f90

@@ -47,6 +47,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
 
         eps = e(p) - e(i) + Om(m)
         num = rhoL(i,p,m)*rhoR(i,p,m)
+!       num = 0.5d0*rhoL(p,i,m)*rhoL(p,i,m)
         Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
         Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
@@ -62,6 +63,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
 
         eps = e(p) - e(a) - Om(m)
         num = rhoL(p,a,m)*rhoR(p,a,m)
+!       num = 0.5d0*rhoL(p,a,m)*rhoL(p,a,m)
         Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
         Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 

src/GT/evGTeh.f90

@@ -69,8 +69,6 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
   double precision,allocatable  :: XmY(:,:)
   double precision,allocatable  :: rhoL(:,:,:)
   double precision,allocatable  :: rhoR(:,:,:)
-  
-  double precision,allocatable  :: eGTlin(:)
 
 ! Hello world
 
@@ -102,7 +100,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 ! Memory allocation
 
   allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),Z(nBas),Sig(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & 
-           rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eGTlin(nBas))
+           rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
 
 ! Initialization
 
@@ -149,23 +147,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 
     ! Solve the quasi-particle equation
 
-    eGTlin(:) = eHF(:) + Sig(:)
+    eGT(:) = eHF(:) + Sig(:)
-
-    ! Linearized or graphical solution?
-
-    if(linearize) then 
- 
-       write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
-       write(*,*)
-
-       eGT(:) = eGTlin(:)
-
-    else 
-
-!      write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
-!      write(*,*)
-
-    end if
 
     ! Convergence criteria
 

src/GW/GW_self_energy.f90

@@ -21,8 +21,7 @@ subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
 ! Local variables
 
   integer                       :: i,j,a,b
-  integer                       :: p,q,r
+  integer                       :: p,q,m
-  integer                       :: jb
   double precision              :: num,eps
 
 ! Output variables
@@ -44,16 +43,16 @@ subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
 
 !$OMP PARALLEL &
 !$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Om) &
-!$OMP PRIVATE(jb,i,q,p,eps,num) &
+!$OMP PRIVATE(m,i,q,p,eps,num) &
 !$OMP DEFAULT(NONE)
 !$OMP DO
   do q=nC+1,nBas-nR
      do p=nC+1,nBas-nR
-        do jb=1,nS
+        do m=1,nS
            do i=nC+1,nO
  
-              eps = e(p) - e(i) + Om(jb)
+              eps = e(p) - e(i) + Om(m)
-              num = 2d0*rho(p,i,jb)*rho(q,i,jb)
+              num = 2d0*rho(p,i,m)*rho(q,i,m)
               Sig(p,q) = Sig(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
  
@@ -68,16 +67,16 @@ subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
 
 !$OMP PARALLEL &
 !$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Om) &
-!$OMP PRIVATE(jb,a,q,p,eps,num) &
+!$OMP PRIVATE(m,a,q,p,eps,num) &
 !$OMP DEFAULT(NONE)
 !$OMP DO  
   do q=nC+1,nBas-nR
      do p=nC+1,nBas-nR
-        do jb=1,nS
+        do m=1,nS
            do a=nO+1,nBas-nR
  
-              eps = e(p) - e(a) - Om(jb)
+              eps = e(p) - e(a) - Om(m)
-              num = 2d0*rho(p,a,jb)*rho(q,a,jb)
+              num = 2d0*rho(p,a,m)*rho(q,a,m)
               Sig(p,q) = Sig(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
  
@@ -91,12 +90,12 @@ subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
 ! Galitskii-Migdal correlation energy
 
   EcGM = 0d0
-  do jb=1,nS
+  do m=1,nS
     do a=nO+1,nBas-nR
       do i=nC+1,nO
 
-        eps = e(a) - e(i) + Om(jb)
+        eps = e(a) - e(i) + Om(m)
-        num = 4d0*rho(a,i,jb)*rho(a,i,jb)
+        num = 4d0*rho(a,i,m)*rho(a,i,m)
         EcGM = EcGM - num*eps/(eps**2 + eta**2)
 
       end do

src/GW/GW_self_energy_diag.f90

@@ -20,7 +20,7 @@ subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
 
 ! Local variables
 
-  integer                       :: i,a,p,q,jb
+  integer                       :: i,a,p,m
   double precision              :: num,eps
 
 ! Output variables
@@ -42,10 +42,10 @@ subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
 
   do p=nC+1,nBas-nR
     do i=nC+1,nO
-      do jb=1,nS
+      do m=1,nS
 
-        eps = e(p) - e(i) + Om(jb)
+        eps = e(p) - e(i) + Om(m)
-        num = 2d0*rho(p,i,jb)**2
+        num = 2d0*rho(p,i,m)**2
         Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
         Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
@@ -57,10 +57,10 @@ subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
 
   do p=nC+1,nBas-nR
     do a=nO+1,nBas-nR
-      do jb=1,nS
+      do m=1,nS
 
-        eps = e(p) - e(a) - Om(jb)
+        eps = e(p) - e(a) - Om(m)
-        num = 2d0*rho(p,a,jb)**2
+        num = 2d0*rho(p,a,m)**2
         Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
         Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
@@ -73,10 +73,10 @@ subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
   EcGM = 0d0
   do i=nC+1,nO
     do a=nO+1,nBas-nR
-      do jb=1,nS
+      do m=1,nS
 
-        eps = e(a) - e(i) + Om(jb)
+        eps = e(a) - e(i) + Om(m)
-        num = 4d0*rho(a,i,jb)**2
+        num = 4d0*rho(a,i,m)**2
         EcGM = EcGM - num*eps/(eps**2 + eta**2)
 
       end do

src/HF/HF.f90

@@ -39,7 +39,6 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues
   double precision              :: start_HF     ,end_HF       ,t_HF
 
   integer                       :: nSCF
-  integer                       :: n_diis
   double precision              :: ET
   double precision              :: EV
   double precision              :: EJ
@@ -83,7 +82,7 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues
     end if
 
     call wall_time(start_HF)
-    call RHF(maxSCF,thresh,n_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
+    call RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
              nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
     call wall_time(end_HF)
 
@@ -103,7 +102,7 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues
     unrestricted = .true.
 
     call wall_time(start_HF)
-    call UHF(maxSCF,thresh,n_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
+    call UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
              nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
     call wall_time(end_HF)
 
@@ -127,7 +126,7 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues
       stop
     end if
 
-!   call RMOM(maxSCF,thresh,n_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, &
+!   call RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, &
 !           nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
 
   end if
@@ -141,7 +140,7 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues
     ! Switch on the unrestricted flag
     unrestricted = .true.
 
-!   call UMOM(maxSCF,thresh,n_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, &
+!   call UMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, &
 !           nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
 
   end if

src/QuAcK/QuAcK.f90

@@ -479,7 +479,7 @@ program QuAcK
     
     call cpu_time(start_GT)
     call GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,  &
-            exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet,spin_conserved,spin_flip, &
+            exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,spin_conserved,spin_flip, &
             linGT,eta_GT,regGT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,                                   &
             ERI_AO,ERI_MO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_AO,dipole_int_MO,dipole_int_aa,dipole_int_bb,   & 
             PHF,cHF,epsHF)