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..eefabdc4f 100644
--- a/tangelo/algorithms/variational/tests/test_vqe_solver.py
+++ b/tangelo/algorithms/variational/tests/test_vqe_solver.py
@@ -137,6 +137,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": False, "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
+        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,
+                                          "h_init": h_init, "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..6f8180696 100644
--- a/tangelo/algorithms/variational/vqe_solver.py
+++ b/tangelo/algorithms/variational/vqe_solver.py
@@ -29,13 +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.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
@@ -50,6 +44,7 @@ class BuiltInAnsatze(Enum):
     UpCCGSD = 4
     QMF = 5
     QCC = 6
+    VSQS = 7
 
 
 class VQESolver:
@@ -186,12 +181,19 @@ 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 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)
         elif not isinstance(self.ansatz, Ansatz):
             raise TypeError(f"Invalid ansatz dataype. Expecting a custom Ansatz (Ansatz class).")
 
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 a9a2a4c16..825f7f587 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,20 @@ 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 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 to generate the exponential.
 
     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 = pauli_order.copy()
+    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 +126,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/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
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..c1d8985e6
--- /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 os
+
+import numpy as np
+from openfermion import load_operator
+
+from tangelo.molecule_library import mol_H2_sto3g
+from tangelo.toolboxes.operators import QubitOperator
+from tangelo.toolboxes.qubit_mappings import jordan_wigner, symmetry_conserving_bravyi_kitaev
+from tangelo.toolboxes.ansatz_generator 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).
+        """
+
+        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, h_init=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_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))
+
+        # 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,
+                           h_nav=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..708ff8f67
--- /dev/null
+++ b/tangelo/toolboxes/ansatz_generator/vsqs.py
@@ -0,0 +1,223 @@
+# 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."""
+
+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 h_init QubitOperator.
+
+    Args:
+        molecule (SecondQuantizedMolecule): The molecular system. Default: None
+        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
+        reference_state (Circuit): The reference state for the propagation as defined by a Circuit. Mandatory if supplying
+            a qubit_hamiltonian. Default: None
+        h_init (QubitOperator): The initial qubit Hamiltonian that corresponds to the reference state. Mandatory if supplying
+            a qubit_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, h_init=None, h_nav=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
+        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(h_init, ofQubitOperator):
+                raise ValueError("When providing a qubit hamiltonian, an initial qubit Hamiltonian must also be provided")
+            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
+        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.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_h_nav = 0
+        else:
+            if isinstance(self.h_nav, ofQubitOperator):
+                self.stride = 3
+                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_h_init + self.n_h_final + self.n_h_nav) * self.trotter_order
+
+        self.var_params = None
+        self.circuit = None
+
+    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()
+        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):
+                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):
+            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)), 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)
+
+    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.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_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_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 = np.zeros(self.intervals+1)
+        b = np.zeros(self.intervals+1)
+        a[0] = 1
+        b[self.intervals] = 1
+        step_size = 1/self.intervals
+        for i in range(1, self.intervals):
+            a[i] = (1 - i * step_size)
+            b[i] = (i * step_size)
+        all_params = np.zeros(self.stride * (self.intervals + 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],...]
+            all_params = np.dstack((a, b, c)).flatten()
+        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- 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.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.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