From a956b4a23f74eb86f12505fa0c92f46d19d48d60 Mon Sep 17 00:00:00 2001
From: James Brown <james.brown@1qbit.com>
Date: Mon, 10 Jan 2022 14:13:42 -0500
Subject: [PATCH 01/10] added vsqs ansatz

---
 .github/workflows/continuous_integration.yml  |   2 +-
 .../variational/tests/test_vqe_solver.py      |  26 +++
 tangelo/algorithms/variational/vqe_solver.py  |  11 +-
 .../ansatz_generator/ansatz_utils.py          |  13 +-
 .../ansatz_generator/tests/test_vsqs.py       | 116 ++++++++++
 tangelo/toolboxes/ansatz_generator/vsqs.py    | 216 ++++++++++++++++++
 6 files changed, 379 insertions(+), 5 deletions(-)
 create mode 100644 tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
 create mode 100644 tangelo/toolboxes/ansatz_generator/vsqs.py

diff --git a/.github/workflows/continuous_integration.yml b/.github/workflows/continuous_integration.yml
index d0edfdbe6..ff7e92122 100755
--- a/.github/workflows/continuous_integration.yml
+++ b/.github/workflows/continuous_integration.yml
@@ -40,7 +40,7 @@ jobs:
         pip install qiskit==0.33.1 # Due to strange behaviour of noise model
         pip install qulacs
         pip install amazon-braket-sdk
-        pip install cirq==0.12.0
+        pip install cirq
         pip install projectq
       if: always()
 
diff --git a/tangelo/algorithms/variational/tests/test_vqe_solver.py b/tangelo/algorithms/variational/tests/test_vqe_solver.py
index c3aaedc09..e910be737 100644
--- a/tangelo/algorithms/variational/tests/test_vqe_solver.py
+++ b/tangelo/algorithms/variational/tests/test_vqe_solver.py
@@ -12,6 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+import re
 import unittest
 import numpy as np
 
@@ -137,6 +138,31 @@ def test_simulate_qcc_h2(self):
         energy = vqe_solver.simulate()
         self.assertAlmostEqual(energy, -1.137270, delta=1e-4)
 
+    def test_simulate_vsqs_h2(self):
+        """Run VQE on H2 molecule, with vsqs ansatz, JW qubit mapping, exact simulator for both molecule input and
+        qubit_hamiltonian/hini/reference_state input
+        """
+        vqe_options = {"molecule": mol_H2_sto3g, "ansatz": BuiltInAnsatze.VSQS, "qubit_mapping": "jw",
+                       "verbose": True, "ansatz_options": {"intervals": 3, "time": 3}}
+        vqe_solver = VQESolver(vqe_options)
+        vqe_solver.build()
+
+        energy = vqe_solver.simulate()
+        self.assertAlmostEqual(energy, -1.137270, delta=1e-4)
+
+        qubit_hamiltonian = vqe_solver.qubit_hamiltonian
+        hini = vqe_solver.ansatz.hini
+        reference_state = vqe_solver.ansatz.prepare_reference_state()
+
+        vqe_options = {"molecule": None, "qubit_hamiltonian": qubit_hamiltonian, "ansatz": BuiltInAnsatze.VSQS, "qubit_mapping": "jw",
+                       "ansatz_options": {"intervals": 3, "time": 3, "qubit_hamiltonian": qubit_hamiltonian,
+                       "hini": hini, "reference_state": reference_state}}
+        vqe_solver = VQESolver(vqe_options)
+        vqe_solver.build()
+
+        energy = vqe_solver.simulate()
+        self.assertAlmostEqual(energy, -1.137270, delta=1e-4)
+
     def test_simulate_h2_qiskit(self):
         """Run VQE on H2 molecule, with UCCSD ansatz, JW qubit mapping, initial
         parameters, exact qiskit simulator.
diff --git a/tangelo/algorithms/variational/vqe_solver.py b/tangelo/algorithms/variational/vqe_solver.py
index 639e3bf3c..e5fe56729 100644
--- a/tangelo/algorithms/variational/vqe_solver.py
+++ b/tangelo/algorithms/variational/vqe_solver.py
@@ -35,6 +35,7 @@
 from tangelo.toolboxes.ansatz_generator.upccgsd import UpCCGSD
 from tangelo.toolboxes.ansatz_generator.qmf import QMF
 from tangelo.toolboxes.ansatz_generator.qcc import QCC
+from tangelo.toolboxes.ansatz_generator.vsqs import VSQS
 from tangelo.toolboxes.ansatz_generator.variational_circuit import VariationalCircuitAnsatz
 from tangelo.toolboxes.ansatz_generator.penalty_terms import combined_penalty
 from tangelo.toolboxes.post_processing.bootstrapping import get_resampled_frequencies
@@ -50,6 +51,7 @@ class BuiltInAnsatze(Enum):
     UpCCGSD = 4
     QMF = 5
     QCC = 6
+    VSQS = 7
 
 
 class VQESolver:
@@ -186,13 +188,20 @@ def build(self):
                     self.ansatz = QMF(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options)
                 elif self.ansatz == BuiltInAnsatze.QCC:
                     self.ansatz = QCC(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options)
+                elif self.ansatz == BuiltInAnsatze.VSQS:
+                    self.ansatz = VSQS(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options)
                 else:
                     raise ValueError(f"Unsupported ansatz. Built-in ansatze:\n\t{self.builtin_ansatze}")
             elif not isinstance(self.ansatz, Ansatz):
                 raise TypeError(f"Invalid ansatz dataype. Expecting instance of Ansatz class, or one of built-in options:\n\t{self.builtin_ansatze}")
 
         # Building with a qubit Hamiltonian.
-        elif not isinstance(self.ansatz, Ansatz):
+        elif not isinstance(self.ansatz, Ansatz) or self.ansatz not in [BuiltInAnsatze.HEA, BuiltInAnsatze.VSQS]:
+            if self.ansatz == BuiltInAnsatze.HEA:
+                self.ansatz = HEA(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options)
+            elif self.ansatz == BuiltInAnsatze.VSQS:
+                self.ansatz = VSQS(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options)
+        else:
             raise TypeError(f"Invalid ansatz dataype. Expecting a custom Ansatz (Ansatz class).")
 
         # Set ansatz initial parameters (default or use input), build corresponding ansatz circuit
diff --git a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
index a9a2a4c16..68cb31630 100644
--- a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
+++ b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
@@ -83,7 +83,7 @@ def exp_pauliword_to_gates(pauli_word, coef, variational=True, control=None):
 
 
 def get_exponentiated_qubit_operator_circuit(qubit_op, time=1., variational=False, trotter_order=1, control=None,
-                                             return_phase=False):
+                                             return_phase=False, pauli_order=None):
     """Generate the exponentiation of a qubit operator in first- or second-order Trotterized form.
     The algorithm is described in Whitfield 2010 https://arxiv.org/pdf/1001.3855.pdf
 
