New release v0.4.1 into main

.github/workflows/continuous_integration.yml

@@ -18,7 +18,7 @@ jobs:
       with:
         python-version: ${{ matrix.python-version }}
 
-    - name: Install pip, wheel, pytest, jupyter
+    - name: Install pip, wheel, pytest, jupyter, pyqsp
       run: |
         python -m pip install --upgrade pip
         pip install wheel
@@ -77,6 +77,12 @@ jobs:
         python -m pip install git+https://github.com/pyscf/semiempirical
       if: always()
 
+    - name: Install rdkit, openbabel-wheel
+      run: |
+        python -m pip install rdkit
+        python -m pip install openbabel-wheel
+      if: always()
+
     - name: tangelo tests
       run: |
         cd tangelo

.github/workflows/run_psi4_test.yml

@@ -31,12 +31,15 @@ jobs:
         conda init
       if: always()
 
-    - name: Install pip, wheel, pytest, jupyter
+    - name: Install pip, wheel, pytest, jupyter, openbabel, rdkit, openmm
       run: |
         python -m pip install --upgrade pip
         pip install wheel
         pip install pytest
         pip install pytest-cov
+        pip install rdkit
+        pip install openbabel-wheel
+        conda install -c conda-forge openmm
 
     - name: Install qulacs
       run: |

CHANGELOG.md

@@ -3,6 +3,32 @@
 This file documents the main changes between versions of the code.
 
 
