diff --git a/tangelo/algorithms/classical/mp2_solver.py b/tangelo/algorithms/classical/mp2_solver.py index 253854bda..0083203fc 100644 --- a/tangelo/algorithms/classical/mp2_solver.py +++ b/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 diff --git a/tangelo/algorithms/classical/tests/test_mp2_solver.py b/tangelo/algorithms/classical/tests/test_mp2_solver.py index e4dd551cd..4c2433755 100644 --- a/tangelo/algorithms/classical/tests/test_mp2_solver.py +++ b/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() diff --git a/tangelo/toolboxes/ansatz_generator/uccsd.py b/tangelo/toolboxes/ansatz_generator/uccsd.py index 0e1ab33ec..5c7094822 100644 --- a/tangelo/toolboxes/ansatz_generator/uccsd.py +++ b/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()