Nuclear repulsion energy (#54)

* add nuclear repulsion energy to reference energy

* test reference energy for both gos and spas system using lih as test system. a molecule is used since it then includes a non-zero repulsion energy.

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):

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):

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