diff --git a/src/eDFT/unrestricted_hartree_energy.f90 b/src/eDFT/unrestricted_hartree_energy.f90
new file mode 100644
index 0000000..9937809
--- /dev/null
+++ b/src/eDFT/unrestricted_hartree_energy.f90
@@ -0,0 +1,29 @@
+subroutine unrestricted_hartree_energy(nBas,P,J,EH)
+
+! Compute the unrestricted version of the Hartree energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: P(nBas,nBas,nspin)
+  double precision,intent(in)   :: J(nBas,nBas,nspin)
+
+! Local variables
+
+  double precision,external     :: trace_matrix
+
+! Output variables
+
+  double precision,intent(out)  :: EH(nsp)
+
+! Compute the components of the Hartree energy
+
+  EH(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
+  EH(2) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) &
+        + 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,1)))
+  EH(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
+
+end subroutine unrestricted_hartree_energy
diff --git a/src/eDFT/unrestricted_hartree_potential.f90 b/src/eDFT/unrestricted_hartree_potential.f90
new file mode 100644
index 0000000..5860e02
--- /dev/null
+++ b/src/eDFT/unrestricted_hartree_potential.f90
@@ -0,0 +1,33 @@
+subroutine unrestricted_hartree_potential(nBas,P,ERI,J)
+
+! Compute the unrestricted version of the Hartree potential 
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: P(nBas,nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer                       :: mu,nu,la,si
+
+! Output variables
+
+  double precision,intent(out)  :: J(nBas,nBas)
+
+  J(:,:) = 0d0
+  do si=1,nBas
+    do la=1,nBas
+      do nu=1,nBas
+        do mu=1,nBas
+          J(mu,nu) = J(mu,nu) + P(la,si)*ERI(mu,la,nu,si)
+        enddo
+      enddo
+    enddo
+  enddo
+
+
+end subroutine unrestricted_hartree_potential