stability analysis

input/methods

@@ -13,7 +13,7 @@
 # G0F2 evGF2 G0F3 evGF3
   F    F     F    F
 # G0W0* evGW* qsGW*
-  T    F    T
+  F    F    F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2

input/options

@@ -1,5 +1,5 @@
-# HF: maxSCF thresh   DIIS n_diis guess_type ortho_type mix_guess
+# HF: maxSCF thresh   DIIS n_diis guess_type ortho_type mix_guess stability
-       128    0.0000001  T    5      1          1         F
+       128    0.0000001  T    5      1          1         T         T
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
@@ -9,7 +9,7 @@
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm 
       256    0.00001 T    5      T   0.0  3      
 # GW/GT: maxSCF thresh  DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 
-         256   0.0000001   T    5     T  0.00367493  F      F     F     F   F  
+         256   0.00001   T    5     T  0.0  F      F     F     F   F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn

mol/h2.xyz

@@ -1,4 +1,4 @@
 2
 
 H 0.0 0.0 0.0
-H 0.0 0.0 0.7408481486
+H 0.0 0.0 1.4

src/HF/RHF_stability.f90

@@ -0,0 +1,147 @@
+subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
+
+! Perform a stability analysis of the RHF solution
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer,parameter             :: maxS = 20
+  integer                       :: ia
+  integer                       :: ispin
+
+  double precision,allocatable  :: A(:,:)
+  double precision,allocatable  :: B(:,:)
+  double precision,allocatable  :: AB(:,:)
+  double precision,allocatable  :: Om(:)
+
+! Memory allocation
+
+  allocate(A(nS,nS),B(nS,nS),AB(nS,nS),Om(nS))
+ 
+!-------------------------------------------------------------!
+! Stability analysis: Real RHF -> Real RHF  
+!-------------------------------------------------------------!
+ 
+  ispin = 1
+
+  call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
+  call linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
+
+  AB(:,:) = A(:,:) + B(:,:)
+
+  call diagonalize_matrix(nS,AB,Om)
+
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)'|       Stability analysis: Real RHF -> Real RHF            |'
+  write(*,*)'-------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
+            '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
+  write(*,*)'-------------------------------------------------------------'
+  do ia=1,min(nS,maxS)
+    write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
+      '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|'
+  enddo
+  write(*,*)'-------------------------------------------------------------'
+
+  if(minval(Om(:)) < 0d0) then
+
+    write(*,'(1X,A40,1X)')        'Too bad, RHF solution is unstable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om(1),' au'
+
+  else
+
+    write(*,'(1X,A40,1X)')        'Well done, RHF solution is stable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au'
+
+  end if
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)
+
+!-------------------------------------------------------------!
+! Stability analysis: Real RHF -> Complex RHF  
+!-------------------------------------------------------------!
+ 
+  AB(:,:) = A(:,:) - B(:,:)
+
+  call diagonalize_matrix(nS,AB,Om)
+
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)'|       Stability analysis: Real RHF -> Complex RHF         |'
+  write(*,*)'-------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
+            '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
+  write(*,*)'-------------------------------------------------------------'
+  do ia=1,min(nS,maxS)
+    write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
+      '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|'
+  enddo
+  write(*,*)'-------------------------------------------------------------'
+
+  if(minval(Om(:)) < 0d0) then
+
+    write(*,'(1X,A40,1X)')        'Too bad, RHF solution is unstable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om(1),' au'
+
+  else
+
+    write(*,'(1X,A40,1X)')        'Well done, RHF solution is stable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au'
+
+  end if
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)
+
+!-------------------------------------------------------------!
+! Stability analysis: Real RHF -> Real UHF  
+!-------------------------------------------------------------!
+ 
+  ispin = 2
+
+  call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
+  call linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
+
+  AB(:,:) = A(:,:) + B(:,:)
+
+  call diagonalize_matrix(nS,AB,Om)
+
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)'|       Stability analysis: Real RHF -> Real UHF            |'
+  write(*,*)'-------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
+            '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
+  write(*,*)'-------------------------------------------------------------'
+  do ia=1,min(nS,maxS)
+    write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
+      '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|'
+  enddo
+  write(*,*)'-------------------------------------------------------------'
+
+  if(minval(Om(:)) < 0d0) then
+
+    write(*,'(1X,A40,1X)')        'Too bad, RHF solution is unstable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om(1),' au'
+
+  else
+
+    write(*,'(1X,A40,1X)')        'Well done, RHF solution is stable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au'
+
+  end if
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)
+    
+end subroutine RHF_stability

src/HF/UHF_stability.f90

