workin on restricted

input/dft

@@ -1,7 +1,7 @@
 # Restricted or unrestricted KS calculation
   GOK-RKS
 # exchange rung: Hartree = 0, LDA = 1 (S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666
-  1 RS51
+  666 HF
 # correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666
   0 HF
 # quadrature grid SG-n

src/eDFT/GOK_RKS.f90

@@ -1,5 +1,5 @@
-subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh,max_diis,guess_type, &
+subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh, & 
-               nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
+                   max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
 
 ! Perform restricted Kohn-Sham calculation for ensembles
 
@@ -206,15 +206,15 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
       rhow(:) = rhow(:) + wEns(iEns)*rho(:,iEns) 
     end do
 
-    if(xc_rung > 1 .and. xc_rung /= 666) then 
+    if(xc_rung > 1) then 
 
-!     Ground state density 
+!     Compute gradient of the one-electron density
 
       do iEns=1,nEns
         call gradient_density(nGrid,nBas,P(:,:,iEns),AO(:,:),dAO(:,:,:),drho(:,:,iEns))
       end do
 
-!     Weight-dependent one-electron density 
+!     Weight-dependent one-electron density gradient
 
       drhow(:,:) = 0d0
       do iEns=1,nEns
@@ -234,11 +234,12 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
 
 !   Compute correlation potential
 
-    call correlation_potential(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fc(:,:))
+    call restricted_correlation_potential(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:), & 
+                                          nBas,AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fc(:,:))
 
 !   Build Fock operator
 
-    F(:,:) = Hc(:,:) + J(:,:) + J(:,:) + Fx(:,:) + Fc(:,:)
+    F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*Fx(:,:) + Fc(:,:)
 
 !   Check convergence 
 
@@ -275,10 +276,11 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
 
     call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas, &
                          Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex)
+    Ex = 0.5d0*Ex
 
 !   Correlation energy
 
-    call correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),Ec)
+    call restricted_correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),Ec)
 
 !   Total energy
 
@@ -321,7 +323,7 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
 ! Compute individual energies from ensemble energy
 !------------------------------------------------------------------------
 
-  call individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas,     &
+  call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas, &
                                     AO(:,:),dAO(:,:,:),nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc,    & 
                                     Pw(:,:),rhow(:),drhow(:,:),J(:,:),Fx(:,:),FxHF(:,:),         & 
                                     Fc(:,:),P(:,:,:),rho(:,:),drho(:,:,:),E(:),Om(:))

src/eDFT/eDFT.f90

@@ -138,7 +138,7 @@ program eDFT
 
     call cpu_time(start_KS)
     call GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & 
-                 nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,EKS)
+                 nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,EKS)
     call cpu_time(end_KS)
 
     t_KS = end_KS - start_KS

src/eDFT/gradient_density.f90

@@ -13,7 +13,7 @@ subroutine gradient_density(nGrid,nBas,P,AO,dAO,drho)
   integer,intent(in)            :: nBas
   double precision,intent(in)   :: P(nBas,nBas)
   double precision,intent(in)   :: AO(nBas,nGrid)
-  double precision,intent(in)   :: dAO(3,nBas,nGrid)
+  double precision,intent(in)   :: dAO(ncart,nBas,nGrid)
 
 ! Local variables
 
@@ -22,13 +22,13 @@ subroutine gradient_density(nGrid,nBas,P,AO,dAO,drho)
 
 ! Output variables
 
-  double precision,intent(out)  :: drho(3,nGrid)
+  double precision,intent(out)  :: drho(ncart,nGrid)
 
   drho(:,:) = 0d0
   do iG=1,nGrid
     do mu=1,nBas
       do nu=1,nBas
-        do ixyz=1,3
+        do ixyz=1,ncart
           drho(ixyz,iG) = drho(ixyz,iG) & 
                         + P(mu,nu)*(dAO(ixyz,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(ixyz,nu,iG))
         enddo
@@ -36,10 +36,10 @@ subroutine gradient_density(nGrid,nBas,P,AO,dAO,drho)
     enddo
   enddo
 
-  do iG=1,nGrid
+! do iG=1,nGrid
-    do ixyz=1,3
+!   do ixyz=1,ncart
-      if(abs(drho(ixyz,iG)) < thresh) drho(ixyz,iG) = thresh
+!     if(abs(drho(ixyz,iG)) < thresh) drho(ixyz,iG) = thresh
-    enddo
+!   enddo
-  enddo
+! enddo
 
 end subroutine gradient_density

src/eDFT/print_RKS.f90

@@ -43,7 +43,6 @@ subroutine print_RKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)
   write(*,'(A40,1X,F16.10,A3)') ' Two-electron    energy: ',EJ + Ex + Ec,' au'
   write(*,'(A40,1X,F16.10,A3)') ' Coulomb         energy: ',EJ,' au'
   write(*,'(A40,1X,F16.10,A3)') ' Exchange        energy: ',Ex,' au'
-  write(*,'(A40,1X,F16.10,A3)') ' Exchange        energy: ',Ex,' au'
   write(*,'(A40,1X,F16.10,A3)') ' Correlation     energy: ',Ec,' au'
   write(*,'(A60)')              '-------------------------------------------------'
   write(*,'(A40,1X,F16.10,A3)') ' Electronic      energy: ',Ew,' au'

src/eDFT/restricted_density_matrix.f90

@@ -29,7 +29,12 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
 ! Doubly-excited state density matrix
 
   iEns = 2 
-  P(:,:,iEns) = 2d0*matmul(c(:,1:nO-1),transpose(c(:,1:nO-1))) &
+  if(nO > 1) then 
-              + 2d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
+
+    P(:,:,iEns) = 2d0*matmul(c(:,1:nO-1),transpose(c(:,1:nO-1))) 
+
+  end if
+
+  P(:,:,iEns) = P(:,:,iEns) + 2d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
 
 end subroutine restricted_density_matrix