fix eDFT for fractional spin

input/dft

@@ -6,7 +6,7 @@
 # GGA          = 2: RB88
 # Hybrid       = 4
 # Hartree-Fock = 666
-  1 S51  
+  666 HF
 # correlation rung: 
 # Hartree      = 0
 # LDA          = 1: RVWN5,RMFL20
@@ -19,20 +19,20 @@
 # Number of states in ensemble (nEns)
   3
 # occupation numbers of orbitals nO and nO+1
- 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
- 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
- 
  1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
  
- 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 
+ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+ 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+ 
+ 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
  1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  0.00 1.00 
+  0.25 0.0
 # N-centered?
- T 
+ F 
 # Parameters for CC weight-dependent exchange functional
-0.445525 0.0901503 -0.286898
+0.0 0.0 0.0
-0.191734 -0.0364788 -0.017035
+0.0 0.0 0.0
 # choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2
 2

input/methods

@@ -1,11 +1,11 @@
 # RHF UHF KS  MOM
-  T   F   F   F
+  F   F   T   F
 # MP2* MP3 MP2-F12
   F   F   F
 # CCD DCD CCSD CCSD(T) 
   F   F   F    F
 # drCCD rCCD lCCD pCCD
-  F      F   F    T
+  F      F   F    F
 # CIS* CIS(D) CID CISD
   F    F      F   F    
 # RPA* RPAx* ppRPA 

input/options

@@ -1,11 +1,11 @@
 # HF: maxSCF thresh   DIIS n_diis guess_type ortho_type mix_guess
-       128    0.0000001  T    5      1          1         T
+       128    0.000001  T    5      1          1         F
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
       64     0.00001  T    5
 # spin: TDA singlet triplet spin_conserved spin_flip 
-         F   T       T       T              F 
+         T   T       T       T              T 
 # 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 
@@ -13,6 +13,6 @@
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn
-        T    F    T    F
+        T    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.740848
+H 0.0 0.0 0.741

src/eDFT/eDFT.f90

@@ -4,7 +4,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,n
 
 ! exchange-correlation density-functional theory calculations
 
-  use xc_f90_lib_m
+! use xc_f90_lib_m
 
   implicit none
   include 'parameters.h'
@@ -92,8 +92,8 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,n
 
 ! Libxc version
 
-  call xc_f90_version(vmajor, vminor, vmicro)
+! call xc_f90_version(vmajor, vminor, vmicro)
-  write(*,'("Libxc version: ",I1,".",I1,".",I1)') vmajor, vminor, vmicro
+! write(*,'("Libxc version: ",I1,".",I1,".",I1)') vmajor, vminor, vmicro
 
 ! call xcinfo()
 

src/eDFT/eDFT_UKS.f90

@@ -337,7 +337,6 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
       nEl(ispin) = electron_number(nGrid,weight,rhow(:,ispin))
     end do
 
-
 !   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,F10.6,1X,A1,1X)') & 
@@ -364,12 +363,11 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
   end if
 
 !!!!!
-    do iEns=1,nEns
+!   do iEns=1,nEns
-      print*,'occnum=',occnum(1,1,iEns),occnum(2,1,iEns),occnum(1,2,iEns),occnum(2,2,iEns)
+!     print*,'occnum=',occnum(1,1,iEns),occnum(2,1,iEns),occnum(1,2,iEns),occnum(2,2,iEns)
-      print*,'nel up and down and total=', electron_number(nGrid,weight,&
+!     print*,'nel up and down and total=', electron_number(nGrid,weight,&
-                                           rho(:,1,iEns)),electron_number(nGrid,weight,rho(:,2,iEns)),sum(nEl(:))
+!                                          rho(:,1,iEns)),electron_number(nGrid,weight,rho(:,2,iEns)),sum(nEl(:))
-
+!   end do
-    end do
 !!!!!
 
 ! Compute final KS energy

src/eDFT/read_options_dft.f90

@@ -114,6 +114,7 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
     write(*,*) (int(occnum(iBas,2,iEns)),iBas=1,nBas)
     write(*,*)
   end do
+
 ! Read ensemble weights for real physical (fractional number of electrons) ensemble (w1,w2)
 
   allocate(nEl(nEns))
@@ -123,8 +124,6 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
       nEl(iEns) = nEl(iEns) + occnum(iBas,1,iEns) + occnum(iBas,2,iEns)
     end do
   end do
-  print*,'nEl'
-  print*,nEl
 
   doNcentered = .false.
 
@@ -164,13 +163,13 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
   read(1,*) Cx_choice
 
   write(*,*)'----------------------------------------------------------'
-  write(*,*)' parameters for w1-dependant exchange functional coefficient '
+  write(*,*)' parameters for w1-dependent exchange functional coefficient '
   write(*,*)'----------------------------------------------------------'
   call matout(3,1,aCC_w1)
   write(*,*)
 
   write(*,*)'----------------------------------------------------------'
-  write(*,*)' parameters for w2-dependant exchange functional coefficient '
+  write(*,*)' parameters for w2-dependent exchange functional coefficient '
   write(*,*)'----------------------------------------------------------'
   call matout(3,1,aCC_w2)
   write(*,*) 

src/eDFT/unrestricted_fock_exchange_energy.f90

@@ -20,6 +20,6 @@ subroutine unrestricted_fock_exchange_energy(nBas,P,Fx,Ex)
 
 ! Compute HF exchange energy
 
-  Ex = trace_matrix(nBas,matmul(P,Fx))
+  Ex = 0.5d0*trace_matrix(nBas,matmul(P,Fx))
 
 end subroutine unrestricted_fock_exchange_energy

src/eDFT/unrestricted_individual_energy.f90

@@ -89,8 +89,6 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
     kappa(iEns) = nEl(iEns)/nEl(1)
   end do
 
-  print*,'test1'
-
 !------------------------------------------------------------------------
 ! Kinetic energy
 !------------------------------------------------------------------------
@@ -104,7 +102,6 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
       end if
     end do
   end do
-  print*,'test2'
 
 !------------------------------------------------------------------------
 ! Potential energy
@@ -120,7 +117,6 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
   end do
   end do
 
-  print*,'test3'
 !------------------------------------------------------------------------
 ! Individual Hartree energy
 !------------------------------------------------------------------------
@@ -145,7 +141,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
     if(doNcentered) EJ(:,iEns) = kappa(iEns)*EJ(:,iEns)
 
   end do
-       print*,'test4'
+
 !------------------------------------------------------------------------
 ! Checking Hartree contributions for each individual states
 !------------------------------------------------------------------------
@@ -174,8 +170,6 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
     end do
   end do
 
-  print*,'test5'
-
 !------------------------------------------------------------------------
 ! Checking exchange contributions for each individual states
 !------------------------------------------------------------------------