@@ -94,12 +94,19 @@ def get_exponentiated_qubit_operator_circuit(qubit_op, time=1., variational=Fals
         variational (bool) : Whether the coefficients are variational
         trotter_order (int): order of trotter approximation, only 1 or 2 are supported.
         return_phase (bool): Return the global-phase generated
+        pauli_order (list): The desired pauli_word order defined as a list with elements (pauli_word, coeff)
+            with corresponding dictionary elements pauli_word: coeff in QubitOperator terms.items()
 
     Returns:
         Circuit: circuit corresponding to exponentiation of qubit operator
         phase : The global phase of the time evolution if return_phase=True else not included
     """
-    pauli_words = list(qubit_op.terms.items())
+    if pauli_order is None:
+        pauli_words = list(qubit_op.terms.items())
+    elif isinstance(pauli_order, list):
+        pauli_words = deepcopy(pauli_order)
+    else:
+        raise ValueError("ordered terms must be a list with elements (keys, values) of qubit_op.terms.items()")
 
     if trotter_order > 2:
         raise ValueError(f"Trotter order of >2 is not supported currently in Tangelo.")
@@ -118,7 +125,7 @@ def get_exponentiated_qubit_operator_circuit(qubit_op, time=1., variational=Fals
     phase = 1.
     exp_pauli_word_gates = list()
     for i in range(trotter_order):
-        if i == 0:
+        if i == 1:
             pauli_words.reverse()
         for pauli_word, coef in pauli_words:
             if pauli_word:  # identity terms do not contribute to evolution outside of a phase
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
new file mode 100644
index 000000000..d53e6bf02
--- /dev/null
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
@@ -0,0 +1,116 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+import numpy as np
+import os
+from openfermion import load_operator
+
+from tangelo.molecule_library import mol_H2_sto3g
+from tangelo.toolboxes.operators.operators import QubitOperator
+from tangelo.toolboxes.qubit_mappings import jordan_wigner, symmetry_conserving_bravyi_kitaev
+from tangelo.toolboxes.ansatz_generator.vsqs import VSQS
+from tangelo.linq import Simulator
+from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit
+
+# For openfermion.load_operator function.
+pwd_this_test = os.path.dirname(os.path.abspath(__file__))
+
+
+class VSQSTest(unittest.TestCase):
+
+    def test_vsqs_set_var_params(self):
+        """Verify behavior of set_var_params for different inputs (list, numpy array).
+        MP2 have their own tests.
+        """
+
+        vsqs_ansatz = VSQS(mol_H2_sto3g)
+
+        two_ones = np.ones((2,))
+
+        vsqs_ansatz.set_var_params([1., 1.])
+        np.testing.assert_array_almost_equal(vsqs_ansatz.var_params, two_ones, decimal=6)
+
+        vsqs_ansatz.set_var_params(np.array([1., 1.]))
+        np.testing.assert_array_almost_equal(vsqs_ansatz.var_params, two_ones, decimal=6)
+
+    def test_vsqs_incorrect_number_var_params(self):
+        """Return an error if user provide incorrect number of variational
+        parameters.
+        """
+
+        vsqs_ansatz = VSQS(mol_H2_sto3g)
+
+        self.assertRaises(ValueError, vsqs_ansatz.set_var_params, np.array([1., 1., 1., 1.]))
+
+    def test_vsqs_H2(self):
+        """Verify closed-shell VSQS functionalities for H2."""
+
+        # Build circuit
+        vsqs_ansatz = VSQS(mol_H2_sto3g, intervals=3, time=3)
+        vsqs_ansatz.build_circuit()
+
+        # Build qubit hamiltonian for energy evaluation
+        qubit_hamiltonian = jordan_wigner(mol_H2_sto3g.fermionic_hamiltonian)
+
+        # Assert energy returned is as expected for given parameters
+        sim = Simulator()
+        vsqs_ansatz.update_var_params([0.66666667, 0.9698286, 0.21132472, 0.6465473])
+        energy = sim.get_expectation_value(qubit_hamiltonian, vsqs_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.1372701255155757, delta=1e-6)
+
+    def test_vsqs_H4_doublecation(self):
+        """Verify closed-shell VSQS functionalities for H4 2+ by using saved qubit hamiltonian and initial hamiltonian"""
+
+        var_params = [-2.53957674, 0.72683888, 1.08799500, 0.49836183,
+                      -0.23020698, 0.93278630, 0.50591026, 0.50486903]
+
+        # Build qubit hamiltonian for energy evaluation
+        qubit_hamiltonian = load_operator("mol_H4_doublecation_minao_qubitham_jw_b.data", data_directory=pwd_this_test+"/data", plain_text=True)
+        initial_hamiltonian = load_operator("mol_H4_doublecation_minao_init_qubitham_jw_b.data", data_directory=pwd_this_test+"/data", plain_text=True)
+        reference_state = get_reference_circuit(8, 2, "jw", up_then_down=True, spin=0)
+
+        # Build circuit
+        vsqs_ansatz = VSQS(qubit_hamiltonian=qubit_hamiltonian, hini=initial_hamiltonian, reference_state=reference_state,
+                           intervals=5, time=5, trotter_order=2)
+        vsqs_ansatz.build_circuit()
+
+        # Assert energy returned is as expected for given parameters
+        sim = Simulator()
+        vsqs_ansatz.update_var_params(var_params)
+        energy = sim.get_expectation_value(qubit_hamiltonian, vsqs_ansatz.circuit)
+        self.assertAlmostEqual(energy, -0.85425, delta=1e-4)
+
+    def test_vsqs_H2_with_hnav(self):
+        """Verify closed-shell VSQS functionalities for H2 with navigator hamiltonian"""
+        navigator_hamiltonian = (QubitOperator('X0 Y1', 0.03632537110234512) + QubitOperator('Y0 X1', 0.03632537110234512)
+                                 + QubitOperator('Y0', 2.e-5) + QubitOperator('Y1', 2.e-5))
+
+        # Build qubit hamiltonian for energy evaluation
+        qubit_hamiltonian = symmetry_conserving_bravyi_kitaev(mol_H2_sto3g.fermionic_hamiltonian, 4, 2, False, 0)
+
+        # Build circuit
+        vsqs_ansatz = VSQS(mol_H2_sto3g, intervals=2, time=1, mapping='scbk', up_then_down=True, trotter_order=2,
+                           hnav=navigator_hamiltonian)
+        vsqs_ansatz.build_circuit()
+
+        # Assert energy returned is as expected for given parameters
+        sim = Simulator()
+        vsqs_ansatz.update_var_params([0.50000001, -0.02494214, 3.15398767])
+        energy = sim.get_expectation_value(qubit_hamiltonian, vsqs_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.1372701255155757, delta=1e-6)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tangelo/toolboxes/ansatz_generator/vsqs.py b/tangelo/toolboxes/ansatz_generator/vsqs.py
new file mode 100644
index 000000000..54bb492dd
--- /dev/null
+++ b/tangelo/toolboxes/ansatz_generator/vsqs.py
@@ -0,0 +1,216 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""This module defines the Variationally Scheduled Quantum Simulation class. It provides an
+Adiabatic State Preparation (ASP) inspired ansatz as described in https://arxiv.org/abs/2003.09913."""
+
+from copy import deepcopy
+
+import numpy as np
+from openfermion import QubitOperator as ofQubitOperator
+
+from .ansatz import Ansatz
+from .ansatz_utils import get_exponentiated_qubit_operator_circuit
+from tangelo.toolboxes.operators import FermionOperator
+from tangelo.linq import Circuit
+from tangelo.toolboxes.qubit_mappings.mapping_transform import get_qubit_number, fermion_to_qubit_mapping
+from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit
+
+
+class VSQS(Ansatz):
+    """This class implements the Variationally Scheduled Quantum Simulator (VSQS) for state preparation as described in
+    https://arxiv.org/abs/2003.09913
+
+    Must supply either a molecule or a qubit_hamiltonian. If supplying a qubit_hamiltonian, must also supply
+    a reference_state Circuit and a hini QubitOperator.
+
+    Args:
+        molecule (SecondQuantizedMolecule): The molecular system. Default: None
+        mapping (str): One of the support fermion to qubit mappings. Default : "JW"
+        up_then_down (bool): change basis ordering putting all spin up orbitals first, followed by all spin down.
+            Default: False (alternating spin up/down ordering)
+        intervals (int): The number of steps in the VSQS process. Must be greater than 1. Default: 2
+        time (float): The propagation time. Default: 1
+        qubit_hamiltonian (QubitOperator): The qubit hamiltonian to evolve. Default: None
+        reference_state (Circuit): The reference state for the propagation as defined by a Circuit. Mandatory if supplying
+            a qubit_hamiltonian. Default: None
+        hini (QubitOperator): The initial qubit Hamiltonian that one corresponds to the reference state. Mandatory if supplying
+            a qubit_hamiltonian. Default: None
+        hnav (QubitOperator): The navigator Hamiltonian. Default: None
+    """
+
+    def __init__(self, molecule=None, mapping="jw", up_then_down=False, intervals=2, time=1, qubit_hamiltonian=None,
+                 reference_state=None, hini=None, hnav=None, trotter_order=1):
+
+        self.up_then_down = up_then_down
+        self.mapping = mapping
+        if intervals > 1:
+            self.intervals = intervals
+        else:
+            raise ValueError("The number of intervals must be greater than 1.")
+        self.time = time
+        self.dt = self.time/self.intervals
+        self.reference_state = reference_state
+        self.trotter_order = trotter_order
+
+        if molecule is None:
+            self.qubit_hamiltonian = qubit_hamiltonian
+            if not isinstance(hini, ofQubitOperator):
+                raise ValueError("When providing a qubit hamiltonian, an initial Hamiltonian must also be provided")
+            self.hini = hini
+            if not isinstance(reference_state, Circuit):
+                raise ValueError("Reference state must be provided when simulating a qubit hamiltonian directly")
+            self.reference_state = reference_state
+        else:
+            self.n_electrons = molecule.n_active_electrons
+            self.n_spinorbitals = int(molecule.n_sos)
+            self.n_qubits = get_qubit_number(mapping, self.n_spinorbitals)
+            self.spin = molecule.spin
+            self.qubit_hamiltonian = fermion_to_qubit_mapping(molecule.fermionic_hamiltonian, n_spinorbitals=self.n_spinorbitals,
+                                                              n_electrons=self.n_electrons, mapping=self.mapping,
+                                                              up_then_down=self.up_then_down, spin=self.spin)
+            self.hini = self.build_hini(molecule) if hini is None else hini
+        self.hfin = self.qubit_hamiltonian
+        self.hnav = hnav
+
+        self.hfin_list, self.n_hfin = self.remove_constant_and_count_terms(self.qubit_hamiltonian)
+        self.hini_list, self.n_hini = self.remove_constant_and_count_terms(self.hini)
+        if self.hnav is None:
+            self.stride = 2
+            self.n_hnav = 0
+        else:
+            if isinstance(self.hnav, ofQubitOperator):
+                self.stride = 3
+                self.hnav_list, self.n_hnav = self.remove_constant_and_count_terms(self.hnav)
+            else:
+                raise ValueError("Navigator Hamiltonian must be a QubitOperator")
+
+        self.n_var_params = (intervals - 1) * self.stride
+        self.n_var_gates = (self.n_hini + self.n_hfin + self.n_hnav) * self.trotter_order
+
+        self.var_params = None
+        self.circuit = None
+
+    def remove_constant_and_count_terms(self, qu_op):
+        """count of non-zero terms in a QubitOperator and return the list of terms.items()"""
+        count = 0
+        new_qu_op = deepcopy(qu_op)
+        for term in new_qu_op.terms.keys():
+            if term:
+                count += 1
+            else:
+                new_qu_op.terms[term] = 0
+        new_qu_op.compress()
+        return list(new_qu_op.terms.items()), count
+
+    def build_hini(self, molecule):
+        """Return the initial Hamiltonian (hini) composed of the one-body terms derived from the diagonal of Fock
+        matrix and one-body off-diagonal terms"""
+        fock = molecule.mean_field.mo_coeff.T @ molecule.mean_field.get_fock() @ molecule.mean_field.mo_coeff
+        h1 = molecule.mean_field.mo_coeff.T @ molecule.mean_field.get_hcore() @ molecule.mean_field.mo_coeff
+        hf_ferm = FermionOperator()
+        for i in range(self.n_spinorbitals//2):
+            for j in range(self.n_spinorbitals//2):
+                if i != j:
+                    hf_ferm += FermionOperator(((i * 2, 1), (j * 2, 0)), h1[i, j])
+                    hf_ferm += FermionOperator(((i * 2 + 1, 1), (j * 2 + 1, 0)), h1[i, j])
+                else:
+                    hf_ferm += FermionOperator(((i * 2, 1), (j * 2, 0)), fock[i, j])
+                    hf_ferm += FermionOperator(((i * 2 + 1, 1), (j * 2 + 1, 0)), fock[i, j])
+        return fermion_to_qubit_mapping(hf_ferm, mapping=self.mapping, n_spinorbitals=self.n_spinorbitals, n_electrons=self.n_electrons,
+                                        up_then_down=self.up_then_down, spin=self.spin)
+
+    def set_var_params(self, var_params=None):
+        """Set values for the variational parameters. Default is linear interpolation."""
+        if var_params is None:
+            var_params = self.init_params()[self.stride:self.n_var_params+self.stride]
+
+        init_var_params = np.array(var_params)
+        if init_var_params.size == self.n_var_params:
+            self.var_params = var_params
+        else:
+            raise ValueError(f"Expected {self.n_var_params} variational parameters but received {init_var_params.size}.")
+        return var_params
+
+    def update_var_params(self, var_params):
+        """Update the variational parameters in the circuit without rebuilding."""
+        for i in range(self.intervals-1):
+            self._update_gate_params_for_qu_op(self.hini_list, self.n_var_gates * i, var_params[self.stride*i], self.n_hini)
+            self._update_gate_params_for_qu_op(self.hfin_list, self.n_var_gates * i + self.n_hini * self.trotter_order,
+                                               var_params[self.stride*i+1], self.n_hfin)
+            if self.hnav is not None:
+                self._update_gate_params_for_qu_op(self.hnav_list, self.n_var_gates * i + (self.n_hini + self.n_hfin) * self.trotter_order,
+                                                   var_params[self.stride*i+2], self.n_hnav)
+
+    def _update_gate_params_for_qu_op(self, qu_op, n_var_start, var_param, num_terms):
+        prefac = 2 / self.trotter_order * self.dt * var_param
+        for i, (_, coeff) in enumerate(qu_op):
+            self.circuit._variational_gates[n_var_start+i].parameter = prefac * coeff
+        if self.trotter_order == 2:
+            for i, (_, coeff) in enumerate(list(reversed(qu_op))):
+                self.circuit._variational_gates[n_var_start+num_terms+i].parameter = prefac * coeff
+
+    def init_params(self):
+        """Generate initial parameters for the VSQS algorithm.
+        a_i = step*i, b_i=1-step*i, c_i= 1-step*i i<n_intervals/2, step*i i>n_intervals/2"""
+        a = np.zeros(self.intervals+1)
+        b = np.zeros(self.intervals+1)
+        c = np.zeros(self.intervals+1)
+        a[0] = 1
+        b[self.intervals] = 1
+        step = 1/self.intervals
+        for i in range(1, self.intervals):
+            a[i] = (1 - i * step)
+            b[i] = (i * step)
+        all_params = np.zeros(self.stride * (self.intervals + 1))
+        if self.hnav is None:
+            for i in range(self.intervals + 1):
+                all_params[2 * i] = a[i]
+                all_params[2 * i + 1] = b[i]
+        else:
+            c[0:self.intervals//2] = b[0:self.intervals//2]
+            c[self.intervals//2:self.intervals+1] = a[self.intervals//2:self.intervals+1]
+            for i in range(self.intervals + 1):
+                all_params[3 * i] = a[i]
+                all_params[3 * i + 1] = b[i]
+                all_params[3 * i + 2] = c[i]
+        return all_params
+
+    def prepare_reference_state(self):
+        """Prepare a circuit generating the HF reference state."""
+        return get_reference_circuit(n_spinorbitals=self.n_spinorbitals, n_electrons=self.n_electrons, mapping=self.mapping,
+                                     up_then_down=self.up_then_down, spin=self.spin)
+
+    def build_circuit(self, var_params=None):
+        """Build the VSQS circuit by successive first-order trotterizations of hini, hfin and possibly hnav"""
+        reference_state_circuit = self.prepare_reference_state() if self.reference_state is None else self.reference_state
+        self.var_params = self.set_var_params(var_params)
+
+        vsqs_circuit = get_exponentiated_qubit_operator_circuit(self.hini, time=self.dt, trotter_order=self.trotter_order, pauli_order=self.hini_list)
+        for i in range(self.intervals-1):
+            vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hini, time=self.var_params[i * self.stride] * self.dt, variational=True,
+                                                                     trotter_order=self.trotter_order, pauli_order=self.hini_list)
+            vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hfin, time=self.var_params[i * self.stride + 1] * self.dt, variational=True,
+                                                                     trotter_order=self.trotter_order, pauli_order=self.hfin_list)
+            if self.hnav is not None:
+                vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hnav, time=self.var_params[i * self.stride + 2] * self.dt, variational=True,
+                                                                         trotter_order=self.trotter_order, pauli_order=self.hnav_list)
+        vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hfin, time=self.dt, trotter_order=self.trotter_order, pauli_order=self.hfin_list)
+
+        if reference_state_circuit.size != 0:
+            self.circuit = reference_state_circuit + vsqs_circuit
+        else:
+            self.circuit = vsqs_circuit
+
+        return self.circuit

From f400be6bed49ab459db0c4b0a1a0cb06476aa364 Mon Sep 17 00:00:00 2001
From: James Brown <james.brown@1qbit.com>
Date: Mon, 10 Jan 2022 14:20:34 -0500
Subject: [PATCH 02/10] underscores in front of internal functions

---
 tangelo/toolboxes/ansatz_generator/vsqs.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/tangelo/toolboxes/ansatz_generator/vsqs.py b/tangelo/toolboxes/ansatz_generator/vsqs.py
index 54bb492dd..8be5ff7bb 100644
--- a/tangelo/toolboxes/ansatz_generator/vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/vsqs.py
@@ -80,19 +80,19 @@ def __init__(self, molecule=None, mapping="jw", up_then_down=False, intervals=2,
             self.qubit_hamiltonian = fermion_to_qubit_mapping(molecule.fermionic_hamiltonian, n_spinorbitals=self.n_spinorbitals,
                                                               n_electrons=self.n_electrons, mapping=self.mapping,
                                                               up_then_down=self.up_then_down, spin=self.spin)
-            self.hini = self.build_hini(molecule) if hini is None else hini
+            self.hini = self._build_hini(molecule) if hini is None else hini
         self.hfin = self.qubit_hamiltonian
         self.hnav = hnav
 
-        self.hfin_list, self.n_hfin = self.remove_constant_and_count_terms(self.qubit_hamiltonian)
-        self.hini_list, self.n_hini = self.remove_constant_and_count_terms(self.hini)
+        self.hfin_list, self.n_hfin = self._remove_constant_and_count_terms(self.qubit_hamiltonian)
+        self.hini_list, self.n_hini = self._remove_constant_and_count_terms(self.hini)
         if self.hnav is None:
             self.stride = 2
             self.n_hnav = 0
         else:
             if isinstance(self.hnav, ofQubitOperator):
                 self.stride = 3
-                self.hnav_list, self.n_hnav = self.remove_constant_and_count_terms(self.hnav)
+                self.hnav_list, self.n_hnav = self._remove_constant_and_count_terms(self.hnav)
             else:
                 raise ValueError("Navigator Hamiltonian must be a QubitOperator")
 
@@ -102,7 +102,7 @@ def __init__(self, molecule=None, mapping="jw", up_then_down=False, intervals=2,
         self.var_params = None
         self.circuit = None
 
-    def remove_constant_and_count_terms(self, qu_op):
+    def _remove_constant_and_count_terms(self, qu_op):
         """count of non-zero terms in a QubitOperator and return the list of terms.items()"""
         count = 0
         new_qu_op = deepcopy(qu_op)
@@ -114,7 +114,7 @@ def remove_constant_and_count_terms(self, qu_op):
         new_qu_op.compress()
         return list(new_qu_op.terms.items()), count
 
-    def build_hini(self, molecule):
+    def _build_hini(self, molecule):
         """Return the initial Hamiltonian (hini) composed of the one-body terms derived from the diagonal of Fock
         matrix and one-body off-diagonal terms"""
         fock = molecule.mean_field.mo_coeff.T @ molecule.mean_field.get_fock() @ molecule.mean_field.mo_coeff
@@ -134,7 +134,7 @@ def build_hini(self, molecule):
     def set_var_params(self, var_params=None):
         """Set values for the variational parameters. Default is linear interpolation."""
         if var_params is None:
-            var_params = self.init_params()[self.stride:self.n_var_params+self.stride]
+            var_params = self._init_params()[self.stride:self.n_var_params+self.stride]
 
         init_var_params = np.array(var_params)
         if init_var_params.size == self.n_var_params:
@@ -161,7 +161,7 @@ def _update_gate_params_for_qu_op(self, qu_op, n_var_start, var_param, num_terms
             for i, (_, coeff) in enumerate(list(reversed(qu_op))):
                 self.circuit._variational_gates[n_var_start+num_terms+i].parameter = prefac * coeff
 
-    def init_params(self):
+    def _init_params(self):
         """Generate initial parameters for the VSQS algorithm.
         a_i = step*i, b_i=1-step*i, c_i= 1-step*i i<n_intervals/2, step*i i>n_intervals/2"""
         a = np.zeros(self.intervals+1)

From 2ff60f5588bcbba3d846735fbb9850f9f0b04535 Mon Sep 17 00:00:00 2001
From: James Brown <james.brown@1qbit.com>
Date: Mon, 10 Jan 2022 14:23:45 -0500
Subject: [PATCH 03/10] fixed description

---
 tangelo/toolboxes/ansatz_generator/vsqs.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tangelo/toolboxes/ansatz_generator/vsqs.py b/tangelo/toolboxes/ansatz_generator/vsqs.py
index 8be5ff7bb..87595cdc4 100644
--- a/tangelo/toolboxes/ansatz_generator/vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/vsqs.py
@@ -193,7 +193,7 @@ def prepare_reference_state(self):
                                      up_then_down=self.up_then_down, spin=self.spin)
 
     def build_circuit(self, var_params=None):
-        """Build the VSQS circuit by successive first-order trotterizations of hini, hfin and possibly hnav"""
+        """Build the VSQS circuit by successive first- or second-order trotterizations of hini, hfin and possibly hnav"""
         reference_state_circuit = self.prepare_reference_state() if self.reference_state is None else self.reference_state
         self.var_params = self.set_var_params(var_params)
 

From 1e9bda68eecda0711cf1af0d7ba30682f1a00ff0 Mon Sep 17 00:00:00 2001
From: James Brown <james.brown@1qbit.com>
Date: Mon, 10 Jan 2022 15:01:13 -0500
Subject: [PATCH 04/10] fix for failing tests

---
 tangelo/algorithms/variational/vqe_solver.py  |   2 +-
 ...doublecation_minao_init_qubitham_jw_b.data |  18 ++
 ...l_H4_doublecation_minao_qubitham_jw_b.data | 186 ++++++++++++++++++
 3 files changed, 205 insertions(+), 1 deletion(-)
 create mode 100644 tangelo/toolboxes/ansatz_generator/tests/data/mol_H4_doublecation_minao_init_qubitham_jw_b.data
 create mode 100644 tangelo/toolboxes/ansatz_generator/tests/data/mol_H4_doublecation_minao_qubitham_jw_b.data

diff --git a/tangelo/algorithms/variational/vqe_solver.py b/tangelo/algorithms/variational/vqe_solver.py
index e5fe56729..d3ca91cd5 100644
--- a/tangelo/algorithms/variational/vqe_solver.py
+++ b/tangelo/algorithms/variational/vqe_solver.py
@@ -196,7 +196,7 @@ def build(self):
                 raise TypeError(f"Invalid ansatz dataype. Expecting instance of Ansatz class, or one of built-in options:\n\t{self.builtin_ansatze}")
 
         # Building with a qubit Hamiltonian.
-        elif not isinstance(self.ansatz, Ansatz) or self.ansatz not in [BuiltInAnsatze.HEA, BuiltInAnsatze.VSQS]:
+        elif not isinstance(self.ansatz, Ansatz) and self.ansatz in [BuiltInAnsatze.HEA, BuiltInAnsatze.VSQS]:
             if self.ansatz == BuiltInAnsatze.HEA:
                 self.ansatz = HEA(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options)
             elif self.ansatz == BuiltInAnsatze.VSQS:
diff --git a/tangelo/toolboxes/ansatz_generator/tests/data/mol_H4_doublecation_minao_init_qubitham_jw_b.data b/tangelo/toolboxes/ansatz_generator/tests/data/mol_H4_doublecation_minao_init_qubitham_jw_b.data
new file mode 100644
index 000000000..3de6e4b08
--- /dev/null
+++ b/tangelo/toolboxes/ansatz_generator/tests/data/mol_H4_doublecation_minao_init_qubitham_jw_b.data
@@ -0,0 +1,18 @@
+QubitOperator:
+(-3.201203136341336+0j) [] +
+(0.018922220448748025+0j) [X0 X1] +
+(0.018922220448748025+0j) [Y0 Y1] +
+(0.6887270305253504+0j) [Z0] +
+(0.37691668685924784+0j) [Z1] +
+(-0.02281117974146273+0j) [X2 X3] +
+(-0.02281117974146273+0j) [Y2 Y3] +
+(0.3753113309246301+0j) [Z2] +
+(0.15964651986143952+0j) [Z3] +
+(0.018922220448748025+0j) [X4 X5] +
+(0.018922220448748025+0j) [Y4 Y5] +
+(0.6887270305253504+0j) [Z4] +
+(0.37691668685924784+0j) [Z5] +
+(-0.02281117974146273+0j) [X6 X7] +
+(-0.02281117974146273+0j) [Y6 Y7] +
+(0.3753113309246301+0j) [Z6] +
+(0.15964651986143952+0j) [Z7]
\ No newline at end of file
diff --git a/tangelo/toolboxes/ansatz_generator/tests/data/mol_H4_doublecation_minao_qubitham_jw_b.data b/tangelo/toolboxes/ansatz_generator/tests/data/mol_H4_doublecation_minao_qubitham_jw_b.data
new file mode 100644
index 000000000..a602becad
--- /dev/null
+++ b/tangelo/toolboxes/ansatz_generator/tests/data/mol_H4_doublecation_minao_qubitham_jw_b.data
@@ -0,0 +1,186 @@
+QubitOperator:
+(-0.845088084959612+0j) [] +
+(-0.0020776113200212075+0j) [X0 X1] +
+(-0.011761015389060145+0j) [X0 X1 X2 X3] +
+(-0.0003365342989052316+0j) [X0 X1 Y2 Y3] +
+(0.002844217831601514+0j) [X0 X1 Z2] +
+(0.005067140284950715+0j) [X0 X1 Z3] +
+(0.02660445410908035+0j) [X0 X1 X4 X5] +
+(0.02660445410908035+0j) [X0 X1 Y4 Y5] +
+(0.009461086978084416+0j) [X0 X1 Z4] +
+(-0.002045316086697597+0j) [X0 X1 Z5] +
+(-0.024395190370696845+0j) [X0 X1 X6 X7] +
+(-0.024395190370696845+0j) [X0 X1 Y6 Y7] +
+(-0.007801932778402895+0j) [X0 X1 Z6] +
+(0.013474635539233077+0j) [X0 X1 Z7] +
+(-0.011424481090154913+0j) [X0 Y1 Y2 X3] +
+(0.024934767638066777+0j) [X0 Z1 X2 X4 Z5 X6] +
+(0.006238164536781785+0j) [X0 Z1 X2 X4 Z5 Z6 X7] +
+(0.024934767638066777+0j) [X0 Z1 X2 Y4 Z5 Y6] +
+(0.006238164536781785+0j) [X0 Z1 X2 Y4 Z5 Z6 Y7] +
+(0.01064615061000441+0j) [X0 Z1 X2 X5 X6] +
+(-0.02405865607179162+0j) [X0 Z1 X2 X5 Z6 X7] +
+(0.01064615061000441+0j) [X0 Z1 X2 Y5 Y6] +
+(-0.02405865607179162+0j) [X0 Z1 X2 Y5 Z6 Y7] +
+(0.006238164536781785+0j) [X0 Z1 Z2 X3 X4 Z5 X6] +
+(0.017676042848046+0j) [X0 Z1 Z2 X3 X4 Z5 Z6 X7] +
+(0.006238164536781785+0j) [X0 Z1 Z2 X3 Y4 Z5 Y6] +
+(0.017676042848046+0j) [X0 Z1 Z2 X3 Y4 Z5 Z6 Y7] +
+(-0.012634174981636702+0j) [X0 Z1 Z2 X3 X5 X6] +
+(-0.008407495254282365+0j) [X0 Z1 Z2 X3 X5 Z6 X7] +
+(-0.012634174981636702+0j) [X0 Z1 Z2 X3 Y5 Y6] +
+(-0.008407495254282365+0j) [X0 Z1 Z2 X3 Y5 Z6 Y7] +
+(-0.011424481090154913+0j) [Y0 X1 X2 Y3] +
+(-0.0020776113200212075+0j) [Y0 Y1] +
+(-0.0003365342989052316+0j) [Y0 Y1 X2 X3] +
+(-0.011761015389060145+0j) [Y0 Y1 Y2 Y3] +
+(0.002844217831601514+0j) [Y0 Y1 Z2] +
+(0.005067140284950715+0j) [Y0 Y1 Z3] +
+(0.02660445410908035+0j) [Y0 Y1 X4 X5] +
+(0.02660445410908035+0j) [Y0 Y1 Y4 Y5] +
+(0.009461086978084416+0j) [Y0 Y1 Z4] +
+(-0.002045316086697597+0j) [Y0 Y1 Z5] +
+(-0.024395190370696845+0j) [Y0 Y1 X6 X7] +
+(-0.024395190370696845+0j) [Y0 Y1 Y6 Y7] +
+(-0.007801932778402895+0j) [Y0 Y1 Z6] +
+(0.013474635539233077+0j) [Y0 Y1 Z7] +
+(0.024934767638066777+0j) [Y0 Z1 Y2 X4 Z5 X6] +
+(0.006238164536781785+0j) [Y0 Z1 Y2 X4 Z5 Z6 X7] +
+(0.024934767638066777+0j) [Y0 Z1 Y2 Y4 Z5 Y6] +
+(0.006238164536781785+0j) [Y0 Z1 Y2 Y4 Z5 Z6 Y7] +
+(0.01064615061000441+0j) [Y0 Z1 Y2 X5 X6] +
+(-0.02405865607179162+0j) [Y0 Z1 Y2 X5 Z6 X7] +
+(0.01064615061000441+0j) [Y0 Z1 Y2 Y5 Y6] +
+(-0.02405865607179162+0j) [Y0 Z1 Y2 Y5 Z6 Y7] +
+(0.006238164536781785+0j) [Y0 Z1 Z2 Y3 X4 Z5 X6] +
+(0.017676042848046+0j) [Y0 Z1 Z2 Y3 X4 Z5 Z6 X7] +
+(0.006238164536781785+0j) [Y0 Z1 Z2 Y3 Y4 Z5 Y6] +
+(0.017676042848046+0j) [Y0 Z1 Z2 Y3 Y4 Z5 Z6 Y7] +
+(-0.012634174981636702+0j) [Y0 Z1 Z2 Y3 X5 X6] +
+(-0.008407495254282365+0j) [Y0 Z1 Z2 Y3 X5 Z6 X7] +
+(-0.012634174981636702+0j) [Y0 Z1 Z2 Y3 Y5 Y6] +
+(-0.008407495254282365+0j) [Y0 Z1 Z2 Y3 Y5 Z6 Y7] +
+(0.23477912416223978+0j) [Z0] +
+(0.07755867751288839+0j) [Z0 Z1] +
+(-0.0025837417313018998+0j) [Z0 X2 X3] +
+(-0.0025837417313018998+0j) [Z0 Y2 Y3] +
+(0.07733247653566544+0j) [Z0 Z2] +
+(0.09333330628587891+0j) [Z0 Z3] +
+(0.009461086978084416+0j) [Z0 X4 X5] +
+(0.009461086978084416+0j) [Z0 Y4 Y5] +
+(0.11171627890094822+0j) [Z0 Z4] +
+(0.10416313162196875+0j) [Z0 Z5] +
+(-0.008821906268083685+0j) [Z0 X6 X7] +
+(-0.008821906268083685+0j) [Z0 Y6 Y7] +
+(0.10226724417373223+0j) [Z0 Z6] +
+(0.1110093491339249+0j) [Z0 Z7] +
+(0.010646150610004408+0j) [X1 X2 X4 Z5 X6] +
+(-0.012634174981636702+0j) [X1 X2 X4 Z5 Z6 X7] +
+(0.010646150610004408+0j) [X1 X2 Y4 Z5 Y6] +
+(-0.012634174981636702+0j) [X1 X2 Y4 Z5 Z6 Y7] +
+(0.01934883030685017+0j) [X1 X2 X5 X6] +
+(-0.008901638567159004+0j) [X1 X2 X5 Z6 X7] +
+(0.01934883030685017+0j) [X1 X2 Y5 Y6] +
+(-0.008901638567159004+0j) [X1 X2 Y5 Z6 Y7] +
+(-0.02405865607179162+0j) [X1 Z2 X3 X4 Z5 X6] +
+(-0.008407495254282365+0j) [X1 Z2 X3 X4 Z5 Z6 X7] +
+(-0.02405865607179162+0j) [X1 Z2 X3 Y4 Z5 Y6] +
+(-0.008407495254282365+0j) [X1 Z2 X3 Y4 Z5 Z6 Y7] +
+(-0.008901638567159004+0j) [X1 Z2 X3 X5 X6] +
+(0.02661467635818076+0j) [X1 Z2 X3 X5 Z6 X7] +
+(-0.008901638567159004+0j) [X1 Z2 X3 Y5 Y6] +
+(0.02661467635818076+0j) [X1 Z2 X3 Y5 Z6 Y7] +
+(0.010646150610004408+0j) [Y1 Y2 X4 Z5 X6] +
+(-0.012634174981636702+0j) [Y1 Y2 X4 Z5 Z6 X7] +
+(0.010646150610004408+0j) [Y1 Y2 Y4 Z5 Y6] +
+(-0.012634174981636702+0j) [Y1 Y2 Y4 Z5 Z6 Y7] +
+(0.01934883030685017+0j) [Y1 Y2 X5 X6] +
+(-0.008901638567159004+0j) [Y1 Y2 X5 Z6 X7] +
+(0.01934883030685017+0j) [Y1 Y2 Y5 Y6] +
+(-0.008901638567159004+0j) [Y1 Y2 Y5 Z6 Y7] +
+(-0.02405865607179162+0j) [Y1 Z2 Y3 X4 Z5 X6] +
+(-0.008407495254282365+0j) [Y1 Z2 Y3 X4 Z5 Z6 X7] +
+(-0.02405865607179162+0j) [Y1 Z2 Y3 Y4 Z5 Y6] +
+(-0.008407495254282365+0j) [Y1 Z2 Y3 Y4 Z5 Z6 Y7] +
+(-0.008901638567159004+0j) [Y1 Z2 Y3 X5 X6] +
+(0.02661467635818076+0j) [Y1 Z2 Y3 X5 Z6 X7] +
+(-0.008901638567159004+0j) [Y1 Z2 Y3 Y5 Y6] +
+(0.02661467635818076+0j) [Y1 Z2 Y3 Y5 Z6 Y7] +
+(0.07873792608304384+0j) [Z1] +
+(-0.00523830590529105+0j) [Z1 X2 X3] +
+(-0.00523830590529105+0j) [Z1 Y2 Y3] +
+(0.08417436673290517+0j) [Z1 Z2] +
+(0.07928157442541509+0j) [Z1 Z3] +
+(-0.002045316086697597+0j) [Z1 X4 X5] +
+(-0.002045316086697597+0j) [Z1 Y4 Y5] +
+(0.10416313162196875+0j) [Z1 Z4] +
+(0.10702518092938991+0j) [Z1 Z5] +
+(0.0036633326618679546+0j) [Z1 X6 X7] +
+(0.0036633326618679546+0j) [Z1 Y6 Y7] +
+(0.10352319703975532+0j) [Z1 Z6] +
+(0.10589625078359582+0j) [Z1 Z7] +
+(-0.004483252597706383+0j) [X2 X3] +
+(-0.024395190370696845+0j) [X2 X3 X4 X5] +
+(-0.024395190370696845+0j) [X2 X3 Y4 Y5] +
+(-0.008821906268083685+0j) [X2 X3 Z4] +
+(0.0036633326618679546+0j) [X2 X3 Z5] +
+(0.025689949571830042+0j) [X2 X3 X6 X7] +
+(0.025689949571830042+0j) [X2 X3 Y6 Y7] +
+(0.007649625311503733+0j) [X2 X3 Z6] +
+(-0.012996931212451403+0j) [X2 X3 Z7] +
+(-0.004483252597706383+0j) [Y2 Y3] +
+(-0.024395190370696845+0j) [Y2 Y3 X4 X5] +
+(-0.024395190370696845+0j) [Y2 Y3 Y4 Y5] +
+(-0.008821906268083685+0j) [Y2 Y3 Z4] +
+(0.0036633326618679546+0j) [Y2 Y3 Z5] +
+(0.025689949571830042+0j) [Y2 Y3 X6 X7] +
+(0.025689949571830042+0j) [Y2 Y3 Y6 Y7] +
+(0.007649625311503733+0j) [Y2 Y3 Z6] +
+(-0.012996931212451403+0j) [Y2 Y3 Z7] +
+(0.07666866090982398+0j) [Z2] +
+(0.07784794028451397+0j) [Z2 Z3] +
+(-0.007801932778402895+0j) [Z2 X4 X5] +
+(-0.007801932778402895+0j) [Z2 Y4 Y5] +
+(0.10226724417373223+0j) [Z2 Z4] +
+(0.10352319703975532+0j) [Z2 Z5] +
+(0.007649625311503733+0j) [Z2 X6 X7] +
+(0.007649625311503733+0j) [Z2 Y6 Y7] +
+(0.1091589968106852+0j) [Z2 Z6] +
+(0.10353788985634402+0j) [Z2 Z7] +
+(-0.12227609382406271+0j) [Z3] +
+(0.013474635539233077+0j) [Z3 X4 X5] +
+(0.013474635539233077+0j) [Z3 Y4 Y5] +
+(0.1110093491339249+0j) [Z3 Z4] +
+(0.10589625078359582+0j) [Z3 Z5] +
+(-0.012996931212451403+0j) [Z3 X6 X7] +
+(-0.012996931212451403+0j) [Z3 Y6 Y7] +
+(0.10353788985634402+0j) [Z3 Z6] +
+(0.11970161375543707+0j) [Z3 Z7] +
+(-0.0020776113200212066+0j) [X4 X5] +
+(-0.011761015389060145+0j) [X4 X5 X6 X7] +
+(-0.0003365342989052316+0j) [X4 X5 Y6 Y7] +
+(0.002844217831601514+0j) [X4 X5 Z6] +
+(0.005067140284950715+0j) [X4 X5 Z7] +
+(-0.011424481090154913+0j) [X4 Y5 Y6 X7] +
+(-0.011424481090154913+0j) [Y4 X5 X6 Y7] +
+(-0.0020776113200212066+0j) [Y4 Y5] +
+(-0.0003365342989052316+0j) [Y4 Y5 X6 X7] +
+(-0.011761015389060145+0j) [Y4 Y5 Y6 Y7] +
+(0.002844217831601514+0j) [Y4 Y5 Z6] +
+(0.005067140284950715+0j) [Y4 Y5 Z7] +
+(0.23477912416223978+0j) [Z4] +
+(0.07755867751288839+0j) [Z4 Z5] +
+(-0.0025837417313018998+0j) [Z4 X6 X7] +
+(-0.0025837417313018998+0j) [Z4 Y6 Y7] +
+(0.07733247653566544+0j) [Z4 Z6] +
+(0.09333330628587891+0j) [Z4 Z7] +
+(0.07873792608304378+0j) [Z5] +
+(-0.00523830590529105+0j) [Z5 X6 X7] +
+(-0.00523830590529105+0j) [Z5 Y6 Y7] +
+(0.08417436673290517+0j) [Z5 Z6] +
+(0.07928157442541509+0j) [Z5 Z7] +
+(-0.004483252597706375+0j) [X6 X7] +
+(-0.004483252597706375+0j) [Y6 Y7] +
+(0.07666866090982404+0j) [Z6] +
+(0.07784794028451397+0j) [Z6 Z7] +
+(-0.12227609382406274+0j) [Z7]
\ No newline at end of file

From f6f6f15987205c0225130efcb1aafd0ce8930ede Mon Sep 17 00:00:00 2001
From: James Brown <james.brown@1qbit.com>
Date: Mon, 10 Jan 2022 15:30:04 -0500
Subject: [PATCH 05/10] fix for failed test

---
 tangelo/algorithms/variational/vqe_solver.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tangelo/algorithms/variational/vqe_solver.py b/tangelo/algorithms/variational/vqe_solver.py
index d3ca91cd5..9b9b83e7a 100644
--- a/tangelo/algorithms/variational/vqe_solver.py
+++ b/tangelo/algorithms/variational/vqe_solver.py
@@ -196,12 +196,12 @@ def build(self):
                 raise TypeError(f"Invalid ansatz dataype. Expecting instance of Ansatz class, or one of built-in options:\n\t{self.builtin_ansatze}")
 
         # Building with a qubit Hamiltonian.
-        elif not isinstance(self.ansatz, Ansatz) and self.ansatz in [BuiltInAnsatze.HEA, BuiltInAnsatze.VSQS]:
+        elif self.ansatz in [BuiltInAnsatze.HEA, BuiltInAnsatze.VSQS]:
             if self.ansatz == BuiltInAnsatze.HEA:
                 self.ansatz = HEA(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options)
             elif self.ansatz == BuiltInAnsatze.VSQS:
                 self.ansatz = VSQS(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options)
-        else:
+        elif not isinstance(self.ansatz, Ansatz):
             raise TypeError(f"Invalid ansatz dataype. Expecting a custom Ansatz (Ansatz class).")
 
         # Set ansatz initial parameters (default or use input), build corresponding ansatz circuit

From ccdd4079562ad00720bd01676c081b4ae05f6c45 Mon Sep 17 00:00:00 2001
From: James Brown <james.brown@1qbit.com>
Date: Wed, 12 Jan 2022 12:10:10 -0500
Subject: [PATCH 06/10] style and docstring improvements

---
 .../variational/tests/test_vqe_solver.py      |  2 +-
 .../ansatz_generator/ansatz_utils.py          |  5 +++--
 .../ansatz_generator/tests/test_vsqs.py       |  1 -
 tangelo/toolboxes/ansatz_generator/vsqs.py    | 20 +++++++++----------
 4 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/tangelo/algorithms/variational/tests/test_vqe_solver.py b/tangelo/algorithms/variational/tests/test_vqe_solver.py
index e910be737..6c117de9d 100644
--- a/tangelo/algorithms/variational/tests/test_vqe_solver.py
+++ b/tangelo/algorithms/variational/tests/test_vqe_solver.py
@@ -156,7 +156,7 @@ def test_simulate_vsqs_h2(self):
 
         vqe_options = {"molecule": None, "qubit_hamiltonian": qubit_hamiltonian, "ansatz": BuiltInAnsatze.VSQS, "qubit_mapping": "jw",
                        "ansatz_options": {"intervals": 3, "time": 3, "qubit_hamiltonian": qubit_hamiltonian,
-                       "hini": hini, "reference_state": reference_state}}
+                                          "hini": hini, "reference_state": reference_state}}
         vqe_solver = VQESolver(vqe_options)
         vqe_solver.build()
 
diff --git a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
index 68cb31630..c74e18550 100644
--- a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
+++ b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
@@ -94,8 +94,9 @@ def get_exponentiated_qubit_operator_circuit(qubit_op, time=1., variational=Fals
         variational (bool) : Whether the coefficients are variational
         trotter_order (int): order of trotter approximation, only 1 or 2 are supported.
         return_phase (bool): Return the global-phase generated
-        pauli_order (list): The desired pauli_word order defined as a list with elements (pauli_word, coeff)
-            with corresponding dictionary elements pauli_word: coeff in QubitOperator terms.items()
+        pauli_order (list): The desired pauli_word order for trotterization defined as a list of (pauli_word, coeff)
+            elements which have matching dictionary elements pauli_word: coeff in QubitOperator terms.items().
+            The coeff in pauli_order is used.
 
     Returns:
         Circuit: circuit corresponding to exponentiation of qubit operator
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
index d53e6bf02..dda9e60c8 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
@@ -32,7 +32,6 @@ class VSQSTest(unittest.TestCase):
 
     def test_vsqs_set_var_params(self):
         """Verify behavior of set_var_params for different inputs (list, numpy array).
-        MP2 have their own tests.
         """
 
         vsqs_ansatz = VSQS(mol_H2_sto3g)
diff --git a/tangelo/toolboxes/ansatz_generator/vsqs.py b/tangelo/toolboxes/ansatz_generator/vsqs.py
index 87595cdc4..efd134979 100644
--- a/tangelo/toolboxes/ansatz_generator/vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/vsqs.py
@@ -37,15 +37,15 @@ class VSQS(Ansatz):
 
     Args:
         molecule (SecondQuantizedMolecule): The molecular system. Default: None
-        mapping (str): One of the support fermion to qubit mappings. Default : "JW"
-        up_then_down (bool): change basis ordering putting all spin up orbitals first, followed by all spin down.
+        mapping (str): One of the supported fermion to qubit mappings. Default : "JW"
+        up_then_down (bool): Change basis ordering, putting all spin up orbitals first, followed by all spin down.
             Default: False (alternating spin up/down ordering)
         intervals (int): The number of steps in the VSQS process. Must be greater than 1. Default: 2
         time (float): The propagation time. Default: 1
-        qubit_hamiltonian (QubitOperator): The qubit hamiltonian to evolve. Default: None
+        qubit_hamiltonian (QubitOperator): The qubit Hamiltonian to evolve. Default: None
         reference_state (Circuit): The reference state for the propagation as defined by a Circuit. Mandatory if supplying
             a qubit_hamiltonian. Default: None
-        hini (QubitOperator): The initial qubit Hamiltonian that one corresponds to the reference state. Mandatory if supplying
+        hini (QubitOperator): The initial qubit Hamiltonian that corresponds to the reference state. Mandatory if supplying
             a qubit_hamiltonian. Default: None
         hnav (QubitOperator): The navigator Hamiltonian. Default: None
     """
@@ -67,10 +67,10 @@ def __init__(self, molecule=None, mapping="jw", up_then_down=False, intervals=2,
         if molecule is None:
             self.qubit_hamiltonian = qubit_hamiltonian
             if not isinstance(hini, ofQubitOperator):
-                raise ValueError("When providing a qubit hamiltonian, an initial Hamiltonian must also be provided")
+                raise ValueError("When providing a qubit hamiltonian, an initial qubit Hamiltonian must also be provided")
             self.hini = hini
             if not isinstance(reference_state, Circuit):
-                raise ValueError("Reference state must be provided when simulating a qubit hamiltonian directly")
+                raise ValueError("Reference state Circuit must be provided when simulating a qubit hamiltonian directly")
             self.reference_state = reference_state
         else:
             self.n_electrons = molecule.n_active_electrons
@@ -166,19 +166,19 @@ def _init_params(self):
         a_i = step*i, b_i=1-step*i, c_i= 1-step*i i<n_intervals/2, step*i i>n_intervals/2"""
         a = np.zeros(self.intervals+1)
         b = np.zeros(self.intervals+1)
-        c = np.zeros(self.intervals+1)
         a[0] = 1
         b[self.intervals] = 1
-        step = 1/self.intervals
+        step_size = 1/self.intervals
         for i in range(1, self.intervals):
-            a[i] = (1 - i * step)
-            b[i] = (i * step)
+            a[i] = (1 - i * step_size)
+            b[i] = (i * step_size)
         all_params = np.zeros(self.stride * (self.intervals + 1))
         if self.hnav is None:
             for i in range(self.intervals + 1):
                 all_params[2 * i] = a[i]
                 all_params[2 * i + 1] = b[i]
         else:
+            c = np.zeros(self.intervals+1)
             c[0:self.intervals//2] = b[0:self.intervals//2]
             c[self.intervals//2:self.intervals+1] = a[self.intervals//2:self.intervals+1]
             for i in range(self.intervals + 1):

From 164837bf1304e65c08c0b19704b3418234b67dc7 Mon Sep 17 00:00:00 2001
From: James Brown <james.brown@1qbit.com>
Date: Thu, 13 Jan 2022 10:52:29 -0500
Subject: [PATCH 07/10] style improvements

---
 .../algorithms/variational/tests/test_vqe_solver.py |  1 -
 tangelo/toolboxes/ansatz_generator/ansatz_utils.py  |  2 +-
 tangelo/toolboxes/ansatz_generator/vsqs.py          | 13 +++++--------
 3 files changed, 6 insertions(+), 10 deletions(-)

diff --git a/tangelo/algorithms/variational/tests/test_vqe_solver.py b/tangelo/algorithms/variational/tests/test_vqe_solver.py
index 6c117de9d..a90b130a1 100644
--- a/tangelo/algorithms/variational/tests/test_vqe_solver.py
+++ b/tangelo/algorithms/variational/tests/test_vqe_solver.py
@@ -12,7 +12,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import re
 import unittest
 import numpy as np
 
diff --git a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
index c74e18550..bf8ad34d7 100644
--- a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
+++ b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
@@ -96,7 +96,7 @@ def get_exponentiated_qubit_operator_circuit(qubit_op, time=1., variational=Fals
         return_phase (bool): Return the global-phase generated
         pauli_order (list): The desired pauli_word order for trotterization defined as a list of (pauli_word, coeff)
             elements which have matching dictionary elements pauli_word: coeff in QubitOperator terms.items().
-            The coeff in pauli_order is used.
+            The coeff in pauli_order is used to generate the exponential.
 
     Returns:
         Circuit: circuit corresponding to exponentiation of qubit operator
diff --git a/tangelo/toolboxes/ansatz_generator/vsqs.py b/tangelo/toolboxes/ansatz_generator/vsqs.py
index efd134979..e694b0ee2 100644
--- a/tangelo/toolboxes/ansatz_generator/vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/vsqs.py
@@ -163,7 +163,7 @@ def _update_gate_params_for_qu_op(self, qu_op, n_var_start, var_param, num_terms
 
     def _init_params(self):
         """Generate initial parameters for the VSQS algorithm.
-        a_i = step*i, b_i=1-step*i, c_i= 1-step*i i<n_intervals/2, step*i i>n_intervals/2"""
+        a_i = step*i, b_i=1-step*i, c_i= 1-step*i i<=n_intervals/2, step*i i>n_intervals/2"""
         a = np.zeros(self.intervals+1)
         b = np.zeros(self.intervals+1)
         a[0] = 1
@@ -174,17 +174,14 @@ def _init_params(self):
             b[i] = (i * step_size)
         all_params = np.zeros(self.stride * (self.intervals + 1))
         if self.hnav is None:
-            for i in range(self.intervals + 1):
-                all_params[2 * i] = a[i]
-                all_params[2 * i + 1] = b[i]
+            # order [a[0],b[0],a[1],b[1],...]
+            all_params = np.dstack((a, b)).flatten()
         else:
             c = np.zeros(self.intervals+1)
             c[0:self.intervals//2] = b[0:self.intervals//2]
             c[self.intervals//2:self.intervals+1] = a[self.intervals//2:self.intervals+1]
-            for i in range(self.intervals + 1):
-                all_params[3 * i] = a[i]
-                all_params[3 * i + 1] = b[i]
-                all_params[3 * i + 2] = c[i]
+            # order [a[0],b[0],c[0],a[1],b[1],c[1],...]
+            all_params = np.dstack((a, b, c)).flatten()
         return all_params
 
     def prepare_reference_state(self):

From ffc77a318f2749416628ac0fb970b27827120565 Mon Sep 17 00:00:00 2001
From: James Brown <james.brown@1qbit.com>
Date: Fri, 14 Jan 2022 13:29:40 -0500
Subject: [PATCH 08/10] PR requested changes

---
 .../variational/tests/test_vqe_solver.py      |   6 +-
 .../toolboxes/ansatz_generator/__init__.py    |  10 ++
 .../ansatz_generator/ansatz_utils.py          |   2 +-
 .../ansatz_generator/tests/test_vsqs.py       |  11 +-
 tangelo/toolboxes/ansatz_generator/vsqs.py    | 154 ++++++++++--------
 .../molecular_computation/molecule.py         |  53 +++---
 6 files changed, 136 insertions(+), 100 deletions(-)

diff --git a/tangelo/algorithms/variational/tests/test_vqe_solver.py b/tangelo/algorithms/variational/tests/test_vqe_solver.py
index a90b130a1..eefabdc4f 100644
--- a/tangelo/algorithms/variational/tests/test_vqe_solver.py
+++ b/tangelo/algorithms/variational/tests/test_vqe_solver.py
@@ -142,7 +142,7 @@ def test_simulate_vsqs_h2(self):
         qubit_hamiltonian/hini/reference_state input
         """
         vqe_options = {"molecule": mol_H2_sto3g, "ansatz": BuiltInAnsatze.VSQS, "qubit_mapping": "jw",
-                       "verbose": True, "ansatz_options": {"intervals": 3, "time": 3}}
+                       "verbose": False, "ansatz_options": {"intervals": 3, "time": 3}}
         vqe_solver = VQESolver(vqe_options)
         vqe_solver.build()
 
@@ -150,12 +150,12 @@ def test_simulate_vsqs_h2(self):
         self.assertAlmostEqual(energy, -1.137270, delta=1e-4)
 
         qubit_hamiltonian = vqe_solver.qubit_hamiltonian
-        hini = vqe_solver.ansatz.hini
+        h_init = vqe_solver.ansatz.h_init
         reference_state = vqe_solver.ansatz.prepare_reference_state()
 
         vqe_options = {"molecule": None, "qubit_hamiltonian": qubit_hamiltonian, "ansatz": BuiltInAnsatze.VSQS, "qubit_mapping": "jw",
                        "ansatz_options": {"intervals": 3, "time": 3, "qubit_hamiltonian": qubit_hamiltonian,
-                                          "hini": hini, "reference_state": reference_state}}
+                                          "h_init": h_init, "reference_state": reference_state}}
         vqe_solver = VQESolver(vqe_options)
         vqe_solver.build()
 
diff --git a/tangelo/toolboxes/ansatz_generator/__init__.py b/tangelo/toolboxes/ansatz_generator/__init__.py
index 532746351..1940ada44 100644
--- a/tangelo/toolboxes/ansatz_generator/__init__.py
+++ b/tangelo/toolboxes/ansatz_generator/__init__.py
@@ -11,3 +11,13 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
+
+from .vsqs import VSQS
+from .qcc import QCC
+from .qmf import QMF
+from .uccsd import UCCSD
+from .adapt_ansatz import ADAPTAnsatz
+from .rucc import RUCC
+from .upccgsd import UpCCGSD
+from .hea import HEA
+from .variational_circuit import VariationalCircuitAnsatz
diff --git a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
index bf8ad34d7..825f7f587 100644
--- a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
+++ b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
@@ -105,7 +105,7 @@ def get_exponentiated_qubit_operator_circuit(qubit_op, time=1., variational=Fals
     if pauli_order is None:
         pauli_words = list(qubit_op.terms.items())
     elif isinstance(pauli_order, list):
-        pauli_words = deepcopy(pauli_order)
+        pauli_words = pauli_order.copy()
     else:
         raise ValueError("ordered terms must be a list with elements (keys, values) of qubit_op.terms.items()")
 
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
index dda9e60c8..b36e19826 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
@@ -13,14 +13,15 @@
 # limitations under the License.
 
 import unittest
+
 import numpy as np
 import os
 from openfermion import load_operator
 
 from tangelo.molecule_library import mol_H2_sto3g
-from tangelo.toolboxes.operators.operators import QubitOperator
+from tangelo.toolboxes.operators import QubitOperator
 from tangelo.toolboxes.qubit_mappings import jordan_wigner, symmetry_conserving_bravyi_kitaev
-from tangelo.toolboxes.ansatz_generator.vsqs import VSQS
+from tangelo.toolboxes.ansatz_generator import VSQS
 from tangelo.linq import Simulator
 from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit
 
@@ -81,7 +82,7 @@ def test_vsqs_H4_doublecation(self):
         reference_state = get_reference_circuit(8, 2, "jw", up_then_down=True, spin=0)
 
         # Build circuit
-        vsqs_ansatz = VSQS(qubit_hamiltonian=qubit_hamiltonian, hini=initial_hamiltonian, reference_state=reference_state,
+        vsqs_ansatz = VSQS(qubit_hamiltonian=qubit_hamiltonian, h_init=initial_hamiltonian, reference_state=reference_state,
                            intervals=5, time=5, trotter_order=2)
         vsqs_ansatz.build_circuit()
 
@@ -91,7 +92,7 @@ def test_vsqs_H4_doublecation(self):
         energy = sim.get_expectation_value(qubit_hamiltonian, vsqs_ansatz.circuit)
         self.assertAlmostEqual(energy, -0.85425, delta=1e-4)
 
-    def test_vsqs_H2_with_hnav(self):
+    def test_vsqs_H2_with_h_nav(self):
         """Verify closed-shell VSQS functionalities for H2 with navigator hamiltonian"""
         navigator_hamiltonian = (QubitOperator('X0 Y1', 0.03632537110234512) + QubitOperator('Y0 X1', 0.03632537110234512)
                                  + QubitOperator('Y0', 2.e-5) + QubitOperator('Y1', 2.e-5))
@@ -101,7 +102,7 @@ def test_vsqs_H2_with_hnav(self):
 
         # Build circuit
         vsqs_ansatz = VSQS(mol_H2_sto3g, intervals=2, time=1, mapping='scbk', up_then_down=True, trotter_order=2,
-                           hnav=navigator_hamiltonian)
+                           h_nav=navigator_hamiltonian)
         vsqs_ansatz.build_circuit()
 
         # Assert energy returned is as expected for given parameters
diff --git a/tangelo/toolboxes/ansatz_generator/vsqs.py b/tangelo/toolboxes/ansatz_generator/vsqs.py
index e694b0ee2..0f4532a6e 100644
--- a/tangelo/toolboxes/ansatz_generator/vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/vsqs.py
@@ -15,8 +15,6 @@
 """This module defines the Variationally Scheduled Quantum Simulation class. It provides an
 Adiabatic State Preparation (ASP) inspired ansatz as described in https://arxiv.org/abs/2003.09913."""
 
-from copy import deepcopy
-
 import numpy as np
 from openfermion import QubitOperator as ofQubitOperator
 
@@ -33,7 +31,7 @@ class VSQS(Ansatz):
     https://arxiv.org/abs/2003.09913
 
     Must supply either a molecule or a qubit_hamiltonian. If supplying a qubit_hamiltonian, must also supply
-    a reference_state Circuit and a hini QubitOperator.
+    a reference_state Circuit and a h_init QubitOperator.
 
     Args:
         molecule (SecondQuantizedMolecule): The molecular system. Default: None
@@ -41,17 +39,18 @@ class VSQS(Ansatz):
         up_then_down (bool): Change basis ordering, putting all spin up orbitals first, followed by all spin down.
             Default: False (alternating spin up/down ordering)
         intervals (int): The number of steps in the VSQS process. Must be greater than 1. Default: 2
-        time (float): The propagation time. Default: 1
+        time (float): The propagation time. Default: 1.
         qubit_hamiltonian (QubitOperator): The qubit Hamiltonian to evolve. Default: None
         reference_state (Circuit): The reference state for the propagation as defined by a Circuit. Mandatory if supplying
             a qubit_hamiltonian. Default: None
-        hini (QubitOperator): The initial qubit Hamiltonian that corresponds to the reference state. Mandatory if supplying
+        h_init (QubitOperator): The initial qubit Hamiltonian that corresponds to the reference state. Mandatory if supplying
             a qubit_hamiltonian. Default: None
-        hnav (QubitOperator): The navigator Hamiltonian. Default: None
+        h_nav (QubitOperator): The navigator Hamiltonian. Default: None
+        trotter_order (int): The order of the Trotterization for each qubit operator. Default: 1
     """
 
-    def __init__(self, molecule=None, mapping="jw", up_then_down=False, intervals=2, time=1, qubit_hamiltonian=None,
-                 reference_state=None, hini=None, hnav=None, trotter_order=1):
+    def __init__(self, molecule=None, mapping="jw", up_then_down=False, intervals=2, time=1., qubit_hamiltonian=None,
+                 reference_state=None, h_init=None, h_nav=None, trotter_order=1):
 
         self.up_then_down = up_then_down
         self.mapping = mapping
@@ -62,13 +61,16 @@ def __init__(self, molecule=None, mapping="jw", up_then_down=False, intervals=2,
         self.time = time
         self.dt = self.time/self.intervals
         self.reference_state = reference_state
-        self.trotter_order = trotter_order
+        if trotter_order in [1, 2]:
+            self.trotter_order = trotter_order
+        else:
+            raise ValueError("Only trotter_order = 1, 2 is supported")
 
         if molecule is None:
             self.qubit_hamiltonian = qubit_hamiltonian
-            if not isinstance(hini, ofQubitOperator):
+            if not isinstance(h_init, ofQubitOperator):
                 raise ValueError("When providing a qubit hamiltonian, an initial qubit Hamiltonian must also be provided")
-            self.hini = hini
+            self.h_init = h_init
             if not isinstance(reference_state, Circuit):
                 raise ValueError("Reference state Circuit must be provided when simulating a qubit hamiltonian directly")
             self.reference_state = reference_state
@@ -80,54 +82,56 @@ def __init__(self, molecule=None, mapping="jw", up_then_down=False, intervals=2,
             self.qubit_hamiltonian = fermion_to_qubit_mapping(molecule.fermionic_hamiltonian, n_spinorbitals=self.n_spinorbitals,
                                                               n_electrons=self.n_electrons, mapping=self.mapping,
                                                               up_then_down=self.up_then_down, spin=self.spin)
-            self.hini = self._build_hini(molecule) if hini is None else hini
-        self.hfin = self.qubit_hamiltonian
-        self.hnav = hnav
-
-        self.hfin_list, self.n_hfin = self._remove_constant_and_count_terms(self.qubit_hamiltonian)
-        self.hini_list, self.n_hini = self._remove_constant_and_count_terms(self.hini)
-        if self.hnav is None:
+            self.h_init = self._build_h_init(molecule) if h_init is None else h_init
+        self.h_final = self.qubit_hamiltonian
+        self.h_nav = h_nav
+
+        def qu_op_to_list(qu_op):
+            """Remove consant term and convert QubitOperator to list of (term, coeff)"""
+            new_qu_op = qu_op - qu_op.constant
+            new_qu_op.compress()
+            return list(new_qu_op.terms.items())
+
+        self.h_final_list = qu_op_to_list(self.h_final)
+        self.n_h_final = len(self.h_final_list)
+        self.h_init_list = qu_op_to_list(self.h_init)
+        self.n_h_init = len(self.h_init_list)
+        if self.h_nav is None:
             self.stride = 2
-            self.n_hnav = 0
+            self.n_h_nav = 0
         else:
-            if isinstance(self.hnav, ofQubitOperator):
+            if isinstance(self.h_nav, ofQubitOperator):
                 self.stride = 3
-                self.hnav_list, self.n_hnav = self._remove_constant_and_count_terms(self.hnav)
+                self.h_nav_list = qu_op_to_list(self.h_nav)
+                self.n_h_nav = len(self.h_nav_list)
             else:
                 raise ValueError("Navigator Hamiltonian must be a QubitOperator")
 
         self.n_var_params = (intervals - 1) * self.stride
-        self.n_var_gates = (self.n_hini + self.n_hfin + self.n_hnav) * self.trotter_order
+        self.n_var_gates = (self.n_h_init + self.n_h_final + self.n_h_nav) * self.trotter_order
 
         self.var_params = None
         self.circuit = None
 
-    def _remove_constant_and_count_terms(self, qu_op):
-        """count of non-zero terms in a QubitOperator and return the list of terms.items()"""
-        count = 0
-        new_qu_op = deepcopy(qu_op)
-        for term in new_qu_op.terms.keys():
-            if term:
-                count += 1
-            else:
-                new_qu_op.terms[term] = 0
-        new_qu_op.compress()
-        return list(new_qu_op.terms.items()), count
-
-    def _build_hini(self, molecule):
-        """Return the initial Hamiltonian (hini) composed of the one-body terms derived from the diagonal of Fock
+    def _build_h_init(self, molecule):
+        """Return the initial Hamiltonian (h_init) composed of the one-body terms derived from the diagonal of Fock
         matrix and one-body off-diagonal terms"""
-        fock = molecule.mean_field.mo_coeff.T @ molecule.mean_field.get_fock() @ molecule.mean_field.mo_coeff
-        h1 = molecule.mean_field.mo_coeff.T @ molecule.mean_field.get_hcore() @ molecule.mean_field.mo_coeff
-        hf_ferm = FermionOperator()
-        for i in range(self.n_spinorbitals//2):
-            for j in range(self.n_spinorbitals//2):
+        core_constant, h1, two_body = molecule.get_active_space_integrals()
+        dfock = np.diag(h1).copy()
+        n_active_occupied = len(molecule.active_occupied)
+        for j in range(molecule.n_active_mos):
+            for i in range(n_active_occupied):
+                dfock[j] += 2*two_body[i, j, j, i] - 1*two_body[i, j, i, j]
+
+        hf_ferm = FermionOperator((), core_constant)
+        for i in range(molecule.n_active_mos):
+            for j in range(molecule.n_active_mos):
                 if i != j:
                     hf_ferm += FermionOperator(((i * 2, 1), (j * 2, 0)), h1[i, j])
                     hf_ferm += FermionOperator(((i * 2 + 1, 1), (j * 2 + 1, 0)), h1[i, j])
                 else:
-                    hf_ferm += FermionOperator(((i * 2, 1), (j * 2, 0)), fock[i, j])
-                    hf_ferm += FermionOperator(((i * 2 + 1, 1), (j * 2 + 1, 0)), fock[i, j])
+                    hf_ferm += FermionOperator(((i * 2, 1), (j * 2, 0)), dfock[i])
+                    hf_ferm += FermionOperator(((i * 2 + 1, 1), (j * 2 + 1, 0)), dfock[j])
         return fermion_to_qubit_mapping(hf_ferm, mapping=self.mapping, n_spinorbitals=self.n_spinorbitals, n_electrons=self.n_electrons,
                                         up_then_down=self.up_then_down, spin=self.spin)
 
@@ -146,24 +150,33 @@ def set_var_params(self, var_params=None):
     def update_var_params(self, var_params):
         """Update the variational parameters in the circuit without rebuilding."""
         for i in range(self.intervals-1):
-            self._update_gate_params_for_qu_op(self.hini_list, self.n_var_gates * i, var_params[self.stride*i], self.n_hini)
-            self._update_gate_params_for_qu_op(self.hfin_list, self.n_var_gates * i + self.n_hini * self.trotter_order,
-                                               var_params[self.stride*i+1], self.n_hfin)
-            if self.hnav is not None:
-                self._update_gate_params_for_qu_op(self.hnav_list, self.n_var_gates * i + (self.n_hini + self.n_hfin) * self.trotter_order,
-                                                   var_params[self.stride*i+2], self.n_hnav)
-
-    def _update_gate_params_for_qu_op(self, qu_op, n_var_start, var_param, num_terms):
+            self._update_gate_params_for_qu_op(self.h_init_list, self.n_var_gates * i, var_params[self.stride*i], self.n_h_init)
+            self._update_gate_params_for_qu_op(self.h_final_list, self.n_var_gates * i + self.n_h_init * self.trotter_order,
+                                               var_params[self.stride*i+1], self.n_h_final)
+            if self.h_nav is not None:
+                self._update_gate_params_for_qu_op(self.h_nav_list, self.n_var_gates * i + (self.n_h_init + self.n_h_final) * self.trotter_order,
+                                                   var_params[self.stride*i+2], self.n_h_nav)
+
+    def _update_gate_params_for_qu_op(self, qu_op_list, n_var_start, var_param, num_terms):
+        """Updates the corresponding circuit variational_gates for a QubitOperator defined by term order qu_op_list
+
+        Args:
+            qu_op_list :: The list with elements (term, coeff) that defines the trotterization of a QubitOperator
+            n_var_start :: The varational_gates position that the trotterization of qu_op starts
+            var_param :: The variational parameter (evolution time) for qu_op. Same for all terms in qu_op
+            num_terms :: The number of terms in qu_op
+        """
         prefac = 2 / self.trotter_order * self.dt * var_param
-        for i, (_, coeff) in enumerate(qu_op):
+        for i, (_, coeff) in enumerate(qu_op_list):
             self.circuit._variational_gates[n_var_start+i].parameter = prefac * coeff
         if self.trotter_order == 2:
-            for i, (_, coeff) in enumerate(list(reversed(qu_op))):
+            for i, (_, coeff) in enumerate(list(reversed(qu_op_list))):
                 self.circuit._variational_gates[n_var_start+num_terms+i].parameter = prefac * coeff
 
     def _init_params(self):
         """Generate initial parameters for the VSQS algorithm.
-        a_i = step*i, b_i=1-step*i, c_i= 1-step*i i<=n_intervals/2, step*i i>n_intervals/2"""
+        a_i = step*i, b_i=1-step*i, c_i= 1-step*i i<=n_intervals/2, step*i i>n_intervals/2
+        """
         a = np.zeros(self.intervals+1)
         b = np.zeros(self.intervals+1)
         a[0] = 1
@@ -173,14 +186,14 @@ def _init_params(self):
             a[i] = (1 - i * step_size)
             b[i] = (i * step_size)
         all_params = np.zeros(self.stride * (self.intervals + 1))
-        if self.hnav is None:
-            # order [a[0],b[0],a[1],b[1],...]
+        if self.h_nav is None:
+            # order [a[0], b[0], a[1], b[1],...]
             all_params = np.dstack((a, b)).flatten()
         else:
             c = np.zeros(self.intervals+1)
             c[0:self.intervals//2] = b[0:self.intervals//2]
             c[self.intervals//2:self.intervals+1] = a[self.intervals//2:self.intervals+1]
-            # order [a[0],b[0],c[0],a[1],b[1],c[1],...]
+            # order [a[0], b[0], c[0], a[1], b[1], c[1],...]
             all_params = np.dstack((a, b, c)).flatten()
         return all_params
 
@@ -190,24 +203,21 @@ def prepare_reference_state(self):
                                      up_then_down=self.up_then_down, spin=self.spin)
 
     def build_circuit(self, var_params=None):
-        """Build the VSQS circuit by successive first- or second-order trotterizations of hini, hfin and possibly hnav"""
+        """Build the VSQS circuit by successive first- or second-order trotterizations of h_init, h_final and possibly h_nav"""
         reference_state_circuit = self.prepare_reference_state() if self.reference_state is None else self.reference_state
         self.var_params = self.set_var_params(var_params)
 
-        vsqs_circuit = get_exponentiated_qubit_operator_circuit(self.hini, time=self.dt, trotter_order=self.trotter_order, pauli_order=self.hini_list)
+        vsqs_circuit = get_exponentiated_qubit_operator_circuit(self.h_init, time=self.dt, trotter_order=self.trotter_order, pauli_order=self.h_init_list)
         for i in range(self.intervals-1):
-            vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hini, time=self.var_params[i * self.stride] * self.dt, variational=True,
-                                                                     trotter_order=self.trotter_order, pauli_order=self.hini_list)
-            vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hfin, time=self.var_params[i * self.stride + 1] * self.dt, variational=True,
-                                                                     trotter_order=self.trotter_order, pauli_order=self.hfin_list)
-            if self.hnav is not None:
-                vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hnav, time=self.var_params[i * self.stride + 2] * self.dt, variational=True,
-                                                                         trotter_order=self.trotter_order, pauli_order=self.hnav_list)
-        vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hfin, time=self.dt, trotter_order=self.trotter_order, pauli_order=self.hfin_list)
-
-        if reference_state_circuit.size != 0:
-            self.circuit = reference_state_circuit + vsqs_circuit
-        else:
-            self.circuit = vsqs_circuit
+            vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.h_init, time=self.var_params[i * self.stride] * self.dt, variational=True,
+                                                                     trotter_order=self.trotter_order, pauli_order=self.h_init_list)
+            vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.h_final, time=self.var_params[i * self.stride + 1] * self.dt, variational=True,
+                                                                     trotter_order=self.trotter_order, pauli_order=self.h_final_list)
+            if self.h_nav is not None:
+                vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.h_nav, time=self.var_params[i * self.stride + 2] * self.dt, variational=True,
+                                                                         trotter_order=self.trotter_order, pauli_order=self.h_nav_list)
+        vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.h_final, time=self.dt, trotter_order=self.trotter_order, pauli_order=self.h_final_list)
+
+        self.circuit = reference_state_circuit + vsqs_circuit if reference_state_circuit.size != 0 else vsqs_circuit
 
         return self.circuit
diff --git a/tangelo/toolboxes/molecular_computation/molecule.py b/tangelo/toolboxes/molecular_computation/molecule.py
index e152145b7..0e00c1389 100644
--- a/tangelo/toolboxes/molecular_computation/molecule.py
+++ b/tangelo/toolboxes/molecular_computation/molecule.py
@@ -21,9 +21,10 @@
 import numpy as np
 from pyscf import gto, scf, ao2mo
 import openfermion
-import openfermionpyscf
+import openfermion.ops.representations as reps
 from openfermionpyscf import run_pyscf
-from openfermion.ops.representations.interaction_operator import get_active_space_integrals
+from openfermion.chem.molecular_data import spinorb_from_spatial
+from openfermion.ops.representations.interaction_operator import get_active_space_integrals as of_get_active_space_integrals
 
 from tangelo.toolboxes.molecular_computation.frozen_orbitals import get_frozen_core
 from tangelo.toolboxes.qubit_mappings.mapping_transform import get_fermion_operator
@@ -275,14 +276,11 @@ def _get_fermionic_hamiltonian(self):
             FermionOperator: Self-explanatory.
         """
 
-        occupied_indices = self.frozen_occupied
-        active_indices = self.active_mos
+        core_constant, one_body_integrals, two_body_integrals = self.get_active_space_integrals()
 
-        of_molecule = self.to_openfermion(self.basis)
-        of_molecule = run_pyscf(of_molecule, run_scf=True, run_mp2=False,
-                                run_cisd=False, run_ccsd=False, run_fci=False)
+        one_body_coefficients, two_body_coefficients = spinorb_from_spatial(one_body_integrals, two_body_integrals)
 
-        molecular_hamiltonian = of_molecule.get_molecular_hamiltonian(occupied_indices, active_indices)
+        molecular_hamiltonian = reps.InteractionOperator(core_constant, one_body_coefficients, 1 / 2 * two_body_coefficients)
 
         return get_fermion_operator(molecular_hamiltonian)
 
@@ -379,6 +377,28 @@ def energy_from_rdms(self, one_rdm, two_rdm):
             float: Molecular energy.
         """
 
+        core_constant, one_electron_integrals, two_electron_integrals = self.get_active_space_integrals()
+
+        # PQRS convention in openfermion:
+        # h[p,q]=\int \phi_p(x)* (T + V_{ext}) \phi_q(x) dx
+        # h[p,q,r,s]=\int \phi_p(x)* \phi_q(y)* V_{elec-elec} \phi_r(y) \phi_s(x) dxdy
+        # The convention is not the same with PySCF integrals. So, a change is
+        # reverse back after performing the truncation for frozen orbitals
+        two_electron_integrals = two_electron_integrals.transpose(0, 3, 1, 2)
+
+        # Computing the total energy from integrals and provided RDMs.
+        e = core_constant + np.sum(one_electron_integrals * one_rdm) + 0.5*np.sum(two_electron_integrals * two_rdm)
+
+        return e.real
+
+    def get_active_space_integrals(self):
+        """Computes core constant, one_body, and two-body coefficients with frozen orbitals folded into one-body coefficients
+        and core constant
+
+        Returns:
+            (float, array, array): (core_constant, one_body coefficients, two_body coefficients)
+        """
+
         # Pyscf molecule to get integrals.
         pyscf_mol = self.to_pyscf(self.basis)
 
@@ -396,19 +416,14 @@ def energy_from_rdms(self, one_rdm, two_rdm):
         # h[p,q]=\int \phi_p(x)* (T + V_{ext}) \phi_q(x) dx
         # h[p,q,r,s]=\int \phi_p(x)* \phi_q(y)* V_{elec-elec} \phi_r(y) \phi_s(x) dxdy
         # The convention is not the same with PySCF integrals. So, a change is
-        # made and reverse back after performing the truncation for frozen
-        # orbitals.
+        # made before performing the truncation for frozen orbitals.
         two_electron_integrals = two_electron_integrals.transpose(0, 2, 3, 1)
-        core_offset, one_electron_integrals, two_electron_integrals = get_active_space_integrals(one_electron_integrals,
-                                                                                                 two_electron_integrals,
-                                                                                                 self.frozen_occupied,
-                                                                                                 self.active_mos)
-        two_electron_integrals = two_electron_integrals.transpose(0, 3, 1, 2)
+        core_offset, one_electron_integrals, two_electron_integrals = of_get_active_space_integrals(one_electron_integrals,
+                                                                                                    two_electron_integrals,
+                                                                                                    self.frozen_occupied,
+                                                                                                    self.active_mos)
 
         # Adding frozen electron contribution to core constant.
         core_constant += core_offset
 
-        # Computing the total energy from integrals and provided RDMs.
-        e = core_constant + np.sum(one_electron_integrals * one_rdm) + 0.5*np.sum(two_electron_integrals * two_rdm)
-
-        return e.real
+        return core_constant, one_electron_integrals, two_electron_integrals

From e7b388d9988c8894ace830cc788e1927a13384fc Mon Sep 17 00:00:00 2001
From: James Brown <james.brown@1qbit.com>
Date: Fri, 14 Jan 2022 13:42:39 -0500
Subject: [PATCH 09/10] style improvements

---
 tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py | 4 ++--
 tangelo/toolboxes/ansatz_generator/vsqs.py            | 8 ++++----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
index b36e19826..c1d8985e6 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
@@ -13,9 +13,9 @@
 # limitations under the License.
 
 import unittest
+import os
 
 import numpy as np
-import os
 from openfermion import load_operator
 
 from tangelo.molecule_library import mol_H2_sto3g
@@ -95,7 +95,7 @@ def test_vsqs_H4_doublecation(self):
     def test_vsqs_H2_with_h_nav(self):
         """Verify closed-shell VSQS functionalities for H2 with navigator hamiltonian"""
         navigator_hamiltonian = (QubitOperator('X0 Y1', 0.03632537110234512) + QubitOperator('Y0 X1', 0.03632537110234512)
-                                 + QubitOperator('Y0', 2.e-5) + QubitOperator('Y1', 2.e-5))
+                              + QubitOperator('Y0', 2.e-5) + QubitOperator('Y1', 2.e-5))
 
         # Build qubit hamiltonian for energy evaluation
         qubit_hamiltonian = symmetry_conserving_bravyi_kitaev(mol_H2_sto3g.fermionic_hamiltonian, 4, 2, False, 0)
diff --git a/tangelo/toolboxes/ansatz_generator/vsqs.py b/tangelo/toolboxes/ansatz_generator/vsqs.py
index 0f4532a6e..708ff8f67 100644
--- a/tangelo/toolboxes/ansatz_generator/vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/vsqs.py
@@ -117,11 +117,11 @@ def _build_h_init(self, molecule):
         """Return the initial Hamiltonian (h_init) composed of the one-body terms derived from the diagonal of Fock
         matrix and one-body off-diagonal terms"""
         core_constant, h1, two_body = molecule.get_active_space_integrals()
-        dfock = np.diag(h1).copy()
+        diag_fock = np.diag(h1).copy()
         n_active_occupied = len(molecule.active_occupied)
         for j in range(molecule.n_active_mos):
             for i in range(n_active_occupied):
-                dfock[j] += 2*two_body[i, j, j, i] - 1*two_body[i, j, i, j]
+                diag_fock[j] += 2*two_body[i, j, j, i] - 1*two_body[i, j, i, j]
 
         hf_ferm = FermionOperator((), core_constant)
         for i in range(molecule.n_active_mos):
@@ -130,8 +130,8 @@ def _build_h_init(self, molecule):
                     hf_ferm += FermionOperator(((i * 2, 1), (j * 2, 0)), h1[i, j])
                     hf_ferm += FermionOperator(((i * 2 + 1, 1), (j * 2 + 1, 0)), h1[i, j])
                 else:
-                    hf_ferm += FermionOperator(((i * 2, 1), (j * 2, 0)), dfock[i])
-                    hf_ferm += FermionOperator(((i * 2 + 1, 1), (j * 2 + 1, 0)), dfock[j])
+                    hf_ferm += FermionOperator(((i * 2, 1), (j * 2, 0)), diag_fock[i])
+                    hf_ferm += FermionOperator(((i * 2 + 1, 1), (j * 2 + 1, 0)), diag_fock[j])
         return fermion_to_qubit_mapping(hf_ferm, mapping=self.mapping, n_spinorbitals=self.n_spinorbitals, n_electrons=self.n_electrons,
                                         up_then_down=self.up_then_down, spin=self.spin)
 

From 44f100b61256eacf8b12f0dc332231096d09d8b1 Mon Sep 17 00:00:00 2001
From: ValentinS4t1qbit <41597680+ValentinS4t1qbit@users.noreply.github.com>
Date: Fri, 14 Jan 2022 15:07:19 -0800
Subject: [PATCH 10/10] Using the shortcut in init

---
 tangelo/algorithms/variational/vqe_solver.py | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/tangelo/algorithms/variational/vqe_solver.py b/tangelo/algorithms/variational/vqe_solver.py
index 9b9b83e7a..6f8180696 100644
--- a/tangelo/algorithms/variational/vqe_solver.py
+++ b/tangelo/algorithms/variational/vqe_solver.py
@@ -29,14 +29,7 @@
 from tangelo.toolboxes.operators import count_qubits, FermionOperator, qubitop_to_qubitham
 from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
 from tangelo.toolboxes.ansatz_generator.ansatz import Ansatz
-from tangelo.toolboxes.ansatz_generator.uccsd import UCCSD
-from tangelo.toolboxes.ansatz_generator.rucc import RUCC
-from tangelo.toolboxes.ansatz_generator.hea import HEA
-from tangelo.toolboxes.ansatz_generator.upccgsd import UpCCGSD
-from tangelo.toolboxes.ansatz_generator.qmf import QMF
-from tangelo.toolboxes.ansatz_generator.qcc import QCC
-from tangelo.toolboxes.ansatz_generator.vsqs import VSQS
-from tangelo.toolboxes.ansatz_generator.variational_circuit import VariationalCircuitAnsatz
+from tangelo.toolboxes.ansatz_generator import UCCSD, RUCC, HEA, UpCCGSD, QMF, QCC, VSQS, VariationalCircuitAnsatz
 from tangelo.toolboxes.ansatz_generator.penalty_terms import combined_penalty
 from tangelo.toolboxes.post_processing.bootstrapping import get_resampled_frequencies
 from tangelo.toolboxes.ansatz_generator.fermionic_operators import number_operator, spinz_operator, spin2_operator