CCSD solver psi4 (#313)

* add ccsd_solver test to psi4 testing
* get_rdm currently not supported, suggestions added in comments

.github/workflows/run_psi4_test.yml

@@ -66,6 +66,7 @@ jobs:
         cd tangelo/algorithms/classical/tests
         conda activate p4env
         pytest --doctest-modules --junitxml=junit/psi4-classical-test-results.xml test_fci_solver.py
+        pytest --doctest-modules --junitxml=junit/psi4-classical-cc-test-results.xml test_ccsd_solver.py
       if: always()
 
     - name: Upload psi4 test results

tangelo/algorithms/classical/ccsd_solver.py

@@ -15,13 +15,25 @@
 """Class performing electronic structure calculation employing the CCSD method.
 """
 
+from typing import Union, Type
+
 import numpy as np
+from sympy.combinatorics.permutations import Permutation
 
-from tangelo.helpers.utils import is_package_installed
+from tangelo.toolboxes.molecular_computation.molecule import SecondQuantizedMolecule
+from tangelo.toolboxes.molecular_computation import IntegralSolverPsi4, IntegralSolverPySCF
 from tangelo.algorithms.electronic_structure_solver import ElectronicStructureSolver
+from tangelo.helpers.utils import installed_chem_backends, is_package_installed
 
+if 'pyscf' in installed_chem_backends:
+    default_ccsd_solver = 'pyscf'
+elif 'psi4' in installed_chem_backends:
+    default_ccsd_solver = 'psi4'
+else:
+    default_ccsd_solver = None
 
-class CCSDSolver(ElectronicStructureSolver):
+
+class CCSDSolverPySCF(ElectronicStructureSolver):
     """Uses the CCSD method to solve the electronic structure problem, through
     pyscf.
 
@@ -110,3 +122,180 @@ def get_rdm(self):
                 two_rdm = np.sum((two_rdm[0], 2*two_rdm[1], two_rdm[2]), axis=0)
 
         return one_rdm, two_rdm
+
+
+class CCSDSolverPsi4(ElectronicStructureSolver):
+    """ Uses the CCSD method to solve the electronic structure problem,
+    through Psi4.
+
+    Args:
+        molecule (SecondQuantizedMolecule): The molecule to simulate.
+
+    Attributes:
+        ccwfn (psi4.core.CCWavefunction): The CCSD wavefunction (float64).
+        backend (psi4): The psi4 module
+        molecule (SecondQuantizedMolecule): The molecule with symmetry=False
+    """
+
+    def __init__(self, molecule: SecondQuantizedMolecule):
+        if not is_package_installed("psi4"):
+            raise ModuleNotFoundError(f"Using {self.__class__.__name__} requires the installation of the Psi4 package")
+
+        import psi4
+        self.backend = psi4
+        self.backend.core.clean_options()
+        self.backend.core.clean()
+        self.backend.core.clean_variables()
+        self.ccwfn = None
+
+        self.n_frozen_vir = len(molecule.frozen_virtual) if not molecule.uhf else len(molecule.frozen_virtual[0])
+        self.n_frozen_occ = len(molecule.frozen_occupied) if not molecule.uhf else len(molecule.frozen_occupied[0])
+        if not molecule.uhf:
+            self.ref = 'rhf' if molecule.spin == 0 else 'rohf'
+        else:
+            self.ref = 'uhf'
+            self.n_frozen_vir_b = len(molecule.frozen_virtual[1])
+            self.n_frozen_occ_b = len(molecule.frozen_occupied[1])
+            if (self.n_frozen_vir, self.n_frozen_occ) != (self.n_frozen_vir_b, self.n_frozen_occ_b):
+                raise ValueError(f"Tangelo does not support unequal number of alpha v. beta frozen or virtual orbitals"
+                                 f"with a UHF reference in {self.__class__.__name__}")
+
+        # Frozen orbitals must be declared before calling compute_mean_field to be saved in ref_wfn for Psi4 ccsd.
+        intsolve = IntegralSolverPsi4()
+        self.backend.set_options({'basis': molecule.basis, 'frozen_docc': [self.n_frozen_occ], 'frozen_uocc': [self.n_frozen_vir],
+                                  'reference': self.ref})
+        self.molecule = SecondQuantizedMolecule(xyz=molecule.xyz, q=molecule.q, spin=molecule.spin,
+                                                solver=intsolve,
+                                                basis=molecule.basis,
+                                                ecp=molecule.ecp,
+                                                symmetry=False,
+                                                uhf=molecule.uhf,
+                                                frozen_orbitals=molecule.frozen_orbitals)
+        self.basis = molecule.basis
+
+    def simulate(self):
+        """Perform the simulation (energy calculation) for the molecule.
+
+        Returns:
+            float: Total CCSD energy.
+        """
+        # Copy reference wavefunction and swap orbitals to obtain correct active space if necessary
+        wfn = self.backend.core.Wavefunction(self.molecule.solver.mol_nosym, self.molecule.solver.wfn.basisset())
+        wfn.deep_copy(self.molecule.solver.wfn)
+        if self.n_frozen_occ or self.n_frozen_vir:
+            if not self.molecule.uhf:
+                mo_order = self.molecule.frozen_occupied + self.molecule.active_occupied + self.molecule.active_virtual + self.molecule.frozen_virtual
+                # Obtain swap operations that will take the unordered list back to ordered with the correct active space in the middle.
+                swap_ops = Permutation(mo_order).transpositions()
+                for swap_op in swap_ops:
+                    wfn.Ca().rotate_columns(0, swap_op[0], swap_op[1], np.deg2rad(90))
+
+            else:
+
+                # Obtain swap operations that will take the unordered list back to ordered with the correct active space in the middle.
+                mo_order = (self.molecule.frozen_occupied[0] + self.molecule.active_occupied[0]
+                            + self.molecule.active_virtual[0] + self.molecule.frozen_virtual[0])
+                swap_ops = Permutation(mo_order).transpositions()
+                for swap_op in swap_ops:
+                    wfn.Ca().rotate_columns(0, swap_op[0], swap_op[1], np.deg2rad(90))
+
+                # Repeat for Beta orbitals
+                mo_order_b = (self.molecule.frozen_occupied[1] + self.molecule.active_occupied[1]
+                              + self.molecule.active_virtual[1] + self.molecule.frozen_virtual[1])
+                swap_ops = Permutation(mo_order_b).transpositions()
+                for swap_op in swap_ops:
+                    wfn.Cb().rotate_columns(0, swap_op[0], swap_op[1], np.deg2rad(90))
+
+        self.backend.set_options({'basis': self.basis, 'frozen_docc': [self.n_frozen_occ], 'frozen_uocc': [self.n_frozen_vir],
+                                  'qc_module': 'ccenergy', 'reference': self.ref})
+        energy, self.ccwfn = self.backend.energy('ccsd', molecule=self.molecule.solver.mol,
+                                                 basis=self.basis, return_wfn=True, ref_wfn=wfn)
+        return energy
+
+    def get_rdm(self):
+        """Compute the Full CI 1- and 2-particle reduced density matrices.
+
+        Returning RDMS from a CCSD calculation in Psi4 is not implemented at this time.
+
+        It may be possible to obtain the one-rdm by running a psi4 CCSD gradient calculation
+        (https://forum.psicode.org/t/saving-ccsd-density-for-read-in/2416/2)
+        Another option to obtain the one-rdm is to use pycc (https://github.com/CrawfordGroup/pycc)
+
+        Raises:
+            NotImplementedError: Returning RDMs from Psi4 in Tangelo is not supported at this time.
+        """
+        raise NotImplementedError(f"{self.__class__.__name__} does not currently support returning RDMs")
+
+
+ccsd_solver_dict = {'pyscf': CCSDSolverPySCF, 'psi4': CCSDSolverPsi4}
+
+
+def get_ccsd_solver(molecule: SecondQuantizedMolecule, solver: Union[None, str, Type[ElectronicStructureSolver]] = default_ccsd_solver, **solver_kwargs):
+    """Return requested target CCSDSolverName object.
+
+    Args:
+        molecule (SecondQuantizedMolecule) : Molecule
+        solver (string or Type[ElectronicStructureSolver] or None): Supported string identifiers can be found in
+            ccsd_solver_dict (from tangelo.algorithms.classical.ccsd_solver). Can also provide a user-defined backend (child to ElectronicStructureSolver class)
+        solver_kwargs: Other arguments that could be passed to a target. Examples are solver type (e.g. lambdacc, fnocc), Convergence options etc.
+     """
+
+    if solver is None:
+        if isinstance(molecule.solver, IntegralSolverPySCF):
+            solver = CCSDSolverPySCF
+        elif isinstance(molecule.solver, IntegralSolverPsi4):
+            solver = CCSDSolverPsi4
+        elif default_ccsd_solver is not None:
+            solver = default_ccsd_solver
+        else:
+            raise ModuleNotFoundError(f"One of the backends for {list(ccsd_solver_dict.keys())} needs to be installed to use a CCSDSolver"
+                                      "without providing a user-defined implementation.")
+
+    # If target is a string use target_dict to return built-in backend
+    elif isinstance(solver, str):
+        try:
+            solver = ccsd_solver_dict[solver.lower()]
+        except KeyError:
+            raise ValueError(f"Error: backend {solver} not supported. Available built-in options: {list(ccsd_solver_dict.keys())}")
+    elif not issubclass(solver, ElectronicStructureSolver):
+        raise TypeError(f"Target must be a str or a subclass of ElectronicStructureSolver but received class {type(solver).__name__}")
+
+    return solver(molecule, **solver_kwargs)
+
+
+class CCSDSolver(ElectronicStructureSolver):
+    """Uses the Full CI method to solve the electronic structure problem.
+
+    Args:
+        molecule (SecondQuantizedMolecule) : Molecule
+        solver (string or Type[ElectronicStructureSolver] or None): Supported string identifiers can be found in
+            available_ccsd_solvers (from tangelo.algorithms.classical.ccsd_solver). Can also provide a user-defined CCSD implementation
+            (child to ElectronicStructureSolver class)
+        solver_kwargs: Other arguments that could be passed to a target. Examples are solver type (e.g. lambdacc, fnocc), Convergence options etc.
+
+    Attributes:
+        solver (Type[ElectronicStructureSolver]): The backend specific CCSD solver
+    """
+
+    def __init__(self, molecule: SecondQuantizedMolecule, solver: Union[None, str, Type[ElectronicStructureSolver]] = default_ccsd_solver, **solver_kwargs):
+        self.solver = get_ccsd_solver(molecule, solver, **solver_kwargs)
+
+    def simulate(self):
+        """Perform the simulation (energy calculation) for the molecule.
+
+        Returns:
+            float: Total CCSD energy.
+        """
+        return self.solver.simulate()
+
+    def get_rdm(self):
+        """Compute the Full CI 1- and 2-particle reduced density matrices.
+
+        Returns:
+            numpy.array: One-particle RDM.
+            numpy.array: Two-particle RDM.
+
+        Raises:
+            RuntimeError: If method "simulate" hasn't been run.
+        """
+        return self.solver.get_rdm()

tangelo/algorithms/classical/tests/test_ccsd_solver.py

@@ -14,11 +14,11 @@
 
 import unittest
 
-from tangelo.algorithms.classical.ccsd_solver import CCSDSolver
-from tangelo.molecule_library import mol_H2_321g, mol_Be_321g, mol_H4_sto3g_uhf_a1_frozen
+from tangelo import SecondQuantizedMolecule
+from tangelo.algorithms.classical.ccsd_solver import CCSDSolver, default_ccsd_solver
+from tangelo.molecule_library import mol_H2_321g, mol_Be_321g, mol_H4_sto3g_uhf_a1_frozen, xyz_H4
 
 
-# TODO: Can we test the get_rdm method on H2 ? How do we get our reference? Whole matrix or its properties?
 class CCSDSolverTest(unittest.TestCase):
 
     def test_ccsd_h2(self):
@@ -27,10 +27,11 @@ def test_ccsd_h2(self):
         solver = CCSDSolver(mol_H2_321g)
         energy = solver.simulate()
 
-        self.assertAlmostEqual(energy, -1.1478300596229851, places=6)
+        self.assertAlmostEqual(energy, -1.1478300, places=5)
 
+    @unittest.skipIf("pyscf" != default_ccsd_solver, "Test Skipped: Only functions for pyscf \n")
     def test_ccsd_h4_uhf_a1_frozen(self):
-        """Test CCSDSolver against result from reference implementation."""
+        """Test CCSDSolver against result from reference implementation for single alpha frozen orbital and rdms returned."""
 
         solver = CCSDSolver(mol_H4_sto3g_uhf_a1_frozen)
         energy = solver.simulate()
@@ -41,13 +42,22 @@ def test_ccsd_h4_uhf_a1_frozen(self):
 
         self.assertAlmostEqual(mol_H4_sto3g_uhf_a1_frozen.energy_from_rdms(one_rdms, two_rdms), -1.95831052, places=6)
 
+    def test_ccsd_h4_uhf_different_alpha_beta_frozen(self):
+        """Test energy for case when different but equal number of alpha/beta orbitals are frozen."""
+
+        mol = SecondQuantizedMolecule(xyz_H4, q=0, spin=0, basis="3-21g", frozen_orbitals=[[2, 3, 4, 5], [2, 3, 6, 5]], symmetry=False, uhf=True)
+        solver = CCSDSolver(mol)
+        energy = solver.simulate()
+
+        self.assertAlmostEqual(energy, -2.0409800, places=5)
+
     def test_ccsd_be(self):
         """Test CCSDSolver against result from reference implementation."""
 
         solver = CCSDSolver(mol_Be_321g)
         energy = solver.simulate()
 
-        self.assertAlmostEqual(energy, -14.531416589890926, places=6)
+        self.assertAlmostEqual(energy, -14.531416, places=5)
 
     def test_ccsd_be_frozen_core(self):
         """ Test CCSDSolver against result from reference implementation, with
@@ -59,7 +69,7 @@ def test_ccsd_be_frozen_core(self):
         solver = CCSDSolver(mol_Be_321g_freeze1)
         energy = solver.simulate()
 
-        self.assertAlmostEqual(energy, -14.530687987160581, places=6)
+        self.assertAlmostEqual(energy, -14.5306879, places=5)
 
     def test_ccsd_be_as_two_levels(self):
         """ Test CCSDSolver against result from reference implementation, with
@@ -72,7 +82,7 @@ def test_ccsd_be_as_two_levels(self):
         solver = CCSDSolver(mol_Be_321g_freeze_list)
         energy = solver.simulate()
 
-        self.assertAlmostEqual(energy, -14.498104489160106, places=6)
+        self.assertAlmostEqual(energy, -14.498104, places=5)
 
     def test_ccsd_get_rdm_without_simulate(self):
         """Test that the runtime error is raised when user calls get RDM without