RMFL20 done

src/eDFT/GOK_RKS.f90

@@ -236,7 +236,7 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
 
     call restricted_correlation_potential(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:), & 
                                           nBas,AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fc(:,:))
-
+!   print*,'Done with restricted_correlation_potential'
 !   Build Fock operator
 
     F(:,:) = Hc(:,:) + J(:,:) + Fx(:,:) + Fc(:,:)
@@ -280,6 +280,7 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
 !   Correlation energy
 
     call restricted_correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),Ec)
+!   print*,'Done with restricted_correlation_energy'
 
 !   Total energy
 

src/eDFT/RMFL20_lda_correlation_energy.f90

@@ -49,8 +49,6 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec)
 
 ! LDA-centered functional
 
-  EcLDA = 0d0
-
   call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),EcLDA)
 
   do iEns=1,nEns

src/eDFT/RMFL20_lda_correlation_individual_energy.f90

@@ -16,7 +16,6 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,
 
 ! Local variables
 
-  logical                       :: LDA_centered = .true.
   integer                       :: iEns
   double precision              :: EcLDA
   double precision,allocatable  :: aMFL(:,:)
@@ -40,27 +39,21 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,
   aMFL(2,2) = -0.0506019d0
   aMFL(3,2) = +0.0331417d0
 
-! Compute correlation energy for ground, singly-excited and doubly-excited states
+! Compute correlation energy for ground- and doubly-excited states
 
   do iEns=1,nEns
-
     call restricted_elda_correlation_individual_energy(nEns,aMFL(:,iEns),nGrid,weight(:),rhow(:),rho(:),EceLDA(iEns))
-
   end do
 
 ! LDA-centered functional
 
   EcLDA = 0d0
 
-  if(LDA_centered) then
+  call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA)
 
-    call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA)
+  do iEns=1,nEns
-
+    EceLDA(iEns) = EceLDA(iEns) + EcLDA - EceLDA(1)
-    do iEns=1,nEns
+  end do
-      EceLDA(iEns) = EceLDA(iEns) + EcLDA - EceLDA(1)
-    end do
-
-  end if
 
 ! Weight-denpendent functional for ensembles
 

src/eDFT/RVWN5_lda_correlation_energy.f90

@@ -22,7 +22,7 @@ subroutine RVWN5_lda_correlation_energy(nGrid,weight,rho,Ec)
 
 ! Output variables
 
-  double precision              :: Ec(nsp)
+  double precision              :: Ec
 
 ! Parameters of the functional
 

src/eDFT/print_restricted_individual_energy.f90

@@ -96,12 +96,13 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,
   write(*,'(A50)')           ' Derivative discontinuities (DD)    '
   write(*,'(A60)')           '-------------------------------------------------'
   do iEns=1,nEns
+    write(*,*)
     write(*,'(A40,I2,A2,F16.10,A3)') '  x ensemble derivative state ',iEns,': ',ExDD(iEns), ' au'
     write(*,'(A40,I2,A2,F16.10,A3)') '  c ensemble derivative state ',iEns,': ',EcDD(iEns), ' au'
     write(*,'(A40,I2,A2,F16.10,A3)') ' xc ensemble derivative state ',iEns,': ',ExcDD(iEns),' au'
-    write(*,*)
   end do
   write(*,'(A60)')           '-------------------------------------------------'
+  write(*,*)
 
 !------------------------------------------------------------------------
 ! Total and Excitation energies

src/eDFT/restricted_elda_correlation_energy.f90

@@ -16,7 +16,7 @@ subroutine restricted_elda_correlation_energy(nEns,aMFL,nGrid,weight,rho,Ec)
 ! Local variables
 
   integer                       :: iG
-  double precision              :: r,ec_p
+  double precision              :: r,e
 
 ! Output variables
 
@@ -33,9 +33,9 @@ subroutine restricted_elda_correlation_energy(nEns,aMFL,nGrid,weight,rho,Ec)
 
     if(r > threshold) then
 
-      ec_p = aMFL(1)/(1d0 + aMFL(2)*r**(-1d0/6d0) + aMFL(3)*r**(-1d0/3d0))
+      e = aMFL(1)/(1d0 + aMFL(2)*r**(-1d0/6d0) + aMFL(3)*r**(-1d0/3d0))
 
-      Ec = Ec + weight(iG)*ec_p*r
+      Ec = Ec + weight(iG)*e*r
 
     end if
 

src/eDFT/restricted_elda_correlation_individual_energy.f90

@@ -43,7 +43,7 @@ subroutine restricted_elda_correlation_individual_energy(nEns,aLF,nGrid,weight,r
       dFcdr = ec_p*dFcdr/(6d0*r)
       dFcdr = ec_p + dFcdr*r
 
-      Ec = Ec + weight(iG)*rI*dFcdr
+      Ec = Ec + weight(iG)*(ec_p*rI + dFcdr*r*rI -  dFcdr*r*r)
 
      end if
 

src/eDFT/restricted_individual_energy.f90

@@ -118,9 +118,9 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
 !------------------------------------------------------------------------
 
   do iEns=1,nEns
-    Exc(iEns) = Ex(iEns) + Ec(iEns)
+    Exc(iEns) =        Ex(iEns) + Ec(iEns)
-    E(iEns) = ENuc + ET(iEns) + EV(iEns) + EJ(iEns) &
+    E(iEns)   = ENuc + ET(iEns) + EV(iEns) + EJ(iEns) &
-                   + Ex(iEns) + Ec(iEns) + ExcDD(iEns)
+                     + Ex(iEns) + Ec(iEns) + ExcDD(iEns)
   end do
 
 !------------------------------------------------------------------------

src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90

@@ -16,8 +16,6 @@ subroutine restricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGr
 
 ! Local variables
 
-  double precision              :: aC
-
 ! Output variables
 
   double precision,intent(out)  :: Ec(nEns)