Converter for SHCI wavefunction

doc/releases/changelog-dev.md

@@ -24,8 +24,12 @@
 
 <h3>Improvements 🛠</h3>
 
+* Extended ``qml.qchem.import_state`` to import wavefunctions from SHCI classical calculations
+  performed with the Dice library.
+  [#4524](https://github.com/PennyLaneAI/pennylane/pull/4524)
+
 * Tensor-network template `qml.MPS` now supports changing `offset` between subsequent blocks for more flexibility.
- [(#4531)](https://github.com/PennyLaneAI/pennylane/pull/4531)
+  [(#4531)](https://github.com/PennyLaneAI/pennylane/pull/4531)
 
 * The qchem ``fermionic_dipole`` and ``particle_number`` functions are updated to use a
   ``FermiSentence``. The deprecated features for using tuples to represent fermionic operations are
@@ -259,6 +263,7 @@
 This release contains contributions from (in alphabetical order):
 
 Utkarsh Azad,
+Stepan Fomichev,
 Diego Guala,
 Soran Jahangiri,
 Christina Lee,

pennylane/qchem/convert.py

@@ -1007,3 +1007,112 @@ def _uccsd_state(ccsd_solver, tol=1e-15):
     dict_fcimatr = {key: value for key, value in dict_fcimatr.items() if abs(value) > tol}
 
     return dict_fcimatr
+
+
+def _shci_state(wavefunction, tol=1e-15):
+    r"""Construct a wavefunction from the SHCI wavefunction obtained from the Dice library.
+
+    The generated wavefunction is a dictionary where the keys represent a configuration, which
+    corresponds to a Slater determinant, and the values are the CI coefficients of the Slater
+    determinant. Each dictionary key is a tuple of two integers. The binary representation of these
+    integers correspond to a specific configuration: the first number represents the
+    configuration of the alpha electrons and the second number represents the configuration of the
+    beta electrons. For instance, the Hartree-Fock state :math:`|1 1 0 0 \rangle` will be
+    represented by the flipped binary string ``0011`` which is split to ``01`` and ``01`` for
+    the alpha and beta electrons. The integer corresponding to ``01`` is ``1`` and the dictionary
+    representation of the Hartree-Fock state will be ``{(1, 1): 1.0}``. The dictionary
+    representation of a state with ``0.99`` contribution from the Hartree-Fock state and ``0.01``
+    contribution from the doubly-excited state, i.e., :math:`|0 0 1 1 \rangle`, will be
+    ``{(1, 1): 0.99, (2, 2): 0.01}``.
+
+    The determinants and coefficients should be supplied externally. They are typically stored under
+    SHCI.outputfile.
+
+    Args:
+        wavefunction tuple(array[int], array[str]): determinants and coefficients in physicist notation
+        tol (float): the tolerance for discarding Slater determinants with small coefficients
+    Returns:
+        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, mcscf
+    >>> from pyscf.shciscf import shci
+    >>> import numpy as np
+    >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g', symmetry="d2h")
+    >>> myhf = scf.RHF(mol).run()
+    >>> ncas, nelecas_a, nelecas_b = mol.nao, mol.nelectron // 2, mol.nelectron // 2
+    >>> myshci = mcscf.CASCI(myhf, ncas, (nelecas_a, nelecas_b))
+    >>> initialStates = myshci.getinitialStateSHCI(myhf, (nelecas_a, nelecas_b))
+    >>> output_file = f"shci_output.out"
+    >>> myshci.fcisolver = shci.SHCI(myhf.mol)
+    >>> myshci.internal_rotation = False
+    >>> myshci.fcisolver.initialStates = initialStates
+    >>> myshci.fcisolver.outputFile = output_file
+    >>> e_shci = np.atleast_1d(myshci.kernel(verbose=5))
+    >>> wf_shci = _shci_state(myshci, tol=1e-1)
+    >>> print(wf_shci)
+    {(7, 7): 0.8874167069, (11, 11): -0.3075774156, (19, 19): -0.3075774156, (35, 35): -0.1450474361}
+    """
+
+    dets, coefs = wavefunction
+
+    xa = []
+    xb = []
+    dat = []
+
+    for coef, det in zip(coefs, dets):
+        if abs(coef) > tol:
+            bin_a, bin_b = _sitevec_to_fock(list(det), "shci")
+
+            xa.append(bin_a)
+            xb.append(bin_b)
+            dat.append(coef)
+
+    ## create the FCI matrix as a dict
+    dict_fcimatr = dict(zip(list(zip(xa, xb)), dat))
+
+    # filter based on tolerance cutoff
+    dict_fcimatr = {key: value for key, value in dict_fcimatr.items() if abs(value) > tol}
+
+    return dict_fcimatr
+
+
+def _sitevec_to_fock(det, format):
+    r"""Covert a Slater determinant from site vector to occupation number vector representation.
+
+    Args:
+        det (list): determinant in site vector representation
+        format (str): the format of the determinant
+
+    Returns:
+        tuple: tuple of integers representing binaries that correspond to occupation vectors in
+            alpha and beta spin sectors
+
+    **Example**
+
+    >>> det = [1, 2, 1, 0, 0, 2]
+    >>> _sitevec_to_fock(det, format = 'dmrg')
+    (5, 34)
+
+    >>> det = ["a", "b", "a", "0", "0", "b"]
+    >>> _sitevec_to_fock(det, format = 'shci')
+    (5, 34)
+    """
+
+    if format == "dmrg":
+        format_map = {0: "00", 1: "10", 2: "01", 3: "11"}
+    elif format == "shci":
+        format_map = {"0": "00", "a": "10", "b": "01", "2": "11"}
+
+    strab = [format_map[key] for key in det]
+
+    stra = "".join(i[0] for i in strab)
+    strb = "".join(i[1] for i in strab)
+
+    inta = int(stra[::-1], 2)
+    intb = int(strb[::-1], 2)
+
+    return inta, intb

tests/qchem/of_tests/test_convert.py

@@ -1076,3 +1076,43 @@ def test_rccsd_state(molecule, basis, symm, tol, wf_ref):
 
     assert wf_ccsd.keys() == wf_ref.keys()
     assert np.allclose(abs(np.array(list(wf_ccsd.values()))), abs(np.array(list(wf_ref.values()))))
+
+
+@pytest.mark.parametrize(
+    ("sitevec", "format", "state_ref"),
+    [([1, 2, 1, 0, 0, 2], "dmrg", (5, 34)), (["a", "b", "a", "0", "0", "b"], "shci", (5, 34))],
+)
+def test_sitevec_to_fock(sitevec, format, state_ref):
+    r"""Test that _sitevec_to_fock returns the correct state."""
+
+    state = qml.qchem.convert._sitevec_to_fock(sitevec, format)
+
+    assert state == state_ref
+
+
+@pytest.mark.parametrize(
+    ("wavefunction", "state_ref"),
+    [
+        (
+            (
+                ["02", "20"],
+                np.array([-0.10660077, 0.9943019]),
+            ),
+            {(2, 2): np.array([-0.10660077]), (1, 1): np.array([0.9943019])},
+        ),
+        (
+            (["02", "ab", "20"], np.array([0.69958765, 0.70211014, 0.1327346])),
+            {
+                (2, 2): np.array([0.69958765]),
+                (1, 2): np.array([0.70211014]),
+                (1, 1): np.array([0.1327346]),
+            },
+        ),
+    ],
+)
+def test_shci_state(wavefunction, state_ref):
+    r"""Test that _shci_state returns the correct state."""
+
+    state = qml.qchem.convert._shci_state(wavefunction)
+
+    assert state == state_ref