diff --git a/quantum_systems/general_orbital_system.py b/quantum_systems/general_orbital_system.py
index 971c65e..6faba7c 100644
--- a/quantum_systems/general_orbital_system.py
+++ b/quantum_systems/general_orbital_system.py
@@ -91,8 +91,11 @@ def compute_reference_energy(self):
 
         o, v = self.o, self.v
 
-        return self.np.trace(self.h[o, o]) + 0.5 * self.np.trace(
-            self.np.trace(self.u[o, o, o, o], axis1=1, axis2=3)
+        return (
+            self.np.trace(self.h[o, o])
+            + 0.5
+            * self.np.trace(self.np.trace(self.u[o, o, o, o], axis1=1, axis2=3))
+            + self.nuclear_repulsion_energy
         )
 
     def construct_fock_matrix(self, h, u, f=None):
diff --git a/quantum_systems/spatial_orbital_system.py b/quantum_systems/spatial_orbital_system.py
index dab7624..5e8bbb9 100644
--- a/quantum_systems/spatial_orbital_system.py
+++ b/quantum_systems/spatial_orbital_system.py
@@ -126,6 +126,7 @@ def compute_reference_energy(self):
             + 2
             * self.np.trace(self.np.trace(self.u[o, o, o, o], axis1=1, axis2=3))
             - self.np.trace(self.np.trace(self.u[o, o, o, o], axis1=1, axis2=2))
+            + self.nuclear_repulsion_energy
         )
 
     def construct_fock_matrix(self, h, u, f=None):
diff --git a/tests/test_custom_system.py b/tests/test_custom_system.py
index 40217b6..b77a865 100644
--- a/tests/test_custom_system.py
+++ b/tests/test_custom_system.py
@@ -89,9 +89,15 @@ def test_change_of_basis():
 
 
 def test_reference_energy():
-    system = construct_pyscf_system_rhf("he", basis="cc-pvdz")
+    gos_system = construct_pyscf_system_rhf(
+        "li 0.0 0.0 0.0; h 0.0 0.0 3.08", basis="cc-pvdz"
+    )
+    spas_system = construct_pyscf_system_rhf(
+        "li 0.0 0.0 0.0; h 0.0 0.0 3.08", basis="cc-pvdz", add_spin=False
+    )
 
     # This energy is found from PySCF's RHF solver
-    he_energy = -2.85516047724274
+    lih_energy = -7.98367215457454
 
-    assert abs(system.compute_reference_energy() - he_energy) < 1e-8
+    assert abs(gos_system.compute_reference_energy() - lih_energy) < 1e-8
+    assert abs(spas_system.compute_reference_energy() - lih_energy) < 1e-8