From a4ed7372fbfac98e09296654b7375ed9babaede0 Mon Sep 17 00:00:00 2001 From: haakoek Date: Wed, 14 Apr 2021 14:14:05 +0200 Subject: [PATCH] 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 | 7 +++++-- quantum_systems/spatial_orbital_system.py | 1 + tests/test_custom_system.py | 12 +++++++++--- 3 files changed, 15 insertions(+), 5 deletions(-) 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