Add converter function for PySCF's RCISD and CCSD/UCCSD to PennyLane statevector

pennylane/qchem/convert.py

@@ -624,20 +624,26 @@ def _wfdict_to_statevector(wf_dict, norbs):
 
 
 def _rcisd_state(cisd_solver, tol=1e-15):
-    r"""
+    r"""Construct a wavefunction from PySCF's `RCISD` solver object.
 
     Args:
-        cisd_solver (PySCF UCISD Class instance): the class object representing
-            the CISD calculation in PySCF. Must have already carried out the
-            calculation, e.g. by calling .kernel() or .run().
+        cisd_solver (PySCF UCISD Class instance): the class object representing the CISD
+        calculation in PySCF.
 
     Returns:
-        cisd_solver (PySCF UCISD Class instance): the class object representing
-            the CISD calculation in PySCF. Must have already carried out the
-            calculation, e.g. by calling .kernel() or .run().
+        dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b`
+        having binary representation corresponding to the Fock occupation vector in alpha and beta
+        spin sectors, respectively, and coeff being the CI coefficients of those configurations.
 
     **Example**
 
+    >>> from pyscf import gto, scf, ci
+    >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g')
+    >>> myhf = scf.RHF(mol).run()
+    >>> myci = ci.CISD(myhf).run()
+    >>> wf_cisd = _rcisd_state(myci, tol=1e-1)
+    >>> print(wf_cisd)
+    {(1, 1): -0.9942969785398775, (2, 2): 0.10664669927602162}
     """
     mol = cisd_solver.mol
     cisdvec = cisd_solver.ci
@@ -715,16 +721,16 @@ def _rccsd_state(ccsd_solver, tol=1e-15):
     state and `0.01` contribution from the doubly-excited state will be
     `{(1, 1): 0.99, (2, 2): 0.01}`.
 
-    In the case of the coupled cluster singles and doubles calculation, in the current version, the exponential
-    ansatz :math:`\exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}}` is expanded to second order, with only
-    single and double excitation terms collected and kept. In the future this will be amended to collect also
-    terms from higher order. The expansion gives
+    In the current version, the exponential ansatz :math:`\exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}}`
+    is expanded to second order, with only single and double excitation terms collected and kept.
+    In the future this will be amended to collect also terms from higher order. The expansion gives
 
     .. math::
         \exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}} = \left[ 1 + \hat{T}_1 +
-                                    \left( \hat{T}_2 + 0.5 * \hat{T}_1^2 \right) \right] \ket{\text{HF}}
+        \left( \hat{T}_2 + 0.5 * \hat{T}_1^2 \right) \right] \ket{\text{HF}}
 
-    The coefficients in this expansion are the CI coefficients used to build the wavefunction representation.
+    The coefficients in this expansion are the CI coefficients used to build the wavefunction
+    representation.
 
     Args:
         ccsd_solver (PySCF CCSD Class instance): the class object representing the CCSD calculation in PySCF
@@ -742,10 +748,10 @@ def _rccsd_state(ccsd_solver, tol=1e-15):
     >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g')
     >>> myhf = scf.RHF(mol).run()
     >>> mycc = cc.CCSD(myhf).run()
-    >>> wf_ccsd = ccsd_state(mycc, tol=1e-1)
+    >>> wf_ccsd = _rccsd_state(mycc, tol=1e-1)
     >>> print(wf_ccsd)
     {(7, 7): 0.8886969878256522, (11, 11): -0.30584590248164206,
-                    (19, 19): -0.30584590248164145, (35, 35): -0.14507552651170982}
+    (19, 19): -0.30584590248164145, (35, 35): -0.14507552651170982}
     """
 
     mol = ccsd_solver.mol
@@ -775,7 +781,7 @@ def _rccsd_state(ccsd_solver, tol=1e-15):
         - 0.5
         * np.kron(t1b, t1b).reshape(nelec_b, nvir_b, nelec_b, nvir_b).transpose(0, 2, 1, 3).numpy()
     )
-    # this just aligns the entries with how the excitations are ordered when generated by _excited_configurations()
+    # align the entries with how the excitations are ordered when generated by _excited_configurations()
     t2ab = t2ab - 0.5 * np.kron(t1a, t1b).reshape(nelec_a, nvir_a, nelec_b, nvir_b).numpy()
 
     # numbers representing the Hartree-Fock vector, e.g., bin(ref_a)[::-1] = 1111...10...0
@@ -852,21 +858,21 @@ def _uccsd_state(ccsd_solver, tol=1e-15):
     state and `0.01` contribution from the doubly-excited state will be
     `{(1, 1): 0.99, (2, 2): 0.01}`.
 
-    In the case of the coupled cluster singles and doubles calculation, in the current version, the exponential
-    ansatz :math:`\exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}}` is expanded to second order, with only
-    single and double excitation terms collected and kept. In the future this will be amended to collect also
-    terms from higher order. The expansion gives
+    In the current version, the exponential ansatz :math:`\exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}}`
+    is expanded to second order, with only single and double excitation terms collected and kept.
+    In the future this will be amended to collect also terms from higher order. The expansion gives
 
     .. math::
         \exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}} = \left[ 1 + \hat{T}_1 +
                                     \left( \hat{T}_2 + 0.5 * \hat{T}_1^2 \right) \right] \ket{\text{HF}}
 
-    The coefficients in this expansion are the CI coefficients used to build the wavefunction representation.
+    The coefficients in this expansion are the CI coefficients used to build the wavefunction
+    representation.
 
     Args:
         ccsd_solver (PySCF UCCSD Class instance): the class object representing the UCCSD calculation in PySCF
         tol (float):  the tolerance to which the wavefunction is being built -- Slater determinants
-         with coefficients smaller than this are discarded. Default is 1e-15 (all coefficients are included).
+        with coefficients smaller than this are discarded. Default is 1e-15 (all coefficients are included).
 
     Returns:
         dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b`
