Fixing CCSDSolver.get_rdm() with frozen orbitals and energy from RDMs

qsdk/algorithms/classical/ccsd_solver.py

@@ -15,7 +15,8 @@
 """Class performing electronic structure calculation employing the CCSD method.
 """
 
-from pyscf import cc
+from pyscf import cc, lib
+from pyscf.cc.ccsd_rdm import _make_rdm1, _make_rdm2, _gamma1_intermediates, _gamma2_outcore
 
 from qsdk.algorithms.electronic_structure_solver import ElectronicStructureSolver
 
@@ -73,7 +74,16 @@ def get_rdm(self):
 
         # Solve the lambda equation and obtain the reduced density matrix from CC calculation
         self.cc_fragment.solve_lambda()
-        one_rdm = self.cc_fragment.make_rdm1()
-        two_rdm = self.cc_fragment.make_rdm2()
+        t1 = self.cc_fragment.t1
+        t2 = self.cc_fragment.t2
+        l1 = self.cc_fragment.l1
+        l2 = self.cc_fragment.l2
+
+        d1 = _gamma1_intermediates(self.cc_fragment, t1, t2, l1, l2)
+        f = lib.H5TmpFile()
+        d2 = _gamma2_outcore(self.cc_fragment, t1, t2, l1, l2, f, False)
+
+        one_rdm = _make_rdm1(self.cc_fragment, d1, with_frozen=False)
+        two_rdm = _make_rdm2(self.cc_fragment, d1, d2, with_dm1=True, with_frozen=False)
 
         return one_rdm, two_rdm

qsdk/toolboxes/molecular_computation/molecule.py

@@ -19,9 +19,10 @@
 import copy
 from dataclasses import dataclass, field
 import numpy as np
-from pyscf import gto, scf
+from pyscf import gto, scf, ao2mo
 import openfermion
 import openfermionpyscf
+from openfermion.ops.representations.interaction_operator import get_active_space_integrals
 
 from qsdk.toolboxes.molecular_computation.frozen_orbitals import get_frozen_core
 from qsdk.toolboxes.qubit_mappings.mapping_transform import get_fermion_operator
@@ -364,3 +365,48 @@ def freeze_mos(self, frozen_orbitals, inplace=True):
             copy_self.frozen_virtual = list_of_active_frozen[3]
 
             return copy_self
+
+    def energy_from_rdms(self, one_rdm, two_rdm):
+        """Computes the molecular energy from one- and two-particle reduced
+        density matrices (RDMs). Coefficients (integrals) are computed
+        on-the-fly using a pyscf object and the mean-field. Frozen orbitals
+        are supported with this method.
+
+        Args:
+            one_rdm (numpy.array): One-particle density matrix in MO basis.
+            two_rdm (numpy.array): Two-particle density matrix in MO basis.
+
+        Returns:
+            float: Molecular energy.
+        """
+
+        # Pyscf molecule to get integrals.
+        pyscf_mol = self.to_pyscf(self.basis)
+
+        # Corresponding to nuclear repulsion energy and static coulomb energy.
+        core_constant = float(pyscf_mol.energy_nuc())
+
+        # get_hcore is equivalent to int1e_kin + int1e_nuc.
+        one_electron_integrals = self.mean_field.mo_coeff.T @ self.mean_field.get_hcore() @ self.mean_field.mo_coeff
+
+        # Getting 2-body integrals in atomic and converting to molecular basis.
+        two_electron_integrals = ao2mo.kernel(pyscf_mol.intor("int2e"), self.mean_field.mo_coeff)
+        two_electron_integrals  = ao2mo.restore(1, two_electron_integrals, len(self.mean_field.mo_coeff))
+
+        # PQRS convention in openfermion:
+        # h[p,q]=\int \phi_p(x)* (T + V_{ext}) \phi_q(x) dx
+        # h[p,q,r,s]=\int \phi_p(x)* \phi_q(y)* V_{elec-elec} \phi_r(y) \phi_s(x) dxdy
+        # The convention is not the same with PySCF integrals. So, a change is
+        # made and reverse back after performing the truncation for frozen
+        # orbitals.
+        two_electron_integrals  = two_electron_integrals.transpose(0, 2, 3, 1)
+        core_offset, one_electron_integrals, two_electron_integrals = get_active_space_integrals(one_electron_integrals, two_electron_integrals, self.frozen_occupied, self.active_mos)
+        two_electron_integrals = two_electron_integrals.transpose(0, 3, 1, 2)
+
+        # Adding frozen electron contribution to core constant.
+        core_constant += core_offset
+
+        # Computing the total energy from integrals and provided RDMs.
+        e = core_constant + np.sum(one_electron_integrals * one_rdm) + 0.5*np.sum(two_electron_integrals * two_rdm)
+
+        return e.real

qsdk/toolboxes/molecular_computation/tests/test_molecule.py

@@ -15,6 +15,7 @@
 import unittest
 
 from qsdk import SecondQuantizedMolecule
+from qsdk.molecule_library import mol_H2_sto3g
 from qsdk.toolboxes.molecular_computation.molecule import atom_string_to_list
 
 
@@ -123,6 +124,26 @@ def test_get_fermionic_hamiltonian(self):
         shift = molecule.fermionic_hamiltonian.constant
         self.assertAlmostEqual(shift, -51.47120372466002, delta=1e-6)
 
+    def test_energy_from_rdms(self):
+        """Verify energy computation from RDMs."""
+
+        rdm1 = [[ 1.97453997e+00, -7.05987336e-17],
+                [-7.05987336e-17,  2.54600303e-02]]
+
+        rdm2 = [
+            [[[ 1.97453997e+00, -7.96423130e-17],
+              [-7.96423130e-17,  3.21234218e-33]],
+             [[-7.96423130e-17, -2.24213843e-01],
+              [ 0.00000000e+00,  9.04357944e-18]]],
+            [[[-7.96423130e-17,  0.00000000e+00],
+              [-2.24213843e-01,  9.04357944e-18]],
+             [[ 3.21234218e-33,  9.04357944e-18],
+              [ 9.04357944e-18,  2.54600303e-02]]]
+            ]
+
+        e_rdms = mol_H2_sto3g.energy_from_rdms(rdm1, rdm2)
+        self.assertAlmostEqual(e_rdms, -1.1372701, delta=1e-5)
+
 
 if __name__ == "__main__":
     unittest.main()