debug from clotilde for N-centered

src/eDFT/US51_lda_exchange_individual_energy.f90

@@ -54,8 +54,13 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered
 
   enddo
 
-  Exrr  = kappa*Exrr
+  if(doNcentered) then
-  ExrI  = kappa*ExrI
+  
+    Exrr  = kappa*Exrr
+    ExrI  = kappa*ExrI
+
+  endif
+
 
   Ex = Exrr + ExrI + ExrrI
 

src/eDFT/read_options_dft.f90

@@ -137,14 +137,15 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
   if (doNcentered) then
 
     wEns(1) = 1d0 - wEns(2) - wEns(3)
+    wEns(2) = (nEl(1)/nEl(2))*wEns(2)
+    wEns(3) = (nEl(1)/nEl(3))*wEns(3)
+
 
   else
 
 !   wEns(1) = 1d0 -  nEl(2)/nEl(1)*wEns(2) - nEl(3)/nEl(1)*wEns(3) 
 
     wEns(1) = 1d0 - wEns(2) - wEns(3)
-    wEns(2) = nEl(1)/nEl(2)*wEns(2)
-    wEns(3) = nEl(1)/nEl(3)*wEns(3)
 
   end if
 

src/eDFT/unrestricted_individual_energy.f90

@@ -97,11 +97,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 
   do ispin=1,nspin
     do iEns=1,nEns
-      if (doNcentered) then
-        ET(ispin,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
-      else 
         ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) 
-      end if
     end do
   end do
 
@@ -111,11 +107,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
 
   do iEns=1,nEns
     do ispin=1,nspin
-      if (doNcentered) then
-        EV(ispin,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
-      else 
         EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
-      end if  
   end do
   end do
 
@@ -129,18 +121,32 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
       call hartree_coulomb(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin))
     end do
 
-    EJ(1,iEns) =       trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) &
+      if(doNcentered) then 
-               - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1)))
  
-    EJ(2,iEns) =       trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) &
+        EJ(1,iEns) =   kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) &
-               +       trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) &
+                   -0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1)))
-               - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2)))     &
-               - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1)))
    
-    EJ(3,iEns) =       trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) &
+        EJ(2,iEns) =     kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) &
-               - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
+                   +     kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) &
+                   -0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) &
+                   - 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1)))
    
-    if(doNcentered) EJ(:,iEns) = kappa(iEns)*EJ(:,iEns)
+        EJ(3,iEns) =    kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) &
+                   -0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
+      else
+
+
+        EJ(1,iEns) =       trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) &
+                   - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1)))
+
+        EJ(2,iEns) =       trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) &
+                   +       trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) &
+                   - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2)))     &
+                   - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1)))
+
+        EJ(3,iEns) =       trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) &
+                   - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
+      end if
 
   end do