@@ -879,10 +885,10 @@ def _uccsd_state(ccsd_solver, tol=1e-15):
     >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g')
     >>> myhf = scf.UHF(mol).run()
     >>> mycc = cc.UCCSD(myhf).run()
-    >>> wf_ccsd = ccsd_state(mycc, tol=1e-1)
+    >>> wf_ccsd = _uccsd_state(mycc, tol=1e-1)
     >>> print(wf_ccsd)
     {(7, 7): 0.8886970081919591, (11, 11): -0.3058459002168582,
-                    (19, 19): -0.30584590021685887, (35, 35): -0.14507552387854625}
+    (19, 19): -0.30584590021685887, (35, 35): -0.14507552387854625}
     """
 
     mol = ccsd_solver.mol
@@ -907,7 +913,7 @@ def _uccsd_state(ccsd_solver, tol=1e-15):
         - 0.5
         * np.kron(t1b, t1b).reshape(nelec_b, nvir_b, nelec_b, nvir_b).transpose(0, 2, 1, 3).numpy()
     )
-    # this just aligns the entries with how the excitations are ordered when generated by _excited_configurations()
+    # align the entries with how the excitations are ordered when generated by _excited_configurations()
     t2ab = (
         t2ab.transpose(0, 2, 1, 3)
         - 0.5 * np.kron(t1a, t1b).reshape(nelec_a, nvir_a, nelec_b, nvir_b).numpy()

tests/qchem/of_tests/test_convert.py

@@ -1000,6 +1000,65 @@ def test_import_state_error():
         _ = qchem.convert.import_state(myci, method)
 
 
+@pytest.mark.parametrize(
+    ("molecule", "basis", "symm", "tol", "wf_ref"),
+    [
+        (
+            [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]],
+            "sto6g",
+            "d2h",
+            1e-1,
+            {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602179},
+        ),
+        (
+            [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]],
+            "cc-pvdz",
+            "d2h",
+            4e-2,
+            {
+                (1, 1): 0.9919704795977625,
+                (2, 2): -0.048530356564386895,
+                (2, 8): 0.044523330850078625,
+                (4, 4): -0.050035945684911876,
+                (8, 2): 0.04452333085007864,
+                (8, 8): -0.052262303220437775,
+                (16, 16): -0.040475973747662694,
+                (32, 32): -0.040475973747662694,
+            },
+        ),
+        (
+            [["Be", (0, 0, 0)]],
+            "sto6g",
+            "d2h",
+            1e-3,
+            {
+                (3, 3): 0.9446343496981953,
+                (6, 5): 0.003359774446779245,
+                (10, 9): 0.003359774446779244,
+                (18, 17): 0.003359774446779245,
+                (5, 6): 0.003359774446779244,
+                (5, 5): -0.18938190575578503,
+                (9, 10): 0.003359774446779243,
+                (9, 9): -0.18938190575578523,
+                (17, 18): 0.003359774446779244,
+                (17, 17): -0.18938190575578503,
+            },
+        ),
+    ],
+)
+def test_rcisd_state(molecule, basis, symm, tol, wf_ref):
+    r"""Test that _rcisd_state returns the correct wavefunction."""
+
+    mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm)
+    myhf = pyscf.scf.RHF(mol).run()
+    myci = pyscf.ci.CISD(myhf).run()
+
+    wf_cisd = qchem.convert._rcisd_state(myci, tol=tol)
+
+    assert wf_cisd.keys() == wf_ref.keys()
+    assert np.allclose(np.array(list(wf_cisd.values())), np.array(list(wf_ref.values())))
+
+
 @pytest.mark.parametrize(
     ("molecule", "basis", "symm", "tol", "wf_ref"),
     [