From 5d4d8ba58b3ea248084aed6a744b2b32523909a9 Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 31 Jul 2023 18:53:10 -0400 Subject: [PATCH 01/60] add excitation function --- pennylane/qchem/convert.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index e75e804b4d0..48fac426f18 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -354,3 +354,36 @@ def import_operator(qubit_observable, format="openfermion", wires=None, tol=1e01 return qml.dot(*_openfermion_to_pennylane(qubit_observable, wires=wires)) return qml.Hamiltonian(*_openfermion_to_pennylane(qubit_observable, wires=wires)) + + +def _excitations(electrons, orbitals): + r"""Generate all possible single and double excitations from a Hartree-Fock reference state. + + Args: + electrons (int): Number of electrons + orbitals (int) – Number of spin orbitals + + Returns: + tuple(list, list): lists with the indices of the spin orbitals involved in the excitations + + **Example** + + >>> electrons = 2 + >>> orbitals = 4 + >>> _excitations(electrons, orbitals) + ([[[0, 2]], [[0, 3]], [[1, 2]], [[1, 3]]], [[0, 1, 2, 3]]) + """ + singles_p, singles_q = [], [] + doubles_pq, doubles_rs = [], [] + + for i in range(electrons): + singles_p += [i] + doubles_pq += [[k, i] for k in range(i)] + for j in range(electrons, orbitals): + singles_q += [j] + doubles_rs += [[k, j] for k in range(electrons, j)] + + singles = [[[p] + [q]] for p in singles_p for q in singles_q] + doubles = [pq + rs for pq in doubles_pq for rs in doubles_rs] + + return singles, doubles From d75a41530ddc007041317e7ef9e821b72acbb41f Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 31 Jul 2023 19:16:05 -0400 Subject: [PATCH 02/60] add excited states --- pennylane/qchem/convert.py | 54 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 52 insertions(+), 2 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 48fac426f18..a0a4de3be7e 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -360,8 +360,8 @@ def _excitations(electrons, orbitals): r"""Generate all possible single and double excitations from a Hartree-Fock reference state. Args: - electrons (int): Number of electrons - orbitals (int) – Number of spin orbitals + electrons (int): number of electrons + orbitals (int): number of spin orbitals Returns: tuple(list, list): lists with the indices of the spin orbitals involved in the excitations @@ -387,3 +387,53 @@ def _excitations(electrons, orbitals): doubles = [pq + rs for pq in doubles_pq for rs in doubles_rs] return singles, doubles + + +def _excitated_states(electrons, orbitals, excitation): + r"""Generate excited states from a Hartree-Fock reference state. + + Args: + electrons (int): number of electrons + orbitals (int): number of spin orbitals + excitation (int): number of excited electrons + + Returns: + tuple(list, list): lists of excited states and signs obtained by reordering of the creation operators + + **Example** + + >>> electrons = 2 + >>> orbitals = 5 + >>> excitation = 2 + >>> _excitated_states(electrons, orbitals, excitation) + ([28, 26, 25], [ 1, -1, 1]) + """ + hf_state = [1] * electrons + [0] * (orbitals - electrons) + + singles, doubles = _excitations(electrons, orbitals) + + states, signs = [], [] + + if excitation == 1: + for s in singles: + state = hf_state.copy() + state[s[0]], state[s[1]] = state[s[1]], state[s[0]] + states += [state] + order = len(range(s[0], electrons - 1)) + signs.append((-1) ** order) + + if excitation == 2: + for d in doubles: + state = hf_state.copy() + state[d[0]], state[d[2]] = state[d[2]], state[d[0]] + state[d[1]], state[d[3]] = state[d[3]], state[d[1]] + states += [state] + + order_pq = len(range(d[0], electrons - 1)) + order_rs = len(range(d[1], electrons - 1)) + signs.append((-1) ** (order_pq + order_rs + 1)) + + states_str = ["".join([str(i) for i in item]) for item in states] + state_int = [int(bb[::-1], 2) for bb in states_str] + + return np.array(state_int), np.array(signs) From a33369a53d37acf2153ecd03e94cb3fe158e7761 Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 31 Jul 2023 19:25:03 -0400 Subject: [PATCH 03/60] modify var names --- pennylane/qchem/convert.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index a0a4de3be7e..ddab25041cb 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -433,7 +433,7 @@ def _excitated_states(electrons, orbitals, excitation): order_rs = len(range(d[1], electrons - 1)) signs.append((-1) ** (order_pq + order_rs + 1)) - states_str = ["".join([str(i) for i in item]) for item in states] - state_int = [int(bb[::-1], 2) for bb in states_str] + states_str = ["".join([str(i) for i in state]) for state in states] + state_int = [int(state[::-1], 2) for state in states_str] return np.array(state_int), np.array(signs) From 3d2aa331c528a5b3dd1b7898f0a9dad627c5114e Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 31 Jul 2023 19:25:21 -0400 Subject: [PATCH 04/60] modify var names --- pennylane/qchem/convert.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index ddab25041cb..62a40bc2e03 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -434,6 +434,6 @@ def _excitated_states(electrons, orbitals, excitation): signs.append((-1) ** (order_pq + order_rs + 1)) states_str = ["".join([str(i) for i in state]) for state in states] - state_int = [int(state[::-1], 2) for state in states_str] + states_int = [int(state[::-1], 2) for state in states_str] - return np.array(state_int), np.array(signs) + return np.array(states_int), np.array(signs) From d41071bc7598eb94f883676a00129a879d9fada9 Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 31 Jul 2023 22:36:43 -0400 Subject: [PATCH 05/60] add tests --- pennylane/qchem/convert.py | 2 +- tests/qchem/of_tests/test_convert.py | 48 ++++++++++++++++++++++++++-- 2 files changed, 47 insertions(+), 3 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 62a40bc2e03..1753f5a6188 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -383,7 +383,7 @@ def _excitations(electrons, orbitals): singles_q += [j] doubles_rs += [[k, j] for k in range(electrons, j)] - singles = [[[p] + [q]] for p in singles_p for q in singles_q] + singles = [[p] + [q] for p in singles_p for q in singles_q] doubles = [pq + rs for pq in doubles_pq for rs in doubles_rs] return singles, doubles diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 6c2f0f008ab..6b3ce829065 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -12,8 +12,8 @@ # See the License for the specific language governing permissions and # limitations under the License. """ -Unit tests for functions needed for converting observables obtained from external libraries to a -PennyLane observable. +Unit tests for functions needed for converting objects obtained from external libraries to a +PennyLane object. """ # pylint: disable=too-many-arguments,protected-access import os @@ -813,3 +813,47 @@ def dummy_ansatz(phis, wires): res = dummy_cost(phis) assert np.abs(res - expected_cost) < tol["atol"] + + +@pytest.mark.parametrize( + ("electrons", "orbitals", "singles_ref", "doubles_ref"), + [ + # trivial case + (2, 4, [[0, 2], [0, 3], [1, 2], [1, 3]], [[0, 1, 2, 3]]), + ], +) +def test_excitations(electrons, orbitals, singles_ref, doubles_ref): + r"""Test if the _excitations function returns correct single and double excitations.""" + singles, doubles = qchem.convert._excitations(electrons, orbitals) + assert singles == singles_ref + assert doubles == doubles_ref + + +@pytest.mark.parametrize( + ("electrons", "orbitals", "excitation", "states_ref", "signs_ref"), + [ + # reference data computed with pyscf: + # pyscf_addrs, pyscf_signs = tn_addrs_signs(orbitals, electrons, excitation) + # pyscf_state = pyscf.fci.cistring.addrs2str(orbitals, electrons, pyscf_addrs) + # pyscf_state, pyscf_signs + ( + 3, + 8, + 1, + np.array([14, 22, 38, 70, 134, 13, 21, 37, 69, 133, 11, 19, 35, 67, 131]), + np.array([1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1]), + ), + ( + 3, + 6, + 2, + np.array([28, 44, 52, 26, 42, 50, 25, 41, 49]), + np.array([1, 1, 1, -1, -1, -1, 1, 1, 1]), + ), + ], +) +def test_excitated_states(electrons, orbitals, excitation, states_ref, signs_ref): + r"""Test if the _excitated_states function returns correct states and signs.""" + states, signs = qchem.convert._excitated_states(electrons, orbitals, excitation) + assert np.allclose(states, states_ref) + assert np.allclose(signs, signs_ref) From 0a8c78d5e93b1579ea1c896f828657e8ea833b30 Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 31 Jul 2023 22:54:52 -0400 Subject: [PATCH 06/60] make code compact --- pennylane/qchem/convert.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 1753f5a6188..5c464b2f1cd 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -408,7 +408,7 @@ def _excitated_states(electrons, orbitals, excitation): >>> _excitated_states(electrons, orbitals, excitation) ([28, 26, 25], [ 1, -1, 1]) """ - hf_state = [1] * electrons + [0] * (orbitals - electrons) + hf_state = qml.qchem.hf_state(electrons, orbitals) singles, doubles = _excitations(electrons, orbitals) @@ -417,18 +417,15 @@ def _excitated_states(electrons, orbitals, excitation): if excitation == 1: for s in singles: state = hf_state.copy() - state[s[0]], state[s[1]] = state[s[1]], state[s[0]] + state[s] = state[s[::-1]] states += [state] - order = len(range(s[0], electrons - 1)) - signs.append((-1) ** order) + signs.append((-1) ** len(range(s[0], electrons - 1))) if excitation == 2: for d in doubles: state = hf_state.copy() - state[d[0]], state[d[2]] = state[d[2]], state[d[0]] - state[d[1]], state[d[3]] = state[d[3]], state[d[1]] + state[d] = state[[d[2], d[3], d[0], d[1]]] states += [state] - order_pq = len(range(d[0], electrons - 1)) order_rs = len(range(d[1], electrons - 1)) signs.append((-1) ** (order_pq + order_rs + 1)) @@ -436,4 +433,4 @@ def _excitated_states(electrons, orbitals, excitation): states_str = ["".join([str(i) for i in state]) for state in states] states_int = [int(state[::-1], 2) for state in states_str] - return np.array(states_int), np.array(signs) + return states_int, signs From 52d8ed4cc8b7b712a7902de18371f580f1ac7f5a Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 31 Jul 2023 22:57:49 -0400 Subject: [PATCH 07/60] make code compact --- pennylane/qchem/convert.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 5c464b2f1cd..dbc384a624e 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -398,7 +398,7 @@ def _excitated_states(electrons, orbitals, excitation): excitation (int): number of excited electrons Returns: - tuple(list, list): lists of excited states and signs obtained by reordering of the creation operators + tuple(array, array): arrays of excited states and signs obtained by reordering of the creation operators **Example** @@ -406,12 +406,10 @@ def _excitated_states(electrons, orbitals, excitation): >>> orbitals = 5 >>> excitation = 2 >>> _excitated_states(electrons, orbitals, excitation) - ([28, 26, 25], [ 1, -1, 1]) + (array([28, 26, 25]), array([ 1, -1, 1])) """ hf_state = qml.qchem.hf_state(electrons, orbitals) - singles, doubles = _excitations(electrons, orbitals) - states, signs = [], [] if excitation == 1: From e06c7a535cd51cd6f2e8ace74e79154b183aa9b8 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Tue, 1 Aug 2023 14:20:49 -0400 Subject: [PATCH 08/60] add ucisd wavefunction constructor --- pennylane/qchem/convert.py | 208 +++++++++++++++++++++++++++++++++++++ 1 file changed, 208 insertions(+) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index dbc384a624e..55acad25d1b 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -432,3 +432,211 @@ def _excitated_states(electrons, orbitals, excitation): states_int = [int(state[::-1], 2) for state in states_str] return states_int, signs + + +def _ucisd_state(cisd_solver, state=0, tol=1e-15): + r"""Construct a wavefunction in the wf_dict format `{ (int_a, int_b): coeff}` from + PySCF's `UCISD` Solver object. + + The format is described in detail in the docstring of :func:`~.cisd_state`. In short, + the wavefunction is returned as a Python `dict` object with the format + `{(int_a, int_b) : coeff}`, where the binary representation of integers `int_a, int_b` + gives the Fock occupation vector over the molecular orbitals of the alpha and beta + electrons, respectively and `coeff` is the value of the CI coefficient. + + The wavefunction is constructed from the configuration interaction coefficients that are + computed by PySCF's `UCISD` object upon calling its `.run()` or `'kernel()` method. In that + computation, the wavefunction is assumed to have the form + + .. math:: + + \ket{\Psi_{CISD}} = c^{(0)} \ket{\text{HF}} + \sum_{i} c^{(1)}_{i} \ket{S_i} + \sum_{i} c^{(2)}_{i} \ket_{D_{i}} + + where :math:`c^{(0,1,2)}` are the + + Parameters + ---------- + PySCF CISD solver + + state + which state to do the conversion for, if multiple were solved for + + Returns + ------- + dict + Dictionary with keys (stra,strb) being binary strings representing + states occupied by alpha and beta electrons, and values being the CI + coefficients. + """ + + norb = mol.nao + nelec = mol.nelectron + nelec_a = int((nelec + mol.spin) / 2) + nelec_b = int((nelec - mol.spin) / 2) + nocc_a, nocc_b = nelec_a, nelec_b + nvir_a, nvir_b = norb - nelec_a, norb - nelec_b + + ## extract the CI coeffs for the chosen state, if multiple were solved for by CISD + assert state in range(cisd_solver.nroots), ( + f"State requested has not " f"been solved for. Re-run with larger nroots." + ) + if cisd_solver.nroots > 1: + cisdvec = cisd_solver.ci[state] + else: + cisdvec = cisd_solver.ci + + # get number of single/double excitations, to partition the cisdvec accordingly + # simple combinatorics (i.e. N orbitals choose n electrons ) gives the answer + size_a = nocc_a * nvir_a + size_b = nocc_b * nvir_b + size_aa = int(nocc_a * (nocc_a - 1) / 2) * int(nvir_a * (nvir_a - 1) / 2) + size_bb = int(nocc_b * (nocc_b - 1) / 2) * int(nvir_b * (nvir_b - 1) / 2) + size_ab = nocc_a * nocc_b * nvir_a * nvir_b + + # cisdvec in PySCF stores excitation coefficients as a flat array [c0, c1a, c1b, c2aa, c2ab, c2bb] + sizes = [1, size_a, size_b, size_aa, size_ab, size_bb] + cumul = np.cumsum(sizes) + idxs = [0] + [slice(cumul[ii], cumul[ii + 1]) for ii in range(len(cumul) - 1)] + + # quick check that the combinatorics is correct + assert np.sum(sizes) == len(cisdvec), ( + f"Number of expected excitations {np.sum(sizes)} does not " + f"correspond to cisdvec size {len(cisdvec)}.\nCheck you combinatorics!" + ) + + # np.cumsum() up to an appropriate size_xx will give the appropriate cisdvec slice + c0, c1a, c1b, c2aa, c2ab, c2bb = [cisdvec[idx] for idx in idxs] + + # with the coefficients c0/c1/c2, can build the wavefunction through helper functions + # use row to store integers representing the Fock vector of alpha electrons, and col for beta + row = [] + col = [] + coeff = [] + + ### Hartree-Fock coefficient ### + # HF vector given by bin(ref_a)[::-1] = 1111...10...0 + ref_a = int(2**nelec_a - 1) + ref_b = int(2**nelec_b - 1) + row.extend([ref_a]) + col.extend([ref_b]) + coeff.extend([c0]) + + ### Single excitations ### + # generate arrays of integers c1x_configs, whose binary form corresponds to + # allowed single excitations, and signs from acting with the excitation operator c^+ c on the HF reference + + # a -> a + # for alpha excitations, beta configs (cols) are unchanged + c1a_configs, c1a_signs = _excitated_states(norb, nelec_a, 1) + row.extend(c1a_configs) + col.extend([ref_b] * size_a) + coeff.extend(c1a * c1a_signs) + + # b -> b + # for beta excitations, alpha configs (rows) are unchanged + c1b_configs, c1b_signs = _excitated_states(norb, nelec_b, 1) + row.extend(c1b_configs) + col.extend([ref_a] * size_b) + coeff.extend(c1b * c1b_signs) + + ### Double excitations ### + # generate arrays of integers c2xx_configs, whose binary form corresponds to + # allowed doubl excitations, and signs from acting with the excitation operator c^+ c c^+ c on the HF reference + + ## aa -> aa + # for alpha excitations, beta configs (cols) are unchanged + c2aa_configs, c2aa_signs = _excitated_states(norb, nelec_a, 2) + row.extend(c2aa_configs) + col.extend([ref_b] * size_aa) + coeff.extend(c2aa * c2aa_signs) + + ## ab -> ab + # generate all possible pairwise combinations of _single_ excitations of alpha and beta sectors + rowvals, colvals = np.array(list(product(c1a_configs, c1b_configs))).T + row.extend(rowvals) + col.extend(colvals) + coeff.extend(c2ab * np.kron(c1a_signs, c1b_signs)) + + ## bb -> bb + # for beta excitations, alpha configs (rows) are unchanged + c2bb_configs, c2bb_signs = _excitated_states(norb, nelec_b, 2) + row.extend([ref_a] * size_bb) + row.extend(c2bb_configs) + coeff.extend(c2bb * c2bb_signs) + + # zip into a Slater det dictionary of the form { (int_a, int_b): coeff } + # where binary form of int_a / int_b gives the Fock occupation vector for alpha/beta sectors + dict_fcimatr = dict(zip(list(zip(row, col)), coeff)) + + # filter based on tolerance cutoff + dict_fcimatr = {key: value for key, value in dict_fcimatr.items() if abs(value) > tol} + + return dict_fcimatr + + +def cisd_state(cisd_solver, hftype, state=0, tol=1e-15): + r"""Wrapper for constructing a wavefunction in the ``wf_dict`` format + ``{ (int_a, int_b): coeff}`` from PySCF's CISD class instance. + + The ``wf_dict`` format is a condensed, unified format for representing + wavefunctions written as a sum of Slater determinants. Each entry in + the dictionary corresponds to a single Slater determinant with coefficient + ``coeff``: the determinant is represented using the integers ``int_a, int_b``, + chosen such that such that their binary representation (i.e. ``bin(int_a)``) + corresponds to strings representing the Fock molecular orbital occupation + vector for the alpha (spin-up) sector and beta (spin-down) sector. + + For example, the Hartree-Fock state of the :math:`\text{H}_2` molecule in + the minimal basis with two electrons, split between the alpha and beta sector, + is written :math:`[1,0] x [1,0]` (high-energy orbitals on the right). In the + `wf_dict` format, it would correspond to integers (1, 1), since + `bin(1) = '0b1' -> [1,0]`. Note that the representation is reversed + (high-energy orbitals are on the left). The doubly excited state + :math:`[0,1] x [0,1]` would correspond to (2, 2), since `bin(2) = '0b10' -> [0,1]`. + + - The wrapper re-directs wavefunction construction to the appropriate method + depending on whether the restricted CISD or unrestricted CISD was used. + + - + + Args: + cisd_solver (PySCF CISD or UCISD Class instance): the class object representing + the CISD / UCISD calculation in PySCF. Must have already carried out the + calculation, e.g. by calling .kernel() or .run(). + + hftype (str): User specifies whether restricted or unrestricted CISD was performed. + The options are 'rhf' for restricted, and 'uhf' for unrestricted. + + state (int): which state to do the conversion for, if several were solved for + when running CISD / UCISD. Default is 0 (the ground state). + + 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). + + Returns: + dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` + having binary represention 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,1.4)]]) + >>> myhf = scf.UHF(mol).run() + >>> myci = ci.UCISD(myhf).run() + >>> wf_cisd = cisd_state(myci, tol=1e-2) + >>> print(wf_cisd) + {(1, 1): -0.9486830343747924, (2, 2): -0.31622855704290276} + + """ + + if hftype == "uhf": + wf = _ucisd_state(cisd_solver, state=state, tol=tol) + # elif hftype == 'rhf': + # wf = _rcisd_state(cisd_solver, state=state, tol=tol) + else: + print(f"Unknown HF reference character. Exiting") + exit() + + return wf From cebd3488a31d07191726ee631f276ec18505656d Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Tue, 1 Aug 2023 14:49:33 -0400 Subject: [PATCH 09/60] added the converter from wf_dict format to pennylane statevector --- pennylane/qchem/convert.py | 47 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 55acad25d1b..ada0ce6dbca 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -469,10 +469,13 @@ def _ucisd_state(cisd_solver, state=0, tol=1e-15): coefficients. """ + mol = cisd_solver.mol + norb = mol.nao nelec = mol.nelectron nelec_a = int((nelec + mol.spin) / 2) nelec_b = int((nelec - mol.spin) / 2) + nocc_a, nocc_b = nelec_a, nelec_b nvir_a, nvir_b = norb - nelec_a, norb - nelec_b @@ -640,3 +643,47 @@ def cisd_state(cisd_solver, hftype, state=0, tol=1e-15): exit() return wf + +def wfdict_to_statevector(wf_dict, norbs): + ''' + Intermediate converter between Overlapper's wavefunction_dict format + and Pennylane's 2^n statevector. + + Parameters + ---------- + + wf_dict : Overlapper wf_dict format + Dictionary with keys (int_a,int_b) being integers whose binary representation + shows states occupied by alpha and beta electrons; and values being the CI + coefficients. + + For example, int_a = 3 --> bin(3) = "0b11" -> [1100] means the + alpha electrons occuppy the lowest two orbitals. If int_b = 3 also, + then the total state can be written [2200], or [11110000] if showing + spin-orbitals explicitly. + + Returns + ------- + list + normalized statevector of length 2^(number_of_spinorbitals), just as in Pennylane + + ''' + + statevector = np.zeros(2**(2*norbs), dtype=complex) + + for (int_a, int_b), coeff in wf_dict.items(): + # convert to binary + bin_a, bin_b = bin(int_a)[2:], bin(int_b)[2:] + # and re-order (different PySCF vs Pennylane conventions) + bin_a, bin_b = bin_a[::-1], bin_b[::-1] + # interleave to get a spin-orbital form + bin_tot = "".join(i + j for i, j in zip(bin_a, bin_b)) + # pad to get the total number of orbitals + bin_tot_full = bin_tot + "0" * (2*norbs - len(bin_tot)) + idx = int(bin_tot_full, 2) + statevector[idx] += coeff + + # normalize + statevector = statevector / np.sqrt( np.sum( statevector**2 ) ) + + return statevector \ No newline at end of file From 1d6bca9065715d0d5caa145eaddea2a27db67128 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Tue, 1 Aug 2023 17:10:11 -0400 Subject: [PATCH 10/60] update cisd_state docstring --- pennylane/qchem/convert.py | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index ada0ce6dbca..068fca4a1cb 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -22,6 +22,7 @@ from pennylane.wires import Wires from pennylane.operation import Tensor, active_new_opmath from pennylane.pauli import pauli_sentence +from itertools import product def _process_wires(wires, n_wires=None): @@ -530,14 +531,14 @@ def _ucisd_state(cisd_solver, state=0, tol=1e-15): # a -> a # for alpha excitations, beta configs (cols) are unchanged - c1a_configs, c1a_signs = _excitated_states(norb, nelec_a, 1) + c1a_configs, c1a_signs = _excitated_states(nelec_a, norb, 1) row.extend(c1a_configs) col.extend([ref_b] * size_a) coeff.extend(c1a * c1a_signs) # b -> b # for beta excitations, alpha configs (rows) are unchanged - c1b_configs, c1b_signs = _excitated_states(norb, nelec_b, 1) + c1b_configs, c1b_signs = _excitated_states(nelec_b, norb, 1) row.extend(c1b_configs) col.extend([ref_a] * size_b) coeff.extend(c1b * c1b_signs) @@ -548,21 +549,21 @@ def _ucisd_state(cisd_solver, state=0, tol=1e-15): ## aa -> aa # for alpha excitations, beta configs (cols) are unchanged - c2aa_configs, c2aa_signs = _excitated_states(norb, nelec_a, 2) + c2aa_configs, c2aa_signs = _excitated_states(nelec_a, norb, 2) row.extend(c2aa_configs) col.extend([ref_b] * size_aa) coeff.extend(c2aa * c2aa_signs) ## ab -> ab # generate all possible pairwise combinations of _single_ excitations of alpha and beta sectors - rowvals, colvals = np.array(list(product(c1a_configs, c1b_configs))).T + rowvals, colvals = np.array( list(product(c1a_configs, c1b_configs)), dtype=int ).T.numpy() row.extend(rowvals) col.extend(colvals) - coeff.extend(c2ab * np.kron(c1a_signs, c1b_signs)) + coeff.extend( c2ab * np.kron(c1a_signs, c1b_signs).numpy() ) ## bb -> bb # for beta excitations, alpha configs (rows) are unchanged - c2bb_configs, c2bb_signs = _excitated_states(norb, nelec_b, 2) + c2bb_configs, c2bb_signs = _excitated_states(nelec_b, norb, 2) row.extend([ref_a] * size_bb) row.extend(c2bb_configs) coeff.extend(c2bb * c2bb_signs) @@ -600,7 +601,11 @@ def cisd_state(cisd_solver, hftype, state=0, tol=1e-15): - The wrapper re-directs wavefunction construction to the appropriate method depending on whether the restricted CISD or unrestricted CISD was used. - - + - The (private) method converts the CISD object's attribute `.ci`, which stores the + Slater determinant coefficients, into the `wf_dict` described above. Optionally, + if multiple states were computed with PySCF CISD, an excited state's wavefunction + can be obtained by passing `state = K`, where :math:`K` points to the :math:`K`'th + excited state. Args: cisd_solver (PySCF CISD or UCISD Class instance): the class object representing @@ -625,7 +630,7 @@ def cisd_state(cisd_solver, hftype, state=0, tol=1e-15): **Example** >>> from pyscf import gto, scf, ci - >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,1.4)]]) + >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]]) >>> myhf = scf.UHF(mol).run() >>> myci = ci.UCISD(myhf).run() >>> wf_cisd = cisd_state(myci, tol=1e-2) From 16de4e4e45e706e5ebcd42809a2a888d945fc248 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Tue, 1 Aug 2023 17:26:04 -0400 Subject: [PATCH 11/60] update _ucisd_state docstring --- pennylane/qchem/convert.py | 54 +++++++++++++++++++++++++++++--------- 1 file changed, 42 insertions(+), 12 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 068fca4a1cb..ac05cb83caf 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -453,23 +453,53 @@ def _ucisd_state(cisd_solver, state=0, tol=1e-15): \ket{\Psi_{CISD}} = c^{(0)} \ket{\text{HF}} + \sum_{i} c^{(1)}_{i} \ket{S_i} + \sum_{i} c^{(2)}_{i} \ket_{D_{i}} - where :math:`c^{(0,1,2)}` are the + where :math:`c^{(0,1,2)}` are the configuration interaction coefficients stored + as the `.ci` attribute of PySCF's `UCISD` object, and :math:`\ket{S_i}, \ket{D_i}` + are the singly and doubly excited Slater determinants, respectively. - Parameters - ---------- - PySCF CISD solver + - The first step is to partition the `.ci` attribute storing the CI coefficients into the :math:`c^{(0,1,2)}` + coefficients. To do this, the number of allowed single and double excitations is computed through combinatorics. + For example, the number of allowed single excitations of :math:`n_e` alpha electrons in :math:`N` orbitals is + :math:`\binom{N}{n_e}`. - state - which state to do the conversion for, if multiple were solved for + - The :func:`~._excited_states` function is then called to generate the configurations corresponding to allowed + single and double excitations, as well as their corresponding signs from permuting the creation operators. These + configurations are again represented by integers whose binary form gives the Fock occupation vector. + + - The configurations for alpha electrons (`row`) and beta electrons (`col`) are collected for all excitations, + and the corresponding coefficients are collected (`coeff`). These three arrays are then used to build the final + `wf_dict`, which is filtered to drop the smallest coefficients according to the `tol` argument. + + Args: + cisd_solver (PySCF UCISD Class instance): the class object representing + the UCISD calculation in PySCF. Must have already carried out the + calculation, e.g. by calling .kernel() or .run(). + + state (int): which state to do the conversion for, if several were solved for + when running UCISD. Default is 0 (the ground state). + + 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). + + Returns: + dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` + having binary represention 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)]]) + >>> myhf = scf.UHF(mol).run() + >>> myci = ci.UCISD(myhf).run() + >>> wf_cisd = cisd_state(myci, tol=1e-2) + >>> print(wf_cisd) + {(1, 1): -0.9486830343747924, (2, 2): -0.31622855704290276} - Returns - ------- - dict - Dictionary with keys (stra,strb) being binary strings representing - states occupied by alpha and beta electrons, and values being the CI - coefficients. """ + mol = cisd_solver.mol norb = mol.nao From 1d29b1f73f86b931449df817d1705bb2b5708933 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Tue, 1 Aug 2023 17:30:55 -0400 Subject: [PATCH 12/60] update docstring of final converter to statevector --- pennylane/qchem/convert.py | 32 +++++++++++--------------------- 1 file changed, 11 insertions(+), 21 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index ac05cb83caf..7d02221a89f 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -680,29 +680,19 @@ def cisd_state(cisd_solver, hftype, state=0, tol=1e-15): return wf def wfdict_to_statevector(wf_dict, norbs): - ''' - Intermediate converter between Overlapper's wavefunction_dict format - and Pennylane's 2^n statevector. + r"""Final converter between wf_dict format and Pennylane's statevector. - Parameters - ---------- - - wf_dict : Overlapper wf_dict format - Dictionary with keys (int_a,int_b) being integers whose binary representation - shows states occupied by alpha and beta electrons; and values being the CI - coefficients. + Args: + wf_dict (wf_dict format): dictionary with keys (int_a,int_b) being integers whose binary representation + shows the Fock occupation vector for alpha and beta electrons; and values being the CI + coefficients. - For example, int_a = 3 --> bin(3) = "0b11" -> [1100] means the - alpha electrons occuppy the lowest two orbitals. If int_b = 3 also, - then the total state can be written [2200], or [11110000] if showing - spin-orbitals explicitly. - - Returns - ------- - list - normalized statevector of length 2^(number_of_spinorbitals), just as in Pennylane - - ''' + norbs (int): number of molecular orbitals in the system + + Returns: + statevector (list): normalized statevector of length 2^(number_of_spinorbitals), standard in Pennylane + + """ statevector = np.zeros(2**(2*norbs), dtype=complex) From f76de224c604a1278f7e80cec6940f3eb6eb9f6d Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Tue, 1 Aug 2023 18:08:30 -0400 Subject: [PATCH 13/60] fixed the example output for cisd_state --- pennylane/qchem/convert.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 7d02221a89f..0320b96b8f7 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -490,12 +490,12 @@ def _ucisd_state(cisd_solver, state=0, tol=1e-15): **Example** >>> from pyscf import gto, scf, ci - >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]]) + >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g') >>> myhf = scf.UHF(mol).run() >>> myci = ci.UCISD(myhf).run() - >>> wf_cisd = cisd_state(myci, tol=1e-2) + >>> wf_cisd = cisd_state(myci, tol=1e-1) >>> print(wf_cisd) - {(1, 1): -0.9486830343747924, (2, 2): -0.31622855704290276} + {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} """ @@ -660,13 +660,12 @@ def cisd_state(cisd_solver, hftype, state=0, tol=1e-15): **Example** >>> from pyscf import gto, scf, ci - >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]]) + >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g') >>> myhf = scf.UHF(mol).run() >>> myci = ci.UCISD(myhf).run() - >>> wf_cisd = cisd_state(myci, tol=1e-2) + >>> wf_cisd = cisd_state(myci, tol=1e-1) >>> print(wf_cisd) - {(1, 1): -0.9486830343747924, (2, 2): -0.31622855704290276} - + {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} """ if hftype == "uhf": @@ -691,7 +690,7 @@ def wfdict_to_statevector(wf_dict, norbs): Returns: statevector (list): normalized statevector of length 2^(number_of_spinorbitals), standard in Pennylane - + """ statevector = np.zeros(2**(2*norbs), dtype=complex) From 4db87786fb6b57ca186b6fbc67430448e1ae5938 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Tue, 1 Aug 2023 21:35:31 -0400 Subject: [PATCH 14/60] added test for private _ucisd func based on h2 molecule --- tests/qchem/of_tests/test_convert.py | 38 ++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 6b3ce829065..89a262e147c 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -857,3 +857,41 @@ def test_excitated_states(electrons, orbitals, excitation, states_ref, signs_ref states, signs = qchem.convert._excitated_states(electrons, orbitals, excitation) assert np.allclose(states, states_ref) assert np.allclose(signs, signs_ref) + +from pyscf import gto +@pytest.mark.parametrize( + ("mol", "hftype", "state", "tol", "wf_ref"), + [ + # reference data computed with pyscf -- identify the FCI matrix entries + # with the correct Fock occupation vector entries + # mycasci = mcscf.CASCI(hf, ncas = orbitals, nelecas = electrons).run() + # print(mycasci.ci) + ( + gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g'), + 'uhf', + 0, + 1e-1, + {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} + ), + # ( + # gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='cc-pvdz'), + # 'uhf', + # 0, + # 4e-2, + # {(1, 1): 0.9919704779918753, (2, 2): 0.048530387857156514, (2, 8): 0.04452334895452825,\ + # (4, 4): -0.050035993275715764, (8, 2): -0.04452335710141664, (8, 8): -0.052262296228815036,\ + # (16, 16): -0.04045330124484208, (32, 32): -0.04045330124484209} + # ), + ], +) +def test_private_ucisd_state(mol, hftype, state, tol, wf_ref): + r"""Test the UCISD wavefunction constructor.""" + + from pyscf import scf, ci + + myhf = scf.UHF(mol).run() + myci = ci.UCISD(myhf).run(state=state) + + wf_cisd = qchem.convert.cisd_state(myci, hftype=hftype, state=state, tol=tol) + + assert { key: round(val, 4) for (key, val) in wf_cisd.items() } == { key: round(val, 4) for (key, val) in wf_ref.items() } \ No newline at end of file From 611715853c2e7826982f14a725ed96a3babbc3f9 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Wed, 2 Aug 2023 09:05:54 -0400 Subject: [PATCH 15/60] added cc-pvdz test --- tests/qchem/of_tests/test_convert.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 89a262e147c..c30ec1717c8 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -873,15 +873,15 @@ def test_excitated_states(electrons, orbitals, excitation, states_ref, signs_ref 1e-1, {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} ), - # ( - # gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='cc-pvdz'), - # 'uhf', - # 0, - # 4e-2, - # {(1, 1): 0.9919704779918753, (2, 2): 0.048530387857156514, (2, 8): 0.04452334895452825,\ - # (4, 4): -0.050035993275715764, (8, 2): -0.04452335710141664, (8, 8): -0.052262296228815036,\ - # (16, 16): -0.04045330124484208, (32, 32): -0.04045330124484209} - # ), + ( + gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='cc-pvdz'), + 'uhf', + 0, + 4e-2, + {(1, 1): 0.9919704779918753, (2, 2): 0.048530387857156514, (2, 8): 0.04452334895452825,\ + (4, 4): -0.050035993275715764, (8, 2): -0.04452335710141664, (8, 8): -0.052262296228815036,\ + (16, 16): -0.04045330124484208, (32, 32): -0.04045330124484209} + ), ], ) def test_private_ucisd_state(mol, hftype, state, tol, wf_ref): From 9f26454bf994dde6a484f8a018ab04c386cd4e0b Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Wed, 2 Aug 2023 09:55:16 -0400 Subject: [PATCH 16/60] added how to generate ref data for test_prive_ucisd --- pennylane/qchem/convert.py | 6 +++--- tests/qchem/of_tests/test_convert.py | 30 ++++++++++++++++------------ 2 files changed, 20 insertions(+), 16 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 0320b96b8f7..f63a390b25e 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -569,8 +569,8 @@ def _ucisd_state(cisd_solver, state=0, tol=1e-15): # b -> b # for beta excitations, alpha configs (rows) are unchanged c1b_configs, c1b_signs = _excitated_states(nelec_b, norb, 1) - row.extend(c1b_configs) - col.extend([ref_a] * size_b) + row.extend([ref_a] * size_b) + col.extend(c1b_configs) coeff.extend(c1b * c1b_signs) ### Double excitations ### @@ -595,7 +595,7 @@ def _ucisd_state(cisd_solver, state=0, tol=1e-15): # for beta excitations, alpha configs (rows) are unchanged c2bb_configs, c2bb_signs = _excitated_states(nelec_b, norb, 2) row.extend([ref_a] * size_bb) - row.extend(c2bb_configs) + col.extend(c2bb_configs) coeff.extend(c2bb * c2bb_signs) # zip into a Slater det dictionary of the form { (int_a, int_b): coeff } diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index c30ec1717c8..b85d89f4691 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -858,40 +858,44 @@ def test_excitated_states(electrons, orbitals, excitation, states_ref, signs_ref assert np.allclose(states, states_ref) assert np.allclose(signs, signs_ref) -from pyscf import gto @pytest.mark.parametrize( - ("mol", "hftype", "state", "tol", "wf_ref"), + ("atom", "basis", "hftype", "state", "tol", "wf_ref"), [ # reference data computed with pyscf -- identify the FCI matrix entries # with the correct Fock occupation vector entries - # mycasci = mcscf.CASCI(hf, ncas = orbitals, nelecas = electrons).run() - # print(mycasci.ci) + # mycasci = mcscf.UCASCI(myhf, norbs, (nelec_a, nelec_b)).run() + # sparse_cascimatr = coo_matrix(mycasci.ci, shape=np.shape(cascivec), dtype=float) + # row, col, dat = sparse_cascimatr.row, sparse_cascimatr.col, sparse_cascimatr.data + # strs_row, strs_col = addrs2str(norbs, nelec_a, row), addrs2str(norbs, nelec_b, col) + # wf_ref = dict(zip(list(zip(strs_row, strs_col)), dat)) ( - gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g'), + [['H', (0, 0, 0)], ['H', (0,0,0.71)]], + 'sto6g', 'uhf', 0, 1e-1, {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} ), ( - gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='cc-pvdz'), + [['H', (0, 0, 0)], ['H', (0, 0, 0.71)]], + 'cc-pvdz', 'uhf', 0, 4e-2, - {(1, 1): 0.9919704779918753, (2, 2): 0.048530387857156514, (2, 8): 0.04452334895452825,\ - (4, 4): -0.050035993275715764, (8, 2): -0.04452335710141664, (8, 8): -0.052262296228815036,\ - (16, 16): -0.04045330124484208, (32, 32): -0.04045330124484209} - ), + {(1, 1): 0.9919704801055201, (2, 2): 0.04853035090478132, (2, 8): -0.04452332640264183, \ + (4, 4): -0.05003594588609278, (8, 2): 0.044523334555907366, (8, 8): -0.05226230845395026}), ], ) -def test_private_ucisd_state(mol, hftype, state, tol, wf_ref): +def test_private_ucisd_state(atom, basis, hftype, state, tol, wf_ref): r"""Test the UCISD wavefunction constructor.""" - from pyscf import scf, ci + from pyscf import gto, scf, ci + mol = gto.M(atom=atom, basis=basis) myhf = scf.UHF(mol).run() myci = ci.UCISD(myhf).run(state=state) wf_cisd = qchem.convert.cisd_state(myci, hftype=hftype, state=state, tol=tol) - assert { key: round(val, 4) for (key, val) in wf_cisd.items() } == { key: round(val, 4) for (key, val) in wf_ref.items() } \ No newline at end of file + assert { key: round(val, 4) for (key, val) in wf_cisd.items() } ==\ + { key: round(val, 4) for (key, val) in wf_ref.items() } \ No newline at end of file From d7168fd8e0bf61420a5b2cff79bfe57fcf7e3827 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Wed, 2 Aug 2023 09:56:02 -0400 Subject: [PATCH 17/60] [skip ci] fixed typos --- tests/qchem/of_tests/test_convert.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index b85d89f4691..c7c8c1735ff 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -897,5 +897,5 @@ def test_private_ucisd_state(atom, basis, hftype, state, tol, wf_ref): wf_cisd = qchem.convert.cisd_state(myci, hftype=hftype, state=state, tol=tol) - assert { key: round(val, 4) for (key, val) in wf_cisd.items() } ==\ - { key: round(val, 4) for (key, val) in wf_ref.items() } \ No newline at end of file + assert { key: round(val, 4) for (key, val) in wf_cisd.items() } == \ + { key: round(val, 4) for (key, val) in wf_ref.items() } \ No newline at end of file From 414965a98a8093671e2886480d7cd1d9830d2ae1 Mon Sep 17 00:00:00 2001 From: soranjh Date: Wed, 2 Aug 2023 14:17:51 -0400 Subject: [PATCH 18/60] modify test to account for sign change --- tests/qchem/of_tests/test_convert.py | 30 ++++++++++++++++++---------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index c7c8c1735ff..0f100d346ce 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -858,6 +858,7 @@ def test_excitated_states(electrons, orbitals, excitation, states_ref, signs_ref assert np.allclose(states, states_ref) assert np.allclose(signs, signs_ref) + @pytest.mark.parametrize( ("atom", "basis", "hftype", "state", "tol", "wf_ref"), [ @@ -869,21 +870,28 @@ def test_excitated_states(electrons, orbitals, excitation, states_ref, signs_ref # strs_row, strs_col = addrs2str(norbs, nelec_a, row), addrs2str(norbs, nelec_b, col) # wf_ref = dict(zip(list(zip(strs_row, strs_col)), dat)) ( - [['H', (0, 0, 0)], ['H', (0,0,0.71)]], - 'sto6g', - 'uhf', + [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], + "sto6g", + "uhf", 0, 1e-1, - {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} + {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159}, ), ( - [['H', (0, 0, 0)], ['H', (0, 0, 0.71)]], - 'cc-pvdz', - 'uhf', + [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], + "cc-pvdz", + "uhf", 0, 4e-2, - {(1, 1): 0.9919704801055201, (2, 2): 0.04853035090478132, (2, 8): -0.04452332640264183, \ - (4, 4): -0.05003594588609278, (8, 2): 0.044523334555907366, (8, 8): -0.05226230845395026}), + { + (1, 1): 0.9919704801055201, + (2, 2): 0.04853035090478132, + (2, 8): -0.04452332640264183, + (4, 4): -0.05003594588609278, + (8, 2): 0.044523334555907366, + (8, 8): -0.05226230845395026, + }, + ), ], ) def test_private_ucisd_state(atom, basis, hftype, state, tol, wf_ref): @@ -897,5 +905,5 @@ def test_private_ucisd_state(atom, basis, hftype, state, tol, wf_ref): wf_cisd = qchem.convert.cisd_state(myci, hftype=hftype, state=state, tol=tol) - assert { key: round(val, 4) for (key, val) in wf_cisd.items() } == \ - { key: round(val, 4) for (key, val) in wf_ref.items() } \ No newline at end of file + assert wf_cisd.keys() == wf_ref.keys() + assert np.allclose(abs(np.array(list(wf_cisd.values()))), abs(np.array(list(wf_ref.values())))) From 782768c4212fa45ce2cd7a4a17b78b3bb31e5665 Mon Sep 17 00:00:00 2001 From: soranjh Date: Wed, 2 Aug 2023 15:18:00 -0400 Subject: [PATCH 19/60] remove comments and modify docstring --- pennylane/qchem/convert.py | 230 +++++++++---------------------------- 1 file changed, 55 insertions(+), 175 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index f63a390b25e..9d6c1564fe5 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -22,7 +22,7 @@ from pennylane.wires import Wires from pennylane.operation import Tensor, active_new_opmath from pennylane.pauli import pauli_sentence -from itertools import product +from itertools import product def _process_wires(wires, n_wires=None): @@ -435,52 +435,25 @@ def _excitated_states(electrons, orbitals, excitation): return states_int, signs -def _ucisd_state(cisd_solver, state=0, tol=1e-15): - r"""Construct a wavefunction in the wf_dict format `{ (int_a, int_b): coeff}` from - PySCF's `UCISD` Solver object. +def _ucisd_state(cisd_solver, tol=1e-15): + r"""Construct a wavefunction from PySCF's `UCISD` Solver object. - The format is described in detail in the docstring of :func:`~.cisd_state`. In short, - the wavefunction is returned as a Python `dict` object with the format - `{(int_a, int_b) : coeff}`, where the binary representation of integers `int_a, int_b` - gives the Fock occupation vector over the molecular orbitals of the alpha and beta - electrons, respectively and `coeff` is the value of the CI coefficient. - - The wavefunction is constructed from the configuration interaction coefficients that are - computed by PySCF's `UCISD` object upon calling its `.run()` or `'kernel()` method. In that - computation, the wavefunction is assumed to have the form - - .. math:: - - \ket{\Psi_{CISD}} = c^{(0)} \ket{\text{HF}} + \sum_{i} c^{(1)}_{i} \ket{S_i} + \sum_{i} c^{(2)}_{i} \ket_{D_{i}} - - where :math:`c^{(0,1,2)}` are the configuration interaction coefficients stored - as the `.ci` attribute of PySCF's `UCISD` object, and :math:`\ket{S_i}, \ket{D_i}` - are the singly and doubly excited Slater determinants, respectively. - - - The first step is to partition the `.ci` attribute storing the CI coefficients into the :math:`c^{(0,1,2)}` - coefficients. To do this, the number of allowed single and double excitations is computed through combinatorics. - For example, the number of allowed single excitations of :math:`n_e` alpha electrons in :math:`N` orbitals is - :math:`\binom{N}{n_e}`. - - - The :func:`~._excited_states` function is then called to generate the configurations corresponding to allowed - single and double excitations, as well as their corresponding signs from permuting the creation operators. These - configurations are again represented by integers whose binary form gives the Fock occupation vector. - - - The configurations for alpha electrons (`row`) and beta electrons (`col`) are collected for all excitations, - and the corresponding coefficients are collected (`coeff`). These three arrays are then used to build the final - `wf_dict`, which is filtered to drop the smallest coefficients according to the `tol` argument. + The wavefunction is represented as a dictionary where the keys are tuples representing a + configuration, which corresponds to a Slater determinant, and the values are the CI coefficient + corresponding to that Slater determinant. Each dictionary key is a tuple of two integers. The + binary representation of these integers correspond to a specific configuration, or Slater + determinant. The first number represents the configuration of alpha electrons and the second + number represents the beta electrons. For instance, the Hartree-Fock state + :math:`|1 1 0 0 \rangle` will be represented by the binary string `0011`. This string can be + splited to `01` and `01` for the alpha and beta electrons. The integer corresponding to `01` is + `1`. Then the dictionary representation of a state with `0.99` contribution of the Hartree-Fock + state and `0.01` contribution from the doubly-excited state will be + `{(1, 1): 0.99, (2, 2): 0.01}`. Args: - cisd_solver (PySCF UCISD Class instance): the class object representing - the UCISD calculation in PySCF. Must have already carried out the - calculation, e.g. by calling .kernel() or .run(). - - state (int): which state to do the conversion for, if several were solved for - when running UCISD. Default is 0 (the ground state). - - 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). + cisd_solver (PySCF UCISD Class instance): the class object representing the UCISD 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). Returns: dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` @@ -496,111 +469,54 @@ def _ucisd_state(cisd_solver, state=0, tol=1e-15): >>> wf_cisd = cisd_state(myci, tol=1e-1) >>> print(wf_cisd) {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} - """ - - mol = cisd_solver.mol + cisdvec = cisd_solver.ci norb = mol.nao nelec = mol.nelectron nelec_a = int((nelec + mol.spin) / 2) nelec_b = int((nelec - mol.spin) / 2) - nocc_a, nocc_b = nelec_a, nelec_b nvir_a, nvir_b = norb - nelec_a, norb - nelec_b - ## extract the CI coeffs for the chosen state, if multiple were solved for by CISD - assert state in range(cisd_solver.nroots), ( - f"State requested has not " f"been solved for. Re-run with larger nroots." - ) - if cisd_solver.nroots > 1: - cisdvec = cisd_solver.ci[state] - else: - cisdvec = cisd_solver.ci - - # get number of single/double excitations, to partition the cisdvec accordingly - # simple combinatorics (i.e. N orbitals choose n electrons ) gives the answer - size_a = nocc_a * nvir_a - size_b = nocc_b * nvir_b - size_aa = int(nocc_a * (nocc_a - 1) / 2) * int(nvir_a * (nvir_a - 1) / 2) - size_bb = int(nocc_b * (nocc_b - 1) / 2) * int(nvir_b * (nvir_b - 1) / 2) - size_ab = nocc_a * nocc_b * nvir_a * nvir_b + size_a, size_b = nelec_a * nvir_a, nelec_b * nvir_b + size_aa = int(nelec_a * (nelec_a - 1) / 2) * int(nvir_a * (nvir_a - 1) / 2) + size_bb = int(nelec_b * (nelec_b - 1) / 2) * int(nvir_b * (nvir_b - 1) / 2) + size_ab = nelec_a * nelec_b * nvir_a * nvir_b - # cisdvec in PySCF stores excitation coefficients as a flat array [c0, c1a, c1b, c2aa, c2ab, c2bb] sizes = [1, size_a, size_b, size_aa, size_ab, size_bb] cumul = np.cumsum(sizes) idxs = [0] + [slice(cumul[ii], cumul[ii + 1]) for ii in range(len(cumul) - 1)] - - # quick check that the combinatorics is correct - assert np.sum(sizes) == len(cisdvec), ( - f"Number of expected excitations {np.sum(sizes)} does not " - f"correspond to cisdvec size {len(cisdvec)}.\nCheck you combinatorics!" - ) - - # np.cumsum() up to an appropriate size_xx will give the appropriate cisdvec slice c0, c1a, c1b, c2aa, c2ab, c2bb = [cisdvec[idx] for idx in idxs] - # with the coefficients c0/c1/c2, can build the wavefunction through helper functions - # use row to store integers representing the Fock vector of alpha electrons, and col for beta - row = [] - col = [] - coeff = [] - - ### Hartree-Fock coefficient ### - # HF vector given by bin(ref_a)[::-1] = 1111...10...0 + # numbers representing the Hartree-Fock vector, e.g., bin(ref_a)[::-1] = 1111...10...0 ref_a = int(2**nelec_a - 1) ref_b = int(2**nelec_b - 1) - row.extend([ref_a]) - col.extend([ref_b]) - coeff.extend([c0]) - ### Single excitations ### - # generate arrays of integers c1x_configs, whose binary form corresponds to - # allowed single excitations, and signs from acting with the excitation operator c^+ c on the HF reference + dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [c0])) - # a -> a - # for alpha excitations, beta configs (cols) are unchanged + # alpha -> alpha excitations c1a_configs, c1a_signs = _excitated_states(nelec_a, norb, 1) - row.extend(c1a_configs) - col.extend([ref_b] * size_a) - coeff.extend(c1a * c1a_signs) + dict_fcimatr.update(dict(zip(list(zip(c1a_configs, [ref_b] * size_a)), c1a * c1a_signs))) - # b -> b - # for beta excitations, alpha configs (rows) are unchanged + # beta -> beta excitations c1b_configs, c1b_signs = _excitated_states(nelec_b, norb, 1) - row.extend([ref_a] * size_b) - col.extend(c1b_configs) - coeff.extend(c1b * c1b_signs) + dict_fcimatr.update(dict(zip(list(zip(c1b_configs, [ref_a] * size_b)), c1b * c1b_signs))) - ### Double excitations ### - # generate arrays of integers c2xx_configs, whose binary form corresponds to - # allowed doubl excitations, and signs from acting with the excitation operator c^+ c c^+ c on the HF reference - - ## aa -> aa - # for alpha excitations, beta configs (cols) are unchanged + # alpha, alpha -> alpha, alpha excitations c2aa_configs, c2aa_signs = _excitated_states(nelec_a, norb, 2) - row.extend(c2aa_configs) - col.extend([ref_b] * size_aa) - coeff.extend(c2aa * c2aa_signs) - - ## ab -> ab - # generate all possible pairwise combinations of _single_ excitations of alpha and beta sectors - rowvals, colvals = np.array( list(product(c1a_configs, c1b_configs)), dtype=int ).T.numpy() - row.extend(rowvals) - col.extend(colvals) - coeff.extend( c2ab * np.kron(c1a_signs, c1b_signs).numpy() ) - - ## bb -> bb - # for beta excitations, alpha configs (rows) are unchanged - c2bb_configs, c2bb_signs = _excitated_states(nelec_b, norb, 2) - row.extend([ref_a] * size_bb) - col.extend(c2bb_configs) - coeff.extend(c2bb * c2bb_signs) + dict_fcimatr.update(dict(zip(list(zip(c2aa_configs, [ref_b] * size_aa)), c2aa * c2aa_signs))) - # zip into a Slater det dictionary of the form { (int_a, int_b): coeff } - # where binary form of int_a / int_b gives the Fock occupation vector for alpha/beta sectors - dict_fcimatr = dict(zip(list(zip(row, col)), coeff)) + # alpha, beta -> alpha, beta excitations + rowvals, colvals = np.array(list(product(c1a_configs, c1b_configs)), dtype=int).T.numpy() + dict_fcimatr.update( + dict(zip(list(zip(rowvals, colvals)), c2ab * np.kron(c1a_signs, c1b_signs).numpy())) + ) + + # beta, beta -> beta, beta excitations + c2bb_configs, c2bb_signs = _excitated_states(nelec_b, norb, 2) + dict_fcimatr.update(dict(zip(list(zip([ref_a] * size_bb, c2bb_configs)), c2bb * c2bb_signs))) # filter based on tolerance cutoff dict_fcimatr = {key: value for key, value in dict_fcimatr.items() if abs(value) > tol} @@ -609,45 +525,15 @@ def _ucisd_state(cisd_solver, state=0, tol=1e-15): def cisd_state(cisd_solver, hftype, state=0, tol=1e-15): - r"""Wrapper for constructing a wavefunction in the ``wf_dict`` format - ``{ (int_a, int_b): coeff}`` from PySCF's CISD class instance. - - The ``wf_dict`` format is a condensed, unified format for representing - wavefunctions written as a sum of Slater determinants. Each entry in - the dictionary corresponds to a single Slater determinant with coefficient - ``coeff``: the determinant is represented using the integers ``int_a, int_b``, - chosen such that such that their binary representation (i.e. ``bin(int_a)``) - corresponds to strings representing the Fock molecular orbital occupation - vector for the alpha (spin-up) sector and beta (spin-down) sector. - - For example, the Hartree-Fock state of the :math:`\text{H}_2` molecule in - the minimal basis with two electrons, split between the alpha and beta sector, - is written :math:`[1,0] x [1,0]` (high-energy orbitals on the right). In the - `wf_dict` format, it would correspond to integers (1, 1), since - `bin(1) = '0b1' -> [1,0]`. Note that the representation is reversed - (high-energy orbitals are on the left). The doubly excited state - :math:`[0,1] x [0,1]` would correspond to (2, 2), since `bin(2) = '0b10' -> [0,1]`. - - - The wrapper re-directs wavefunction construction to the appropriate method - depending on whether the restricted CISD or unrestricted CISD was used. - - - The (private) method converts the CISD object's attribute `.ci`, which stores the - Slater determinant coefficients, into the `wf_dict` described above. Optionally, - if multiple states were computed with PySCF CISD, an excited state's wavefunction - can be obtained by passing `state = K`, where :math:`K` points to the :math:`K`'th - excited state. + r"""Construct a wavefunction from PySCF output. Args: - cisd_solver (PySCF CISD or UCISD Class instance): the class object representing - the CISD / UCISD calculation in PySCF. Must have already carried out the - calculation, e.g. by calling .kernel() or .run(). + cisd_solver (PySCF CISD or UCISD Class instance): the class object representing the + CISD / UCISD calculation in PySCF - hftype (str): User specifies whether restricted or unrestricted CISD was performed. + hftype (str): keyword specifying whether restricted or unrestricted CISD was performed. The options are 'rhf' for restricted, and 'uhf' for unrestricted. - state (int): which state to do the conversion for, if several were solved for - when running CISD / UCISD. Default is 0 (the ground state). - 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). @@ -669,45 +555,39 @@ def cisd_state(cisd_solver, hftype, state=0, tol=1e-15): """ if hftype == "uhf": - wf = _ucisd_state(cisd_solver, state=state, tol=tol) - # elif hftype == 'rhf': - # wf = _rcisd_state(cisd_solver, state=state, tol=tol) + wf_dict = _ucisd_state(cisd_solver, state=state, tol=tol) else: - print(f"Unknown HF reference character. Exiting") - exit() + raise ValueError("Only restricted, 'rhf', and unrestricted, 'uhf', are supported.") + wf = wfdict_to_statevector(wf_dict, cisd_solver.mol.nao) return wf + def wfdict_to_statevector(wf_dict, norbs): - r"""Final converter between wf_dict format and Pennylane's statevector. + r"""Convert a wavefunction in sparce dictionary format to a Pennylane's statevector. Args: wf_dict (wf_dict format): dictionary with keys (int_a,int_b) being integers whose binary representation shows the Fock occupation vector for alpha and beta electrons; and values being the CI - coefficients. - + coefficients. + norbs (int): number of molecular orbitals in the system - + Returns: statevector (list): normalized statevector of length 2^(number_of_spinorbitals), standard in Pennylane """ - statevector = np.zeros(2**(2*norbs), dtype=complex) + statevector = np.zeros(2 ** (2 * norbs), dtype=complex) for (int_a, int_b), coeff in wf_dict.items(): - # convert to binary bin_a, bin_b = bin(int_a)[2:], bin(int_b)[2:] - # and re-order (different PySCF vs Pennylane conventions) bin_a, bin_b = bin_a[::-1], bin_b[::-1] - # interleave to get a spin-orbital form bin_tot = "".join(i + j for i, j in zip(bin_a, bin_b)) - # pad to get the total number of orbitals - bin_tot_full = bin_tot + "0" * (2*norbs - len(bin_tot)) + bin_tot_full = bin_tot + "0" * (2 * norbs - len(bin_tot)) idx = int(bin_tot_full, 2) statevector[idx] += coeff - # normalize - statevector = statevector / np.sqrt( np.sum( statevector**2 ) ) + statevector = statevector / np.sqrt(np.sum(statevector**2)) - return statevector \ No newline at end of file + return statevector From 4acf9c8e7741ba61c1508c0dd8e0870458aa8464 Mon Sep 17 00:00:00 2001 From: soranjh <40344468+soranjh@users.noreply.github.com> Date: Wed, 2 Aug 2023 15:28:15 -0400 Subject: [PATCH 20/60] Add functions to return excited basis states (#4417) --- pennylane/qchem/convert.py | 30 +++++++++++++++++++++++++++- tests/qchem/of_tests/test_convert.py | 10 +++++----- 2 files changed, 34 insertions(+), 6 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 9d6c1564fe5..91b4827f85c 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -360,6 +360,25 @@ def import_operator(qubit_observable, format="openfermion", wires=None, tol=1e01 def _excitations(electrons, orbitals): r"""Generate all possible single and double excitations from a Hartree-Fock reference state. + This function is a more performant version of `qchem.excitations`, where the order of the + generated excitations is consistent with PySCF. + + Single and double excitations can be generated by acting with the operators + :math:`\hat T_1` and :math:`\hat T_2` on the Hartree-Fock reference state: + + .. math:: + + && \hat{T}_1 = \sum_{p \in \mathrm{occ} \\ q \in \mathrm{unocc}} + \hat{c}_q^\dagger \hat{c}_p \\ + && \hat{T}_2 = \sum_{p>q \in \mathrm{occ} \\ r>s \in + \mathrm{unocc}} \hat{c}_r^\dagger \hat{c}_s^\dagger \hat{c}_p \hat{c}_q. + + + In the equations above the indices :math:`p, q, r, s` run over the + occupied (occ) and unoccupied (unocc) *spin* orbitals and :math:`\hat c` and + :math:`\hat c^\dagger` are the electron annihilation and creation operators, + respectively. + Args: electrons (int): number of electrons orbitals (int): number of spin orbitals @@ -390,9 +409,18 @@ def _excitations(electrons, orbitals): return singles, doubles -def _excitated_states(electrons, orbitals, excitation): +def _excited_configurations(electrons, orbitals, excitation): r"""Generate excited states from a Hartree-Fock reference state. + This function generates excited states in the form of binary strings or integers representing + them, from a Hartree-Fock (HF) reference state. The HF state is assumed to be + :math:`|1 1 ...1 0 ... 0 0 \rangle` where the number of :math:`1` and :math:`0` elements + represents the number of occupied and unoccupied spin orbitals, respectively. The string + representation of the state is obtained by converting the occupation-number vector to a string, + e.g., ``111000`` represents :math:`|1 1 1 0 0 0 \rangle. The number representation of the state + is the integer corresonding to the binary number, e.g., ``111000`` is represented by + :math:`int('111000', 2) = 56`. + Args: electrons (int): number of electrons orbitals (int): number of spin orbitals diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 0f100d346ce..1e27f09c56c 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -818,7 +818,7 @@ def dummy_ansatz(phis, wires): @pytest.mark.parametrize( ("electrons", "orbitals", "singles_ref", "doubles_ref"), [ - # trivial case + # trivial case, e.g., H2/STO-3G (2, 4, [[0, 2], [0, 3], [1, 2], [1, 3]], [[0, 1, 2, 3]]), ], ) @@ -833,7 +833,7 @@ def test_excitations(electrons, orbitals, singles_ref, doubles_ref): ("electrons", "orbitals", "excitation", "states_ref", "signs_ref"), [ # reference data computed with pyscf: - # pyscf_addrs, pyscf_signs = tn_addrs_signs(orbitals, electrons, excitation) + # pyscf_addrs, pyscf_signs = pyscf.ci.cisd.tn_addrs_signs(orbitals, electrons, excitation) # pyscf_state = pyscf.fci.cistring.addrs2str(orbitals, electrons, pyscf_addrs) # pyscf_state, pyscf_signs ( @@ -852,9 +852,9 @@ def test_excitations(electrons, orbitals, singles_ref, doubles_ref): ), ], ) -def test_excitated_states(electrons, orbitals, excitation, states_ref, signs_ref): - r"""Test if the _excitated_states function returns correct states and signs.""" - states, signs = qchem.convert._excitated_states(electrons, orbitals, excitation) +def test_excited_configurations(electrons, orbitals, excitation, states_ref, signs_ref): + r"""Test if the _excited_configurations function returns correct states and signs.""" + states, signs = qchem.convert._excited_configurations(electrons, orbitals, excitation) assert np.allclose(states, states_ref) assert np.allclose(signs, signs_ref) From 9be65fa9305bfbce0df666d3f80ec1968fcc6a4b Mon Sep 17 00:00:00 2001 From: soranjh Date: Wed, 2 Aug 2023 15:37:50 -0400 Subject: [PATCH 21/60] correct test --- pennylane/qchem/convert.py | 16 ++++++++-------- tests/qchem/of_tests/test_convert.py | 16 ++++++---------- 2 files changed, 14 insertions(+), 18 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 91b4827f85c..654b3516b68 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -15,14 +15,14 @@ This module contains the functions for converting an external operator to a Pennylane operator. """ import warnings +from itertools import product # pylint: disable=import-outside-toplevel import pennylane as qml from pennylane import numpy as np -from pennylane.wires import Wires from pennylane.operation import Tensor, active_new_opmath from pennylane.pauli import pauli_sentence -from itertools import product +from pennylane.wires import Wires def _process_wires(wires, n_wires=None): @@ -525,15 +525,15 @@ def _ucisd_state(cisd_solver, tol=1e-15): dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [c0])) # alpha -> alpha excitations - c1a_configs, c1a_signs = _excitated_states(nelec_a, norb, 1) + c1a_configs, c1a_signs = _excited_configurations(nelec_a, norb, 1) dict_fcimatr.update(dict(zip(list(zip(c1a_configs, [ref_b] * size_a)), c1a * c1a_signs))) # beta -> beta excitations - c1b_configs, c1b_signs = _excitated_states(nelec_b, norb, 1) + c1b_configs, c1b_signs = _excited_configurations(nelec_b, norb, 1) dict_fcimatr.update(dict(zip(list(zip(c1b_configs, [ref_a] * size_b)), c1b * c1b_signs))) # alpha, alpha -> alpha, alpha excitations - c2aa_configs, c2aa_signs = _excitated_states(nelec_a, norb, 2) + c2aa_configs, c2aa_signs = _excited_configurations(nelec_a, norb, 2) dict_fcimatr.update(dict(zip(list(zip(c2aa_configs, [ref_b] * size_aa)), c2aa * c2aa_signs))) # alpha, beta -> alpha, beta excitations @@ -543,7 +543,7 @@ def _ucisd_state(cisd_solver, tol=1e-15): ) # beta, beta -> beta, beta excitations - c2bb_configs, c2bb_signs = _excitated_states(nelec_b, norb, 2) + c2bb_configs, c2bb_signs = _excited_configurations(nelec_b, norb, 2) dict_fcimatr.update(dict(zip(list(zip([ref_a] * size_bb, c2bb_configs)), c2bb * c2bb_signs))) # filter based on tolerance cutoff @@ -552,7 +552,7 @@ def _ucisd_state(cisd_solver, tol=1e-15): return dict_fcimatr -def cisd_state(cisd_solver, hftype, state=0, tol=1e-15): +def cisd_state(cisd_solver, hftype, tol=1e-15): r"""Construct a wavefunction from PySCF output. Args: @@ -583,7 +583,7 @@ def cisd_state(cisd_solver, hftype, state=0, tol=1e-15): """ if hftype == "uhf": - wf_dict = _ucisd_state(cisd_solver, state=state, tol=tol) + wf_dict = _ucisd_state(cisd_solver, tol=tol) else: raise ValueError("Only restricted, 'rhf', and unrestricted, 'uhf', are supported.") wf = wfdict_to_statevector(wf_dict, cisd_solver.mol.nao) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 1e27f09c56c..0cb4ca0cc9a 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -24,7 +24,7 @@ import pennylane as qml from pennylane import numpy as np from pennylane import qchem -from pennylane.operation import enable_new_opmath, disable_new_opmath +from pennylane.operation import disable_new_opmath, enable_new_opmath openfermion = pytest.importorskip("openfermion") openfermionpyscf = pytest.importorskip("openfermionpyscf") @@ -860,7 +860,7 @@ def test_excited_configurations(electrons, orbitals, excitation, states_ref, sig @pytest.mark.parametrize( - ("atom", "basis", "hftype", "state", "tol", "wf_ref"), + ("molecule", "basis", "tol", "wf_ref"), [ # reference data computed with pyscf -- identify the FCI matrix entries # with the correct Fock occupation vector entries @@ -872,16 +872,12 @@ def test_excited_configurations(electrons, orbitals, excitation, states_ref, sig ( [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], "sto6g", - "uhf", - 0, 1e-1, {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159}, ), ( [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], "cc-pvdz", - "uhf", - 0, 4e-2, { (1, 1): 0.9919704801055201, @@ -894,16 +890,16 @@ def test_excited_configurations(electrons, orbitals, excitation, states_ref, sig ), ], ) -def test_private_ucisd_state(atom, basis, hftype, state, tol, wf_ref): +def test_ucisd_state(molecule, basis, tol, wf_ref): r"""Test the UCISD wavefunction constructor.""" from pyscf import gto, scf, ci - mol = gto.M(atom=atom, basis=basis) + mol = gto.M(atom=molecule, basis=basis) myhf = scf.UHF(mol).run() - myci = ci.UCISD(myhf).run(state=state) + myci = ci.UCISD(myhf).run() - wf_cisd = qchem.convert.cisd_state(myci, hftype=hftype, state=state, tol=tol) + wf_cisd = qchem.convert._ucisd_state(myci, tol=tol) assert wf_cisd.keys() == wf_ref.keys() assert np.allclose(abs(np.array(list(wf_cisd.values()))), abs(np.array(list(wf_ref.values())))) From a4eceafbfb5f4eb1f54076e699b69642b3b4e355 Mon Sep 17 00:00:00 2001 From: soranjh Date: Wed, 2 Aug 2023 15:39:48 -0400 Subject: [PATCH 22/60] correct import --- tests/qchem/of_tests/test_convert.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 0cb4ca0cc9a..d33f9a9bb8a 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -28,7 +28,7 @@ openfermion = pytest.importorskip("openfermion") openfermionpyscf = pytest.importorskip("openfermionpyscf") - +pyscf = pytest.importorskip("pyscf") pauli_ops_and_prod = (qml.PauliX, qml.PauliY, qml.PauliZ, qml.Identity, qml.ops.Prod) pauli_ops_and_tensor = (qml.PauliX, qml.PauliY, qml.PauliZ, qml.Identity, qml.operation.Tensor) @@ -893,11 +893,9 @@ def test_excited_configurations(electrons, orbitals, excitation, states_ref, sig def test_ucisd_state(molecule, basis, tol, wf_ref): r"""Test the UCISD wavefunction constructor.""" - from pyscf import gto, scf, ci - - mol = gto.M(atom=molecule, basis=basis) - myhf = scf.UHF(mol).run() - myci = ci.UCISD(myhf).run() + mol = pyscf.gto.M(atom=molecule, basis=basis) + myhf = pyscf.scf.UHF(mol).run() + myci = pyscf.ci.UCISD(myhf).run() wf_cisd = qchem.convert._ucisd_state(myci, tol=tol) From 99b0fee2f2e4a7376efd0ce1240d6ae3e295785d Mon Sep 17 00:00:00 2001 From: soranjh Date: Wed, 2 Aug 2023 16:02:50 -0400 Subject: [PATCH 23/60] add test for statevector --- tests/qchem/of_tests/test_convert.py | 39 +++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index d33f9a9bb8a..39a131214ac 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -891,7 +891,7 @@ def test_excited_configurations(electrons, orbitals, excitation, states_ref, sig ], ) def test_ucisd_state(molecule, basis, tol, wf_ref): - r"""Test the UCISD wavefunction constructor.""" + r"""Test that _ucisd_state returns the correct wavefunction.""" mol = pyscf.gto.M(atom=molecule, basis=basis) myhf = pyscf.scf.UHF(mol).run() @@ -901,3 +901,40 @@ def test_ucisd_state(molecule, basis, tol, wf_ref): assert wf_cisd.keys() == wf_ref.keys() assert np.allclose(abs(np.array(list(wf_cisd.values()))), abs(np.array(list(wf_ref.values())))) + + +@pytest.mark.parametrize( + ("wf_dict", "n_orbitals", "statevector"), + [ + ( # -0.99 |1100> + 0.11 |0011> + {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159}, + 2, + np.array( + [ + 0.0, + 0.0, + 0.0, + 0.1066467, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + -0.99429698, + 0.0, + 0.0, + 0.0, + ] + ), + ), + ], +) +def test_wfdict_to_statevector(wf_dict, n_orbitals, statevector): + r"""Test that wfdict_to_statevector returns the correct statevector.""" + + wf_comp = qchem.convert.wfdict_to_statevector(wf_dict, n_orbitals) + + assert np.allclose(wf_comp, statevector) From 41238e618fe1730ca542e51968f2fde835bf4e10 Mon Sep 17 00:00:00 2001 From: soranjh Date: Wed, 2 Aug 2023 16:14:58 -0400 Subject: [PATCH 24/60] add test for cisd --- tests/qchem/of_tests/test_convert.py | 59 ++++++++++++++++++++++------ 1 file changed, 46 insertions(+), 13 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 39a131214ac..6d79b9b1916 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -862,13 +862,6 @@ def test_excited_configurations(electrons, orbitals, excitation, states_ref, sig @pytest.mark.parametrize( ("molecule", "basis", "tol", "wf_ref"), [ - # reference data computed with pyscf -- identify the FCI matrix entries - # with the correct Fock occupation vector entries - # mycasci = mcscf.UCASCI(myhf, norbs, (nelec_a, nelec_b)).run() - # sparse_cascimatr = coo_matrix(mycasci.ci, shape=np.shape(cascivec), dtype=float) - # row, col, dat = sparse_cascimatr.row, sparse_cascimatr.col, sparse_cascimatr.data - # strs_row, strs_col = addrs2str(norbs, nelec_a, row), addrs2str(norbs, nelec_b, col) - # wf_ref = dict(zip(list(zip(strs_row, strs_col)), dat)) ( [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], "sto6g", @@ -904,10 +897,10 @@ def test_ucisd_state(molecule, basis, tol, wf_ref): @pytest.mark.parametrize( - ("wf_dict", "n_orbitals", "statevector"), + ("wf_dict", "n_orbitals", "wf_ref"), [ ( # -0.99 |1100> + 0.11 |0011> - {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159}, + {(1, 1): 0.9942969785398778, (2, 2): 0.10664669927602159}, 2, np.array( [ @@ -923,7 +916,7 @@ def test_ucisd_state(molecule, basis, tol, wf_ref): 0.0, 0.0, 0.0, - -0.99429698, + 0.99429698, 0.0, 0.0, 0.0, @@ -932,9 +925,49 @@ def test_ucisd_state(molecule, basis, tol, wf_ref): ), ], ) -def test_wfdict_to_statevector(wf_dict, n_orbitals, statevector): +def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): r"""Test that wfdict_to_statevector returns the correct statevector.""" - wf_comp = qchem.convert.wfdict_to_statevector(wf_dict, n_orbitals) + assert np.allclose(wf_comp, wf_ref) + + +@pytest.mark.parametrize( + ("molecule", "basis", "hftype", "wf_ref"), + [ + ( + [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], + "sto6g", + "uhf", + np.array( + [ + 0.0, + 0.0, + 0.0, + 0.1066467, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.99429698, + 0.0, + 0.0, + 0.0, + ] + ), + ), + ], +) +def test_cisd_state(molecule, basis, hftype, wf_ref): + r"""Test that cisd_state returns the correct wavefunction.""" + + mol = pyscf.gto.M(atom=molecule, basis=basis) + myhf = pyscf.scf.UHF(mol).run() + myci = pyscf.ci.UCISD(myhf).run() + + wf_comp = qchem.convert.cisd_state(myci, hftype) - assert np.allclose(wf_comp, statevector) + assert np.allclose(abs(wf_comp), wf_ref) From 7e9e1c1df0312dd31248fa4abd397dadd3f96b9d Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Wed, 2 Aug 2023 18:05:10 -0400 Subject: [PATCH 25/60] added the wavefunction constructor for restricted cisd, matched ucisd syntax --- pennylane/qchem/convert.py | 68 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 66 insertions(+), 2 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 654b3516b68..111afa7ba76 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -434,7 +434,7 @@ def _excited_configurations(electrons, orbitals, excitation): >>> electrons = 2 >>> orbitals = 5 >>> excitation = 2 - >>> _excitated_states(electrons, orbitals, excitation) + >>> _excited_configurations(electrons, orbitals, excitation) (array([28, 26, 25]), array([ 1, -1, 1])) """ hf_state = qml.qchem.hf_state(electrons, orbitals) @@ -462,6 +462,70 @@ def _excited_configurations(electrons, orbitals, excitation): return states_int, signs +def _rcisd_state(cisd_solver, state=0, tol=1e-15): + r""" + + 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(). + + 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(). + + **Example** + + """ + mol = cisd_solver.mol + cisdvec = cisd_solver.ci + + norb = mol.nao + nelec = mol.nelectron + nocc, nvir = nelec // 2, norb - nelec // 2 + + c0, c1, c2 = cisdvec[0], cisdvec[1:nocc*nvir+1], cisdvec[nocc*nvir+1:].reshape(nocc,nocc,nvir,nvir) + + # numbers representing the Hartree-Fock vector, e.g., bin(ref_a)[::-1] = 1111...10...0 + ref_a = int(2**nelec - 1) + ref_b = ref_a + + dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [c0])) + + # alpha -> alpha excitations + c1a_configs, c1a_signs = _excited_configurations(nocc, norb, 1) + dict_fcimatr.update(dict(zip(list(zip(c1a_configs, [ref_b] * len(c1a_configs))), c1 * c1a_signs))) + # beta -> beta excitations + dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(c1a_configs), c1a_configs)), c1 * c1a_signs))) + + # check that double excitations within one spin sector are possible + if nocc > 1 and nvir > 1: + # get rid of excitations from same orbitals + c2_tr = c2 - c2.transpose(1, 0, 2, 3) + # select only unqiue excitations, via lower triangle of matrix + ooidx, vvidx = np.tril_indices(nocc, -1), np.tril_indices(nvir, -1) + c2aa = c2_tr[ooidx][:, vvidx[0], vvidx[1]] + + # alpha, alpha -> alpha, alpha excitations + c2aa_configs, c2aa_signs = _excited_configurations(nocc, norb, 2) + dict_fcimatr.update(dict(zip(list(zip(c2aa_configs, [ref_b] * len(c2aa_configs))), c2aa.ravel() * c2aa_signs))) + # beta, beta -> beta, beta excitations + dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(c2aa_configs), c2aa_configs)), c2aa.ravel() * c2aa_signs))) + + # alpha, beta -> alpha, beta excitations + # generate all possible pairwise combinations of _single_ excitations of alpha and beta sectors + rowvals, colvals = np.array( list(product(c1a_configs, c1a_configs)), dtype=int ).T.numpy() + c2ab = c2.transpose(0, 2, 1, 3).reshape(nocc*nvir, -1) + dict_fcimatr.update( + dict(zip(list(zip(rowvals, colvals)), c2ab.ravel() * np.kron(c1a_signs, c1a_signs).numpy())) + ) + + # filter based on tolerance cutoff + dict_fcimatr = {key: value for key, value in dict_fcimatr.items() if abs(value) > tol} + + return dict_fcimatr + def _ucisd_state(cisd_solver, tol=1e-15): r"""Construct a wavefunction from PySCF's `UCISD` Solver object. @@ -530,7 +594,7 @@ def _ucisd_state(cisd_solver, tol=1e-15): # beta -> beta excitations c1b_configs, c1b_signs = _excited_configurations(nelec_b, norb, 1) - dict_fcimatr.update(dict(zip(list(zip(c1b_configs, [ref_a] * size_b)), c1b * c1b_signs))) + dict_fcimatr.update(dict(zip(list(zip([ref_a] * size_b), c1b_configs), c1b * c1b_signs))) # alpha, alpha -> alpha, alpha excitations c2aa_configs, c2aa_signs = _excited_configurations(nelec_a, norb, 2) From 0b8566bfb7a5fceecae221c74a8f906da2fc64ca Mon Sep 17 00:00:00 2001 From: soranjh Date: Wed, 2 Aug 2023 18:06:30 -0400 Subject: [PATCH 26/60] add error test and fix bug --- pennylane/qchem/convert.py | 6 ++++-- tests/qchem/of_tests/test_convert.py | 6 ++++++ 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 654b3516b68..c676f399999 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -530,7 +530,7 @@ def _ucisd_state(cisd_solver, tol=1e-15): # beta -> beta excitations c1b_configs, c1b_signs = _excited_configurations(nelec_b, norb, 1) - dict_fcimatr.update(dict(zip(list(zip(c1b_configs, [ref_a] * size_b)), c1b * c1b_signs))) + dict_fcimatr.update(dict(zip(list(zip([ref_a] * size_b, c1b_configs)), c1b * c1b_signs))) # alpha, alpha -> alpha, alpha excitations c2aa_configs, c2aa_signs = _excited_configurations(nelec_a, norb, 2) @@ -585,7 +585,9 @@ def cisd_state(cisd_solver, hftype, tol=1e-15): if hftype == "uhf": wf_dict = _ucisd_state(cisd_solver, tol=tol) else: - raise ValueError("Only restricted, 'rhf', and unrestricted, 'uhf', are supported.") + raise ValueError( + "The supported hftype options are 'rhf' for restricted," " and 'uhf' for unrestricted." + ) wf = wfdict_to_statevector(wf_dict, cisd_solver.mol.nao) return wf diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 6d79b9b1916..6b128561b86 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -971,3 +971,9 @@ def test_cisd_state(molecule, basis, hftype, wf_ref): wf_comp = qchem.convert.cisd_state(myci, hftype) assert np.allclose(abs(wf_comp), wf_ref) + + +def test_cisd_state_error(): + r"""Test that an error is raised if a wrong/not-supported hftype symbol is entered.""" + with pytest.raises(ValueError, match="The supported hftype options are"): + _ = qchem.convert.cisd_state("myci", "hf") From 541318a9bcd4e3098247fdb4c7277b07b273ccf4 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Wed, 2 Aug 2023 18:08:37 -0400 Subject: [PATCH 27/60] minor comments for clarification added --- pennylane/qchem/convert.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 111afa7ba76..432fa4493ec 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -499,26 +499,26 @@ def _rcisd_state(cisd_solver, state=0, tol=1e-15): # beta -> beta excitations dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(c1a_configs), c1a_configs)), c1 * c1a_signs))) - # check that double excitations within one spin sector are possible + # check if double excitations within one spin sector (aa->aa and bb->bb) are possible if nocc > 1 and nvir > 1: - # get rid of excitations from same orbitals + # get rid of excitations from identical orbitals, double-count the allowed ones c2_tr = c2 - c2.transpose(1, 0, 2, 3) - # select only unqiue excitations, via lower triangle of matrix + # select only unique excitations, via lower triangle of matrix (already double-counted) ooidx, vvidx = np.tril_indices(nocc, -1), np.tril_indices(nvir, -1) - c2aa = c2_tr[ooidx][:, vvidx[0], vvidx[1]] + c2aa = c2_tr[ooidx][:, vvidx[0], vvidx[1]].ravel() # alpha, alpha -> alpha, alpha excitations c2aa_configs, c2aa_signs = _excited_configurations(nocc, norb, 2) - dict_fcimatr.update(dict(zip(list(zip(c2aa_configs, [ref_b] * len(c2aa_configs))), c2aa.ravel() * c2aa_signs))) + dict_fcimatr.update(dict(zip(list(zip(c2aa_configs, [ref_b] * len(c2aa_configs))), c2aa * c2aa_signs))) # beta, beta -> beta, beta excitations - dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(c2aa_configs), c2aa_configs)), c2aa.ravel() * c2aa_signs))) + dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(c2aa_configs), c2aa_configs)), c2aa * c2aa_signs))) # alpha, beta -> alpha, beta excitations # generate all possible pairwise combinations of _single_ excitations of alpha and beta sectors rowvals, colvals = np.array( list(product(c1a_configs, c1a_configs)), dtype=int ).T.numpy() - c2ab = c2.transpose(0, 2, 1, 3).reshape(nocc*nvir, -1) + c2ab = (c2.transpose(0, 2, 1, 3).reshape(nocc*nvir, -1)).ravel() dict_fcimatr.update( - dict(zip(list(zip(rowvals, colvals)), c2ab.ravel() * np.kron(c1a_signs, c1a_signs).numpy())) + dict(zip(list(zip(rowvals, colvals)), c2ab * np.kron(c1a_signs, c1a_signs).numpy())) ) # filter based on tolerance cutoff From ae363b3c9b40750b85cebda4456962806e207cf2 Mon Sep 17 00:00:00 2001 From: soranjh Date: Wed, 2 Aug 2023 18:12:29 -0400 Subject: [PATCH 28/60] add correct input args to test --- tests/qchem/of_tests/test_convert.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 6b128561b86..cdc0eaec40d 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -975,5 +975,9 @@ def test_cisd_state(molecule, basis, hftype, wf_ref): def test_cisd_state_error(): r"""Test that an error is raised if a wrong/not-supported hftype symbol is entered.""" + + myci = pyscf.ci.UCISD + hftype = "wrongtype" + with pytest.raises(ValueError, match="The supported hftype options are"): - _ = qchem.convert.cisd_state("myci", "hf") + _ = qchem.convert.cisd_state(myci, hftype) From 0ddab7bd5281526ae1315dcbee5eaee99ff3c836 Mon Sep 17 00:00:00 2001 From: soranjh Date: Wed, 2 Aug 2023 18:16:22 -0400 Subject: [PATCH 29/60] fix typo --- pennylane/qchem/convert.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index c676f399999..3e17fb710ea 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -586,7 +586,7 @@ def cisd_state(cisd_solver, hftype, tol=1e-15): wf_dict = _ucisd_state(cisd_solver, tol=tol) else: raise ValueError( - "The supported hftype options are 'rhf' for restricted," " and 'uhf' for unrestricted." + "The supported hftype options are 'rhf' for restricted and 'uhf' for unrestricted." ) wf = wfdict_to_statevector(wf_dict, cisd_solver.mol.nao) From c80f7e5f776b67df9297ac3e9d1bfd78c2de17c3 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Wed, 2 Aug 2023 18:58:12 -0400 Subject: [PATCH 30/60] added tests for rcisd and improved ucisd tests --- tests/qchem/of_tests/test_convert.py | 114 ++++++++++++++++++++++----- 1 file changed, 95 insertions(+), 19 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 6d79b9b1916..e1f6cb57e14 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -860,33 +860,37 @@ def test_excited_configurations(electrons, orbitals, excitation, states_ref, sig @pytest.mark.parametrize( - ("molecule", "basis", "tol", "wf_ref"), + ("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.10664669927602159}, + {(1, 1): 0.9942969785398776, (2, 2): -0.10664669927602176}, ), ( [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], "cc-pvdz", + "d2h", 4e-2, { - (1, 1): 0.9919704801055201, - (2, 2): 0.04853035090478132, - (2, 8): -0.04452332640264183, - (4, 4): -0.05003594588609278, - (8, 2): 0.044523334555907366, - (8, 8): -0.05226230845395026, + (1, 1): 0.9919704795977625, + (2, 2): -0.048530356564387034, + (2, 8): 0.0445233308500785, + (4, 4): -0.05003594568491194, + (8, 2): 0.04452333085007853, + (8, 8): -0.05226230322043741, + (16, 16): -0.0404759737476627, + (32, 32): -0.0404759737476627 }, ), ], ) -def test_ucisd_state(molecule, basis, tol, wf_ref): +def test_ucisd_state(molecule, basis, symm, tol, wf_ref): r"""Test that _ucisd_state returns the correct wavefunction.""" - mol = pyscf.gto.M(atom=molecule, basis=basis) + mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) myhf = pyscf.scf.UHF(mol).run() myci = pyscf.ci.UCISD(myhf).run() @@ -895,12 +899,52 @@ def test_ucisd_state(molecule, basis, tol, wf_ref): assert wf_cisd.keys() == wf_ref.keys() assert np.allclose(abs(np.array(list(wf_cisd.values()))), abs(np.array(list(wf_ref.values())))) +@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 + }, + ), + ], +) +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( ("wf_dict", "n_orbitals", "wf_ref"), [ ( # -0.99 |1100> + 0.11 |0011> - {(1, 1): 0.9942969785398778, (2, 2): 0.10664669927602159}, + {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602179}, 2, np.array( [ @@ -916,7 +960,7 @@ def test_ucisd_state(molecule, basis, tol, wf_ref): 0.0, 0.0, 0.0, - 0.99429698, + -0.99429698, 0.0, 0.0, 0.0, @@ -932,18 +976,19 @@ def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): @pytest.mark.parametrize( - ("molecule", "basis", "hftype", "wf_ref"), + ("molecule", "basis", "symm", "hftype", "wf_ref"), [ ( [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], "sto6g", + "d2h", "uhf", np.array( [ 0.0, 0.0, 0.0, - 0.1066467, + -0.1066467, 0.0, 0.0, 0.0, @@ -959,15 +1004,46 @@ def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): ] ), ), + ( + [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], + "sto6g", + "d2h", + "rhf", + np.array( + [ + 0.0, + 0.0, + 0.0, + 0.10664669927602179, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + -0.9942969785398778, + 0.0, + 0.0, + 0.0, + ] + ), + ), ], ) -def test_cisd_state(molecule, basis, hftype, wf_ref): +def test_cisd_state(molecule, basis, symm, hftype, wf_ref): r"""Test that cisd_state returns the correct wavefunction.""" - mol = pyscf.gto.M(atom=molecule, basis=basis) - myhf = pyscf.scf.UHF(mol).run() - myci = pyscf.ci.UCISD(myhf).run() + mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) + if hftype == "rhf": + myhf = pyscf.scf.RHF(mol).run() + myci = pyscf.ci.CISD(myhf).run() + elif hftype =="uhf": + myhf = pyscf.scf.UHF(mol).run() + myci = pyscf.ci.UCISD(myhf).run() wf_comp = qchem.convert.cisd_state(myci, hftype) - assert np.allclose(abs(wf_comp), wf_ref) + # overall sign could be +/-1 different + assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) From 84067eba501e1da7abbce109b0d4e6a7ff10e80f Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Wed, 2 Aug 2023 18:59:13 -0400 Subject: [PATCH 31/60] small fixes to _rcisd func --- pennylane/qchem/convert.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 432fa4493ec..321276ba277 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -462,7 +462,7 @@ def _excited_configurations(electrons, orbitals, excitation): return states_int, signs -def _rcisd_state(cisd_solver, state=0, tol=1e-15): +def _rcisd_state(cisd_solver, tol=1e-15): r""" Args: @@ -488,7 +488,7 @@ def _rcisd_state(cisd_solver, state=0, tol=1e-15): c0, c1, c2 = cisdvec[0], cisdvec[1:nocc*nvir+1], cisdvec[nocc*nvir+1:].reshape(nocc,nocc,nvir,nvir) # numbers representing the Hartree-Fock vector, e.g., bin(ref_a)[::-1] = 1111...10...0 - ref_a = int(2**nelec - 1) + ref_a = int(2**nocc - 1) ref_b = ref_a dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [c0])) @@ -594,7 +594,7 @@ def _ucisd_state(cisd_solver, tol=1e-15): # beta -> beta excitations c1b_configs, c1b_signs = _excited_configurations(nelec_b, norb, 1) - dict_fcimatr.update(dict(zip(list(zip([ref_a] * size_b), c1b_configs), c1b * c1b_signs))) + dict_fcimatr.update(dict(zip(list(zip([ref_a] * size_b, c1b_configs)), c1b * c1b_signs))) # alpha, alpha -> alpha, alpha excitations c2aa_configs, c2aa_signs = _excited_configurations(nelec_a, norb, 2) @@ -646,7 +646,9 @@ def cisd_state(cisd_solver, hftype, tol=1e-15): {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} """ - if hftype == "uhf": + if hftype == "rhf": + wf_dict = _rcisd_state(cisd_solver, tol=tol) + elif hftype == "uhf": wf_dict = _ucisd_state(cisd_solver, tol=tol) else: raise ValueError("Only restricted, 'rhf', and unrestricted, 'uhf', are supported.") From 0ceffc85a1d56a3b96ddfd518f21c8e285886015 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Wed, 2 Aug 2023 19:04:14 -0400 Subject: [PATCH 32/60] updating tests with symm argument --- tests/qchem/of_tests/test_convert.py | 135 +++++++++++++++++++++------ 1 file changed, 106 insertions(+), 29 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index c7c8c1735ff..c30276950b9 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -28,7 +28,7 @@ openfermion = pytest.importorskip("openfermion") openfermionpyscf = pytest.importorskip("openfermionpyscf") - +pyscf = pytest.importorskip("pyscf") pauli_ops_and_prod = (qml.PauliX, qml.PauliY, qml.PauliZ, qml.Identity, qml.ops.Prod) pauli_ops_and_tensor = (qml.PauliX, qml.PauliY, qml.PauliZ, qml.Identity, qml.operation.Tensor) @@ -859,43 +859,120 @@ def test_excitated_states(electrons, orbitals, excitation, states_ref, signs_ref assert np.allclose(signs, signs_ref) @pytest.mark.parametrize( - ("atom", "basis", "hftype", "state", "tol", "wf_ref"), + ("molecule", "basis", "symm", "tol", "wf_ref"), [ - # reference data computed with pyscf -- identify the FCI matrix entries - # with the correct Fock occupation vector entries - # mycasci = mcscf.UCASCI(myhf, norbs, (nelec_a, nelec_b)).run() - # sparse_cascimatr = coo_matrix(mycasci.ci, shape=np.shape(cascivec), dtype=float) - # row, col, dat = sparse_cascimatr.row, sparse_cascimatr.col, sparse_cascimatr.data - # strs_row, strs_col = addrs2str(norbs, nelec_a, row), addrs2str(norbs, nelec_b, col) - # wf_ref = dict(zip(list(zip(strs_row, strs_col)), dat)) ( - [['H', (0, 0, 0)], ['H', (0,0,0.71)]], - 'sto6g', - 'uhf', - 0, + [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], + "sto6g", + "d2h", 1e-1, - {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} + {(1, 1): 0.9942969785398776, (2, 2): -0.10664669927602176}, ), ( - [['H', (0, 0, 0)], ['H', (0, 0, 0.71)]], - 'cc-pvdz', - 'uhf', - 0, + [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], + "cc-pvdz", + "d2h", 4e-2, - {(1, 1): 0.9919704801055201, (2, 2): 0.04853035090478132, (2, 8): -0.04452332640264183, \ - (4, 4): -0.05003594588609278, (8, 2): 0.044523334555907366, (8, 8): -0.05226230845395026}), + { + (1, 1): 0.9919704795977625, + (2, 2): -0.048530356564387034, + (2, 8): 0.0445233308500785, + (4, 4): -0.05003594568491194, + (8, 2): 0.04452333085007853, + (8, 8): -0.05226230322043741, + (16, 16): -0.0404759737476627, + (32, 32): -0.0404759737476627 + }, + ), + ], +) +def test_ucisd_state(molecule, basis, symm, tol, wf_ref): + r"""Test that _ucisd_state returns the correct wavefunction.""" + + mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) + myhf = pyscf.scf.UHF(mol).run() + myci = pyscf.ci.UCISD(myhf).run() + + wf_cisd = qchem.convert._ucisd_state(myci, tol=tol) + + assert wf_cisd.keys() == wf_ref.keys() + assert np.allclose(abs(np.array(list(wf_cisd.values()))), abs(np.array(list(wf_ref.values())))) + + +@pytest.mark.parametrize( + ("wf_dict", "n_orbitals", "wf_ref"), + [ + ( # -0.99 |1100> + 0.11 |0011> + {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602179}, + 2, + np.array( + [ + 0.0, + 0.0, + 0.0, + 0.1066467, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + -0.99429698, + 0.0, + 0.0, + 0.0, + ] + ), + ), ], ) -def test_private_ucisd_state(atom, basis, hftype, state, tol, wf_ref): - r"""Test the UCISD wavefunction constructor.""" +def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): + r"""Test that wfdict_to_statevector returns the correct statevector.""" + wf_comp = qchem.convert.wfdict_to_statevector(wf_dict, n_orbitals) + assert np.allclose(wf_comp, wf_ref) + - from pyscf import gto, scf, ci +@pytest.mark.parametrize( + ("molecule", "basis", "symm", "hftype", "wf_ref"), + [ + ( + [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], + "sto6g", + "d2h", + "uhf", + np.array( + [ + 0.0, + 0.0, + 0.0, + -0.1066467, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.99429698, + 0.0, + 0.0, + 0.0, + ] + ), + ), + ], +) +def test_cisd_state(molecule, basis, symm, hftype, wf_ref): + r"""Test that cisd_state returns the correct wavefunction.""" - mol = gto.M(atom=atom, basis=basis) - myhf = scf.UHF(mol).run() - myci = ci.UCISD(myhf).run(state=state) + mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) + myhf = pyscf.scf.UHF(mol).run() + myci = pyscf.ci.UCISD(myhf).run() - wf_cisd = qchem.convert.cisd_state(myci, hftype=hftype, state=state, tol=tol) + wf_comp = qchem.convert.cisd_state(myci, hftype) - assert { key: round(val, 4) for (key, val) in wf_cisd.items() } == \ - { key: round(val, 4) for (key, val) in wf_ref.items() } \ No newline at end of file + # overall sign could be +/-1 different + assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) \ No newline at end of file From eee04106820200757e0fff80611c87f587ec50b1 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Wed, 2 Aug 2023 19:08:37 -0400 Subject: [PATCH 33/60] reformatted with black --- tests/qchem/of_tests/test_convert.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index a72a6a9f239..844a75b33a7 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -875,14 +875,14 @@ def test_excited_configurations(electrons, orbitals, excitation, states_ref, sig "d2h", 4e-2, { - (1, 1): 0.9919704795977625, - (2, 2): -0.048530356564387034, - (2, 8): 0.0445233308500785, - (4, 4): -0.05003594568491194, - (8, 2): 0.04452333085007853, - (8, 8): -0.05226230322043741, - (16, 16): -0.0404759737476627, - (32, 32): -0.0404759737476627 + (1, 1): 0.9919704795977625, + (2, 2): -0.048530356564387034, + (2, 8): 0.0445233308500785, + (4, 4): -0.05003594568491194, + (8, 2): 0.04452333085007853, + (8, 8): -0.05226230322043741, + (16, 16): -0.0404759737476627, + (32, 32): -0.0404759737476627, }, ), ], @@ -978,6 +978,7 @@ def test_cisd_state(molecule, basis, symm, hftype, wf_ref): # overall sign could be +/-1 different assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) + def test_cisd_state_error(): r"""Test that an error is raised if a wrong/not-supported hftype symbol is entered.""" From 6afd6be65557124c6c1dd3524fcca60b0b8efd01 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Wed, 2 Aug 2023 19:11:36 -0400 Subject: [PATCH 34/60] fixed black formatting --- pennylane/qchem/convert.py | 29 ++++++++++++++------ tests/qchem/of_tests/test_convert.py | 41 ++++++++++++++-------------- 2 files changed, 42 insertions(+), 28 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 321276ba277..b3bfbf1147f 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -462,6 +462,7 @@ def _excited_configurations(electrons, orbitals, excitation): return states_int, signs + def _rcisd_state(cisd_solver, tol=1e-15): r""" @@ -484,8 +485,12 @@ def _rcisd_state(cisd_solver, tol=1e-15): norb = mol.nao nelec = mol.nelectron nocc, nvir = nelec // 2, norb - nelec // 2 - - c0, c1, c2 = cisdvec[0], cisdvec[1:nocc*nvir+1], cisdvec[nocc*nvir+1:].reshape(nocc,nocc,nvir,nvir) + + c0, c1, c2 = ( + cisdvec[0], + cisdvec[1 : nocc * nvir + 1], + cisdvec[nocc * nvir + 1 :].reshape(nocc, nocc, nvir, nvir), + ) # numbers representing the Hartree-Fock vector, e.g., bin(ref_a)[::-1] = 1111...10...0 ref_a = int(2**nocc - 1) @@ -495,9 +500,13 @@ def _rcisd_state(cisd_solver, tol=1e-15): # alpha -> alpha excitations c1a_configs, c1a_signs = _excited_configurations(nocc, norb, 1) - dict_fcimatr.update(dict(zip(list(zip(c1a_configs, [ref_b] * len(c1a_configs))), c1 * c1a_signs))) + dict_fcimatr.update( + dict(zip(list(zip(c1a_configs, [ref_b] * len(c1a_configs))), c1 * c1a_signs)) + ) # beta -> beta excitations - dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(c1a_configs), c1a_configs)), c1 * c1a_signs))) + dict_fcimatr.update( + dict(zip(list(zip([ref_a] * len(c1a_configs), c1a_configs)), c1 * c1a_signs)) + ) # check if double excitations within one spin sector (aa->aa and bb->bb) are possible if nocc > 1 and nvir > 1: @@ -509,14 +518,18 @@ def _rcisd_state(cisd_solver, tol=1e-15): # alpha, alpha -> alpha, alpha excitations c2aa_configs, c2aa_signs = _excited_configurations(nocc, norb, 2) - dict_fcimatr.update(dict(zip(list(zip(c2aa_configs, [ref_b] * len(c2aa_configs))), c2aa * c2aa_signs))) + dict_fcimatr.update( + dict(zip(list(zip(c2aa_configs, [ref_b] * len(c2aa_configs))), c2aa * c2aa_signs)) + ) # beta, beta -> beta, beta excitations - dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(c2aa_configs), c2aa_configs)), c2aa * c2aa_signs))) + dict_fcimatr.update( + dict(zip(list(zip([ref_a] * len(c2aa_configs), c2aa_configs)), c2aa * c2aa_signs)) + ) # alpha, beta -> alpha, beta excitations # generate all possible pairwise combinations of _single_ excitations of alpha and beta sectors - rowvals, colvals = np.array( list(product(c1a_configs, c1a_configs)), dtype=int ).T.numpy() - c2ab = (c2.transpose(0, 2, 1, 3).reshape(nocc*nvir, -1)).ravel() + rowvals, colvals = np.array(list(product(c1a_configs, c1a_configs)), dtype=int).T.numpy() + c2ab = (c2.transpose(0, 2, 1, 3).reshape(nocc * nvir, -1)).ravel() dict_fcimatr.update( dict(zip(list(zip(rowvals, colvals)), c2ab * np.kron(c1a_signs, c1a_signs).numpy())) ) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index e1f6cb57e14..89dcd14a026 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -875,14 +875,14 @@ def test_excited_configurations(electrons, orbitals, excitation, states_ref, sig "d2h", 4e-2, { - (1, 1): 0.9919704795977625, - (2, 2): -0.048530356564387034, - (2, 8): 0.0445233308500785, - (4, 4): -0.05003594568491194, - (8, 2): 0.04452333085007853, - (8, 8): -0.05226230322043741, - (16, 16): -0.0404759737476627, - (32, 32): -0.0404759737476627 + (1, 1): 0.9919704795977625, + (2, 2): -0.048530356564387034, + (2, 8): 0.0445233308500785, + (4, 4): -0.05003594568491194, + (8, 2): 0.04452333085007853, + (8, 8): -0.05226230322043741, + (16, 16): -0.0404759737476627, + (32, 32): -0.0404759737476627, }, ), ], @@ -899,30 +899,31 @@ def test_ucisd_state(molecule, basis, symm, tol, wf_ref): assert wf_cisd.keys() == wf_ref.keys() assert np.allclose(abs(np.array(list(wf_cisd.values()))), abs(np.array(list(wf_ref.values())))) + @pytest.mark.parametrize( ("molecule", "basis", "symm", "tol", "wf_ref"), [ ( [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], "sto6g", - 'd2h', + "d2h", 1e-1, {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602179}, ), ( [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], "cc-pvdz", - 'd2h', + "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 + (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, }, ), ], @@ -937,7 +938,7 @@ def test_rcisd_state(molecule, basis, symm, tol, wf_ref): 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())) ) + assert np.allclose(np.array(list(wf_cisd.values())), np.array(list(wf_ref.values()))) @pytest.mark.parametrize( @@ -1039,7 +1040,7 @@ def test_cisd_state(molecule, basis, symm, hftype, wf_ref): if hftype == "rhf": myhf = pyscf.scf.RHF(mol).run() myci = pyscf.ci.CISD(myhf).run() - elif hftype =="uhf": + elif hftype == "uhf": myhf = pyscf.scf.UHF(mol).run() myci = pyscf.ci.UCISD(myhf).run() From 3f17546b5cdcbed6641e98e274b4955a45f05f30 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Thu, 3 Aug 2023 08:33:56 -0400 Subject: [PATCH 35/60] added Be atom test for rcisd aa->aa and bb->bb check --- tests/qchem/of_tests/test_convert.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 89dcd14a026..7640474055c 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -926,6 +926,23 @@ def test_ucisd_state(molecule, basis, symm, tol, wf_ref): (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): From 9981f155c914bed17254c1810d6b9adfaeebcb73 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Thu, 3 Aug 2023 08:37:20 -0400 Subject: [PATCH 36/60] black format fix --- tests/qchem/of_tests/test_convert.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 7640474055c..3e28c868d90 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -932,16 +932,17 @@ def test_ucisd_state(molecule, basis, symm, tol, wf_ref): "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}, + (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, + }, ), ], ) From 8320ceca1fee73e516e6eeb6e38041a4dbef11ec Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Thu, 3 Aug 2023 09:16:38 -0400 Subject: [PATCH 37/60] first version of uccsd converter --- pennylane/qchem/convert.py | 64 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index b3bfbf1147f..72392ffd1cd 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -670,6 +670,70 @@ def cisd_state(cisd_solver, hftype, tol=1e-15): return wf +def _uccsd_state(ccsd_solver, tol=1e-15): + + mol = ccsd_solver.mol + + norb = mol.nao + nelec = mol.nelectron + nelec_a = int((nelec + mol.spin) / 2) + nelec_b = int((nelec - mol.spin) / 2) + + nvir_a, nvir_b = norb - nelec_a, norb - nelec_b + + t1a, t1b = ccsd_solver.t1 + t2aa, t2ab, t2bb = ccsd_solver.t2 + # compute the disconnected part of double excitations, + + # numbers representing the Hartree-Fock vector, e.g., bin(ref_a)[::-1] = 1111...10...0 + ref_a = int(2**nelec_a - 1) + ref_b = int(2**nelec_b - 1) + + dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [1.])) + + # alpha -> alpha excitations + t1a_configs, t1a_signs = _excited_configurations(nelec_a, norb, 1) + dict_fcimatr.update(dict(zip(list(zip(t1a_configs, [ref_b] * len(t1a_configs))), t1a.ravel() * t1a_signs))) + + # beta -> beta excitations + t1b_configs, t1b_signs = _excited_configurations(nelec_b, norb, 1) + dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(t1b_configs), t1b_configs)), t1b.ravel() * t1b_signs))) + + # alpha, alpha -> alpha, alpha excitations + if nelec_a > 1 and nvir_a > 1: + t2aa_configs, t2aa_signs = _excited_configurations(nelec_a, norb, 2) + # select only unique excitations, via lower triangle of matrix + ooidx = np.tril_indices(nelec_a, -1) + vvidx = np.tril_indices(nvir_a, -1) + t2aa = t2aa[ooidx][:, vvidx[0], vvidx[1]] + #TODO: add the (1/2) T_1^2! + dict_fcimatr.update(dict(zip(list(zip(t2aa_configs, [ref_b] * len(t2aa_configs))), t2aa.ravel() * t2aa_signs))) + + if nelec_b > 1 and nvir_b > 1: + t2bb_configs, t2bb_signs = _excited_configurations(nelec_b, norb, 2) + # select only unique excitations, via lower triangle of matrix + ooidx = np.tril_indices(nelec_b, -1) + vvidx = np.tril_indices(nvir_b, -1) + t2bb = t2bb[ooidx][:, vvidx[0], vvidx[1]] + #TODO: add the (1/2) T_1^2! + dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(t2bb_configs), t2bb_configs)), t2bb.ravel() * t2bb_signs))) + + # alpha, beta -> alpha, beta excitations + rowvals, colvals = np.array(list(product(t1a_configs, t1b_configs)), dtype=int).T.numpy() + #TODO: add the (1/2) T_1^2! + dict_fcimatr.update( + dict(zip(list(zip(rowvals, colvals)), t2ab.ravel() * np.kron(t1a_signs, t1b_signs).numpy())) + ) + + # renormalize, to get the HF coefficient (CC wavefunction not normalized) + norm = np.sqrt(np.sum(dict_fcimatr.values()**2)) + dict_fcimatr = {key: value / norm for (key, value) in dict_fcimatr.items()} + + # filter based on tolerance cutoff + dict_fcimatr = {key: value for key, value in dict_fcimatr.items() if abs(value) > tol} + + return dict_fcimatr + def wfdict_to_statevector(wf_dict, norbs): r"""Convert a wavefunction in sparce dictionary format to a Pennylane's statevector. From 38011470ca1090490b837d9ac8d3ff0461f45040 Mon Sep 17 00:00:00 2001 From: soranjh Date: Thu, 3 Aug 2023 09:43:10 -0400 Subject: [PATCH 38/60] make code compact --- pennylane/qchem/convert.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 3e17fb710ea..67cb6aad944 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -594,7 +594,7 @@ def cisd_state(cisd_solver, hftype, tol=1e-15): def wfdict_to_statevector(wf_dict, norbs): - r"""Convert a wavefunction in sparce dictionary format to a Pennylane's statevector. + r"""Convert a wavefunction in sparce dictionary format to a PennyLane statevector. Args: wf_dict (wf_dict format): dictionary with keys (int_a,int_b) being integers whose binary representation @@ -611,12 +611,10 @@ def wfdict_to_statevector(wf_dict, norbs): statevector = np.zeros(2 ** (2 * norbs), dtype=complex) for (int_a, int_b), coeff in wf_dict.items(): - bin_a, bin_b = bin(int_a)[2:], bin(int_b)[2:] - bin_a, bin_b = bin_a[::-1], bin_b[::-1] - bin_tot = "".join(i + j for i, j in zip(bin_a, bin_b)) - bin_tot_full = bin_tot + "0" * (2 * norbs - len(bin_tot)) - idx = int(bin_tot_full, 2) - statevector[idx] += coeff + bin_a, bin_b = bin(int_a)[2:][::-1], bin(int_b)[2:][::-1] + bin_ab = "".join(i + j for i, j in zip(bin_a, bin_b)) + bin_ab_full = bin_ab + "0" * (2 * norbs - len(bin_ab)) + statevector[int(bin_ab_full, 2)] += coeff statevector = statevector / np.sqrt(np.sum(statevector**2)) From c8e69837f4ba4caa7b405f8dbd327d96bdb38c37 Mon Sep 17 00:00:00 2001 From: soranjh Date: Thu, 3 Aug 2023 10:33:55 -0400 Subject: [PATCH 39/60] rename cisd_state to import_state --- pennylane/qchem/convert.py | 20 +++++++++----------- tests/qchem/of_tests/test_convert.py | 22 +++++++++++----------- 2 files changed, 20 insertions(+), 22 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 67cb6aad944..208590c29e7 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -552,15 +552,12 @@ def _ucisd_state(cisd_solver, tol=1e-15): return dict_fcimatr -def cisd_state(cisd_solver, hftype, tol=1e-15): +def import_state(solver, method, tol=1e-15): r"""Construct a wavefunction from PySCF output. Args: - cisd_solver (PySCF CISD or UCISD Class instance): the class object representing the - CISD / UCISD calculation in PySCF - - hftype (str): keyword specifying whether restricted or unrestricted CISD was performed. - The options are 'rhf' for restricted, and 'uhf' for unrestricted. + solver: external wavefunction object that will be converted + method (str): keyword specifying the method of calculation tol (float): the tolerance to which the wavefunction is being built -- Slater determinants with coefficients smaller than this are discarded. Default @@ -582,18 +579,19 @@ def cisd_state(cisd_solver, hftype, tol=1e-15): {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} """ - if hftype == "uhf": - wf_dict = _ucisd_state(cisd_solver, tol=tol) + if method == "ucisd": + wf_dict = _ucisd_state(solver, tol=tol) else: raise ValueError( - "The supported hftype options are 'rhf' for restricted and 'uhf' for unrestricted." + "The supported method options are 'rcisd' for restricted and 'ucisd' for unrestricted" + " CISD calculations." ) - wf = wfdict_to_statevector(wf_dict, cisd_solver.mol.nao) + wf = _wfdict_to_statevector(wf_dict, solver.mol.nao) return wf -def wfdict_to_statevector(wf_dict, norbs): +def _wfdict_to_statevector(wf_dict, norbs): r"""Convert a wavefunction in sparce dictionary format to a PennyLane statevector. Args: diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 844a75b33a7..8b678c5f7f6 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -930,19 +930,19 @@ def test_ucisd_state(molecule, basis, symm, tol, wf_ref): ], ) def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): - r"""Test that wfdict_to_statevector returns the correct statevector.""" - wf_comp = qchem.convert.wfdict_to_statevector(wf_dict, n_orbitals) + r"""Test that _wfdict_to_statevector returns the correct statevector.""" + wf_comp = qchem.convert._wfdict_to_statevector(wf_dict, n_orbitals) assert np.allclose(wf_comp, wf_ref) @pytest.mark.parametrize( - ("molecule", "basis", "symm", "hftype", "wf_ref"), + ("molecule", "basis", "symm", "method", "wf_ref"), [ ( [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], "sto6g", "d2h", - "uhf", + "ucisd", np.array( [ 0.0, @@ -966,24 +966,24 @@ def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): ), ], ) -def test_cisd_state(molecule, basis, symm, hftype, wf_ref): +def test_import_state(molecule, basis, symm, method, wf_ref): r"""Test that cisd_state returns the correct wavefunction.""" mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) myhf = pyscf.scf.UHF(mol).run() myci = pyscf.ci.UCISD(myhf).run() - wf_comp = qchem.convert.cisd_state(myci, hftype) + wf_comp = qchem.convert.import_state(myci, method) # overall sign could be +/-1 different assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) -def test_cisd_state_error(): - r"""Test that an error is raised if a wrong/not-supported hftype symbol is entered.""" +def test_import_state_error(): + r"""Test that an error is raised if a wrong/not-supported method symbol is entered.""" myci = pyscf.ci.UCISD - hftype = "wrongtype" + method = "wrongmethod" - with pytest.raises(ValueError, match="The supported hftype options are"): - _ = qchem.convert.cisd_state(myci, hftype) + with pytest.raises(ValueError, match="The supported method options are"): + _ = qchem.convert.import_state(myci, method) From fafef6a9f52ecb5f42c45f3f5c021bf2d112216e Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Thu, 3 Aug 2023 11:44:26 -0400 Subject: [PATCH 40/60] polished uccsd converter and added test for it --- pennylane/qchem/convert.py | 97 ++++++++++++++++++++++++++-- tests/qchem/of_tests/test_convert.py | 39 +++++++++++ 2 files changed, 131 insertions(+), 5 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 72392ffd1cd..a272d7c9a70 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -671,6 +671,50 @@ def cisd_state(cisd_solver, hftype, tol=1e-15): def _uccsd_state(ccsd_solver, tol=1e-15): + r"""Construct a wavefunction from PySCF's `UCCSD` Solver object. + + The wavefunction is represented as a dictionary where the keys are tuples representing a + configuration, which corresponds to a Slater determinant, and the values are the CI coefficient + corresponding to that Slater determinant. Each dictionary key is a tuple of two integers. The + binary representation of these integers correspond to a specific configuration, or Slater + determinant. The first number represents the configuration of alpha electrons and the second + number represents the beta electrons. For instance, the Hartree-Fock state + :math:`|1 1 0 0 \rangle` will be represented by the binary string `0011`. This string can be + splited to `01` and `01` for the alpha and beta electrons. The integer corresponding to `01` is + `1`. Then the dictionary representation of a state with `0.99` contribution of the Hartree-Fock + 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 + + .. 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 \right) \right] \ket{\text{HF}} + + 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). + + Returns: + dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` + having binary represention 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, cc + >>> 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) + >>> print(wf_ccsd) + {(7, 7): 0.8886988662500356, (11, 11): -0.30584653966926967, + (19, 19): -0.3058465396692694, (35, 35): -0.14507582719761972} + """ mol = ccsd_solver.mol @@ -683,7 +727,11 @@ def _uccsd_state(ccsd_solver, tol=1e-15): t1a, t1b = ccsd_solver.t1 t2aa, t2ab, t2bb = ccsd_solver.t2 - # compute the disconnected part of double excitations, + # add in the disconnected part ( + 0.5 T_1^2) of double excitations + t2aa = t2aa - 0*0.5 * np.kron(t1a, t1a).reshape(nelec_a, nvir_a, nelec_a, nvir_a).transpose(0,2,1,3).numpy() + t2bb = t2bb - 0*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() + t2ab = t2ab.transpose(0,2,1,3) - 0*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 ref_a = int(2**nelec_a - 1) @@ -706,7 +754,6 @@ def _uccsd_state(ccsd_solver, tol=1e-15): ooidx = np.tril_indices(nelec_a, -1) vvidx = np.tril_indices(nvir_a, -1) t2aa = t2aa[ooidx][:, vvidx[0], vvidx[1]] - #TODO: add the (1/2) T_1^2! dict_fcimatr.update(dict(zip(list(zip(t2aa_configs, [ref_b] * len(t2aa_configs))), t2aa.ravel() * t2aa_signs))) if nelec_b > 1 and nvir_b > 1: @@ -715,18 +762,16 @@ def _uccsd_state(ccsd_solver, tol=1e-15): ooidx = np.tril_indices(nelec_b, -1) vvidx = np.tril_indices(nvir_b, -1) t2bb = t2bb[ooidx][:, vvidx[0], vvidx[1]] - #TODO: add the (1/2) T_1^2! dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(t2bb_configs), t2bb_configs)), t2bb.ravel() * t2bb_signs))) # alpha, beta -> alpha, beta excitations rowvals, colvals = np.array(list(product(t1a_configs, t1b_configs)), dtype=int).T.numpy() - #TODO: add the (1/2) T_1^2! dict_fcimatr.update( dict(zip(list(zip(rowvals, colvals)), t2ab.ravel() * np.kron(t1a_signs, t1b_signs).numpy())) ) # renormalize, to get the HF coefficient (CC wavefunction not normalized) - norm = np.sqrt(np.sum(dict_fcimatr.values()**2)) + norm = np.sqrt(np.sum(np.array(list(dict_fcimatr.values()))**2)).numpy() dict_fcimatr = {key: value / norm for (key, value) in dict_fcimatr.items()} # filter based on tolerance cutoff @@ -734,6 +779,48 @@ def _uccsd_state(ccsd_solver, tol=1e-15): return dict_fcimatr +def ccsd_state(ccsd_solver, hftype, tol=1e-15): + r"""Construct a wavefunction from PySCF CCSD output. + + Args: + ccsd_solver (PySCF CCSD or UCCSD Class instance): the class object representing the + CCSD / UCCSD calculation in PySCF + + hftype (str): keyword specifying whether restricted or unrestricted CCSD was performed. + The options are 'rhf' for restricted, and 'uhf' for unrestricted. + + 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). + + Returns: + dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` + having binary represention 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.UHF(mol).run() + >>> mycc = ci.UCCSD(myhf).run() + >>> wf_ccsd = ccsd_state(mycc, tol=1e-1) + >>> print(wf_cisd) + [ 0. 0. 0. -0.1066467 0. 0. 0. + 0. 0. 0. 0. 0. 0.99429698 + 0. 0. 0. ] + """ + + if hftype == "rhf": + wf_dict = _rccsd_state(ccsd_solver, tol=tol) + elif hftype == "uhf": + wf_dict = _uccsd_state(ccsd_solver, tol=tol) + else: + raise ValueError("Only restricted, 'rhf', and unrestricted, 'uhf', are supported.") + wf = wfdict_to_statevector(wf_dict, ccsd_solver.mol.nao) + + return wf + def wfdict_to_statevector(wf_dict, norbs): r"""Convert a wavefunction in sparce dictionary format to a Pennylane's statevector. diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 3e28c868d90..09eb5ea8a05 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -1066,3 +1066,42 @@ def test_cisd_state(molecule, basis, symm, hftype, wf_ref): # overall sign could be +/-1 different assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) + +@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.9942969785398776, + (2, 2): -0.10664669927602176 + } + ), + ( + [["Li", (0, 0, 0)], ["Li", (0, 0, 0.71)]], + "sto6g", + "d2h", + 1e-1, + { + (7, 7): 0.8886988662500364, + (11, 11): -0.3058465396692694, + (19, 19): -0.3058465396692695, + (35, 35): -0.14507582719761614 + }, + ), + ], +) +def test_uccsd_state(molecule, basis, symm, tol, wf_ref): + r"""Test that _uccsd_state returns the correct wavefunction.""" + + mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) + myhf = pyscf.scf.UHF(mol).run() + mycc = pyscf.cc.UCCSD(myhf).run() + + wf_ccsd = qchem.convert._uccsd_state(mycc, tol=tol) + + assert wf_ccsd.keys() == wf_ref.keys() + assert np.allclose(abs(np.array(list(wf_ccsd.values()))), abs(np.array(list(wf_ref.values())))) \ No newline at end of file From 6079d551c88c32d7e5a8df0c6cf9d39cb5eacd2a Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Thu, 3 Aug 2023 11:57:19 -0400 Subject: [PATCH 41/60] added the rccsd converter --- pennylane/qchem/convert.py | 129 +++++++++++++++++++++++++++++++++++-- 1 file changed, 123 insertions(+), 6 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index a272d7c9a70..c031e17de9f 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -669,6 +669,121 @@ def cisd_state(cisd_solver, hftype, tol=1e-15): return wf +def _rccsd_state(ccsd_solver, tol=1e-15): + r"""Construct a wavefunction from PySCF's `CCSD` Solver object. + + The wavefunction is represented as a dictionary where the keys are tuples representing a + configuration, which corresponds to a Slater determinant, and the values are the CI coefficient + corresponding to that Slater determinant. Each dictionary key is a tuple of two integers. The + binary representation of these integers correspond to a specific configuration, or Slater + determinant. The first number represents the configuration of alpha electrons and the second + number represents the beta electrons. For instance, the Hartree-Fock state + :math:`|1 1 0 0 \rangle` will be represented by the binary string `0011`. This string can be + splited to `01` and `01` for the alpha and beta electrons. The integer corresponding to `01` is + `1`. Then the dictionary representation of a state with `0.99` contribution of the Hartree-Fock + 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 + + .. 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. + + Args: + ccsd_solver (PySCF CCSD Class instance): the class object representing the CCSD 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). + + Returns: + dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` + having binary represention 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, cc + >>> 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) + >>> print(wf_ccsd) + {(7, 7): 0.8886969878256522, (11, 11): -0.30584590248164206, + (19, 19): -0.30584590248164145, (35, 35): -0.14507552651170982} + """ + + mol = ccsd_solver.mol + + norb = mol.nao + nelec = mol.nelectron + nelec_a = int((nelec + mol.spin) / 2) + nelec_b = int((nelec - mol.spin) / 2) + + nvir_a, nvir_b = norb - nelec_a, norb - nelec_b + + # build the full, unrestricted representation of the coupled cluster amplitudes + t1a = ccsd_solver.t1 + t1b = t1a + t2aa = ccsd_solver.t2 - ccsd_solver.t2.transpose(1,0,2,3) + t2ab = ccsd_solver.t2.transpose(0, 2, 1, 3) + t2bb = t2aa + + # add in the disconnected part ( + 0.5 T_1^2) of double excitations + t2aa = t2aa - 0.5 * np.kron(t1a, t1a).reshape(nelec_a, nvir_a, nelec_a, nvir_a).transpose(0,2,1,3).numpy() + t2bb = t2bb - 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() + 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 + ref_a = int(2**nelec_a - 1) + ref_b = int(2**nelec_b - 1) + + dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [1.])) + + # alpha -> alpha excitations + t1a_configs, t1a_signs = _excited_configurations(nelec_a, norb, 1) + dict_fcimatr.update(dict(zip(list(zip(t1a_configs, [ref_b] * len(t1a_configs))), t1a.ravel() * t1a_signs))) + + # beta -> beta excitations + t1b_configs, t1b_signs = _excited_configurations(nelec_b, norb, 1) + dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(t1b_configs), t1b_configs)), t1b.ravel() * t1b_signs))) + + # alpha, alpha -> alpha, alpha excitations + if nelec_a > 1 and nvir_a > 1: + t2aa_configs, t2aa_signs = _excited_configurations(nelec_a, norb, 2) + # select only unique excitations, via lower triangle of matrix + ooidx = np.tril_indices(nelec_a, -1) + vvidx = np.tril_indices(nvir_a, -1) + t2aa = t2aa[ooidx][:, vvidx[0], vvidx[1]] + dict_fcimatr.update(dict(zip(list(zip(t2aa_configs, [ref_b] * len(t2aa_configs))), t2aa.ravel() * t2aa_signs))) + + if nelec_b > 1 and nvir_b > 1: + t2bb_configs, t2bb_signs = _excited_configurations(nelec_b, norb, 2) + # select only unique excitations, via lower triangle of matrix + ooidx = np.tril_indices(nelec_b, -1) + vvidx = np.tril_indices(nvir_b, -1) + t2bb = t2bb[ooidx][:, vvidx[0], vvidx[1]] + dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(t2bb_configs), t2bb_configs)), t2bb.ravel() * t2bb_signs))) + + # alpha, beta -> alpha, beta excitations + rowvals, colvals = np.array(list(product(t1a_configs, t1b_configs)), dtype=int).T.numpy() + dict_fcimatr.update( + dict(zip(list(zip(rowvals, colvals)), t2ab.ravel() * np.kron(t1a_signs, t1b_signs).numpy())) + ) + + # renormalize, to get the HF coefficient (CC wavefunction not normalized) + norm = np.sqrt(np.sum(np.array(list(dict_fcimatr.values()))**2)).numpy() + dict_fcimatr = {key: value / norm for (key, value) in dict_fcimatr.items()} + + # filter based on tolerance cutoff + dict_fcimatr = {key: value for key, value in dict_fcimatr.items() if abs(value) > tol} + + return dict_fcimatr def _uccsd_state(ccsd_solver, tol=1e-15): r"""Construct a wavefunction from PySCF's `UCCSD` Solver object. @@ -692,7 +807,9 @@ def _uccsd_state(ccsd_solver, tol=1e-15): .. 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 \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. Args: ccsd_solver (PySCF UCCSD Class instance): the class object representing the UCCSD calculation in PySCF @@ -712,8 +829,8 @@ def _uccsd_state(ccsd_solver, tol=1e-15): >>> mycc = cc.UCCSD(myhf).run() >>> wf_ccsd = ccsd_state(mycc, tol=1e-1) >>> print(wf_ccsd) - {(7, 7): 0.8886988662500356, (11, 11): -0.30584653966926967, - (19, 19): -0.3058465396692694, (35, 35): -0.14507582719761972} + {(7, 7): 0.8886970081919591, (11, 11): -0.3058459002168582, + (19, 19): -0.30584590021685887, (35, 35): -0.14507552387854625} """ mol = ccsd_solver.mol @@ -728,10 +845,10 @@ def _uccsd_state(ccsd_solver, tol=1e-15): t1a, t1b = ccsd_solver.t1 t2aa, t2ab, t2bb = ccsd_solver.t2 # add in the disconnected part ( + 0.5 T_1^2) of double excitations - t2aa = t2aa - 0*0.5 * np.kron(t1a, t1a).reshape(nelec_a, nvir_a, nelec_a, nvir_a).transpose(0,2,1,3).numpy() - t2bb = t2bb - 0*0.5 * np.kron(t1b, t1b).reshape(nelec_b, nvir_b, nelec_b, nvir_b).transpose(0,2,1,3).numpy() + t2aa = t2aa - 0.5 * np.kron(t1a, t1a).reshape(nelec_a, nvir_a, nelec_a, nvir_a).transpose(0,2,1,3).numpy() + t2bb = t2bb - 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() - t2ab = t2ab.transpose(0,2,1,3) - 0*0.5 * np.kron(t1a, t1b).reshape(nelec_a, nvir_a, nelec_b, nvir_b).numpy() + t2ab = t2ab.transpose(0,2,1,3) - 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 ref_a = int(2**nelec_a - 1) From 32cd0b1e650b41ee2be0e51af248f5e9dc6a40ed Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Thu, 3 Aug 2023 11:57:40 -0400 Subject: [PATCH 42/60] added tests for rccsd and ccsd_state --- tests/qchem/of_tests/test_convert.py | 127 +++++++++++++++++++++++++-- 1 file changed, 120 insertions(+), 7 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 09eb5ea8a05..0e77774a83c 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -1076,8 +1076,8 @@ def test_cisd_state(molecule, basis, symm, hftype, wf_ref): "d2h", 1e-1, { - (1, 1): 0.9942969785398776, - (2, 2): -0.10664669927602176 + (1, 1): 0.9942969030671576, + (2, 2): -0.10664740292693174 } ), ( @@ -1086,10 +1086,10 @@ def test_cisd_state(molecule, basis, symm, hftype, wf_ref): "d2h", 1e-1, { - (7, 7): 0.8886988662500364, - (11, 11): -0.3058465396692694, - (19, 19): -0.3058465396692695, - (35, 35): -0.14507582719761614 + (7, 7): 0.8886970081919591, + (11, 11): -0.3058459002168582, + (19, 19): -0.30584590021685887, + (35, 35): -0.14507552387854625 }, ), ], @@ -1104,4 +1104,117 @@ def test_uccsd_state(molecule, basis, symm, tol, wf_ref): wf_ccsd = qchem.convert._uccsd_state(mycc, tol=tol) assert wf_ccsd.keys() == wf_ref.keys() - assert np.allclose(abs(np.array(list(wf_ccsd.values()))), abs(np.array(list(wf_ref.values())))) \ No newline at end of file + assert np.allclose(abs(np.array(list(wf_ccsd.values()))), abs(np.array(list(wf_ref.values())))) + +@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.9942969030671576, + (2, 2): -0.10664740292693242 + }, + ), + ( + [["Li", (0, 0, 0)], ["Li", (0, 0, 0.71)]], + "sto6g", + "d2h", + 1e-1, + { + (7, 7): 0.8886969878256522, + (11, 11): -0.30584590248164206, + (19, 19): -0.30584590248164145, + (35, 35): -0.14507552651170982 + }, + ), + ], +) +def test_rccsd_state(molecule, basis, symm, tol, wf_ref): + r"""Test that _rccsd_state returns the correct wavefunction.""" + + mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) + myhf = pyscf.scf.RHF(mol).run() + mycc = pyscf.cc.CCSD(myhf).run() + + wf_ccsd = qchem.convert._rccsd_state(mycc, tol=tol) + + 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( + ("molecule", "basis", "symm", "hftype", "wf_ref"), + [ + ( + [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], + "sto6g", + "d2h", + "uhf", + np.array( + [ + 0.0, + 0.0, + 0.0, + -0.1066467, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.99429698, + 0.0, + 0.0, + 0.0, + ] + ), + ), + ( + [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], + "sto6g", + "d2h", + "rhf", + np.array( + [ + 0.0, + 0.0, + 0.0, + 0.10664669927602179, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + -0.9942969785398778, + 0.0, + 0.0, + 0.0, + ] + ), + ), + ], +) +def test_ccsd_state(molecule, basis, symm, hftype, wf_ref): + r"""Test that ccsd_state returns the correct wavefunction.""" + + mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) + if hftype == "rhf": + myhf = pyscf.scf.RHF(mol).run() + mycc = pyscf.cc.CCSD(myhf).run() + elif hftype == "uhf": + myhf = pyscf.scf.UHF(mol).run() + mycc = pyscf.cc.UCCSD(myhf).run() + + wf_comp = qchem.convert.ccsd_state(mycc, hftype) + + # overall sign could be +/-1 different + assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) From c05ddae687adcafcb12d5b1f996a7d3cff573549 Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Thu, 3 Aug 2023 11:59:02 -0400 Subject: [PATCH 43/60] fix black formatting --- pennylane/qchem/convert.py | 119 +++++++++++++++++++-------- tests/qchem/of_tests/test_convert.py | 26 +++--- 2 files changed, 94 insertions(+), 51 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index c031e17de9f..56bc678b22a 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -669,6 +669,7 @@ def cisd_state(cisd_solver, hftype, tol=1e-15): return wf + def _rccsd_state(ccsd_solver, tol=1e-15): r"""Construct a wavefunction from PySCF's `CCSD` Solver object. @@ -684,13 +685,13 @@ 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 + 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 .. math:: - \exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}} = \left[ 1 + \hat{T}_1 + + \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. @@ -713,7 +714,7 @@ def _rccsd_state(ccsd_solver, tol=1e-15): >>> mycc = cc.CCSD(myhf).run() >>> wf_ccsd = ccsd_state(mycc, tol=1e-1) >>> print(wf_ccsd) - {(7, 7): 0.8886969878256522, (11, 11): -0.30584590248164206, + {(7, 7): 0.8886969878256522, (11, 11): -0.30584590248164206, (19, 19): -0.30584590248164145, (35, 35): -0.14507552651170982} """ @@ -729,13 +730,21 @@ def _rccsd_state(ccsd_solver, tol=1e-15): # build the full, unrestricted representation of the coupled cluster amplitudes t1a = ccsd_solver.t1 t1b = t1a - t2aa = ccsd_solver.t2 - ccsd_solver.t2.transpose(1,0,2,3) + t2aa = ccsd_solver.t2 - ccsd_solver.t2.transpose(1, 0, 2, 3) t2ab = ccsd_solver.t2.transpose(0, 2, 1, 3) t2bb = t2aa - # add in the disconnected part ( + 0.5 T_1^2) of double excitations - t2aa = t2aa - 0.5 * np.kron(t1a, t1a).reshape(nelec_a, nvir_a, nelec_a, nvir_a).transpose(0,2,1,3).numpy() - t2bb = t2bb - 0.5 * np.kron(t1b, t1b).reshape(nelec_b, nvir_b, nelec_b, nvir_b).transpose(0,2,1,3).numpy() + # add in the disconnected part ( + 0.5 T_1^2) of double excitations + t2aa = ( + t2aa + - 0.5 + * np.kron(t1a, t1a).reshape(nelec_a, nvir_a, nelec_a, nvir_a).transpose(0, 2, 1, 3).numpy() + ) + t2bb = ( + t2bb + - 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() t2ab = t2ab - 0.5 * np.kron(t1a, t1b).reshape(nelec_a, nvir_a, nelec_b, nvir_b).numpy() @@ -743,15 +752,19 @@ def _rccsd_state(ccsd_solver, tol=1e-15): ref_a = int(2**nelec_a - 1) ref_b = int(2**nelec_b - 1) - dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [1.])) + dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [1.0])) # alpha -> alpha excitations t1a_configs, t1a_signs = _excited_configurations(nelec_a, norb, 1) - dict_fcimatr.update(dict(zip(list(zip(t1a_configs, [ref_b] * len(t1a_configs))), t1a.ravel() * t1a_signs))) + dict_fcimatr.update( + dict(zip(list(zip(t1a_configs, [ref_b] * len(t1a_configs))), t1a.ravel() * t1a_signs)) + ) # beta -> beta excitations t1b_configs, t1b_signs = _excited_configurations(nelec_b, norb, 1) - dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(t1b_configs), t1b_configs)), t1b.ravel() * t1b_signs))) + dict_fcimatr.update( + dict(zip(list(zip([ref_a] * len(t1b_configs), t1b_configs)), t1b.ravel() * t1b_signs)) + ) # alpha, alpha -> alpha, alpha excitations if nelec_a > 1 and nvir_a > 1: @@ -759,16 +772,24 @@ def _rccsd_state(ccsd_solver, tol=1e-15): # select only unique excitations, via lower triangle of matrix ooidx = np.tril_indices(nelec_a, -1) vvidx = np.tril_indices(nvir_a, -1) - t2aa = t2aa[ooidx][:, vvidx[0], vvidx[1]] - dict_fcimatr.update(dict(zip(list(zip(t2aa_configs, [ref_b] * len(t2aa_configs))), t2aa.ravel() * t2aa_signs))) + t2aa = t2aa[ooidx][:, vvidx[0], vvidx[1]] + dict_fcimatr.update( + dict( + zip(list(zip(t2aa_configs, [ref_b] * len(t2aa_configs))), t2aa.ravel() * t2aa_signs) + ) + ) if nelec_b > 1 and nvir_b > 1: t2bb_configs, t2bb_signs = _excited_configurations(nelec_b, norb, 2) # select only unique excitations, via lower triangle of matrix ooidx = np.tril_indices(nelec_b, -1) vvidx = np.tril_indices(nvir_b, -1) - t2bb = t2bb[ooidx][:, vvidx[0], vvidx[1]] - dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(t2bb_configs), t2bb_configs)), t2bb.ravel() * t2bb_signs))) + t2bb = t2bb[ooidx][:, vvidx[0], vvidx[1]] + dict_fcimatr.update( + dict( + zip(list(zip([ref_a] * len(t2bb_configs), t2bb_configs)), t2bb.ravel() * t2bb_signs) + ) + ) # alpha, beta -> alpha, beta excitations rowvals, colvals = np.array(list(product(t1a_configs, t1b_configs)), dtype=int).T.numpy() @@ -777,7 +798,7 @@ def _rccsd_state(ccsd_solver, tol=1e-15): ) # renormalize, to get the HF coefficient (CC wavefunction not normalized) - norm = np.sqrt(np.sum(np.array(list(dict_fcimatr.values()))**2)).numpy() + norm = np.sqrt(np.sum(np.array(list(dict_fcimatr.values())) ** 2)).numpy() dict_fcimatr = {key: value / norm for (key, value) in dict_fcimatr.items()} # filter based on tolerance cutoff @@ -785,6 +806,7 @@ def _rccsd_state(ccsd_solver, tol=1e-15): return dict_fcimatr + def _uccsd_state(ccsd_solver, tol=1e-15): r"""Construct a wavefunction from PySCF's `UCCSD` Solver object. @@ -800,13 +822,13 @@ 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 + 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 .. math:: - \exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}} = \left[ 1 + \hat{T}_1 + + \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. @@ -829,7 +851,7 @@ def _uccsd_state(ccsd_solver, tol=1e-15): >>> mycc = cc.UCCSD(myhf).run() >>> wf_ccsd = ccsd_state(mycc, tol=1e-1) >>> print(wf_ccsd) - {(7, 7): 0.8886970081919591, (11, 11): -0.3058459002168582, + {(7, 7): 0.8886970081919591, (11, 11): -0.3058459002168582, (19, 19): -0.30584590021685887, (35, 35): -0.14507552387854625} """ @@ -843,26 +865,41 @@ def _uccsd_state(ccsd_solver, tol=1e-15): nvir_a, nvir_b = norb - nelec_a, norb - nelec_b t1a, t1b = ccsd_solver.t1 - t2aa, t2ab, t2bb = ccsd_solver.t2 - # add in the disconnected part ( + 0.5 T_1^2) of double excitations - t2aa = t2aa - 0.5 * np.kron(t1a, t1a).reshape(nelec_a, nvir_a, nelec_a, nvir_a).transpose(0,2,1,3).numpy() - t2bb = t2bb - 0.5 * np.kron(t1b, t1b).reshape(nelec_b, nvir_b, nelec_b, nvir_b).transpose(0,2,1,3).numpy() + t2aa, t2ab, t2bb = ccsd_solver.t2 + # add in the disconnected part ( + 0.5 T_1^2) of double excitations + t2aa = ( + t2aa + - 0.5 + * np.kron(t1a, t1a).reshape(nelec_a, nvir_a, nelec_a, nvir_a).transpose(0, 2, 1, 3).numpy() + ) + t2bb = ( + t2bb + - 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() - t2ab = t2ab.transpose(0,2,1,3) - 0.5 * np.kron(t1a, t1b).reshape(nelec_a, nvir_a, nelec_b, nvir_b).numpy() + t2ab = ( + t2ab.transpose(0, 2, 1, 3) + - 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 ref_a = int(2**nelec_a - 1) ref_b = int(2**nelec_b - 1) - dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [1.])) + dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [1.0])) # alpha -> alpha excitations t1a_configs, t1a_signs = _excited_configurations(nelec_a, norb, 1) - dict_fcimatr.update(dict(zip(list(zip(t1a_configs, [ref_b] * len(t1a_configs))), t1a.ravel() * t1a_signs))) + dict_fcimatr.update( + dict(zip(list(zip(t1a_configs, [ref_b] * len(t1a_configs))), t1a.ravel() * t1a_signs)) + ) # beta -> beta excitations t1b_configs, t1b_signs = _excited_configurations(nelec_b, norb, 1) - dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(t1b_configs), t1b_configs)), t1b.ravel() * t1b_signs))) + dict_fcimatr.update( + dict(zip(list(zip([ref_a] * len(t1b_configs), t1b_configs)), t1b.ravel() * t1b_signs)) + ) # alpha, alpha -> alpha, alpha excitations if nelec_a > 1 and nvir_a > 1: @@ -870,16 +907,24 @@ def _uccsd_state(ccsd_solver, tol=1e-15): # select only unique excitations, via lower triangle of matrix ooidx = np.tril_indices(nelec_a, -1) vvidx = np.tril_indices(nvir_a, -1) - t2aa = t2aa[ooidx][:, vvidx[0], vvidx[1]] - dict_fcimatr.update(dict(zip(list(zip(t2aa_configs, [ref_b] * len(t2aa_configs))), t2aa.ravel() * t2aa_signs))) + t2aa = t2aa[ooidx][:, vvidx[0], vvidx[1]] + dict_fcimatr.update( + dict( + zip(list(zip(t2aa_configs, [ref_b] * len(t2aa_configs))), t2aa.ravel() * t2aa_signs) + ) + ) if nelec_b > 1 and nvir_b > 1: t2bb_configs, t2bb_signs = _excited_configurations(nelec_b, norb, 2) # select only unique excitations, via lower triangle of matrix ooidx = np.tril_indices(nelec_b, -1) vvidx = np.tril_indices(nvir_b, -1) - t2bb = t2bb[ooidx][:, vvidx[0], vvidx[1]] - dict_fcimatr.update(dict(zip(list(zip([ref_a] * len(t2bb_configs), t2bb_configs)), t2bb.ravel() * t2bb_signs))) + t2bb = t2bb[ooidx][:, vvidx[0], vvidx[1]] + dict_fcimatr.update( + dict( + zip(list(zip([ref_a] * len(t2bb_configs), t2bb_configs)), t2bb.ravel() * t2bb_signs) + ) + ) # alpha, beta -> alpha, beta excitations rowvals, colvals = np.array(list(product(t1a_configs, t1b_configs)), dtype=int).T.numpy() @@ -888,7 +933,7 @@ def _uccsd_state(ccsd_solver, tol=1e-15): ) # renormalize, to get the HF coefficient (CC wavefunction not normalized) - norm = np.sqrt(np.sum(np.array(list(dict_fcimatr.values()))**2)).numpy() + norm = np.sqrt(np.sum(np.array(list(dict_fcimatr.values())) ** 2)).numpy() dict_fcimatr = {key: value / norm for (key, value) in dict_fcimatr.items()} # filter based on tolerance cutoff @@ -896,6 +941,7 @@ def _uccsd_state(ccsd_solver, tol=1e-15): return dict_fcimatr + def ccsd_state(ccsd_solver, hftype, tol=1e-15): r"""Construct a wavefunction from PySCF CCSD output. @@ -924,7 +970,7 @@ def ccsd_state(ccsd_solver, hftype, tol=1e-15): >>> wf_ccsd = ccsd_state(mycc, tol=1e-1) >>> print(wf_cisd) [ 0. 0. 0. -0.1066467 0. 0. 0. - 0. 0. 0. 0. 0. 0.99429698 + 0. 0. 0. 0. 0. 0.99429698 0. 0. 0. ] """ @@ -938,6 +984,7 @@ def ccsd_state(ccsd_solver, hftype, tol=1e-15): return wf + def wfdict_to_statevector(wf_dict, norbs): r"""Convert a wavefunction in sparce dictionary format to a Pennylane's statevector. diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 0e77774a83c..384675513fa 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -1067,6 +1067,7 @@ def test_cisd_state(molecule, basis, symm, hftype, wf_ref): # overall sign could be +/-1 different assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) + @pytest.mark.parametrize( ("molecule", "basis", "symm", "tol", "wf_ref"), [ @@ -1075,10 +1076,7 @@ def test_cisd_state(molecule, basis, symm, hftype, wf_ref): "sto6g", "d2h", 1e-1, - { - (1, 1): 0.9942969030671576, - (2, 2): -0.10664740292693174 - } + {(1, 1): 0.9942969030671576, (2, 2): -0.10664740292693174}, ), ( [["Li", (0, 0, 0)], ["Li", (0, 0, 0.71)]], @@ -1087,9 +1085,9 @@ def test_cisd_state(molecule, basis, symm, hftype, wf_ref): 1e-1, { (7, 7): 0.8886970081919591, - (11, 11): -0.3058459002168582, - (19, 19): -0.30584590021685887, - (35, 35): -0.14507552387854625 + (11, 11): -0.3058459002168582, + (19, 19): -0.30584590021685887, + (35, 35): -0.14507552387854625, }, ), ], @@ -1106,6 +1104,7 @@ def test_uccsd_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( ("molecule", "basis", "symm", "tol", "wf_ref"), [ @@ -1114,10 +1113,7 @@ def test_uccsd_state(molecule, basis, symm, tol, wf_ref): "sto6g", "d2h", 1e-1, - { - (1, 1): 0.9942969030671576, - (2, 2): -0.10664740292693242 - }, + {(1, 1): 0.9942969030671576, (2, 2): -0.10664740292693242}, ), ( [["Li", (0, 0, 0)], ["Li", (0, 0, 0.71)]], @@ -1125,10 +1121,10 @@ def test_uccsd_state(molecule, basis, symm, tol, wf_ref): "d2h", 1e-1, { - (7, 7): 0.8886969878256522, - (11, 11): -0.30584590248164206, - (19, 19): -0.30584590248164145, - (35, 35): -0.14507552651170982 + (7, 7): 0.8886969878256522, + (11, 11): -0.30584590248164206, + (19, 19): -0.30584590248164145, + (35, 35): -0.14507552651170982, }, ), ], From 699373d0477b321f26043448cf81e781811f450e Mon Sep 17 00:00:00 2001 From: soranjh Date: Thu, 3 Aug 2023 13:03:03 -0400 Subject: [PATCH 44/60] modify docs --- pennylane/qchem/convert.py | 65 ++++++++++++++-------------- tests/qchem/of_tests/test_convert.py | 2 +- 2 files changed, 33 insertions(+), 34 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 208590c29e7..9fa3c779086 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -410,16 +410,19 @@ def _excitations(electrons, orbitals): def _excited_configurations(electrons, orbitals, excitation): - r"""Generate excited states from a Hartree-Fock reference state. + r"""Generate excited configurations from a Hartree-Fock reference state. - This function generates excited states in the form of binary strings or integers representing - them, from a Hartree-Fock (HF) reference state. The HF state is assumed to be - :math:`|1 1 ...1 0 ... 0 0 \rangle` where the number of :math:`1` and :math:`0` elements - represents the number of occupied and unoccupied spin orbitals, respectively. The string - representation of the state is obtained by converting the occupation-number vector to a string, - e.g., ``111000`` represents :math:`|1 1 1 0 0 0 \rangle. The number representation of the state - is the integer corresonding to the binary number, e.g., ``111000`` is represented by - :math:`int('111000', 2) = 56`. + This function generates excited configurations in the form of integers representing a binary + string, e.g., :math:`|1 1 0 1 0 0 \rangle is represented by :math:`int('110100', 2) = 52`. + + The excited configurations are generated from a Hartree-Fock (HF) reference state. The HF state + is assumed to have the form :math:`|1 1 ...1 0 ... 0 0 \rangle` where the number of :math:`1` + and :math:`0` elements are the number of occupied and unoccupied spin orbitals, respectively. + The string representation of the state is obtained by converting the occupation-number vector to + a string, e.g., ``111000`` to represent :math:`|1 1 1 0 0 0 \rangle. + + Each excited configuration has a sign, :math:`+1` or :math:`-1`, that is obtained by reordering + the creation operators. Args: electrons (int): number of electrons @@ -427,7 +430,8 @@ def _excited_configurations(electrons, orbitals, excitation): excitation (int): number of excited electrons Returns: - tuple(array, array): arrays of excited states and signs obtained by reordering of the creation operators + tuple(array, array): arrays of excited configurations and signs obtained by reordering the + creation operators **Example** @@ -464,19 +468,18 @@ def _excited_configurations(electrons, orbitals, excitation): def _ucisd_state(cisd_solver, tol=1e-15): - r"""Construct a wavefunction from PySCF's `UCISD` Solver object. - - The wavefunction is represented as a dictionary where the keys are tuples representing a - configuration, which corresponds to a Slater determinant, and the values are the CI coefficient - corresponding to that Slater determinant. Each dictionary key is a tuple of two integers. The - binary representation of these integers correspond to a specific configuration, or Slater - determinant. The first number represents the configuration of alpha electrons and the second - number represents the beta electrons. For instance, the Hartree-Fock state - :math:`|1 1 0 0 \rangle` will be represented by the binary string `0011`. This string can be - splited to `01` and `01` for the alpha and beta electrons. The integer corresponding to `01` is - `1`. Then the dictionary representation of a state with `0.99` contribution of the Hartree-Fock - state and `0.01` contribution from the doubly-excited state will be - `{(1, 1): 0.99, (2, 2): 0.01}`. + r"""Construct a wavefunction from PySCF's `UCISD` solver object. + + 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 that Slater + determinant. Each dictionary key is a tuple of two integers. The binary representation of these + integers correspond to a specific configuration such that 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`. This string can be splitted to `01` and `01` for + the alpha and beta electrons. The integer corresponding to `01` is `1`. Then the dictionary + representation of a state with `0.99` contribution of the Hartree-Fock state and `0.01` + contribution from the doubly-excited state will be `{(1, 1): 0.99, (2, 2): 0.01}`. Args: cisd_solver (PySCF UCISD Class instance): the class object representing the UCISD calculation in PySCF @@ -485,7 +488,7 @@ def _ucisd_state(cisd_solver, tol=1e-15): Returns: dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` - having binary represention corresponding to the Fock occupation vector in alpha and beta + 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** @@ -553,20 +556,17 @@ def _ucisd_state(cisd_solver, tol=1e-15): def import_state(solver, method, tol=1e-15): - r"""Construct a wavefunction from PySCF output. + r"""Convert an external wavefunction to a PennyLane state vector. Args: solver: external wavefunction object that will be converted - method (str): keyword specifying the method of calculation + method (str): keyword specifying the calculation method 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). + determinants with coefficients smaller than this tolerance are discarded. Returns: - dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` - having binary represention corresponding to the Fock occupation vector in alpha and beta - spin sectors, respectively, and coeff being the CI coefficients of those configurations. + statevector (array): normalized state vector of length 2**len(number_of_spinorbitals) **Example** @@ -602,10 +602,9 @@ def _wfdict_to_statevector(wf_dict, norbs): norbs (int): number of molecular orbitals in the system Returns: - statevector (list): normalized statevector of length 2^(number_of_spinorbitals), standard in Pennylane + statevector (array): normalized state vector of length 2^(number_of_spinorbitals) """ - statevector = np.zeros(2 ** (2 * norbs), dtype=complex) for (int_a, int_b), coeff in wf_dict.items(): diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 8b678c5f7f6..c332d3ede45 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -975,7 +975,7 @@ def test_import_state(molecule, basis, symm, method, wf_ref): wf_comp = qchem.convert.import_state(myci, method) - # overall sign could be +/-1 different + # overall sign could be different in each PySCF run assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) From e3e8daefc2343c0720897cd10974e54fc918278d Mon Sep 17 00:00:00 2001 From: soranjh Date: Thu, 3 Aug 2023 13:54:47 -0400 Subject: [PATCH 45/60] resolve conflicts in tests --- tests/qchem/of_tests/test_convert.py | 34 ++-------------------------- 1 file changed, 2 insertions(+), 32 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index fb66359efc2..4aa07598cb9 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -964,44 +964,14 @@ def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): ] ), ), - ( - [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], - "sto6g", - "d2h", - "rcisd", - np.array( - [ - 0.0, - 0.0, - 0.0, - 0.10664669927602179, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - -0.9942969785398778, - 0.0, - 0.0, - 0.0, - ] - ), - ), ], ) def test_import_state(molecule, basis, symm, method, wf_ref): r"""Test that cisd_state returns the correct wavefunction.""" mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) - if method == "rcisd": - myhf = pyscf.scf.RHF(mol).run() - myci = pyscf.ci.CISD(myhf).run() - elif method == "ucisd": - myhf = pyscf.scf.UHF(mol).run() - myci = pyscf.ci.UCISD(myhf).run() + myhf = pyscf.scf.UHF(mol).run() + myci = pyscf.ci.UCISD(myhf).run() wf_comp = qchem.convert.import_state(myci, method) From 3f55198cc21bcaf4edc97af5944fe87359981ff6 Mon Sep 17 00:00:00 2001 From: soranjh Date: Thu, 3 Aug 2023 14:16:51 -0400 Subject: [PATCH 46/60] merge functions --- pennylane/qchem/convert.py | 54 +++------------- tests/qchem/of_tests/test_convert.py | 95 +++++----------------------- 2 files changed, 24 insertions(+), 125 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index ad7ffffbcbe..c1c3b366808 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -578,13 +578,18 @@ def import_state(solver, method, tol=1e-15): >>> print(wf_cisd) {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} """ - if method == "ucisd": wf_dict = _ucisd_state(solver, tol=tol) + elif method == "rcisd": + wf_dict = _rcisd_state(solver, tol=tol) + elif method == "rccsd": + wf_dict = _rccsd_state(solver, tol=tol) + elif method == "uccsd": + wf_dict = _uccsd_state(solver, tol=tol) else: raise ValueError( - "The supported method options are 'rcisd' for restricted and 'ucisd' for unrestricted" - " CISD calculations." + "The supported method options are 'rcisd','ucisd', 'rccsd', and 'uccsd' for restricted" + " and unrestricted CISD and CCSD calculations." ) wf = _wfdict_to_statevector(wf_dict, solver.mol.nao) @@ -965,46 +970,3 @@ 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 ccsd_state(ccsd_solver, hftype, tol=1e-15): - r"""Construct a wavefunction from PySCF CCSD output. - - Args: - ccsd_solver (PySCF CCSD or UCCSD Class instance): the class object representing the - CCSD / UCCSD calculation in PySCF - - hftype (str): keyword specifying whether restricted or unrestricted CCSD was performed. - The options are 'rhf' for restricted, and 'uhf' for unrestricted. - - 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). - - Returns: - dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` - having binary represention 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.UHF(mol).run() - >>> mycc = ci.UCCSD(myhf).run() - >>> wf_ccsd = ccsd_state(mycc, tol=1e-1) - >>> print(wf_cisd) - [ 0. 0. 0. -0.1066467 0. 0. 0. - 0. 0. 0. 0. 0. 0.99429698 - 0. 0. 0. ] - """ - - if hftype == "rhf": - wf_dict = _rccsd_state(ccsd_solver, tol=tol) - elif hftype == "uhf": - wf_dict = _uccsd_state(ccsd_solver, tol=tol) - else: - raise ValueError("Only restricted, 'rhf', and unrestricted, 'uhf', are supported.") - wf = _wfdict_to_statevector(wf_dict, ccsd_solver.mol.nao) - - return wf diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 4aa07598cb9..dd995e08fad 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -936,13 +936,12 @@ def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): @pytest.mark.parametrize( - ("molecule", "basis", "symm", "method", "wf_ref"), + ("molecule", "basis", "symm", "wf_ref"), [ ( [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], "sto6g", "d2h", - "ucisd", np.array( [ 0.0, @@ -966,14 +965,26 @@ def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): ), ], ) +@pytest.mark.parametrize("method", ["rcisd", "ucisd", "rccsd", "uccsd"]) def test_import_state(molecule, basis, symm, method, wf_ref): r"""Test that cisd_state returns the correct wavefunction.""" mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) - myhf = pyscf.scf.UHF(mol).run() - myci = pyscf.ci.UCISD(myhf).run() - wf_comp = qchem.convert.import_state(myci, method) + if method == "rcisd": + myhf = pyscf.scf.RHF(mol).run() + solver = pyscf.ci.CISD(myhf).run() + elif method == "ucisd": + myhf = pyscf.scf.UHF(mol).run() + solver = pyscf.ci.UCISD(myhf).run() + elif method == "rccsd": + myhf = pyscf.scf.RHF(mol).run() + solver = pyscf.cc.CCSD(myhf).run() + elif method == "uccsd": + myhf = pyscf.scf.UHF(mol).run() + solver = pyscf.cc.UCCSD(myhf).run() + + wf_comp = qchem.convert.import_state(solver, method) # overall sign could be different in each PySCF run assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) @@ -1061,77 +1072,3 @@ 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( - ("molecule", "basis", "symm", "hftype", "wf_ref"), - [ - ( - [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], - "sto6g", - "d2h", - "uhf", - np.array( - [ - 0.0, - 0.0, - 0.0, - -0.1066467, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.99429698, - 0.0, - 0.0, - 0.0, - ] - ), - ), - ( - [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], - "sto6g", - "d2h", - "rhf", - np.array( - [ - 0.0, - 0.0, - 0.0, - 0.10664669927602179, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - -0.9942969785398778, - 0.0, - 0.0, - 0.0, - ] - ), - ), - ], -) -def test_ccsd_state(molecule, basis, symm, hftype, wf_ref): - r"""Test that ccsd_state returns the correct wavefunction.""" - - mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) - if hftype == "rhf": - myhf = pyscf.scf.RHF(mol).run() - mycc = pyscf.cc.CCSD(myhf).run() - elif hftype == "uhf": - myhf = pyscf.scf.UHF(mol).run() - mycc = pyscf.cc.UCCSD(myhf).run() - - wf_comp = qchem.convert.ccsd_state(mycc, hftype) - - # overall sign could be +/-1 different - assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) From 9ad777f42ae3586f903ca09d909e0b1d0c583a3b Mon Sep 17 00:00:00 2001 From: soranjh Date: Thu, 3 Aug 2023 15:09:30 -0400 Subject: [PATCH 47/60] add small modifications --- pennylane/qchem/convert.py | 56 ++++++++++++++------------ tests/qchem/of_tests/test_convert.py | 59 ++++++++++++++++++++++++++++ 2 files changed, 90 insertions(+), 25 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index c1c3b366808..faf351f5990 100644 --- a/pennylane/qchem/convert.py +++ b/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() diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index dd995e08fad..34b6bade9a3 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/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"), [ From f2bee885e874d88ec168069cf5b48fb1a6c2d0e1 Mon Sep 17 00:00:00 2001 From: soranjh Date: Fri, 4 Aug 2023 15:16:07 -0400 Subject: [PATCH 48/60] update with master --- pennylane/qchem/convert.py | 76 +++++++++++++---------- tests/qchem/of_tests/test_convert.py | 93 ++++++++++++++-------------- 2 files changed, 89 insertions(+), 80 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index faf351f5990..e504c61d587 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -391,7 +391,7 @@ def _excitations(electrons, orbitals): >>> electrons = 2 >>> orbitals = 4 >>> _excitations(electrons, orbitals) - ([[[0, 2]], [[0, 3]], [[1, 2]], [[1, 3]]], [[0, 1, 2, 3]]) + ([[0, 2], [0, 3], [1, 2], [1, 3]], [[0, 1, 2, 3]]) """ singles_p, singles_q = [], [] doubles_pq, doubles_rs = [], [] @@ -435,12 +435,17 @@ def _excited_configurations(electrons, orbitals, excitation): **Example** - >>> electrons = 2 + >>> electrons = 3 >>> orbitals = 5 >>> excitation = 2 - >>> _excitated_states(electrons, orbitals, excitation) - (array([28, 26, 25]), array([ 1, -1, 1])) + >>> _excited_configurations(electrons, orbitals, excitation) + ([28, 26, 25], [1, -1, 1]) """ + if excitation not in [1, 2]: + raise ValueError( + "Only single (excitation = 1) and double (excitation = 2) excitations are supported." + ) + hf_state = qml.qchem.hf_state(electrons, orbitals) singles, doubles = _excitations(electrons, orbitals) states, signs = [], [] @@ -448,14 +453,14 @@ def _excited_configurations(electrons, orbitals, excitation): if excitation == 1: for s in singles: state = hf_state.copy() - state[s] = state[s[::-1]] + state[s] = state[[s[1], s[0]]] # apply single excitation states += [state] signs.append((-1) ** len(range(s[0], electrons - 1))) if excitation == 2: for d in doubles: state = hf_state.copy() - state[d] = state[[d[2], d[3], d[0], d[1]]] + state[d] = state[[d[2], d[3], d[0], d[1]]] # apply double excitation states += [state] order_pq = len(range(d[0], electrons - 1)) order_rs = len(range(d[1], electrons - 1)) @@ -497,7 +502,7 @@ def _ucisd_state(cisd_solver, tol=1e-15): >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g') >>> myhf = scf.UHF(mol).run() >>> myci = ci.UCISD(myhf).run() - >>> wf_cisd = cisd_state(myci, tol=1e-1) + >>> wf_cisd = _ucisd_state(myci, tol=1e-1) >>> print(wf_cisd) {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} """ @@ -555,15 +560,15 @@ def _ucisd_state(cisd_solver, tol=1e-15): return dict_fcimatr -def import_state(solver, method, tol=1e-15): +def import_state(solver, tol=1e-15): r"""Convert an external wavefunction to a PennyLane state vector. Args: solver: external wavefunction object that will be converted - method (str): keyword specifying the calculation method + tol (float): the tolerance for discarding Slater determinants with small coefficients - tol (float): the tolerance to which the wavefunction is being built -- Slater - determinants with coefficients smaller than this tolerance are discarded. + Raises: + ValueError: if ``method`` is not supported Returns: statevector (array): normalized state vector of length 2**len(number_of_spinorbitals) @@ -574,49 +579,52 @@ def import_state(solver, method, tol=1e-15): >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g') >>> myhf = scf.UHF(mol).run() >>> myci = ci.UCISD(myhf).run() - >>> wf_cisd = cisd_state(myci, tol=1e-1) + >>> wf_cisd = qml.qchem.convert.import_state(myci, tol=1e-1) >>> print(wf_cisd) - {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159} + [ 0. +0.j 0. +0.j 0. +0.j 0.1066467 +0.j + 0. +0.j 0. +0.j 0. +0.j 0. +0.j + 0. +0.j 0. +0.j 0. +0.j 0. +0.j + -0.99429698+0.j 0. +0.j 0. +0.j 0. +0.j] """ - if method == "ucisd": + try: + import pyscf + except ImportError as Error: + raise ImportError( + "This feature requires pyscf. It can be installed with: pip install pyscf" + ) from Error + + if isinstance(solver, pyscf.ci.ucisd.UCISD): wf_dict = _ucisd_state(solver, tol=tol) - elif method == "rcisd": - wf_dict = _rcisd_state(solver, tol=tol) - elif method == "rccsd": - wf_dict = _rccsd_state(solver, tol=tol) - elif method == "uccsd": - wf_dict = _uccsd_state(solver, tol=tol) else: - raise ValueError( - "The supported method options are 'rcisd','ucisd', 'rccsd', and 'uccsd' for restricted" - " and unrestricted CISD and CCSD calculations." - ) + raise ValueError("The supported option is 'ucisd' for unrestricted CISD calculations.") wf = _wfdict_to_statevector(wf_dict, solver.mol.nao) return wf def _wfdict_to_statevector(wf_dict, norbs): - r"""Convert a wavefunction in sparce dictionary format to a PennyLane statevector. + r"""Convert a wavefunction in sparse dictionary format to a PennyLane statevector. - Args: - wf_dict (wf_dict format): dictionary with keys (int_a,int_b) being integers whose binary representation - shows the Fock occupation vector for alpha and beta electrons; and values being the CI - coefficients. + In the sparse dictionary format, the keys (int_a, int_b) are integers whose binary + representation shows the Fock occupation vector for alpha and beta electrons and values are the + CI coefficients. - norbs (int): number of molecular orbitals in the system + Args: + wf_dict (dict): the sparse dictionary format of a wavefunction + norbs (int): number of molecular orbitals Returns: statevector (array): normalized state vector of length 2^(number_of_spinorbitals) - """ statevector = np.zeros(2 ** (2 * norbs), dtype=complex) for (int_a, int_b), coeff in wf_dict.items(): - bin_a, bin_b = bin(int_a)[2:][::-1], bin(int_b)[2:][::-1] + bin_a = bin(int_a)[2:][::-1] + bin_b = bin(int_b)[2:][::-1] + bin_a += "0" * (norbs - len(bin_a)) + bin_b += "0" * (norbs - len(bin_b)) bin_ab = "".join(i + j for i, j in zip(bin_a, bin_b)) - bin_ab_full = bin_ab + "0" * (2 * norbs - len(bin_ab)) - statevector[int(bin_ab_full, 2)] += coeff + statevector[int(bin_ab, 2)] += coeff statevector = statevector / np.sqrt(np.sum(statevector**2)) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 34b6bade9a3..c3985aac48d 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -901,37 +901,30 @@ def test_ucisd_state(molecule, basis, symm, tol, wf_ref): @pytest.mark.parametrize( - ("wf_dict", "n_orbitals", "wf_ref"), - [ - ( # -0.99 |1100> + 0.11 |0011> - {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602179}, + ("wf_dict", "n_orbitals", "string_ref", "coeff_ref"), + [ # reference data were obtained manually + ( # 0.87006284 |1100> + 0.3866946 |1001> + 0.29002095 |0110> + 0.09667365 |0011> + {(1, 1): 0.87006284, (1, 2): 0.3866946, (2, 1): 0.29002095, (2, 2): 0.09667365}, 2, - np.array( - [ - 0.0, - 0.0, - 0.0, - 0.1066467, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - -0.99429698, - 0.0, - 0.0, - 0.0, - ] - ), + ["1100", "1001", "0110", "0011"], + [0.87006284, 0.3866946, 0.29002095, 0.09667365], + ), + ( # 0.80448616 |110000> + 0.53976564 |001100> + 0.22350293 |000011> + 0.10724511 |100100> + {(1, 1): 0.80448616, (2, 2): 0.53976564, (4, 4): 0.22350293, (1, 2): 0.10724511}, + 3, + ["110000", "001100", "000011", "100100"], + [0.80448616, 0.53976564, 0.22350293, 0.10724511], ), ], ) -def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): +def test_wfdict_to_statevector(wf_dict, n_orbitals, string_ref, coeff_ref): r"""Test that _wfdict_to_statevector returns the correct statevector.""" + wf_ref = np.zeros(2 ** (n_orbitals * 2)) + idx_nonzero = [int(s, 2) for s in string_ref] + wf_ref[idx_nonzero] = coeff_ref + wf_comp = qchem.convert._wfdict_to_statevector(wf_dict, n_orbitals) + assert np.allclose(wf_comp, wf_ref) @@ -965,39 +958,47 @@ def test_wfdict_to_statevector(wf_dict, n_orbitals, wf_ref): ), ], ) -@pytest.mark.parametrize("method", ["rcisd", "ucisd", "rccsd", "uccsd"]) -def test_import_state(molecule, basis, symm, method, wf_ref): +def test_import_state(molecule, basis, symm, wf_ref): r"""Test that cisd_state returns the correct wavefunction.""" mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) + myhf = pyscf.scf.UHF(mol).run() + myci = pyscf.ci.UCISD(myhf).run() - if method == "rcisd": - myhf = pyscf.scf.RHF(mol).run() - solver = pyscf.ci.CISD(myhf).run() - elif method == "ucisd": - myhf = pyscf.scf.UHF(mol).run() - solver = pyscf.ci.UCISD(myhf).run() - elif method == "rccsd": - myhf = pyscf.scf.RHF(mol).run() - solver = pyscf.cc.CCSD(myhf).run() - elif method == "uccsd": - myhf = pyscf.scf.UHF(mol).run() - solver = pyscf.cc.UCCSD(myhf).run() - - wf_comp = qchem.convert.import_state(solver, method) + wf_comp = qchem.convert.import_state(myci) # overall sign could be different in each PySCF run assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) def test_import_state_error(): - r"""Test that an error is raised if a wrong/not-supported method symbol is entered.""" + r"""Test that an error is raised by import_state if a wrong object is entered.""" + + myci = "wrongobject" - myci = pyscf.ci.UCISD - method = "wrongmethod" + with pytest.raises(ValueError, match="The supported option"): + _ = qchem.convert.import_state(myci) + + +@pytest.mark.parametrize(("excitation"), [-1, 0, 3]) +def test_excited_configurations_error(excitation): + r"""Test that an error is raised by _excited_configurations if a wrong excitation is entered.""" + with pytest.raises(ValueError, match="excitations are supported"): + _ = qchem.convert._excited_configurations(2, 4, excitation) + + +def test_fail_import_pyscf(monkeypatch): + """Test if an ImportError is raised when pyscf is requested but not installed.""" + + mol = pyscf.gto.M(atom=[["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], basis="sto6g") + myhf = pyscf.scf.UHF(mol).run() + myci = pyscf.ci.UCISD(myhf).run() + + with monkeypatch.context() as m: + m.setitem(sys.modules, "pyscf", None) - with pytest.raises(ValueError, match="The supported method options are"): - _ = qchem.convert.import_state(myci, method) + with pytest.raises(ImportError, match="This feature requires pyscf"): + qml.qchem.convert.import_state(myci) @pytest.mark.parametrize( From 8738750b9a1fd24ca34f47e935b667935433133c Mon Sep 17 00:00:00 2001 From: soranjh Date: Fri, 4 Aug 2023 16:26:28 -0400 Subject: [PATCH 49/60] expand import state test --- pennylane/qchem/convert.py | 13 +++++++++++-- tests/qchem/of_tests/test_convert.py | 24 ++++++++++++++++++------ 2 files changed, 29 insertions(+), 8 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index e504c61d587..7395617169e 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -593,10 +593,19 @@ def import_state(solver, tol=1e-15): "This feature requires pyscf. It can be installed with: pip install pyscf" ) from Error - if isinstance(solver, pyscf.ci.ucisd.UCISD): + if "RCISD" in str(solver.__str__): + wf_dict = _rcisd_state(solver, tol=tol) + elif "UCISD" in str(solver.__str__): wf_dict = _ucisd_state(solver, tol=tol) + elif "RCCSD" in str(solver.__str__): + wf_dict = _rccsd_state(solver, tol=tol) + elif "UCCSD" in str(solver.__str__): + wf_dict = _uccsd_state(solver, tol=tol) else: - raise ValueError("The supported option is 'ucisd' for unrestricted CISD calculations.") + raise ValueError( + "The supported objects are RCISD, UCISD, RCCSD, and UCCSD for restricted and" + " unrestricted configuration interaction and coupled cluster calculations." + ) wf = _wfdict_to_statevector(wf_dict, solver.mol.nao) return wf diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index c3985aac48d..df35184cdfb 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -958,14 +958,26 @@ def test_wfdict_to_statevector(wf_dict, n_orbitals, string_ref, coeff_ref): ), ], ) -def test_import_state(molecule, basis, symm, wf_ref): - r"""Test that cisd_state returns the correct wavefunction.""" +@pytest.mark.parametrize("method", ["rcisd", "ucisd", "rccsd", "uccsd"]) +def test_import_state(molecule, basis, symm, method, wf_ref): + r"""Test that import_state returns the correct wavefunction.""" mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) - myhf = pyscf.scf.UHF(mol).run() - myci = pyscf.ci.UCISD(myhf).run() - wf_comp = qchem.convert.import_state(myci) + if method == "rcisd": + myhf = pyscf.scf.RHF(mol).run() + solver = pyscf.ci.cisd.RCISD(myhf).run() + elif method == "ucisd": + myhf = pyscf.scf.UHF(mol).run() + solver = pyscf.ci.ucisd.UCISD(myhf).run() + elif method == "rccsd": + myhf = pyscf.scf.RHF(mol).run() + solver = pyscf.cc.rccsd.RCCSD(myhf).run() + elif method == "uccsd": + myhf = pyscf.scf.UHF(mol).run() + solver = pyscf.cc.uccsd.UCCSD(myhf).run() + + wf_comp = qchem.convert.import_state(solver) # overall sign could be different in each PySCF run assert np.allclose(wf_comp, wf_ref) or np.allclose(wf_comp, -wf_ref) @@ -976,7 +988,7 @@ def test_import_state_error(): myci = "wrongobject" - with pytest.raises(ValueError, match="The supported option"): + with pytest.raises(ValueError, match="The supported objects"): _ = qchem.convert.import_state(myci) From b17773ef746651a0077b6de39c9b91b1c67f6af4 Mon Sep 17 00:00:00 2001 From: soranjh Date: Fri, 4 Aug 2023 17:57:23 -0400 Subject: [PATCH 50/60] remove import pyscf --- pennylane/qchem/convert.py | 17 +++++------------ tests/qchem/of_tests/test_convert.py | 16 +--------------- 2 files changed, 6 insertions(+), 27 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 7395617169e..c96ecfbb855 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -586,13 +586,6 @@ def import_state(solver, tol=1e-15): 0. +0.j 0. +0.j 0. +0.j 0. +0.j -0.99429698+0.j 0. +0.j 0. +0.j 0. +0.j] """ - try: - import pyscf - except ImportError as Error: - raise ImportError( - "This feature requires pyscf. It can be installed with: pip install pyscf" - ) from Error - if "RCISD" in str(solver.__str__): wf_dict = _rcisd_state(solver, tol=tol) elif "UCISD" in str(solver.__str__): @@ -644,8 +637,8 @@ def _rcisd_state(cisd_solver, tol=1e-15): 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. + cisd_solver (PySCF CISD Class instance): the class object representing the CISD + calculation in PySCF Returns: dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` @@ -655,7 +648,7 @@ def _rcisd_state(cisd_solver, tol=1e-15): **Example** >>> from pyscf import gto, scf, ci - >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g') + >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g', symm="d2h") >>> myhf = scf.RHF(mol).run() >>> myci = ci.CISD(myhf).run() >>> wf_cisd = _rcisd_state(myci, tol=1e-1) @@ -762,7 +755,7 @@ def _rccsd_state(ccsd_solver, tol=1e-15): **Example** >>> from pyscf import gto, scf, cc - >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g') + >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g', symm="d2h") >>> myhf = scf.RHF(mol).run() >>> mycc = cc.CCSD(myhf).run() >>> wf_ccsd = _rccsd_state(mycc, tol=1e-1) @@ -899,7 +892,7 @@ def _uccsd_state(ccsd_solver, tol=1e-15): **Example** >>> from pyscf import gto, scf, cc - >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g') + >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g', symm="d2h") >>> myhf = scf.UHF(mol).run() >>> mycc = cc.UCCSD(myhf).run() >>> wf_ccsd = _uccsd_state(mycc, tol=1e-1) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index df35184cdfb..2d8c00725e7 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -999,20 +999,6 @@ def test_excited_configurations_error(excitation): _ = qchem.convert._excited_configurations(2, 4, excitation) -def test_fail_import_pyscf(monkeypatch): - """Test if an ImportError is raised when pyscf is requested but not installed.""" - - mol = pyscf.gto.M(atom=[["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], basis="sto6g") - myhf = pyscf.scf.UHF(mol).run() - myci = pyscf.ci.UCISD(myhf).run() - - with monkeypatch.context() as m: - m.setitem(sys.modules, "pyscf", None) - - with pytest.raises(ImportError, match="This feature requires pyscf"): - qml.qchem.convert.import_state(myci) - - @pytest.mark.parametrize( ("molecule", "basis", "symm", "tol", "wf_ref"), [ @@ -1069,7 +1055,7 @@ def test_rcisd_state(molecule, basis, symm, tol, wf_ref): 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()))) + assert np.allclose(abs(np.array(list(wf_cisd.values()))), abs(np.array(list(wf_ref.values())))) @pytest.mark.parametrize( From c0dcbd1c2e4f0141b1a0fb0f79bc82a8e0221541 Mon Sep 17 00:00:00 2001 From: soranjh Date: Fri, 4 Aug 2023 18:00:01 -0400 Subject: [PATCH 51/60] update docstring --- pennylane/qchem/convert.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index c96ecfbb855..3bec7211a09 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -561,14 +561,14 @@ def _ucisd_state(cisd_solver, tol=1e-15): def import_state(solver, tol=1e-15): - r"""Convert an external wavefunction to a PennyLane state vector. + r"""Convert an external wavefunction to a state vector. Args: solver: external wavefunction object that will be converted - tol (float): the tolerance for discarding Slater determinants with small coefficients + tol (float): the tolerance for discarding Slater determinants based on their coefficients Raises: - ValueError: if ``method`` is not supported + ValueError: if external object type is not supported Returns: statevector (array): normalized state vector of length 2**len(number_of_spinorbitals) From 69a2b84771ab7f1ec3bd561512766e8a725403c7 Mon Sep 17 00:00:00 2001 From: soranjh Date: Sat, 5 Aug 2023 12:16:13 -0400 Subject: [PATCH 52/60] update changelog --- doc/releases/changelog-dev.md | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index ce09034bb36..4ba2334442c 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -39,6 +39,7 @@ array([False, False]) * Functions added to convert wavefunctions obtained from `PySCF` to a state vector. [(#4427)](https://github.com/PennyLaneAI/pennylane/pull/4427) + [(#4433)](https://github.com/PennyLaneAI/pennylane/pull/4433)

Improvements 🛠

From 876f671fa1ef2e5563d0cb9d871f3507c0f924eb Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 7 Aug 2023 16:10:47 -0400 Subject: [PATCH 53/60] modify docstring --- pennylane/qchem/convert.py | 98 ++++++++++++++++++++++---------------- 1 file changed, 56 insertions(+), 42 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 3bec7211a09..7dff649cfda 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -476,25 +476,26 @@ def _ucisd_state(cisd_solver, tol=1e-15): r"""Construct a wavefunction from PySCF's `UCISD` solver object. 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 that Slater + 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 such that the first number represents the + 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`. This string can be splitted to `01` and `01` for - the alpha and beta electrons. The integer corresponding to `01` is `1`. Then the dictionary - representation of a state with `0.99` contribution of the Hartree-Fock state and `0.01` - contribution from the doubly-excited state will be `{(1, 1): 0.99, (2, 2): 0.01}`. + represented by the flipped binary string `0011` which is splitted 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}`. Args: cisd_solver (PySCF UCISD Class instance): the class object representing the UCISD 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). + tol (float): the tolerance for discarding Slater determinants based on their coefficients Returns: - dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` + 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. + spin sectors, respectively, and coeff being the CI coefficients of those configurations **Example** @@ -636,14 +637,27 @@ def _wfdict_to_statevector(wf_dict, norbs): def _rcisd_state(cisd_solver, tol=1e-15): r"""Construct a wavefunction from PySCF's `RCISD` solver object. + 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 splitted 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}`. + Args: - cisd_solver (PySCF CISD Class instance): the class object representing the CISD - calculation in PySCF + cisd_solver (PySCF CISD Class instance): the class object representing the CISD calculation in PySCF + tol (float): the tolerance for discarding Slater determinants based on their coefficients Returns: - dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` + 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. + spin sectors, respectively, and coeff being the CI coefficients of those configurations **Example** @@ -719,16 +733,17 @@ def _rcisd_state(cisd_solver, tol=1e-15): def _rccsd_state(ccsd_solver, tol=1e-15): r"""Construct a wavefunction from PySCF's `CCSD` Solver object. - The wavefunction is represented as a dictionary where the keys are tuples representing a - configuration, which corresponds to a Slater determinant, and the values are the CI coefficient - corresponding to that Slater determinant. Each dictionary key is a tuple of two integers. The - binary representation of these integers correspond to a specific configuration, or Slater - determinant. The first number represents the configuration of alpha electrons and the second - number represents the beta electrons. For instance, the Hartree-Fock state - :math:`|1 1 0 0 \rangle` will be represented by the binary string `0011`. This string can be - splited to `01` and `01` for the alpha and beta electrons. The integer corresponding to `01` is - `1`. Then the dictionary representation of a state with `0.99` contribution of the Hartree-Fock - state and `0.01` contribution from the doubly-excited state will be + 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 splitted 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}`. In the current version, the exponential ansatz :math:`\exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}}` @@ -744,13 +759,12 @@ def _rccsd_state(ccsd_solver, tol=1e-15): Args: ccsd_solver (PySCF CCSD Class instance): the class object representing the CCSD 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). + 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` + dict: dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` having binary represention corresponding to the Fock occupation vector in alpha and beta - spin sectors, respectively, and coeff being the CI coefficients of those configurations. + spin sectors, respectively, and coeff being the CI coefficients of those configurations **Example** @@ -856,16 +870,17 @@ def _rccsd_state(ccsd_solver, tol=1e-15): def _uccsd_state(ccsd_solver, tol=1e-15): r"""Construct a wavefunction from PySCF's `UCCSD` Solver object. - The wavefunction is represented as a dictionary where the keys are tuples representing a - configuration, which corresponds to a Slater determinant, and the values are the CI coefficient - corresponding to that Slater determinant. Each dictionary key is a tuple of two integers. The - binary representation of these integers correspond to a specific configuration, or Slater - determinant. The first number represents the configuration of alpha electrons and the second - number represents the beta electrons. For instance, the Hartree-Fock state - :math:`|1 1 0 0 \rangle` will be represented by the binary string `0011`. This string can be - splited to `01` and `01` for the alpha and beta electrons. The integer corresponding to `01` is - `1`. Then the dictionary representation of a state with `0.99` contribution of the Hartree-Fock - state and `0.01` contribution from the doubly-excited state will be + 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 splitted 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}`. In the current version, the exponential ansatz :math:`\exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}}` @@ -881,13 +896,12 @@ def _uccsd_state(ccsd_solver, tol=1e-15): 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). + 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` + dict: dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b` having binary represention corresponding to the Fock occupation vector in alpha and beta - spin sectors, respectively, and coeff being the CI coefficients of those configurations. + spin sectors, respectively, and coeff being the CI coefficients of those configurations **Example** From a7ffcdd9d6d3503f485a37b0ea23c77671625fe1 Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 7 Aug 2023 16:26:46 -0400 Subject: [PATCH 54/60] modify docstring --- tests/qchem/of_tests/test_convert.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index 2d8c00725e7..a05465f73c2 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -918,7 +918,7 @@ def test_ucisd_state(molecule, basis, symm, tol, wf_ref): ], ) def test_wfdict_to_statevector(wf_dict, n_orbitals, string_ref, coeff_ref): - r"""Test that _wfdict_to_statevector returns the correct statevector.""" + r"""Test that _wfdict_to_statevector returns the correct state vector.""" wf_ref = np.zeros(2 ** (n_orbitals * 2)) idx_nonzero = [int(s, 2) for s in string_ref] wf_ref[idx_nonzero] = coeff_ref @@ -960,7 +960,7 @@ def test_wfdict_to_statevector(wf_dict, n_orbitals, string_ref, coeff_ref): ) @pytest.mark.parametrize("method", ["rcisd", "ucisd", "rccsd", "uccsd"]) def test_import_state(molecule, basis, symm, method, wf_ref): - r"""Test that import_state returns the correct wavefunction.""" + r"""Test that import_state returns the correct state vector.""" mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) From 8b393620d6670fe24dd0546e35e00ca2fca54c8c Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 7 Aug 2023 17:03:17 -0400 Subject: [PATCH 55/60] update changelog --- doc/releases/changelog-dev.md | 20 ++++++++++++++++++++ pennylane/qchem/convert.py | 8 ++++---- 2 files changed, 24 insertions(+), 4 deletions(-) diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index 4ba2334442c..79c26b71e7e 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -41,6 +41,26 @@ array([False, False]) [(#4427)](https://github.com/PennyLaneAI/pennylane/pull/4427) [(#4433)](https://github.com/PennyLaneAI/pennylane/pull/4433) + The `qml.qchem.import_state` function can be used to import a `PySCF` solver object and return the + corresponding state vector. + + ```pycon + >>> from pyscf import gto, scf, ci + >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g') + >>> myhf = scf.UHF(mol).run() + >>> myci = ci.UCISD(myhf).run() + >>> wf_cisd = qml.qchem.import_state(myci, tol=1e-1) + >>> print(wf_cisd) + [ 0. +0.j 0. +0.j 0. +0.j 0.1066467 +0.j + 0. +0.j 0. +0.j 0. +0.j 0. +0.j + 0. +0.j 0. +0.j 0. +0.j 0. +0.j + -0.99429698+0.j 0. +0.j 0. +0.j 0. +0.j] + ``` + + The currently supported objects are RCISD, UCISD, RCCSD, and UCCSD which correspond to + restricted (R) and unrestricted (U) configuration interaction (CI )and coupled cluster (CC) + calculations with single and double (SD) excitations. +

Improvements 🛠

* Transform Programs, `qml.transforms.core.TransformProgram`, can now be called on a batch of circuits diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 7dff649cfda..8bacb5770d8 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -580,7 +580,7 @@ def import_state(solver, tol=1e-15): >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g') >>> myhf = scf.UHF(mol).run() >>> myci = ci.UCISD(myhf).run() - >>> wf_cisd = qml.qchem.convert.import_state(myci, tol=1e-1) + >>> wf_cisd = qml.qchem.import_state(myci, tol=1e-1) >>> print(wf_cisd) [ 0. +0.j 0. +0.j 0. +0.j 0.1066467 +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j @@ -662,7 +662,7 @@ def _rcisd_state(cisd_solver, tol=1e-15): **Example** >>> from pyscf import gto, scf, ci - >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g', symm="d2h") + >>> 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) @@ -769,7 +769,7 @@ def _rccsd_state(ccsd_solver, tol=1e-15): **Example** >>> from pyscf import gto, scf, cc - >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g', symm="d2h") + >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g', symmetry="d2h") >>> myhf = scf.RHF(mol).run() >>> mycc = cc.CCSD(myhf).run() >>> wf_ccsd = _rccsd_state(mycc, tol=1e-1) @@ -906,7 +906,7 @@ def _uccsd_state(ccsd_solver, tol=1e-15): **Example** >>> from pyscf import gto, scf, cc - >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g', symm="d2h") + >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g', symmetry="d2h") >>> myhf = scf.UHF(mol).run() >>> mycc = cc.UCCSD(myhf).run() >>> wf_ccsd = _uccsd_state(mycc, tol=1e-1) From 55da678970be062b437ae0ffe0ee4634fd961a1d Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 7 Aug 2023 17:13:02 -0400 Subject: [PATCH 56/60] modify docstring --- pennylane/qchem/convert.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 8bacb5770d8..42f38f99ba4 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -430,7 +430,7 @@ def _excited_configurations(electrons, orbitals, excitation): excitation (int): number of excited electrons Returns: - tuple(array, array): arrays of excited configurations and signs obtained by reordering the + tuple(list, list): lists of excited configurations and signs obtained by reordering the creation operators **Example** @@ -572,7 +572,7 @@ def import_state(solver, tol=1e-15): ValueError: if external object type is not supported Returns: - statevector (array): normalized state vector of length 2**len(number_of_spinorbitals) + array: normalized state vector of length 2**len(number_of_spinorbitals) **Example** @@ -617,7 +617,7 @@ def _wfdict_to_statevector(wf_dict, norbs): norbs (int): number of molecular orbitals Returns: - statevector (array): normalized state vector of length 2^(number_of_spinorbitals) + array: normalized state vector of length 2^(number_of_spinorbitals) """ statevector = np.zeros(2 ** (2 * norbs), dtype=complex) From ae19005aaaee69f6505a6e7726bd1283ffc29686 Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 7 Aug 2023 17:20:48 -0400 Subject: [PATCH 57/60] update changelog --- doc/releases/changelog-dev.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index 79c26b71e7e..c9ffc96f8f2 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -37,7 +37,7 @@ def circuit(): >>> circuit(shots=1) array([False, False]) -* Functions added to convert wavefunctions obtained from `PySCF` to a state vector. +* Functions are available to obtain a state vector from `PySCF` solver objects. [(#4427)](https://github.com/PennyLaneAI/pennylane/pull/4427) [(#4433)](https://github.com/PennyLaneAI/pennylane/pull/4433) From 9f3d1d4735953d0a46f2503072ca014c440109ee Mon Sep 17 00:00:00 2001 From: soranjh Date: Mon, 7 Aug 2023 18:13:46 -0400 Subject: [PATCH 58/60] extend object support --- pennylane/qchem/convert.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 42f38f99ba4..23ed0141045 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -587,13 +587,16 @@ def import_state(solver, tol=1e-15): 0. +0.j 0. +0.j 0. +0.j 0. +0.j -0.99429698+0.j 0. +0.j 0. +0.j 0. +0.j] """ - if "RCISD" in str(solver.__str__): + + method = str(solver.__str__) + + if "CISD" in method and "UCISD" not in method: wf_dict = _rcisd_state(solver, tol=tol) - elif "UCISD" in str(solver.__str__): + elif "UCISD" in method: wf_dict = _ucisd_state(solver, tol=tol) - elif "RCCSD" in str(solver.__str__): + elif "CCSD" in method and "UCCSD" not in method: wf_dict = _rccsd_state(solver, tol=tol) - elif "UCCSD" in str(solver.__str__): + elif "UCCSD" in method: wf_dict = _uccsd_state(solver, tol=tol) else: raise ValueError( From c9f7ac30907a55193897bec863d05c151964f2ac Mon Sep 17 00:00:00 2001 From: Stepan Fomichev Date: Tue, 8 Aug 2023 09:30:26 -0400 Subject: [PATCH 59/60] added symmetry=d2h to examples in rcisd / ucisd converters --- pennylane/qchem/convert.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py index 23ed0141045..557074a69cb 100644 --- a/pennylane/qchem/convert.py +++ b/pennylane/qchem/convert.py @@ -500,7 +500,7 @@ def _ucisd_state(cisd_solver, tol=1e-15): **Example** >>> from pyscf import gto, scf, ci - >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g') + >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g', symmetry='d2h') >>> myhf = scf.UHF(mol).run() >>> myci = ci.UCISD(myhf).run() >>> wf_cisd = _ucisd_state(myci, tol=1e-1) @@ -665,7 +665,7 @@ def _rcisd_state(cisd_solver, tol=1e-15): **Example** >>> from pyscf import gto, scf, ci - >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g') + >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g', symmetry='d2h') >>> myhf = scf.RHF(mol).run() >>> myci = ci.CISD(myhf).run() >>> wf_cisd = _rcisd_state(myci, tol=1e-1) From f4d0abcb20d43968dd386fb943ceacd8053f8811 Mon Sep 17 00:00:00 2001 From: soranjh Date: Tue, 8 Aug 2023 14:27:33 -0400 Subject: [PATCH 60/60] define test data --- tests/qchem/of_tests/test_convert.py | 150 +++++++++------------------ 1 file changed, 48 insertions(+), 102 deletions(-) diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py index a05465f73c2..302c87c40f1 100644 --- a/tests/qchem/of_tests/test_convert.py +++ b/tests/qchem/of_tests/test_convert.py @@ -859,47 +859,6 @@ def test_excited_configurations(electrons, orbitals, excitation, states_ref, sig assert np.allclose(signs, signs_ref) -@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.9942969785398776, (2, 2): -0.10664669927602176}, - ), - ( - [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]], - "cc-pvdz", - "d2h", - 4e-2, - { - (1, 1): 0.9919704795977625, - (2, 2): -0.048530356564387034, - (2, 8): 0.0445233308500785, - (4, 4): -0.05003594568491194, - (8, 2): 0.04452333085007853, - (8, 8): -0.05226230322043741, - (16, 16): -0.0404759737476627, - (32, 32): -0.0404759737476627, - }, - ), - ], -) -def test_ucisd_state(molecule, basis, symm, tol, wf_ref): - r"""Test that _ucisd_state returns the correct wavefunction.""" - - mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) - myhf = pyscf.scf.UHF(mol).run() - myci = pyscf.ci.UCISD(myhf).run() - - wf_cisd = qchem.convert._ucisd_state(myci, tol=tol) - - assert wf_cisd.keys() == wf_ref.keys() - assert np.allclose(abs(np.array(list(wf_cisd.values()))), abs(np.array(list(wf_ref.values())))) - - @pytest.mark.parametrize( ("wf_dict", "n_orbitals", "string_ref", "coeff_ref"), [ # reference data were obtained manually @@ -999,32 +958,53 @@ def test_excited_configurations_error(excitation): _ = qchem.convert._excited_configurations(2, 4, excitation) +h2_molecule = [["H", (0, 0, 0)], ["H", (0, 0, 0.71)]] +h2_wf_sto6g = {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602179} # tol = 1e-1 +h2_wf_ccpvdz = { # tol = 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, +} + +li2_molecule = [["Li", (0, 0, 0)], ["Li", (0, 0, 0.71)]] +li2_wf_sto6g = { # tol = 1e-1 + (7, 7): 0.8886970081919591, + (11, 11): -0.3058459002168582, + (19, 19): -0.30584590021685887, + (35, 35): -0.14507552387854625, +} + + @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, - }, - ), + (h2_molecule, "sto6g", "d2h", 1e-1, h2_wf_sto6g), + (h2_molecule, "cc-pvdz", "d2h", 4e-2, h2_wf_ccpvdz), + ], +) +def test_ucisd_state(molecule, basis, symm, tol, wf_ref): + r"""Test that _ucisd_state returns the correct wavefunction.""" + + mol = pyscf.gto.M(atom=molecule, basis=basis, symmetry=symm) + myhf = pyscf.scf.UHF(mol).run() + myci = pyscf.ci.UCISD(myhf).run() + + wf_cisd = qchem.convert._ucisd_state(myci, tol=tol) + + assert wf_cisd.keys() == wf_ref.keys() + assert np.allclose(abs(np.array(list(wf_cisd.values()))), abs(np.array(list(wf_ref.values())))) + + +@pytest.mark.parametrize( + ("molecule", "basis", "symm", "tol", "wf_ref"), + [ + (h2_molecule, "sto6g", "d2h", 1e-1, h2_wf_sto6g), + (h2_molecule, "cc-pvdz", "d2h", 4e-2, h2_wf_ccpvdz), ( [["Be", (0, 0, 0)]], "sto6g", @@ -1061,25 +1041,8 @@ def test_rcisd_state(molecule, basis, symm, tol, wf_ref): @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.9942969030671576, (2, 2): -0.10664740292693174}, - ), - ( - [["Li", (0, 0, 0)], ["Li", (0, 0, 0.71)]], - "sto6g", - "d2h", - 1e-1, - { - (7, 7): 0.8886970081919591, - (11, 11): -0.3058459002168582, - (19, 19): -0.30584590021685887, - (35, 35): -0.14507552387854625, - }, - ), + (h2_molecule, "sto6g", "d2h", 1e-1, h2_wf_sto6g), + (li2_molecule, "sto6g", "d2h", 1e-1, li2_wf_sto6g), ], ) def test_uccsd_state(molecule, basis, symm, tol, wf_ref): @@ -1098,25 +1061,8 @@ def test_uccsd_state(molecule, basis, symm, tol, wf_ref): @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.9942969030671576, (2, 2): -0.10664740292693242}, - ), - ( - [["Li", (0, 0, 0)], ["Li", (0, 0, 0.71)]], - "sto6g", - "d2h", - 1e-1, - { - (7, 7): 0.8886969878256522, - (11, 11): -0.30584590248164206, - (19, 19): -0.30584590248164145, - (35, 35): -0.14507552651170982, - }, - ), + (h2_molecule, "sto6g", "d2h", 1e-1, h2_wf_sto6g), + (li2_molecule, "sto6g", "d2h", 1e-1, li2_wf_sto6g), ], ) def test_rccsd_state(molecule, basis, symm, tol, wf_ref):