Include nuclear repulsion energy in energy functions

[Schoyen]

If the specified quantum system has a non-zero nuclear repulsion energy (or, a constant term in the Hamiltonian which as of now is named "nuclear repulsion energy") this has to be added to the energy outside the energy calculations by the user if they want to get the full expression. It would be nice to include the constant term in all the energy expressions in the coupled cluster classes so the user does not have to constantly remember to add this term. Perhaps the simplest would be to include this in the reference energy calculations (e.g., here), viz.

class CoupledCluster(metaclass=abc.ABCMeta):

    # ...

    def compute_reference_energy(self):
        """Computes reference energy
        Returns
        -------
        np.array
            Reference energy
        """

        return compute_reference_energy(
            self.f, self.u, self.o, self.v, np=self.np
        ) + self.system.nuclear_repulsion_energy

    # ....

[haakoek]

Agree

[Schoyen]

This should be seen in relation with #34. The nuclear repulsion energy should be computed in the reference energy and will therefore be moved to quantum-systems as a constant addition to the appropriate reference energy function.

[haakoek]

I think the nuclear repulsion energy should be added in quantum systems compute_reference energy. However, keep in mind that for orbital adaptive methods such as OMP, OCC and OACC, the reference energy changes as the basis changes.

Thus, we have to keep in mind, espacially in time-dependent calculations (if we want to check conservation of energy in the field-free case or after a pulse) that IF we compute the reference energy it is with respect to the updated basis.

If I recall correctly, the way we have set up the time-dependence for the orbital adaptive methods we do not necessarily change the integrals in system.

[Schoyen]

I agree! Now the functions compute_reference_energy in SpatialOrbitalSystem and GeneralOrbitalSystem compute the reference energy for the basis set in the system. That is, there is no option to pass in a h and u. I think it would make sense to make this function accept an h and a u from the caller similar to the construct_fock_matrix-functions (in SpatialOrbitalSystem and in GeneralOrbitalSystem). This would make it easy to use for both solvers with and without static orbitals.

[Schoyen]

Once merged, this will be fixed in #64.