print GHF

input/methods

@@ -1,19 +1,18 @@
 # RHF UHF GHF ROHF 
-  F   T   F   F   
+  F   F   T   F   
-# MP2* MP3 
+# MP2 MP3 
   F   F   
 # CCD pCCD DCD CCSD CCSD(T) 
   F   F    F   F    F
 # drCCD rCCD crCCD lCCD
   F     F    F     F
-# CIS* CIS(D) CID CISD FCI
+# CIS CIS(D) CID CISD FCI
-  F    F      F   F    F
+  F   F      F   F    F
-# phRPA* phRPAx* crRPA ppRPA 
+# phRPA phRPAx crRPA ppRPA 
-  F      F       F     F 
+  F     F      F     F 
-# G0F2* evGF2* qsGF2* G0F3 evGF3
+# G0F2 evGF2 qsGF2 G0F3 evGF3
-  F     F      F      F    F
+  F    F     F     F    F
-# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
+# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
-  T     F     F     F        F      F
+  F    F    F    F        F      F
-# G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh
+# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
-  F       F       F       F      F      F
+  F      F      F      F      F      F
-# * unrestricted version available

input/options

@@ -1,5 +1,5 @@
 # HF: maxSCF thresh    DIIS guess mix_guess level_shift stability
-      10000   0.000001  5     2    0.0       0.0         T
+      1000  0.0000001  5     3    0.0       0.0         T
 # MP: reg
       F
 # CC: maxSCF thresh     DIIS

src/HF/GHF.f90

@@ -42,7 +42,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
   double precision              :: ET,ETaa,ETbb
   double precision              :: EV,EVaa,EVbb
   double precision              :: EJ,EJaaaa,EJaabb,EJbbaa,EJbbbb
-  double precision              :: Ex,Exaaaa,Exabba,Exbaab,Exbbbb
+  double precision              :: EK,EKaaaa,EKabba,EKbaab,EKbbbb
   double precision              :: dipole(ncart)
 
   double precision,allocatable  :: Jaa(:,:),Jbb(:,:)
@@ -139,7 +139,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
   write(*,*)
   write(*,*)'-----------------------------------------------------------------------------'
   write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & 
-            '|','#','|','E(GHF)','|','EJ(GHF)','|','Ex(GHF)','|','Conv','|'
+            '|','#','|','E(GHF)','|','EJ(GHF)','|','EK(GHF)','|','Conv','|'
   write(*,*)'-----------------------------------------------------------------------------'
   
   do while(Conv > thresh .and. nSCF < maxSCF)
@@ -239,7 +239,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 
    EV = EVaa + EVbb
 
-!  Hartree energy: 16 terms?
+!  Hartree energy
 
     EJaaaa = 0.5d0*trace_matrix(nBas,matmul(Paa,Jaa))
     EJaabb = 0.5d0*trace_matrix(nBas,matmul(Paa,Jbb))
@@ -250,21 +250,21 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
   
 !   Exchange energy
 
-    Exaaaa = 0.5d0*trace_matrix(nBas,matmul(Paa,Kaa))
+    EKaaaa = 0.5d0*trace_matrix(nBas,matmul(Paa,Kaa))
-    Exabba = 0.5d0*trace_matrix(nBas,matmul(Pab,Kba))
+    EKabba = 0.5d0*trace_matrix(nBas,matmul(Pab,Kba))
-    Exbaab = 0.5d0*trace_matrix(nBas,matmul(Pba,Kab))
+    EKbaab = 0.5d0*trace_matrix(nBas,matmul(Pba,Kab))
-    Exbbbb = 0.5d0*trace_matrix(nBas,matmul(Pbb,Kbb))
+    EKbbbb = 0.5d0*trace_matrix(nBas,matmul(Pbb,Kbb))
 
-    Ex = Exaaaa + Exabba + Exbaab + Exbbbb
+    EK = EKaaaa + EKabba + EKbaab + EKbbbb
 
 !   Total energy
 
-    EHF = ET + EV + EJ + Ex
+    EHF = ET + EV + EJ + EK
 
 !   Dump results
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & 
-      '|',nSCF,'|',EHF + ENuc,'|',EJ,'|',Ex,'|',Conv,'|'
+      '|',nSCF,'|',EHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|'
  
   end do
   write(*,*)'-----------------------------------------------------------------------------'
@@ -292,6 +292,6 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 
 ! Compute final GHF energy
 
-! call print_GHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
+  call print_GHF(nBas2,nO,e,C,ENuc,ET,EV,EJ,EK,EHF,dipole)
 
 end subroutine 

src/HF/print_GHF.f90

@@ -0,0 +1,77 @@
+subroutine print_GHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,EGHF,dipole)
+
+! Print one-electron energies and other stuff for GHF
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)                 :: nBas
+  integer,intent(in)                 :: nO
+  double precision,intent(in)        :: eHF(nBas)
+  double precision,intent(in)        :: cHF(nBas,nBas)
+  double precision,intent(in)        :: ENuc
+  double precision,intent(in)        :: ET
+  double precision,intent(in)        :: EV
+  double precision,intent(in)        :: EJ
+  double precision,intent(in)        :: EK
+  double precision,intent(in)        :: EGHF
+  double precision,intent(in)        :: dipole(ncart)
+
+! Local variables
+
+  integer                            :: ixyz
+  integer                            :: HOMO
+  integer                            :: LUMO
+  double precision                   :: Gap
+
+! HOMO and LUMO
+
+  HOMO = nO
+  LUMO = HOMO + 1
+  Gap = eHF(LUMO)-eHF(HOMO)
+
+! Dump results
+
+
+  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 + EK,' au'
+  write(*,'(A32,1X,F16.10,A3)') ' Hartree      energy: ',EJ,' au'
+  write(*,'(A32,1X,F16.10,A3)') ' Exchange     energy: ',EK,' au'
+  write(*,'(A50)')           '-----------------------------------------'
+  write(*,'(A32,1X,F16.10,A3)') ' Electronic   energy: ',EGHF,' au'
+  write(*,'(A32,1X,F16.10,A3)') ' Nuclear   repulsion: ',ENuc,' au'
+  write(*,'(A32,1X,F16.10,A3)') ' GHF          energy: ',EGHF + ENuc,' au'
+  write(*,'(A50)')           '-----------------------------------------'
+  write(*,'(A32,1X,F16.6,A3)')  ' HF HOMO      energy: ',eHF(HOMO)*HaToeV,' eV'
+  write(*,'(A32,1X,F16.6,A3)')  ' HF LUMO      energy: ',eHF(LUMO)*HaToeV,' eV'
+  write(*,'(A32,1X,F16.6,A3)')  ' HF HOMO-LUMO gap   : ',Gap*HaToeV,' eV'
+  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(*,*)
+
+! Print results
+
+  write(*,'(A50)')  '---------------------------------------'
+  write(*,'(A32)') 'MO coefficients'
+  write(*,'(A50)')  '---------------------------------------'
+  call matout(nBas,nBas,cHF)
+  write(*,*)
+  write(*,'(A50)') '---------------------------------------'
+  write(*,'(A32)') 'MO energies'
+  write(*,'(A50)') '---------------------------------------'
+  call matout(nBas,1,eHF)
+  write(*,*)
+
+end subroutine