+## [0.4.1] - 2023-10-18
+
+### Added
+
+- QM/MM problem decomposition
+- QPE framework
+- Truncated taylor series function returning qubits
+- simplify method on Circuit
+
+### Changed
+
+- Automated testing currently covering python 3.8, 3.9, 3.10, 3.11 (#333)
+- Installation: pyscf removed from requirements
+- Performance improvement: combinatorial mapping
+- Feature: ILC iteration now implements exact expansion / parameters
+- Feature: VQE-like algorithms verbose mode now prints and tracks energies, for users interested in the convergence of the algorithm.
+- Bugfix: Combinatorial mapping now handles spin != 0 (#330)
+- Bugfix: get_expectation_value takes into account n_shots for all backends supporting the option.
+- Bugfix: Fix corner case of FCI and CCSD solvers calculations (mo coefficients were occasionally recomputed differently).
+- Bugfix: Updates in IBMConnection to keep up with changes in qiskit-runtime.
+
+
+### Deprecatedv / Removed
+
+
+
 ## [0.4.0] - 2023-06-29
 
 ### Added

README.rst

@@ -124,8 +124,19 @@ The contribution process is detailed in the `contributions <./CONTRIBUTIONS.rst>
 Citations
 ---------
 
-If you use Tangelo in your research, please cite:
-
-   Valentin Senicourt, James Brown, Alexandre Fleury, Ryan Day, Erika Lloyd, Marc P. Coons, Krzysztof Bieniasz, Lee Huntington, Alejandro J. Garza, Shunji Matsuura, Rudi Plesch, Takeshi Yamazaki, and Arman Zaribafiyan Tangelo: An Open-source Python Package for End-to-end Chemistry Workflows on Quantum Computers 2022 arXiv:2206.12424
+If you use Tangelo in your research, please cite the `Tangelo release paper <https://arxiv.org/abs/2206.12424>`_ and consider mentioning Tangelo in your talks.
+
+.. code-block:: latex
+
+   @article{tangelo,
+      author = {Valentin Senicourt and James Brown and Alexandre Fleury and Ryan Day and Erika Lloyd and Marc P. Coons and Krzysztof Bieniasz and Lee Huntington and Alejandro J. Garza and Shunji Matsuura and Rudi Plesch and Takeshi Yamazaki and Arman Zaribafiyan},
+      title = {Tangelo: An Open-source Python Package for End-to-end Chemistry Workflows on Quantum Computers},
+      year = {2022},
+      url= {https://arxiv.org/abs/2206.12424},
+      number = {arXiv:2206.12424},
+      eprint = {arXiv:2206.12424},
+      publisher = {{arXiv}},
+      doi = {10.48550/arXiv.2206.12424}
+   }
 
 © Good Chemistry Company 2023. This software is released under the Apache Software License version 2.0.

tangelo/_version.py

@@ -14,4 +14,4 @@
 
 """ Define version number here. It is read in setup.py, and bumped automatically
 when using the new release Github action. """
-__version__ = "0.4.0"
+__version__ = "0.4.1"

tangelo/algorithms/classical/ccsd_solver.py

@@ -21,7 +21,7 @@
 from sympy.combinatorics.permutations import Permutation
 
 from tangelo.toolboxes.molecular_computation.molecule import SecondQuantizedMolecule
-from tangelo.toolboxes.molecular_computation import IntegralSolverPsi4, IntegralSolverPySCF
+from tangelo.toolboxes.molecular_computation import IntegralSolverPsi4, IntegralSolverPySCF, IntegralSolverPsi4QMMM
 from tangelo.algorithms.electronic_structure_solver import ElectronicStructureSolver
 from tangelo.helpers.utils import installed_chem_backends, is_package_installed
 
@@ -161,16 +161,20 @@ def __init__(self, molecule: SecondQuantizedMolecule):
                                  f"with a UHF reference in {self.__class__.__name__}")
 
         # Frozen orbitals must be declared before calling compute_mean_field to be saved in ref_wfn for Psi4 ccsd.
-        intsolve = IntegralSolverPsi4()
+        intsolve = IntegralSolverPsi4() if not hasattr(molecule.solver, "charges") else IntegralSolverPsi4QMMM(molecule.solver.combinedcharges)
         self.backend.set_options({'basis': molecule.basis, 'frozen_docc': [self.n_frozen_occ], 'frozen_uocc': [self.n_frozen_vir],
                                   'reference': self.ref})
+
         self.molecule = SecondQuantizedMolecule(xyz=molecule.xyz, q=molecule.q, spin=molecule.spin,
                                                 solver=intsolve,
                                                 basis=molecule.basis,
                                                 ecp=molecule.ecp,
                                                 symmetry=False,
                                                 uhf=molecule.uhf,
                                                 frozen_orbitals=molecule.frozen_orbitals)
+
+        self.init_mo_coeff = molecule.mo_coeff
+        self.extra_nuc_energy = 0. if not hasattr(molecule.solver, "charges") else self.molecule.solver.ext_pot.computeNuclearEnergy(self.molecule.solver.mol)
         self.basis = molecule.basis
 
     def simulate(self):
@@ -185,13 +189,15 @@ def simulate(self):
         if self.n_frozen_occ or self.n_frozen_vir:
             if not self.molecule.uhf:
                 mo_order = self.molecule.frozen_occupied + self.molecule.active_occupied + self.molecule.active_virtual + self.molecule.frozen_virtual
+                self.molecule.solver.modify_c(wfn, self.init_mo_coeff)
                 # Obtain swap operations that will take the unordered list back to ordered with the correct active space in the middle.
                 swap_ops = Permutation(mo_order).transpositions()
                 for swap_op in swap_ops:
                     wfn.Ca().rotate_columns(0, swap_op[0], swap_op[1], np.deg2rad(90))
 
             else:
-
+                # Modify Ca mo_coeff in reference wavefunction to initial mo_coeff
+                self.molecule.solver.modify_c(wfn, self.init_mo_coeff[0])
                 # Obtain swap operations that will take the unordered list back to ordered with the correct active space in the middle.
                 mo_order = (self.molecule.frozen_occupied[0] + self.molecule.active_occupied[0]
                             + self.molecule.active_virtual[0] + self.molecule.frozen_virtual[0])
@@ -200,6 +206,7 @@ def simulate(self):
                     wfn.Ca().rotate_columns(0, swap_op[0], swap_op[1], np.deg2rad(90))
 
                 # Repeat for Beta orbitals
+                self.molecule.solver.modify_c(wfn, self.init_mo_coeff[1], False)
                 mo_order_b = (self.molecule.frozen_occupied[1] + self.molecule.active_occupied[1]
                               + self.molecule.active_virtual[1] + self.molecule.frozen_virtual[1])
                 swap_ops = Permutation(mo_order_b).transpositions()
@@ -208,9 +215,14 @@ def simulate(self):
 
         self.backend.set_options({'basis': self.basis, 'frozen_docc': [self.n_frozen_occ], 'frozen_uocc': [self.n_frozen_vir],
                                   'qc_module': 'ccenergy', 'reference': self.ref})
-        energy, self.ccwfn = self.backend.energy('ccsd', molecule=self.molecule.solver.mol,
-                                                 basis=self.basis, return_wfn=True, ref_wfn=wfn)
-        return energy
+        if hasattr(self.molecule, "charges") and self.backend.__version__ >= "1.6":
+            energy, self.ccwfn = self.backend.energy('ccsd', molecule=self.molecule.solver.mol,
+                                                     basis=self.basis, return_wfn=True, ref_wfn=wfn,
+                                                     external_potentials=self.molecule.solver.external_potentials)
+        else:
+            energy, self.ccwfn = self.backend.energy('ccsd', molecule=self.molecule.solver.mol,
+                                                     basis=self.basis, return_wfn=True, ref_wfn=wfn)
+        return energy + self.extra_nuc_energy
 
     def get_rdm(self):
         """Compute the Full CI 1- and 2-particle reduced density matrices.

tangelo/algorithms/classical/fci_solver.py

@@ -174,8 +174,12 @@ def __init__(self, molecule: SecondQuantizedMolecule):
         self.backend.core.clean_variables()
         self.ciwfn = None
         self.degenerate_mo_energies = False
+        self.extra_nuc_energy = 0
         if isinstance(molecule.solver, IntegralSolverPsi4) and not molecule.symmetry:
             self.molecule = molecule
+            if hasattr(molecule.solver, "chrgfield"):
+                self.chrgfield = molecule.solver.chrgfield
+                self.extra_nuc_energy = self.molecule.solver.ext_pot.computeNuclearEnergy(self.molecule.solver.mol)
         else:
             self.degenerate_mo_energies = np.any(np.isclose(molecule.mo_energies[1:], molecule.mo_energies[:-1]))
             self.molecule = SecondQuantizedMolecule(xyz=molecule.xyz, q=molecule.q, spin=molecule.spin,
@@ -210,10 +214,20 @@ def simulate(self):
             for swap_op in swap_ops:
                 wfn.Ca().rotate_columns(0, swap_op[0], swap_op[1], np.deg2rad(90))
 
+        if hasattr(self, "chrgfield"):
+            # TODO: Remove support for older version.
+            if self.backend.__version__ < "1.6":
+                self.backend.core.set_global_option_python('EXTERN', self.chrgfield.extern)
+            else:
+                energy, self.ciwfn = self.backend.energy('fci', molecule=self.molecule.solver.mol,
+                                                         basis=self.basis, return_wfn=True,
+                                                         ref_wfn=wfn, external_potentials=self.molecule.solver.external_potentials)
+                return energy + self.extra_nuc_energy
+
         energy, self.ciwfn = self.backend.energy('fci', molecule=self.molecule.solver.mol,
                                                  basis=self.basis, return_wfn=True,
                                                  ref_wfn=wfn)
-        return energy
+        return energy + self.extra_nuc_energy
 
     def get_rdm(self):
         """Compute the Full CI 1- and 2-particle reduced density matrices.

tangelo/algorithms/classical/tests/test_ccsd_solver.py

@@ -14,9 +14,11 @@
 
 import unittest
 
+import numpy as np
+
 from tangelo import SecondQuantizedMolecule
 from tangelo.algorithms.classical.ccsd_solver import CCSDSolver, default_ccsd_solver
-from tangelo.molecule_library import mol_H2_321g, mol_Be_321g, mol_H4_sto3g_uhf_a1_frozen, xyz_H4
+from tangelo.molecule_library import mol_H2_321g, mol_Be_321g, mol_H4_sto3g_uhf_a1_frozen, xyz_H4, xyz_Be, xyz_H2O
 
 
 class CCSDSolverTest(unittest.TestCase):
@@ -59,6 +61,33 @@ def test_ccsd_be(self):
 
         self.assertAlmostEqual(energy, -14.531416, places=5)
 
+    def test_ccsd_h2o_change_mo_coeff(self):
+        """Test CCSDSolver against result from a reference implementation when rotating the molecular coefficients"""
+        mol = SecondQuantizedMolecule(xyz_H2O, q=0, spin=0, basis="sto-3g", frozen_orbitals=[6], symmetry=False)
+        mol_5 = SecondQuantizedMolecule(xyz_H2O, q=0, spin=0, basis="sto-3g", frozen_orbitals=[5], symmetry=False)
+        swap_mat = np.eye(mol.mo_coeff.shape[1])
+        swap_mat[6, 6] = 0.
+        swap_mat[5, 5] = 0.
+        swap_mat[6, 5] = 1.
+        swap_mat[5, 6] = 1.
+
+        solver = CCSDSolver(mol)
+        energy_before = solver.simulate()
+
+        solver_5 = CCSDSolver(mol_5)
+        energy_5 = solver_5.simulate()
+
+        # Swap orbital 6 and 5
+        mol.mo_coeff = mol.mo_coeff@swap_mat
+        energy = solver.simulate()
+        # Compare to SecondQunatizedMolecule with orbital 5 frozen
+        self.assertAlmostEqual(energy, energy_5, places=5)
+
+        # Swap back orbital 1 and 5
+        mol.mo_coeff = mol.mo_coeff@swap_mat
+        energy = solver.simulate()
+        self.assertAlmostEqual(energy, energy_before, places=5)
+
     def test_ccsd_be_frozen_core(self):
         """ Test CCSDSolver against result from reference implementation, with
         no mean-field provided as input. Frozen core is considered.

tangelo/algorithms/classical/tests/test_fci_solver.py

@@ -14,8 +14,13 @@
 
 import unittest
 
+import numpy as np
+from openfermion import get_sparse_operator
+
+from tangelo import SecondQuantizedMolecule
 from tangelo.algorithms import FCISolver
-from tangelo.molecule_library import mol_H2_321g, mol_Be_321g, mol_H4_cation_sto3g, mol_H4_sto3g
+from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
+from tangelo.molecule_library import mol_H2_321g, mol_Be_321g, mol_H4_cation_sto3g, mol_H4_sto3g, xyz_H4
 
 
 class FCISolverTest(unittest.TestCase):
@@ -81,6 +86,28 @@ def test_fci_H4_interior_frozen_orbitals(self):
         one_rdm, two_rdm = solver.get_rdm()
         self.assertAlmostEqual(energy, mol_H4_sto3g_freeze3.energy_from_rdms(one_rdm, two_rdm), places=5)
 
+    def test_fci_H4_modified_mo_coeff(self):
+        """ Test FCISolver against result from reference implementation with interior frozen orbitals by manually defining
+        a swap operation.
+        """
+
+        mol = SecondQuantizedMolecule(xyz_H4, frozen_orbitals=[1, 3, 4])
+
+        unitary_mat = np.array([[ 0.99799630, -0.05866044,  0.02358562,  0.00246221],
+                                [ 0.05673690,  0.99556613,  0.07416274,  0.01135272],
+                                [-0.02767097, -0.07212301,  0.99605136, -0.04375241],
+                                [-0.00431649, -0.01432819,  0.04272342,  0.99897486]])
+        # Modify orbitals
+        mol.mo_coeff = [email protected]_coeff
+
+        op = get_sparse_operator(fermion_to_qubit_mapping(mol.fermionic_hamiltonian, "JW")).toarray()
+        energy_op = np.linalg.eigh(op)[0][0]
+
+        solver = FCISolver(mol)
+        energy = solver.simulate()
+        self.assertAlmostEqual(energy, -1.8057990, places=5)
+        self.assertAlmostEqual(energy_op, -1.8057990, places=5)
+
 
 if __name__ == "__main__":
     unittest.main()

tangelo/algorithms/projective/quantum_imaginary_time.py

@@ -26,7 +26,6 @@
 from tangelo.toolboxes.operators.operators import FermionOperator, QubitOperator
 from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
 from tangelo.toolboxes.ansatz_generator._general_unitary_cc import uccgsd_generator as uccgsd_pool
-from tangelo.toolboxes.operators import qubitop_to_qubitham
 from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit
 from tangelo.linq import Circuit, get_backend
 
@@ -128,7 +127,7 @@ def build(self):
                                                 up_then_down=self.up_then_down,
                                                 spin=0)
 
-            self.qubit_hamiltonian = qubitop_to_qubitham(qubit_op, self.qubit_mapping, self.up_then_down)
+            self.qubit_hamiltonian = qubit_op
 
         # Getting the pool of operators for the ansatz. If more functionalities
         # are added, this part must be modified and generalized.
@@ -179,10 +178,8 @@ def build(self):
         for term in self.pool_operators:
             self.pool_qubit_op += term
 
-        self.qubit_operator = self.qubit_hamiltonian.to_qubitoperator()
-
         # Obtain all qubit terms that need to be measured
-        self.pool_h = [element*self.qubit_operator for element in self.pool_operators]
+        self.pool_h = [element*self.qubit_hamiltonian for element in self.pool_operators]
         self.pool_pool = [[element1*element2 for element2 in self.pool_operators] for element1 in self.pool_operators]
 
         # Obtain initial state preparation circuit
@@ -294,7 +291,7 @@ def energy_expectation(self, backend):
             float: energy computed by the backend
         """
         circuit = Circuit(n_qubits=self.final_circuit.width) if self.use_statevector else self.final_circuit
-        energy = backend.get_expectation_value(self.qubit_hamiltonian.to_qubitoperator(),
+        energy = backend.get_expectation_value(self.qubit_hamiltonian,
                                                circuit,
                                                initial_statevector=self.final_statevector)
         return energy

tangelo/algorithms/variational/iqcc_ilc_solver.py

@@ -256,21 +256,21 @@ def _update_ilc_solver(self, e_ilc):
         if self.compress_qubit_ham:
             self.ilc_ansatz.qubit_ham.frobenius_norm_compression(self.compress_eps, n_qubits)
 
-        # set dis, acs, and var_params to none to rebuild the dis & acs and initialize new parameters
-        self.ilc_ansatz.dis = None
-        self.ilc_ansatz.acs = None
-        self.ilc_ansatz.var_params = None
-        self.ilc_ansatz.build_circuit()
-
-        self.vqe_solver.qubit_hamiltonian = self.ilc_ansatz.qubit_ham
-        self.vqe_solver.initial_var_params = self.ilc_ansatz.var_params
-
         if self.iteration == self.max_ilc_iter:
             self.terminate_ilc = True
             self.final_optimal_qmf_params = optimal_qmf_var_params
             self.final_optimal_ilc_params = optimal_ilc_var_params
             if self.verbose:
                 print(f"Terminating the ILC-VQE solver after {self.max_ilc_iter} iterations.")
+        else:
+            # Set dis, acs, and var_params to None to rebuild the dis & acs and initialize new parameters
+            self.ilc_ansatz.dis = None
+            self.ilc_ansatz.acs = None
+            self.ilc_ansatz.var_params = None
+            self.ilc_ansatz.build_circuit()
+
+            self.vqe_solver.qubit_hamiltonian = self.ilc_ansatz.qubit_ham
+            self.vqe_solver.initial_var_params = self.ilc_ansatz.var_params
 
     def _update_qcc_solver(self, e_iqcc_ilc):
         """Updates after the final QCC energy evaluation with the final ILC dressed