Support for UMP2 initial parameters (sandbox-quantum#310)

* Support for UMP2 initial parameters.

tangelo/algorithms/classical/mp2_solver.py

@@ -15,8 +15,12 @@
 """Class performing electronic structure calculation employing the Moller-Plesset perturbation theory (MP2) method.
 """
 
+from itertools import combinations, product
+from math import ceil
+
 from tangelo.algorithms.electronic_structure_solver import ElectronicStructureSolver
 from tangelo.helpers.utils import is_package_installed
+from tangelo.toolboxes.ansatz_generator._unitary_cc_openshell import uccsd_openshell_get_packed_amplitudes
 
 
 class MP2Solver(ElectronicStructureSolver):
@@ -37,9 +41,6 @@ def __init__(self, molecule):
             raise ModuleNotFoundError(f"Using {self.__class__.__name__} requires the installation of the pyscf package")
         from pyscf import mp
 
-        if molecule.uhf:
-            raise NotImplementedError(f"SecondQuantizedMolecule that use UHF are not currently supported in {self.__class__.__name__}. Use CCSDSolver")
-
         self.mp = mp
         self.mp2_fragment = None
 
@@ -49,14 +50,28 @@ def __init__(self, molecule):
         self.frozen = molecule.frozen_mos
         self.uhf = molecule.uhf
 
+        # Define variables used to transform the MP2 parameters into an ordered
+        # list of parameters with single and double excitations.
+        if self.spin != 0 or self.uhf:
+            self.n_alpha, self.n_beta = molecule.n_active_ab_electrons
+            self.n_active_moa, self.n_active_mob = molecule.n_active_mos if self.uhf else (molecule.n_active_mos,)*2
+        else:
+            self.n_occupied = ceil(molecule.n_active_electrons / 2)
+            self.n_virtual = molecule.n_active_mos - self.n_occupied
+
     def simulate(self):
         """Perform the simulation (energy calculation) for the molecule.
 
         Returns:
             float: MP2 energy.
         """
+
         # Execute MP2 calculation
-        self.mp2_fragment = self.mp.MP2(self.mean_field, frozen=self.frozen)
+        if self.uhf:
+            self.mp2_fragment = self.mp.UMP2(self.mean_field, frozen=self.frozen)
+        else:
+            self.mp2_fragment = self.mp.RMP2(self.mean_field, frozen=self.frozen)
+
         self.mp2_fragment.verbose = 0
         _, self.mp2_t2 = self.mp2_fragment.kernel()
 
@@ -75,13 +90,54 @@ def get_rdm(self):
             RuntimeError: If no simulation has been run.
         """
 
-        # Check if MP2 is performed
+        # Check if MP2 has been performed
         if self.mp2_fragment is None:
-            raise RuntimeError("MP2Solver: Cannot retrieve RDM. Please run the 'simulate' method first")
+            raise RuntimeError(f"{self.__class__.__name__}: Cannot retrieve RDM. Please run the 'simulate' method first")
         if self.frozen is not None:
-            raise RuntimeError("MP2Solver: RDM calculation is not implemented with frozen orbitals.")
+            raise RuntimeError(f"{self.__class__.__name__}: RDM calculation is not implemented with frozen orbitals.")
 
         one_rdm = self.mp2_fragment.make_rdm1()
         two_rdm = self.mp2_fragment.make_rdm2()
 
         return one_rdm, two_rdm
+
+    def get_mp2_amplitudes(self):
+        """Compute the double amplitudes from the MP2 perturbative method, and
+        then reorder the elements into a dense list. The single (T1) amplitudes
+        are set to a small non-zero value. The ordering is single, double
+        (diagonal), double (non-diagonal).
+
+        Returns:
+            list of float: The electronic excitation amplitudes.
+        """
+
+        # Check if MP2 has been performed.
+        if self.mp2_fragment is None:
+            raise RuntimeError(f"{self.__class__.__name__}: Cannot retrieve MP2 parameters. Please run the 'simulate' method first")
+
+        if self.spin != 0 or self.uhf:
+            # Reorder the T2 amplitudes in a dense list.
+            mp2_params = uccsd_openshell_get_packed_amplitudes(
+                self.mp2_t2[0],  # aa
+                self.mp2_t2[2],  # bb
+                self.mp2_t2[1],  # ab
+                self.n_alpha,
+                self.n_beta,
+                self.n_active_moa,
+                self.n_active_mob
+            )
+        else:
+            # Get singles amplitude. Just get "up" amplitude, since "down" should be the same
+            singles = [2.e-5] * (self.n_virtual * self.n_occupied)
+
+            # Get singles and doubles amplitudes associated with one spatial occupied-virtual pair
+            doubles_1 = [-self.mp2_t2[q, q, p, p]/2. if (abs(-self.mp2_t2[q, q, p, p]/2.) > 1e-15) else 0.
+                        for p, q in product(range(self.n_virtual), range(self.n_occupied))]
+
+            # Get doubles amplitudes associated with two spatial occupied-virtual pairs
+            doubles_2 = [-self.mp2_t2[q, s, p, r] for (p, q), (r, s)
+                        in combinations(product(range(self.n_virtual), range(self.n_occupied)), 2)]
+
+            mp2_params = singles + doubles_1 + doubles_2
+
+        return mp2_params

tangelo/algorithms/classical/tests/test_mp2_solver.py

@@ -14,8 +14,10 @@
 
 import unittest
 
+import numpy as np
+
 from tangelo.algorithms.classical import MP2Solver
-from tangelo.molecule_library import mol_H2_321g, mol_Be_321g, mol_H4_cation_sto3g
+from tangelo.molecule_library import mol_H2_321g, mol_Be_321g, mol_H2_sto3g, mol_H2_sto3g_uhf
 
 
 class MP2SolverTest(unittest.TestCase):
@@ -59,6 +61,29 @@ def test_be_frozen_core(self):
 
         self.assertAlmostEqual(energy, -14.5092873, places=6)
 
+    def test_get_mp2_params_restricted(self):
+        """Test the packing of RMP2 amplitudes as initial parameters for coupled
+        cluster based methods.
+        """
+
+        solver = MP2Solver(mol_H2_sto3g)
+        solver.simulate()
+
+        ref_params = [2.e-05, 3.632537e-02]
+
+        np.testing.assert_array_almost_equal(ref_params, solver.get_mp2_amplitudes())
+
+    def test_get_mp2_params_unrestricted(self):
+        """Test the packing of UMP2 amplitudes as initial parameters for coupled
+        cluster based methods.
+        """
+
+        solver = MP2Solver(mol_H2_sto3g_uhf)
+        solver.simulate()
+        ref_params = [0., 0., 0.030736]
+
+        np.testing.assert_array_almost_equal(ref_params, solver.get_mp2_amplitudes())
+
 
 if __name__ == "__main__":
     unittest.main()

tangelo/toolboxes/ansatz_generator/uccsd.py

@@ -28,13 +28,10 @@
         Physical Review A 95, 020501 (2017).
 """
 
-import itertools
-
 import numpy as np
 from openfermion.circuits import uccsd_singlet_generator
 
 from tangelo import SecondQuantizedMolecule
-from tangelo.algorithms.classical.mp2_solver import MP2Solver
 from tangelo.linq import Circuit
 from tangelo.helpers.utils import is_package_installed
 from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
@@ -99,11 +96,10 @@ def __init__(self, molecule, mapping="JW", up_then_down=False, spin=None, refere
         # TODO: support for others
         self.supported_reference_state = {"HF", "zero"}
         # Supported var param initialization
-        self.supported_initial_var_params = {"ones", "random", "mp2"} if (self.spin == 0 and not self.molecule.uhf) else {"ones", "random"}
+        self.supported_initial_var_params = {"ones", "random", "mp2"}
 
         # Default initial parameters for initialization
-        # TODO: support for openshell MP2 initialization
-        self.var_params_default = "mp2" if "mp2" in self.supported_initial_var_params else "ones"
+        self.var_params_default = "mp2"
         self.reference_state = reference_state
 
         self.var_params = None
@@ -255,22 +251,22 @@ def _get_openshell_qubit_operator(self):
         return qubit_op
 
     def _compute_mp2_params(self):
-        """Computes the MP2 initial variational parameters. Compute the initial
-        variational parameters with PySCF MP2 calculation, and then reorders the
-        elements into the appropriate convention. MP2 only has doubles (T2)
-        amplitudes, thus the single (T1) amplitudes are set to a small non-zero
-        value and added. The ordering is single, double (diagonal), double
-        (non-diagonal).
+        """Compute the MP2 initial variational parameters. Only double
+        excitation amplitudes are retrieved from an MP2 calculation (singles
+        set to a small number close to 0).
 
         Returns:
             list of float: The initial variational parameters.
         """
 
-        if self.molecule.uhf:
-            raise NotImplementedError(f"MP2 initialization is not currently implemented for UHF reference in {self.__class__}")
         if not is_package_installed("pyscf"):
             raise ValueError(f"pyscf is required for MP2 initial parameters in {self.__class__.__name__}.")
 
+        # Import here to solve an AttributeError: partially initialized module
+        # tangelo.toolboxes.ansatz_generator' has no attribute 'UCCSD'
+        # (most likely due to a circular import).
+        from tangelo.algorithms.classical.mp2_solver import MP2Solver
+
         if not isinstance(self.molecule.solver, IntegralSolverPySCF):
             pymol = SecondQuantizedMolecule(self.molecule.xyz, self.molecule.q, self.molecule.spin, basis=self.molecule.basis,
                                             ecp=self.molecule.ecp, symmetry=self.molecule.symmetry, uhf=self.molecule.uhf,
@@ -281,15 +277,4 @@ def _compute_mp2_params(self):
         mp2 = MP2Solver(pymol)
         mp2.simulate()
 
-        # Get singles amplitude. Just get "up" amplitude, since "down" should be the same
-        singles = [2.e-5] * (self.n_virtual * self.n_occupied)
-
-        # Get singles and doubles amplitudes associated with one spatial occupied-virtual pair
-        doubles_1 = [-mp2.mp2_t2[q, q, p, p]/2. if (abs(-mp2.mp2_t2[q, q, p, p]/2.) > 1e-15) else 0.
-                     for p, q in itertools.product(range(self.n_virtual), range(self.n_occupied))]
-
-        # Get doubles amplitudes associated with two spatial occupied-virtual pairs
-        doubles_2 = [-mp2.mp2_t2[q, s, p, r] for (p, q), (r, s)
-                     in itertools.combinations(itertools.product(range(self.n_virtual), range(self.n_occupied)), 2)]
-
-        return singles + doubles_1 + doubles_2
+        return mp2.get_mp2_amplitudes()