@@ -0,0 +1,174 @@
+subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb)
+
+! Perform a stability analysis of the UHF solution
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  integer,intent(in)            :: nS(nspin)
+
+  double precision,intent(in)   :: eHF(nBas,nspin)
+  double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer,parameter             :: maxS = 20
+  integer                       :: ia
+  integer                       :: ispin
+
+  integer                       :: nS_aa,nS_bb,nS_sc
+  double precision,allocatable  :: Om_sc(:)
+  double precision,allocatable  :: A_sc(:,:)
+  double precision,allocatable  :: B_sc(:,:)
+  double precision,allocatable  :: AB_sc(:,:)
+
+  integer                       :: nS_ab,nS_ba,nS_sf
+  double precision,allocatable  :: Om_sf(:)
+  double precision,allocatable  :: A_sf(:,:)
+  double precision,allocatable  :: B_sf(:,:)
+  double precision,allocatable  :: AB_sf(:,:)
+
+! Menory allocation
+
+  nS_aa = nS(1)
+  nS_bb = nS(2)
+  nS_sc = nS_aa + nS_bb
+
+  allocate(Om_sc(nS_sc),A_sc(nS_sc,nS_sc),B_sc(nS_sc,nS_sc),AB_sc(nS_sc,nS_sc))
+ 
+!-------------------------------------------------------------!
+! Stability analysis: Real UHF -> Real UHF  
+!-------------------------------------------------------------!
+ 
+  ispin = 1
+
+  call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eHF, &
+                                             ERI_aaaa,ERI_aabb,ERI_bbbb,A_sc)
+  call unrestricted_linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, &
+                                               ERI_aaaa,ERI_aabb,ERI_bbbb,B_sc)
+
+  AB_sc(:,:) = A_sc(:,:) + B_sc(:,:)
+
+  call diagonalize_matrix(nS_sc,AB_sc,Om_sc)
+
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)'|       Stability analysis: Real UHF -> Real UHF            |'
+  write(*,*)'-------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
+            '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
+  write(*,*)'-------------------------------------------------------------'
+  do ia=1,min(nS_sc,maxS)
+    write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
+      '|',ia,'|',Om_sc(ia),'|',Om_sc(ia)*HaToeV,'|'
+  enddo
+  write(*,*)'-------------------------------------------------------------'
+
+  if(minval(Om_sc(:)) < 0d0) then
+
+    write(*,'(1X,A40,1X)')        'Too bad, UHF solution is unstable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om_sc(1),' au'
+
+  else
+
+    write(*,'(1X,A40,1X)')        'Well done, UHF solution is stable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om_sc(1),' au'
+
+  end if
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)
+
+!-------------------------------------------------------------!
+! Stability analysis: Real UHF -> Complex UHF  
+!-------------------------------------------------------------!
+ 
+  AB_sc(:,:) = A_sc(:,:) - B_sc(:,:)
+
+  call diagonalize_matrix(nS_sc,AB_sc,Om_sc)
+
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)'|       Stability analysis: Real UHF -> Complex UHF         |'
+  write(*,*)'-------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
+            '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
+  write(*,*)'-------------------------------------------------------------'
+  do ia=1,min(nS_sc,maxS)
+    write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
+      '|',ia,'|',Om_sc(ia),'|',Om_sc(ia)*HaToeV,'|'
+  enddo
+  write(*,*)'-------------------------------------------------------------'
+
+  if(minval(Om_sc(:)) < 0d0) then
+
+    write(*,'(1X,A40,1X)')        'Too bad, UHF solution is unstable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om_sc(1),' au'
+
+  else
+
+    write(*,'(1X,A40,1X)')        'Well done, UHF solution is stable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om_sc(1),' au'
+
+  end if
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)
+
+! Menory (de)allocation
+
+
+  nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2))
+  nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1))
+  nS_sf = nS_ab + nS_ba
+
+  deallocate(Om_sc,A_sc,B_sc,AB_sc)
+  allocate(Om_sf(nS_sf),A_sf(nS_sf,nS_sf),B_sf(nS_sf,nS_sf),AB_sf(nS_sf,nS_sf))
+
+!-------------------------------------------------------------!
+! Stability analysis: Real UHF -> Real GHF  
+!-------------------------------------------------------------!
+ 
+  ispin = 2
+
+  call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,1d0,eHF, &
+                                             ERI_aaaa,ERI_aabb,ERI_bbbb,A_sf)
+  call unrestricted_linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,1d0, &
+                                               ERI_aaaa,ERI_aabb,ERI_bbbb,B_sf)
+
+  AB_sf(:,:) = A_sf(:,:) + B_sf(:,:)
+
+  call diagonalize_matrix(nS_sf,AB_sf,Om_sf)
+
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)'|       Stability analysis: Real UHF -> Real GHF            |'
+  write(*,*)'-------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
+            '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
+  write(*,*)'-------------------------------------------------------------'
+  do ia=1,min(nS_sf,maxS)
+    write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
+      '|',ia,'|',Om_sf(ia),'|',Om_sf(ia)*HaToeV,'|'
+  enddo
+  write(*,*)'-------------------------------------------------------------'
+
+  if(minval(Om_sf(:)) < 0d0) then
+
+    write(*,'(1X,A40,1X)')        'Too bad, UHF solution is unstable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om_sf(1),' au'
+
+  else
+
+    write(*,'(1X,A40,1X)')        'Well done, UHF solution is stable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om_sf(1),' au'
+
+  end if
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)
+    
+end subroutine UHF_stability

