guess mix in eDFT

input/options

@@ -1,5 +1,5 @@
 # HF: maxSCF thresh   DIIS n_diis guess_type ortho_type mix_guess
-       128    0.000001  T    5      1          1         F
+       128    0.000001  T    5      1          1         T
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
@@ -13,6 +13,6 @@
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn
-        T    T    T    F
+        F    T    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T

mol/h2.xyz

@@ -1,4 +1,4 @@
 2
 
 H 0.0 0.0 0.0
-H 0.0 0.0 0.741
+H 0.0 0.0 2.645875

src/LR/linear_response_B_pp.f90

@@ -86,4 +86,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp)
 
   end if
 
+  print*,'B pp-matrix'
+  call matout(nVV,nOO,B_pp)
+
 end subroutine linear_response_B_pp

src/LR/linear_response_C_pp.f90

@@ -96,4 +96,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
 
   end if
 
+  print*,'C pp-matrix'
+  call matout(nVV,nVV,C_pp)
+
 end subroutine linear_response_C_pp

src/LR/linear_response_D_pp.f90

@@ -96,4 +96,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
 
   end if
 
+  print*,'D pp-matrix'
+  call matout(nOO,nOO,D_pp)
+
 end subroutine linear_response_D_pp

src/QuAcK/QuAcK.f90

@@ -316,7 +316,7 @@ program QuAcK
   if(doKS) then
 
     call cpu_time(start_KS)
-    call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, & 
+    call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, & 
               nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
               max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO)
              

src/eDFT/eDFT.f90

@@ -1,4 +1,4 @@
-subroutine eDFT(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, & 
+subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, & 
                 nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
                 max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI,dipole_int)
 
@@ -14,6 +14,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,n
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   integer,intent(in)            :: guess_type
+  logical,intent(in)            :: mix
   double precision,intent(in)   :: thresh
 
   integer,intent(in)            :: nNuc
@@ -201,13 +202,13 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,n
   end if
 
 !------------------------------------------------------------------------
-! Compute N-centered UKS energy (UNBROKEN)
+! Compute N-centered UKS energy 
 !------------------------------------------------------------------------
 
   if(method == 'eDFT-UKS') then
 
     call cpu_time(start_KS)
-    call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & 
+    call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type,mix, & 
                   nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice,doNcentered)
     call cpu_time(end_KS)
 

src/eDFT/eDFT_UKS.f90

@@ -1,4 +1,4 @@
-subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh,max_diis,guess_type, &
+subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
                    nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice,doNcentered)
 
 ! Perform unrestricted Kohn-Sham calculation for ensembles
@@ -17,6 +17,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   integer,intent(in)            :: maxSCF,max_diis,guess_type
+  logical,intent(in)            :: mix
   double precision,intent(in)   :: thresh
   integer,intent(in)            :: nBas
   double precision,intent(in)   :: AO(nBas,nGrid)
@@ -39,6 +40,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   logical                       :: LDA_centered = .false.
   integer                       :: nSCF,nBasSq
   integer                       :: n_diis
+  integer                       :: nO(nspin)
   double precision              :: conv
   double precision              :: rcond(nspin)
   double precision              :: ET(nspin)
@@ -75,7 +77,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   double precision              :: E(nEns)
   double precision              :: Om(nEns)
 
-  integer                       :: ispin,iEns
+  integer                       :: ispin,iEns,iBas
 
 ! Hello world
 
@@ -136,6 +138,15 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
       c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
     end do
 
+    ! Mix guess to enforce symmetry breaking
+
+    nO(:) = 0
+    do ispin=1,nspin
+      nO(ispin) = int(sum(occnum(:,ispin,1)))
+    end do
+  
+    if(mix) call mix_guess(nBas,nO,c)
+
   else if(guess_type == 2) then
 
     do ispin=1,nspin
@@ -260,7 +271,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
       err(:,:,ispin) = matmul(F(:,:,ispin),matmul(Pw(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),Pw(:,:,ispin)),F(:,:,ispin))
     end do
 
-    conv = maxval(abs(err(:,:,:)))
+    if(nSCF > 1) conv = maxval(abs(err(:,:,:)))
     
 !   DIIS extrapolation