src/MBPT/print_qsUGW.f90

@@ -113,6 +113,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
 ! Dump results for final iteration
 
   if(Conv < thresh) then
+! if(.true.) then
 
     write(*,*)
     write(*,'(A60)')              '-------------------------------------------------'

src/MP/MP2.f90

@@ -6,9 +6,15 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
 
 ! Input variables
 
-  integer,intent(in)            :: nBas,nC,nO,nV,nR
+  integer,intent(in)            :: nBas
-  double precision,intent(in)   :: ENuc,EHF
+  integer,intent(in)            :: nC
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas),e(nBas)
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: EHF
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
 
 ! Local variables
 

src/QuAcK/QuAcK.f90

@@ -6,6 +6,7 @@ program QuAcK
   logical                       :: doSph
   logical                       :: unrestricted = .false.
   logical                       :: doRHF,doUHF,doMOM 
+  logical                       :: dostab
   logical                       :: doKS
   logical                       :: doMP2,doMP3,doMP2F12
   logical                       :: doCCD,doDCD,doCCSD,doCCSDT
@@ -78,6 +79,7 @@ program QuAcK
   double precision              :: start_QuAcK  ,end_QuAcK    ,t_QuAcK
   double precision              :: start_int    ,end_int      ,t_int  
   double precision              :: start_HF     ,end_HF       ,t_HF
+  double precision              :: start_stab   ,end_stab     ,t_stab
   double precision              :: start_KS     ,end_KS       ,t_KS
   double precision              :: start_MOM    ,end_MOM      ,t_MOM
   double precision              :: start_AOtoMO ,end_AOtoMO   ,t_AOtoMO
@@ -169,14 +171,14 @@ program QuAcK
 
 ! Read options for methods
 
-  call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix, &
+  call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, &
-                    maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC,                           &
+                    maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC,                                  &
-                    TDA,singlet,triplet,spin_conserved,spin_flip,                    &
+                    TDA,singlet,triplet,spin_conserved,spin_flip,                           &
-                    maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,     &
+                    maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,            &
-                    maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,              & 
+                    maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,                     & 
-                    COHSEX,SOSEX,TDA_W,G0W,GW0,                                      &  
+                    COHSEX,SOSEX,TDA_W,G0W,GW0,                                             &  
-                    doACFDT,exchange_kernel,doXBS,                                   &
+                    doACFDT,exchange_kernel,doXBS,                                          &
-                    BSE,dBSE,dTDA,evDyn,                                             &
+                    BSE,dBSE,dTDA,evDyn,                                                    &
                     nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
 
 ! Weird stuff
@@ -439,6 +441,32 @@ program QuAcK
   write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds'
   write(*,*)
 
+!------------------------------------------------------------------------
+! Stability analysis of HF solution
+!------------------------------------------------------------------------
+
+  if(dostab) then
+
+    call cpu_time(start_stab)
+
+    if(unrestricted) then
+
+      call UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb)
+
+    else
+
+      call RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_MO)
+
+    end if
+
+    call cpu_time(end_stab)
+
+    t_stab = end_stab - start_stab
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds'
+    write(*,*)
+
+  end if
+
 !------------------------------------------------------------------------
 ! Compute MP2 energy
 !------------------------------------------------------------------------

src/QuAcK/read_options.f90

@@ -1,11 +1,11 @@
-subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix, &
+subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, &
-                        maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC,                           &
+                        maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC,                                  &
-                        TDA,singlet,triplet,spin_conserved,spin_flip,                    &
+                        TDA,singlet,triplet,spin_conserved,spin_flip,                           &
-                        maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,     &
+                        maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,            &
-                        maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,              &
+                        maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,                     &
-                        COHSEX,SOSEX,TDA_W,G0W,GW0,                                      &
+                        COHSEX,SOSEX,TDA_W,G0W,GW0,                                             &
-                        doACFDT,exchange_kernel,doXBS,                                   &
+                        doACFDT,exchange_kernel,doXBS,                                          &
-                        BSE,dBSE,dTDA,evDyn,                                             &
+                        BSE,dBSE,dTDA,evDyn,                                                    &
                         nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
 
 ! Read desired methods 
@@ -21,6 +21,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
   integer,intent(out)           :: guess_type
   integer,intent(out)           :: ortho_type
   logical,intent(out)           :: mix
+  logical,intent(out)           :: dostab
 
   integer,intent(out)           :: maxSCF_CC
   double precision,intent(out)  :: thresh_CC
@@ -87,12 +88,14 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
   guess_type   = 1
   ortho_type   = 1
   mix          = .false.
+  dostab       = .false.
 
   read(1,*) 
-  read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2
+  read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,answer3
 
   if(answer1 == 'T') DIIS_HF = .true.
   if(answer2 == 'T') mix     = .true.
+  if(answer3 == 'T') dostab  = .true.
 
   if(.not.DIIS_HF) n_diis_HF = 1