diff --git a/.github/workflows/continuous_integration.yml b/.github/workflows/continuous_integration.yml
index dc34bc80..9b16e306 100755
--- a/.github/workflows/continuous_integration.yml
+++ b/.github/workflows/continuous_integration.yml
@@ -73,7 +73,7 @@ jobs:
 
     - name: Install pyscf
       run: |
-        python -m pip install pyscf==2.4.0
+        python -m pip install pyscf
         python -m pip install git+https://github.com/pyscf/semiempirical
       if: always()
 
diff --git a/.github/workflows/deploy_docs.yml b/.github/workflows/deploy_docs.yml
new file mode 100644
index 00000000..aedaf16e
--- /dev/null
+++ b/.github/workflows/deploy_docs.yml
@@ -0,0 +1,42 @@
+name: Build & deploy sphinx document
+
+on:
+  workflow_dispatch:
+  push:
+    branches:
+      - main
+
+jobs:
+  build-deploy:
+    name: Build & deploy sphinx document
+    runs-on: ubuntu-latest
+    permissions:
+      contents: write
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+
+      - name: Set up python
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+
+      - name: Install Dependencies
+        run: |
+          python -m pip install sphinx sphinx_rtd_theme nbsphinx
+          python -m pip install .
+
+      - name: Run build command
+        run: |
+          sphinx-apidoc -o ./docs/source ./tangelo
+          sphinx-build -M html ./docs/source ./docs/build -E
+
+      - name: Deploy to GitHub Pages
+        uses: peaceiris/actions-gh-pages@v4
+        with:
+          publish_branch: gh-pages
+          github_token: ${{ secrets.GITHUB_TOKEN }}
+          publish_dir: ./docs/build/html
+          force_orphan: true
+          user_name: 'github-actions[bot]'
+          user_email: 'github-actions[bot]@users.noreply.github.com'
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 00000000..8cd6c408
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+**/__pycache__
+.pytest_cache
+**egg-info
+**/.DS_Store
+**/build
diff --git a/tangelo/algorithms/classical/mp2_solver.py b/tangelo/algorithms/classical/mp2_solver.py
index 90f0fecd..2e9566b1 100644
--- a/tangelo/algorithms/classical/mp2_solver.py
+++ b/tangelo/algorithms/classical/mp2_solver.py
@@ -83,7 +83,14 @@ def simulate(self):
         if self.uhf:
             self.mp2_fragment = self.mp.UMP2(self.mean_field, frozen=self.frozen)
         else:
-            self.mp2_fragment = self.mp.RMP2(self.mean_field, frozen=self.frozen)
+            import pyscf
+            if pyscf.__version__ == '2.5.0' and self.mean_field.istype('ROHF'):
+                mf = self.mean_field
+                mf = mf.remove_soscf()
+                mf = mf.to_uhf()
+                self.mp2_fragment = self.mp.UMP2(mf, frozen=self.frozen, mo_coeff=mf.mo_coeff, mo_occ=None)
+            else:
+                self.mp2_fragment = self.mp.RMP2(self.mean_field, frozen=self.frozen)
 
         self.mp2_fragment.verbose = 0
         _, self.mp2_t2 = self.mp2_fragment.kernel()
diff --git a/tangelo/algorithms/variational/adapt_vqe_solver.py b/tangelo/algorithms/variational/adapt_vqe_solver.py
index d275650f..7e4fe1b4 100644
--- a/tangelo/algorithms/variational/adapt_vqe_solver.py
+++ b/tangelo/algorithms/variational/adapt_vqe_solver.py
@@ -319,7 +319,7 @@ def choose_operator(self, gradients, tolerance=1e-3):
         max_partial = gradients[sorted_op_indices[-1]]
 
         if self.verbose:
-            print(f"LARGEST PARTIAL DERIVATIVE: {max_partial :4E}")
+            print(f"LARGEST PARTIAL DERIVATIVE: {max_partial:4E}")
 
         return [sorted_op_indices[-1]] if max_partial >= tolerance else []
 
diff --git a/tangelo/linq/circuit.py b/tangelo/linq/circuit.py
index c51ea165..51dbc3ab 100644
--- a/tangelo/linq/circuit.py
+++ b/tangelo/linq/circuit.py
@@ -24,6 +24,7 @@
 import warnings
 
 import numpy as np
+import sympy as sp
 from cirq.contrib.svg import SVGCircuit
 
 from tangelo.linq import Gate
@@ -180,14 +181,26 @@ def applied_gates(self):
         return self._applied_gates if "CMEASURE" in self.counts else self._gates
 
     def draw(self):
-        """Method to output a prettier version of the circuit for use in jupyter notebooks that uses cirq SVGCircuit"""
-        # circular import
+        """Method to output a prettier version of the circuit
+        for use in jupyter notebooks that uses cirq SVGCircuit"""
         from tangelo.linq.translator.translate_cirq import translate_c_to_cirq
-        cirq_circ = translate_c_to_cirq(self)
-        # Remove identity gates that are added in translate_c_to_cirq (to ensure all qubits are initialized) before drawing.
+        circuit_copy = self.copy()
+        for gate in circuit_copy._gates:
+            if gate.parameter and isinstance(gate.parameter, str):
+                gate.parameter = self._string_to_sympy(gate)
+
+        cirq_circ = translate_c_to_cirq(circuit_copy)
         cirq_circ.__delitem__(0)
         return SVGCircuit(cirq_circ)
 
+    def _string_to_sympy(self, gate):
+        """Convert a gate parameter (type string) to a sympy symbol"""
+        try:
+            return sp.symbols(gate.parameter)
+        except Exception as e:
+            print(f"Error converting {gate.parameter} to sympy symbol: {e}")
+            return gate.parameter
+
     def copy(self):
         """Return a deepcopy of circuit"""
         return Circuit(copy.deepcopy(self._gates), n_qubits=self._qubits_simulated, name=self.name, cmeasure_control=copy.deepcopy(self._cmeasure_control))
@@ -427,6 +440,13 @@ def finalize_cmeasure_control(self):
         if isinstance(self._cmeasure_control, ClassicalControl):
             self._cmeasure_control.finalize()
 
+    def fix_variational_parameters(self):
+        """Fix all variational parameters in this circuit, making the corresponding gates non-variational."""
+
+        for gate in self._variational_gates:
+            gate.is_variational = False
+        self._variational_gates = []
+
 
 def stack(*circuits):
     """ Take list of circuits as input, and stack them (e.g concatenate them along the
diff --git a/tangelo/linq/tests/test_translator_circuit.py b/tangelo/linq/tests/test_translator_circuit.py
index 7e3f0de8..3c103168 100644
--- a/tangelo/linq/tests/test_translator_circuit.py
+++ b/tangelo/linq/tests/test_translator_circuit.py
@@ -180,6 +180,25 @@ def test_qiskit(self):
         sim_results = qiskit_simulator.run(translated_circuit).result()
         np.testing.assert_array_almost_equal(sim_results.get_statevector(translated_circuit), reference_big_msq, decimal=6)
 
+    @unittest.skipIf("qiskit" not in installed_backends, "Test Skipped: Backend not available \n")
+    def test_qiskit_multi_control(self):
+        from qiskit import QuantumCircuit
+
+        for g in ["CRX", "CRY", "CRZ", "CPHASE"]:
+            c = Circuit([Gate(g, 2, control=[0, 1], parameter=1.)])
+            c_qiskit = translate_c(c, target="qiskit")
+            q = QuantumCircuit(3)
+            if g == "CRX":
+                q.mcrx(1., [0, 1], 2)
+            elif g == "CRY":
+                q.mcry(1., [0, 1], 2)
+            elif g == "CRZ":
+                q.mcrz(1., [0, 1], 2)
+            elif g == "CPHASE":
+                q.mcp(1., [0, 1], 2)
+
+            assert c_qiskit.data == q.data
+
     @unittest.skipIf("cirq" not in installed_backends, "Test Skipped: Backend not available \n")
     def test_cirq(self):
         """
diff --git a/tangelo/linq/translator/translate_qiskit.py b/tangelo/linq/translator/translate_qiskit.py
index dd662f6a..6a7df4a6 100644
--- a/tangelo/linq/translator/translate_qiskit.py
+++ b/tangelo/linq/translator/translate_qiskit.py
@@ -55,6 +55,10 @@ def get_qiskit_gates():
     GATE_QISKIT["CRX"] = qiskit.QuantumCircuit.crx
     GATE_QISKIT["CRY"] = qiskit.QuantumCircuit.cry
     GATE_QISKIT["CRZ"] = qiskit.QuantumCircuit.crz
+    GATE_QISKIT["MCRX"] = qiskit.QuantumCircuit.mcrx
+    GATE_QISKIT["MCRY"] = qiskit.QuantumCircuit.mcry
+    GATE_QISKIT["MCRZ"] = qiskit.QuantumCircuit.mcrz
+    GATE_QISKIT["MCPHASE"] = qiskit.QuantumCircuit.mcp
     GATE_QISKIT["CNOT"] = qiskit.QuantumCircuit.cx
     GATE_QISKIT["SWAP"] = qiskit.QuantumCircuit.swap
     GATE_QISKIT["XX"] = qiskit.QuantumCircuit.rxx
@@ -95,14 +99,17 @@ def translate_c_to_qiskit(source_circuit: Circuit, save_measurements=False, no_c
     # Maps the gate information properly. Different for each backend (order, values)
     for gate in source_circuit._gates:
         if gate.control is not None:
-            if len(gate.control) > 1:
+            if (len(gate.control) > 1) and (gate.name not in {"CRX", "CRY", "CRZ", "CPHASE"}):
                 raise ValueError('Multi-controlled gates not supported with qiskit. Gate {gate.name} with controls {gate.control} is not allowed')
         if gate.name in {"H", "Y", "X", "Z", "S", "T"}:
             (GATE_QISKIT[gate.name])(target_circuit, gate.target[0])
         elif gate.name in {"RX", "RY", "RZ", "PHASE"}:
             (GATE_QISKIT[gate.name])(target_circuit, gate.parameter, gate.target[0])
         elif gate.name in {"CRX", "CRY", "CRZ", "CPHASE"}:
-            (GATE_QISKIT[gate.name])(target_circuit, gate.parameter, gate.control[0], gate.target[0])
+            if len(gate.control) > 1:
+                (GATE_QISKIT["M" + gate.name])(target_circuit, gate.parameter, gate.control, gate.target[0])
+            else:
+                (GATE_QISKIT[gate.name])(target_circuit, gate.parameter, gate.control[0], gate.target[0])
         elif gate.name in {"CNOT", "CH", "CX", "CY", "CZ"}:
             (GATE_QISKIT[gate.name])(target_circuit, gate.control[0], gate.target[0])
         elif gate.name in {"SWAP"}:
diff --git a/tangelo/toolboxes/ansatz_generator/adapt_ansatz.py b/tangelo/toolboxes/ansatz_generator/adapt_ansatz.py
index 31137631..d0cd42f8 100644
--- a/tangelo/toolboxes/ansatz_generator/adapt_ansatz.py
+++ b/tangelo/toolboxes/ansatz_generator/adapt_ansatz.py
@@ -106,7 +106,11 @@ def update_var_params(self, var_params):
 
     def prepare_reference_state(self):
         """Prepare a circuit generating the HF reference state."""
-        if self.reference_state.upper() == "HF":
+        if isinstance(self.reference_state, Circuit):
+            ref_circuit = self.reference_state.copy()
+            ref_circuit.fix_variational_parameters()
+            return ref_circuit
+        elif self.reference_state.upper() == "HF":
             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)
         else:
diff --git a/tangelo/toolboxes/ansatz_generator/hea.py b/tangelo/toolboxes/ansatz_generator/hea.py
index 5520a83c..b8039f49 100644
--- a/tangelo/toolboxes/ansatz_generator/hea.py
+++ b/tangelo/toolboxes/ansatz_generator/hea.py
@@ -44,18 +44,25 @@ class HEA(Ansatz):
         n_qubits (int) : The number of qubits in the ansatz.
             Default, None.
         n_electrons (int) : Self-explanatory.
-        reference_state (str): "HF": Hartree-Fock reference state. "zero": for
-            no reference state. Default: "HF".
+        reference_state (str, Circuit): "HF": Hartree-Fock reference state. "zero": for
+            no reference state. Can also be a Circuit object, in which case a copy of
+            circuit with variational parameters fixed is used. Default: "HF".
         """
 
     def __init__(self, molecule=None, mapping="jw", up_then_down=False,
                  n_layers=2, rot_type="euler", n_qubits=None, n_electrons=None,
                  spin=None, reference_state="HF"):
 
-        if not (bool(molecule) ^ (bool(n_qubits) and (bool(n_electrons) | (reference_state == "zero")))):
-            raise ValueError(f"A molecule OR qubit + electrons number must be "
-                             "provided when instantiating the HEA with the HF reference state. "
-                             "For reference_state='zero', only the number of qubits is needed.")
+        # Ensure sufficient parameters are passed to instantiate this HEA with the given reference state
+        if isinstance(reference_state, Circuit):
+            if not bool(molecule) and not bool(n_qubits):
+                raise ValueError('Either a molecule or a qubit number must be specified to instantiate a HEA with a Circuit reference state.')
+        elif reference_state == 'HF':
+            if not bool(molecule) and not (bool(n_qubits) and bool(n_electrons)):
+                raise ValueError('Either a molecule or a qubit number + electron number must be specified to instantiate a HEA with the "HF" reference state.')
+        elif reference_state == 'zero':
+            if not bool(molecule) and not bool(n_qubits):
+                raise ValueError('Either a molecule or a qubit number must be specified to instantiate a HEA with the "zero" reference state.')
 
         if n_qubits:
             self.n_qubits = n_qubits
@@ -124,6 +131,12 @@ def set_var_params(self, var_params=None):
 
     def prepare_reference_state(self):
         """Prepare a circuit generating the HF reference state."""
+
+        if isinstance(self.reference_state, Circuit):
+            ref_circuit = self.reference_state.copy()
+            ref_circuit.fix_variational_parameters()
+            return ref_circuit
+
         if self.reference_state not in self.supported_reference_state:
             raise ValueError(f"{self.reference_state} not in supported reference state methods of:{self.supported_reference_state}")
 
diff --git a/tangelo/toolboxes/ansatz_generator/ilc.py b/tangelo/toolboxes/ansatz_generator/ilc.py
index c7824c71..5e5ca37e 100644
--- a/tangelo/toolboxes/ansatz_generator/ilc.py
+++ b/tangelo/toolboxes/ansatz_generator/ilc.py
@@ -73,9 +73,10 @@ class ILC(Ansatz):
         max_ilc_gens (int or None): Maximum number of generators allowed in the ansatz. If None,
             one generator from each DIS group is selected. If int, then min(|DIS|, max_ilc_gens)
             generators are selected in order of decreasing |dEILC/dtau|. Default, None.
-        reference_state (string): The reference state id for the ansatz. The
-            supported reference states are stored in the supported_reference_state
-            attributes. Default, "HF".
+        reference_state (string, Circuit): The reference state id for the ansatz. If a Circuit object
+            is passed, then a copy of this circuit overrides the qmf_circuit and its variational
+            parameters override qmf_var_params. The supported string reference states are stored in the
+            supported_reference_state attributes. Default, "HF".
     """
 
     def __init__(self, molecule, mapping="jw", up_then_down=False, acs=None,
@@ -114,14 +115,24 @@ def __init__(self, molecule, mapping="jw", up_then_down=False, acs=None,
                                                       self.n_spinorbitals, self.n_electrons,
                                                       self.up_then_down, self.spin)
 
+        # If a circuit is supplied as the reference state use this as the QMF circuit
+        # while retaining all variational parameters:
+        if isinstance(reference_state, Circuit):
+            self.qmf_circuit = reference_state.copy()
+            self.qmf_var_params = [
+                gate.parameter for gate in reference_state._variational_gates
+            ]
+            self.n_qmf_params = len(self.qmf_var_params)
+        else:
+            self.qmf_circuit = qmf_circuit
+            self.n_qmf_params = 2 * self.n_qubits
+
         self.qmf_var_params = np.array(qmf_var_params) if isinstance(qmf_var_params, list) else qmf_var_params
         if not isinstance(self.qmf_var_params, np.ndarray):
             self.qmf_var_params = init_qmf_from_hf(self.n_spinorbitals, self.n_electrons,
                                                    self.mapping, self.up_then_down, self.spin)
         if self.qmf_var_params.size != 2 * self.n_qubits:
             raise ValueError("The number of QMF variational parameters must be 2 * n_qubits.")
-        self.n_qmf_params = 2 * self.n_qubits
-        self.qmf_circuit = qmf_circuit
 
         self.acs = acs
         self.ilc_tau_guess = ilc_tau_guess
@@ -172,6 +183,7 @@ def set_var_params(self, var_params=None):
             # Initialize ILC parameters by matrix diagonalization (see Appendix B, Refs. 1 & 2).
             elif var_params == "diag":
                 initial_var_params = get_ilc_params_by_diag(self.qubit_ham, self.acs, self.qmf_var_params)
+
             # Insert the QMF variational parameters at the beginning.
             initial_var_params = np.concatenate((self.qmf_var_params, initial_var_params))
         else:
@@ -187,11 +199,18 @@ def prepare_reference_state(self):
         wavefunction with HF, multi-reference state, etc). These preparations must be consistent
         with the transform used to obtain the qubit operator. """
 
+        # Note: because reference state parameters are needed in this ansatz, the reference state circuit
+        # is simply copied and stored in self.qmf_circuit:
+        if isinstance(self.reference_state, Circuit):
+            return self.qmf_circuit
+
         if self.reference_state not in self.supported_reference_state:
             raise ValueError(f"Only supported reference state methods are: "
                              f"{self.supported_reference_state}.")
+
         if self.reference_state == "HF":
             reference_state_circuit = get_qmf_circuit(self.qmf_var_params, True)
+
         return reference_state_circuit
 
     def build_circuit(self, var_params=None):
diff --git a/tangelo/toolboxes/ansatz_generator/puccd.py b/tangelo/toolboxes/ansatz_generator/puccd.py
index a75136e8..93de60f1 100644
--- a/tangelo/toolboxes/ansatz_generator/puccd.py
+++ b/tangelo/toolboxes/ansatz_generator/puccd.py
@@ -38,8 +38,9 @@ class pUCCD(Ansatz):
 
     Args:
         molecule (SecondQuantizedMolecule): Self-explanatory.
-        reference_state (string): String refering to an initial state.
-            Default: "HF".
+        reference_state (string, Circuit): String refering to an initial state.
+            Can also be a Circuit object, in which case a copy of
+            circuit with variational parameters fixed is used. Default: "HF".
     """
 
     def __init__(self, molecule, reference_state="HF"):
@@ -100,6 +101,11 @@ def prepare_reference_state(self):
         the qubit operator.
         """
 
+        if isinstance(self.reference_state, Circuit):
+            reference_state_circuit = self.reference_state.copy()
+            reference_state_circuit.fix_variational_parameters()
+            return reference_state_circuit
+
         if self.reference_state not in self.supported_reference_state:
             raise ValueError(f"Only supported reference state methods are:{self.supported_reference_state}")
 
diff --git a/tangelo/toolboxes/ansatz_generator/qcc.py b/tangelo/toolboxes/ansatz_generator/qcc.py
index c56c0c31..c83fd977 100644
--- a/tangelo/toolboxes/ansatz_generator/qcc.py
+++ b/tangelo/toolboxes/ansatz_generator/qcc.py
@@ -78,9 +78,10 @@ class QCC(Ansatz):
         max_qcc_gens (int or None): Maximum number of generators allowed in the ansatz. If None,
             one generator from each DIS group is selected. If int, then min(|DIS|, max_qcc_gens)
             generators are selected in order of decreasing |dEQCC/dtau|. Default, None.
-        reference_state (string): The reference state id for the ansatz. The
-            supported reference states are stored in the supported_reference_state
-            attributes. Default, "HF".
+        reference_state (string): The reference state id for the ansatz. If a Circuit object
+            is passed, then a copy of this circuit overrides the qmf_circuit and its variational
+            parameters override qmf_var_params. The supported reference states are stored in the
+            supported_reference_state attributes. Default, "HF".
     """
 
     def __init__(self, molecule, mapping="jw", up_then_down=False, dis=None,
@@ -121,14 +122,24 @@ def __init__(self, molecule, mapping="jw", up_then_down=False, dis=None,
                                                       self.n_spinorbitals, self.n_electrons,
                                                       self.up_then_down, self.spin)
 
+        # If a circuit is supplied as the reference state use this as the QMF circuit
+        # while retaining all variational parameters:
+        if isinstance(reference_state, Circuit):
+            self.qmf_circuit = reference_state.copy()
+            self.qmf_var_params = [
+                gate.parameter for gate in reference_state._variational_gates
+            ]
+            self.n_qmf_params = len(self.qmf_var_params)
+        else:
+            self.qmf_circuit = qmf_circuit
+            self.n_qmf_params = 2 * self.n_qubits
+
         self.qmf_var_params = np.array(qmf_var_params) if isinstance(qmf_var_params, list) else qmf_var_params
         if not isinstance(self.qmf_var_params, np.ndarray):
             self.qmf_var_params = init_qmf_from_hf(self.n_spinorbitals, self.n_electrons,
                                                    self.mapping, self.up_then_down, self.spin)
         if self.qmf_var_params.size != 2 * self.n_qubits:
             raise ValueError("The number of QMF variational parameters must be 2 * n_qubits.")
-        self.n_qmf_params = 2 * self.n_qubits
-        self.qmf_circuit = qmf_circuit
 
         self.dis = dis
         self.qcc_tau_guess = qcc_tau_guess
@@ -191,6 +202,11 @@ def prepare_reference_state(self):
         wavefunction with HF, multi-reference state, etc). These preparations must be consistent
         with the transform used to obtain the qubit operator. """
 
+        # Note: because reference state parameters are needed in this ansatz, the reference state circuit
+        # is simply copied and stored in self.qmf_circuit:
+        if isinstance(self.reference_state, Circuit):
+            return self.qmf_circuit
+
         if self.reference_state not in self.supported_reference_state:
             raise ValueError(f"Only supported reference state methods are: "
                              f"{self.supported_reference_state}.")
diff --git a/tangelo/toolboxes/ansatz_generator/qmf.py b/tangelo/toolboxes/ansatz_generator/qmf.py
index f5091d05..2e3bb755 100755
--- a/tangelo/toolboxes/ansatz_generator/qmf.py
+++ b/tangelo/toolboxes/ansatz_generator/qmf.py
@@ -33,6 +33,7 @@
 
 import numpy as np
 
+from tangelo.linq import Circuit
 from tangelo.toolboxes.qubit_mappings.mapping_transform import get_qubit_number,\
                                                                fermion_to_qubit_mapping
 from tangelo.toolboxes.ansatz_generator.ansatz import Ansatz
@@ -69,9 +70,10 @@ class QMF(Ansatz):
             <N> = n_electrons, <S^2> = spin_z * (spin_z + 1), and <Sz> = spin_z, where
             spin_z = spin // 2. Key, value pairs are case sensitive and mu > 0.
             Default, {"init_params": "hf_state"}.
-        reference_state (string): The reference state id for the ansatz. The
-            supported reference states are stored in the supported_reference_state
-            attributes. Default, "HF".
+        reference_state (string, Circuit): The reference state id for the ansatz.  Can also be a
+            Circuit object, in which case a copy of circuit with variational parameters fixed is used.
+            The supported reference states are stored in the supported_reference_state attributes.
+            Default, "HF".
     """
 
     def __init__(self, molecule, mapping="jw", up_then_down=False, init_qmf=None, reference_state="HF"):
@@ -195,6 +197,11 @@ def prepare_reference_state(self):
         wavefunction with HF, multi-reference state, etc). These preparations must be consistent
         with the transform used to obtain the qubit operator. """
 
+        if isinstance(self.reference_state, Circuit):
+            ref_circuit = self.reference_state.copy()
+            ref_circuit.fix_variational_parameters()
+            return ref_circuit
+
         if self.reference_state not in self.supported_reference_state:
             raise ValueError(f"Only supported reference state methods are: "
                              f"{self.supported_reference_state}")
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_adapt_ansatz.py b/tangelo/toolboxes/ansatz_generator/tests/test_adapt_ansatz.py
index 53fb6efe..b34b4b66 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_adapt_ansatz.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_adapt_ansatz.py
@@ -18,6 +18,7 @@
 from tangelo.toolboxes.ansatz_generator.adapt_ansatz import ADAPTAnsatz
 from tangelo.toolboxes.operators import FermionOperator
 from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
+from tangelo.linq import Gate, Circuit
 
 f_op = FermionOperator("2^ 3^ 0 1") - FermionOperator("0^ 1^ 2 3")
 qu_op = fermion_to_qubit_mapping(f_op, "jw")
@@ -57,6 +58,51 @@ def test_adaptansatz_set_var_params(self):
         with self.assertRaises(ValueError):
             ansatz.set_var_params([1.999, 2.999])
 
+    def test_adaptansatz_reference_state_circuit(self):
+        """ Verify variational gate parameters are frozen in a reference state circuit"""
+
+        # Create simple reference circuit:
+        ref_circuit = Circuit([
+            Gate('RY', i, parameter=1.0, is_variational=True)
+            for i in range(4)
+        ])
+
+        # Create ADAPTAnsatz ansatz with reference circuit
+        adapt_ansatz = ADAPTAnsatz(n_spinorbitals=2, n_electrons=2, spin=0,
+                            ansatz_options=dict(reference_state=ref_circuit))
+
+        adapt_ansatz.build_circuit()
+
+        adapt_circ = adapt_ansatz.circuit
+        adapt_circ_gates = list(adapt_circ)
+
+        # Ensure gates are correctly prepended to circuit
+        self.assertEqual(adapt_circ.width, 4)
+        self.assertEqual(adapt_circ_gates[0].name, 'RY')
+
+        # Ensure reference circuit gates were correctly converted to
+        # non-variational gates
+        self.assertFalse(adapt_circ_gates[0].is_variational)
+
+        # Add qu_op to ansatz
+        adapt_ansatz.add_operator(qu_op)
+
+        # Check ansatz parameters
+        self.assertEqual(adapt_ansatz.n_var_params, 1)
+        self.assertEqual(adapt_ansatz._n_terms_operators, [8])
+
+        adapt_circ = adapt_ansatz.circuit
+        self.assertEqual(adapt_circ.width, 4)
+
+        # Ensure reference circuit gates were correctly converted to non-variational gates
+        # with the same name
+        ref_circ_gate_names = [ gate.name for gate in ref_circuit ]
+        adapt_circ_gate_names = [ gate.name for gate in adapt_circ ]
+        adapt_circ_gate_variationals = [ gate.is_variational for gate in adapt_circ ]
+
+        self.assertListEqual(ref_circ_gate_names, adapt_circ_gate_names[:ref_circuit.size])
+        self.assertTrue(not any(adapt_circ_gate_variationals[:ref_circuit.size]))
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_hea.py b/tangelo/toolboxes/ansatz_generator/tests/test_hea.py
index 630bf37d..70509430 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_hea.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_hea.py
@@ -20,8 +20,10 @@
 from tangelo.molecule_library import mol_H2_sto3g, mol_H4_doublecation_minao
 from tangelo.toolboxes.qubit_mappings import jordan_wigner
 from tangelo.toolboxes.ansatz_generator.hea import HEA
+from tangelo.toolboxes.operators import count_qubits
+from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit
 
-from tangelo.linq import get_backend
+from tangelo.linq import Circuit, Gate, get_backend
 
 # Initiate simulator
 sim = get_backend()
@@ -104,6 +106,67 @@ def test_hea_H4_doublecation(self):
         energy = sim.get_expectation_value(qubit_hamiltonian, hea_ansatz.circuit)
         self.assertAlmostEqual(energy, -0.795317, delta=1e-4)
 
+    def test_hea_circuit_variational_reference_state(self):
+        """ Verify variational gate parameters are frozen in a reference state circuit"""
+
+        # Create simple reference circuit:
+        ref_circuit = Circuit([
+            Gate('RY', 0, parameter=1.0, is_variational=True),
+            Gate('RY', 1, parameter=1.0, is_variational=True),
+            Gate('CX', 0, 1),
+        ])
+
+        # Create HEA ansatz with reference circuit
+        hea_ansatz = HEA(n_qubits=2, reference_state=ref_circuit)
+        hea_circ = hea_ansatz.build_circuit()
+
+        # Ensure gates are correctly prepended to circuit
+        self.assertEqual(hea_circ.width, 2)
+
+        # Ensure reference circuit gates were correctly converted
+        # to non-variational gates with the same name
+        hea_circ_gate_names = [ gate.name for gate in hea_circ ]
+        ref_circ_gate_names = [ gate.name for gate in ref_circuit ]
+        hea_circ_gate_variationals = [ gate.is_variational for gate in hea_circ ]
+
+        self.assertListEqual(ref_circ_gate_names, hea_circ_gate_names[:ref_circuit.size])
+        self.assertTrue(not any(hea_circ_gate_variationals[:ref_circuit.size]))
+
+    def test_hea_circuit_reference_state_H2(self):
+        """ Verify construction of H2 ansatz works using a circuit reference state."""
+        sim = get_backend()
+
+        # Build qubit hamiltonian
+        qubit_hamiltonian = jordan_wigner(mol_H2_sto3g.fermionic_hamiltonian)
+        n_qubits = count_qubits(qubit_hamiltonian)
+
+        # Construct reference circuit by hand
+        ref_h2_circuit = get_reference_circuit(
+                n_spinorbitals=n_qubits,
+                n_electrons=mol_H2_sto3g.n_electrons,
+                mapping='jw',
+                up_then_down=False,
+                spin=mol_H2_sto3g.spin)
+
+        # Build ansatz circuit
+        hea_ansatz = HEA(n_qubits=n_qubits, reference_state=ref_h2_circuit)
+        hea_ansatz.build_circuit()
+
+        params = [ 1.96262489e+00, -9.83505909e-01, -1.92659544e+00,  3.68638855e+00,
+                4.71410736e+00,  4.78247991e+00,  4.71258582e+00, -4.79077006e+00,
+                4.60613188e+00, -3.14130503e+00,  4.71232383e+00,  1.35715841e+00,
+                4.30998410e+00, -3.00415626e+00, -1.06784872e+00, -2.05119893e+00,
+                6.44114344e+00, -1.56358255e+00, -6.28254779e+00,  3.14118427e+00,
+                -3.10505551e+00,  3.15123780e+00, -3.64794717e+00,  1.09127829e+00,
+                4.67093656e-01, -1.19912860e-01,  3.99351728e-03, -1.57104046e+00,
+                1.56811666e+00, -3.14050540e+00,  4.71181097e+00,  1.57036595e+00,
+                -2.16414405e+00,  3.40295404e+00, -2.87986715e+00, -1.69054279e+00]
+
+        # Assert energy returned is as expected for given parameters
+        hea_ansatz.update_var_params(params)
+        energy = sim.get_expectation_value(qubit_hamiltonian, hea_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.1372661564779496, delta=1e-6)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_ilc.py b/tangelo/toolboxes/ansatz_generator/tests/test_ilc.py
index 299e34a4..6217abf6 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_ilc.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_ilc.py
@@ -22,6 +22,7 @@
 from tangelo.linq import get_backend
 from tangelo.toolboxes.ansatz_generator.ilc import ILC
 from tangelo.toolboxes.ansatz_generator._qubit_ilc import gauss_elim_over_gf2, get_ilc_params_by_diag
+from tangelo.toolboxes.ansatz_generator._qubit_mf import get_qmf_circuit
 from tangelo.toolboxes.operators.operators import QubitOperator
 from tangelo.molecule_library import mol_H2_sto3g, mol_H4_cation_sto3g
 
@@ -182,6 +183,33 @@ def test_ilc_h4_cation(self):
         energy = sim.get_expectation_value(qubit_hamiltonian, ilc_ansatz.circuit)
         self.assertAlmostEqual(energy_diag, energy, delta=1e-5)
 
+    def test_ilc_circuit_reference_state_H2(self):
+        """ Verify construction of H2 ansatz works using a circuit reference state."""
+
+        # Construct reference circuit
+        qmf_var_params = [3.14159265e+00,  3.14159265e+00, -7.59061327e-12,  0.]
+        ref_qmf_circuit = get_qmf_circuit(np.array(qmf_var_params), True)
+
+        # Specify the qubit operators from the anticommuting set (ACS) of ILC generators.
+        acs = [QubitOperator("Y0 X1")]
+        ilc_ansatz = ILC(mol_H2_sto3g, mapping="scbk", up_then_down=True, acs=acs,
+                            reference_state=ref_qmf_circuit)
+
+        # Build the ILC circuit, which is prepended by the qubit mean field (QMF) circuit.
+        ilc_ansatz.build_circuit()
+
+        # Get qubit hamiltonian for energy evaluation
+        qubit_hamiltonian = ilc_ansatz.qubit_ham
+
+        # The QMF and ILC parameters can both be specified; determined automatically otherwise.
+        ilc_var_params = [-1.12894599e-01]
+        var_params = qmf_var_params + ilc_var_params
+        ilc_ansatz.update_var_params(var_params)
+
+        # Assert energy returned is as expected for given parameters
+        energy = sim.get_expectation_value(qubit_hamiltonian, ilc_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.137270126, delta=1e-6)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_puccd.py b/tangelo/toolboxes/ansatz_generator/tests/test_puccd.py
index ce0f5eb0..a76839d8 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_puccd.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_puccd.py
@@ -18,6 +18,7 @@
 from tangelo.molecule_library import mol_H2_sto3g
 from tangelo.toolboxes.qubit_mappings.hcb import hard_core_boson_operator, boson_to_qubit_mapping
 from tangelo.toolboxes.ansatz_generator.puccd import pUCCD
+from tangelo.toolboxes.qubit_mappings.statevector_mapping import vector_to_circuit
 from tangelo.linq import get_backend
 
 
@@ -68,6 +69,31 @@ def test_puccd_H2(self):
         energy = sim.get_expectation_value(qubit_hamiltonian, puccd_ansatz.circuit)
         self.assertAlmostEqual(energy, -1.13727, delta=1e-4)
 
+    def test_puccd_circuit_reference_state_H2(self):
+        """ Verify construction of H2 ansatz works using a circuit reference state."""
+
+        sim = get_backend()
+
+        # Construct H2 reference state circuit
+        ref_circuit = vector_to_circuit([
+            int(i < (mol_H2_sto3g.n_active_electrons // 2))
+            for i in range(mol_H2_sto3g.n_active_mos)
+        ])
+
+        # Build circuit.
+        puccd_ansatz = pUCCD(mol_H2_sto3g, reference_state=ref_circuit)
+        puccd_ansatz.build_circuit()
+
+        # Build qubit hamiltonian for energy evaluation.
+        qubit_hamiltonian = boson_to_qubit_mapping(
+            hard_core_boson_operator(mol_H2_sto3g.fermionic_hamiltonian)
+        )
+
+        # Assert energy returned is as expected for given parameters.
+        puccd_ansatz.update_var_params([-0.22617753])
+        energy = sim.get_expectation_value(qubit_hamiltonian, puccd_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.13727, delta=1e-4)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_qcc.py b/tangelo/toolboxes/ansatz_generator/tests/test_qcc.py
index 1c2e7b11..87b1b8d7 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_qcc.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_qcc.py
@@ -22,6 +22,7 @@
 
 from tangelo.linq import get_backend
 from tangelo.toolboxes.ansatz_generator.qcc import QCC
+from tangelo.toolboxes.ansatz_generator._qubit_mf import get_qmf_circuit
 from tangelo.toolboxes.operators.operators import QubitOperator
 from tangelo.molecule_library import mol_H2_sto3g, mol_H4_cation_sto3g, mol_H4_doublecation_minao, mol_H4_sto3g_uhf_a1_frozen
 from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
@@ -154,6 +155,32 @@ def test_qmf_qcc_h4_uhf_ref(self):
         energy = sim.get_expectation_value(qu_op, qcc_ansatz.prepare_reference_state())
         self.assertAlmostEqual(energy, mol.mean_field.e_tot, delta=1e-6)
 
+    def test_qcc_circuit_reference_state_h2(self):
+        """ Verify construction of H2 ansatz works using a circuit reference state."""
+
+        qmf_var_params = [ 3.14159265e+00, -2.42743256e-08,  3.14159266e+00, -3.27162543e-08,
+                           3.08514545e-09,  3.08514545e-09,  3.08514545e-09,  3.08514545e-09]
+        qmf_ref_circuit = get_qmf_circuit(np.array(qmf_var_params), True)
+
+        # Build the ansatz with qmf reference circuit:
+        dis = [QubitOperator("Y0 X1 X2 X3")]
+        qcc_ansatz = QCC(mol_H2_sto3g, up_then_down=True, dis=dis, reference_state=qmf_ref_circuit)
+
+        # Build the QCC circuit, which is prepended by the qubit mean field (QMF) circuit.
+        qcc_ansatz.build_circuit()
+
+        # Get qubit hamiltonian for energy evaluation
+        qubit_hamiltonian = qcc_ansatz.qubit_ham
+
+        # The QMF and QCC parameters can both be specified; determined automatically otherwise.
+        qcc_var_params = [-2.26136280e-01]
+        var_params = qmf_var_params + qcc_var_params
+
+        # Assert energy returned is as expected for given parameters
+        qcc_ansatz.update_var_params(var_params)
+        energy = sim.get_expectation_value(qubit_hamiltonian, qcc_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.1372701746609022, delta=1e-6)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_qmf.py b/tangelo/toolboxes/ansatz_generator/tests/test_qmf.py
index c680703a..a5308e19 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_qmf.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_qmf.py
@@ -21,6 +21,7 @@
 from tangelo.linq import get_backend
 from tangelo.toolboxes.ansatz_generator.qmf import QMF
 from tangelo.molecule_library import mol_H2_sto3g, mol_H4_sto3g, mol_H4_cation_sto3g, mol_H4_sto3g_uhf_a1_frozen
+from tangelo.toolboxes.ansatz_generator._qubit_mf import get_qmf_circuit
 
 sim = get_backend()
 
@@ -151,6 +152,23 @@ def test_qmf_uhf_h4_cation(self):
         energy = sim.get_expectation_value(qubit_hamiltonian, qmf_ansatz.circuit)
         self.assertAlmostEqual(energy, mol_H4_sto3g_uhf_a1_frozen.mean_field.e_tot, delta=1e-6)
 
+    def test_qmf_circuit_reference_state_h2(self):
+        """ Verify construction of H2 ansatz works using a circuit reference state."""
+
+        # Build ansatz and circuit
+        qmf_var_params = [np.pi] * 2 + [0.] * 6
+        qmf_circuit = get_qmf_circuit(np.array(qmf_var_params), True)
+        qmf_ansatz = QMF(mol_H2_sto3g, reference_state=qmf_circuit)
+        qmf_ansatz.build_circuit()
+
+        # Build qubit hamiltonian for energy evaluation
+        qubit_hamiltonian = qmf_ansatz.qubit_ham
+
+        # Assert energy returned is as expected for given parameters
+        qmf_ansatz.update_var_params(qmf_var_params)
+        energy = sim.get_expectation_value(qubit_hamiltonian, qmf_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.1166843870853400, delta=1e-6)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_uccgd.py b/tangelo/toolboxes/ansatz_generator/tests/test_uccgd.py
index 723c1ea0..4bb7e70f 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_uccgd.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_uccgd.py
@@ -19,6 +19,8 @@
 from tangelo.molecule_library import mol_H2_sto3g, mol_H4_cation_sto3g
 from tangelo.toolboxes.qubit_mappings import jordan_wigner
 from tangelo.toolboxes.ansatz_generator.uccgd import UCCGD
+from tangelo.toolboxes.operators import count_qubits
+from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit
 from tangelo.linq import get_backend
 
 
@@ -91,6 +93,31 @@ def test_uccgd_H4_open(self):
         energy = sim.get_expectation_value(qubit_hamiltonian, uccgd_ansatz.circuit)
         self.assertAlmostEqual(energy, -1.64190668, delta=1e-6)
 
+    def test_uccgd_circuit_reference_state_H2(self):
+        """ Verify construction of H2 ansatz works using a circuit reference state."""
+
+        # Build qubit hamiltonian
+        qubit_hamiltonian = jordan_wigner(mol_H2_sto3g.fermionic_hamiltonian)
+        n_qubits = count_qubits(qubit_hamiltonian)
+
+        # Construct reference circuit by hand
+        ref_h2_circuit = get_reference_circuit(
+                n_spinorbitals=n_qubits,
+                n_electrons=mol_H2_sto3g.n_electrons,
+                mapping='jw',
+                up_then_down=False,
+                spin=mol_H2_sto3g.spin)
+
+        # Build circuit
+        uccgd_ansatz = UCCGD(mol_H2_sto3g, reference_state=ref_h2_circuit)
+        uccgd_ansatz.build_circuit()
+
+        # Assert energy returned is as expected for given parameters
+        sim = get_backend()
+        uccgd_ansatz.update_var_params([0.78525105, 1.14993361, 1.57070471])
+        energy = sim.get_expectation_value(qubit_hamiltonian, uccgd_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.1372701, delta=1e-6)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_uccsd.py b/tangelo/toolboxes/ansatz_generator/tests/test_uccsd.py
index bf7a6831..44f464e7 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_uccsd.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_uccsd.py
@@ -20,6 +20,8 @@
 from tangelo.molecule_library import mol_H2_sto3g, mol_H4_sto3g, mol_H4_doublecation_minao, mol_H4_cation_sto3g
 from tangelo.toolboxes.qubit_mappings import jordan_wigner
 from tangelo.toolboxes.ansatz_generator.uccsd import UCCSD
+from tangelo.toolboxes.operators import count_qubits
+from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit
 from tangelo.linq import get_backend
 
 # For openfermion.load_operator function.
@@ -146,6 +148,30 @@ def test_uccsd_H4_doublecation(self):
         energy = sim.get_expectation_value(qubit_hamiltonian, uccsd_ansatz.circuit)
         self.assertAlmostEqual(energy, -0.854607, delta=1e-4)
 
+    def test_uccsd_circuit_reference_state_H2(self):
+        """ Verify construction of H2 ansatz works using a circuit reference state."""
+
+        # Build qubit hamiltonian
+        qubit_hamiltonian = jordan_wigner(mol_H2_sto3g.fermionic_hamiltonian)
+
+        # Construct reference circuit by hand
+        ref_h2_circuit = get_reference_circuit(
+                n_spinorbitals=mol_H2_sto3g.n_sos,
+                n_electrons=mol_H2_sto3g.n_electrons,
+                mapping='jw',
+                up_then_down=False,
+                spin=mol_H2_sto3g.spin)
+
+        # Build circuit
+        uccsd_ansatz = UCCSD(mol_H2_sto3g, reference_state=ref_h2_circuit)
+        uccsd_ansatz.build_circuit()
+
+        # Assert energy returned is as expected for given parameters
+        sim = get_backend()
+        uccsd_ansatz.update_var_params([5.86665842e-06, 0.0565317429])
+        energy = sim.get_expectation_value(qubit_hamiltonian, uccsd_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.137270174551959, delta=1e-6)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_upccgsd.py b/tangelo/toolboxes/ansatz_generator/tests/test_upccgsd.py
index 5373aacc..595e3e33 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_upccgsd.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_upccgsd.py
@@ -20,6 +20,7 @@
 from tangelo.molecule_library import mol_H2_sto3g, mol_H4_doublecation_minao, mol_H4_cation_sto3g
 from tangelo.toolboxes.qubit_mappings import jordan_wigner
 from tangelo.toolboxes.ansatz_generator.upccgsd import UpCCGSD
+from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit
 from tangelo.linq import get_backend
 
 # For openfermion.load_operator function.
@@ -118,6 +119,31 @@ def test_upccgsd_H4_doublecation(self):
         energy = sim.get_expectation_value(qubit_hamiltonian, upccgsd_ansatz.circuit)
         self.assertAlmostEqual(energy, -0.854608, delta=1e-4)
 
+    def test_upccgsd_circuit_reference_state_H2(self):
+        """ Verify construction of H2 ansatz works using a circuit reference state."""
+
+        # Build qubit hamiltonian
+        qubit_hamiltonian = jordan_wigner(mol_H2_sto3g.fermionic_hamiltonian)
+
+        # Construct reference circuit by hand
+        ref_h2_circuit = get_reference_circuit(
+                n_spinorbitals=mol_H2_sto3g.n_sos,
+                n_electrons=mol_H2_sto3g.n_electrons,
+                mapping='jw',
+                up_then_down=False,
+                spin=mol_H2_sto3g.spin)
+
+        # Build circuit
+        upccgsd_ansatz = UpCCGSD(mol_H2_sto3g, reference_state=ref_h2_circuit)
+        upccgsd_ansatz.build_circuit()
+
+        # Assert energy returned is as expected for given parameters
+        sim = get_backend()
+        upccgsd_ansatz.update_var_params([0.03518165, -0.02986551,  0.02897598, -0.03632711,
+                                          0.03044071,  0.08252277])
+        energy = sim.get_expectation_value(qubit_hamiltonian, upccgsd_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.1372658, delta=1e-6)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tangelo/toolboxes/ansatz_generator/uccgd.py b/tangelo/toolboxes/ansatz_generator/uccgd.py
index 4e1781d6..1a4fe307 100644
--- a/tangelo/toolboxes/ansatz_generator/uccgd.py
+++ b/tangelo/toolboxes/ansatz_generator/uccgd.py
@@ -45,9 +45,10 @@ class UCCGD(Ansatz):
         up_then_down (bool): change basis ordering putting all spin up orbitals
             first, followed by all spin down. Default, False (i.e. has
             alternating spin up/down ordering).
-        reference_state (string): The reference state id for the ansatz. The
+        reference_state (string, Circuit): The reference state id for the ansatz. The
             supported reference states are stored in the supported_reference_state
-            attributes. Default, "HF".
+            attributes. Can also be a Circuit object, in which case a copy of
+            circuit with variational parameters fixed is used. Default: "HF".
     """
 
     def __init__(self, molecule, mapping="JW", up_then_down=False, reference_state="HF"):
@@ -112,6 +113,11 @@ def prepare_reference_state(self):
         the qubit operator.
         """
 
+        if isinstance(self.reference_state, Circuit):
+            ref_circuit = self.reference_state.copy()
+            ref_circuit.fix_variational_parameters()
+            return ref_circuit
+
         if self.reference_state not in self.supported_reference_state:
             raise ValueError(f"Only supported reference state methods are:{self.supported_reference_state}")
 
diff --git a/tangelo/toolboxes/ansatz_generator/uccsd.py b/tangelo/toolboxes/ansatz_generator/uccsd.py
index 3f87c79e..31b53e3c 100644
--- a/tangelo/toolboxes/ansatz_generator/uccsd.py
+++ b/tangelo/toolboxes/ansatz_generator/uccsd.py
@@ -54,8 +54,10 @@ class UCCSD(Ansatz):
         up_then_down (bool): change basis ordering putting all spin up orbitals
             first, followed by all spin down. Default, False (i.e. has
             alternating spin up/down ordering).
-        reference_state (string): The reference state id for the ansatz. The
-            supported reference states are stored in the supported_reference_state
+        reference_state (string, Circuit): The reference state id for the ansatz.
+            Can also be a Circuit object, in which case a copy of
+            circuit with variational parameters fixed is used. The supported
+            string reference states are stored in the supported_reference_state
             attributes. Default, "HF".
     """
 
@@ -140,6 +142,11 @@ def prepare_reference_state(self):
         the qubit operator.
         """
 
+        if isinstance(self.reference_state, Circuit):
+            ref_circuit = self.reference_state.copy()
+            ref_circuit.fix_variational_parameters()
+            return ref_circuit
+
         if self.reference_state not in self.supported_reference_state:
             raise ValueError(f"Only supported reference state methods are:{self.supported_reference_state}")
 
diff --git a/tangelo/toolboxes/ansatz_generator/upccgsd.py b/tangelo/toolboxes/ansatz_generator/upccgsd.py
index 2b7f7c48..df7391e1 100644
--- a/tangelo/toolboxes/ansatz_generator/upccgsd.py
+++ b/tangelo/toolboxes/ansatz_generator/upccgsd.py
@@ -48,9 +48,10 @@ class UpCCGSD(Ansatz):
         up_then_down (bool): change basis ordering putting all spin up orbitals
             first, followed by all spin down. Default, False (i.e. has
             alternating spin up/down ordering).
-        reference_state (string): The reference state id for the ansatz. The
-            supported reference states are stored in the supported_reference_state
-            attributes. Default, "HF".
+        reference_state (string, Circuit): The reference state id for the ansatz. The
+            supported string reference states are stored in the supported_reference_state
+            attributes.  Can also be a Circuit object, in which case a copy of
+            circuit with variational parameters fixed is used. Default: "HF".
     """
 
     def __init__(self, molecule, mapping="JW", up_then_down=False, k=2, reference_state="HF"):
@@ -117,6 +118,11 @@ def prepare_reference_state(self):
         the qubit operator.
         """
 
+        if isinstance(self.reference_state, Circuit):
+            ref_circuit = self.reference_state.copy()
+            ref_circuit.fix_variational_parameters()
+            return ref_circuit
+
         if self.reference_state not in self.supported_reference_state:
             raise ValueError(f"Only supported reference state methods are:{self.supported_reference_state}")
 
diff --git a/tangelo/toolboxes/optimizers/rotosolve.py b/tangelo/toolboxes/optimizers/rotosolve.py
index 808a57a6..e8cf9b69 100644
--- a/tangelo/toolboxes/optimizers/rotosolve.py
+++ b/tangelo/toolboxes/optimizers/rotosolve.py
@@ -15,10 +15,35 @@
 import numpy as np
 
 
-def rotosolve_step(func, var_params, i, *func_args):
+def extrapolate_expval(theta, m_0, m_minus, m_plus, phi=0.0):
+    """Extrapolates the expectation value of an observable
+    M with respect to a single parameterized rotation (i.e.
+    RX, RY, RZ) with angle theta. The extrapolation uses
+    samples taken at the angles phi, phi+pi/2, and phi-pi/2.
+    This function uses the formula in Appendix A from
+    arXiv:1905.09692, by Mateusz Ostaszewski et al.
+
+    Args:
+        theta (float): Gate rotation angle to extrapolate to
+        m_0 (float): Expectation value of M mat angle phi
+        m_minus (float): Expectation value of M mat angle phi - pi/2
+        m_plus (float): Expectation value of M mat angle phi + pi/2
+        phi (float, optional): Angle of phi. Defaults to 0.0
+
+    Returns:
+        float: The expectation value of M estimated for theta.
+    """
+    a = 0.5*np.sqrt((2 * m_0 - m_plus - m_minus)**2 + (m_plus - m_minus)**2)
+    b = np.arctan2(2 * m_0 - m_plus - m_minus, m_plus - m_minus) - phi
+    c = 0.5*(m_plus + m_minus)
+
+    return a*np.sin(theta + b) + c
+
+
+def rotosolve_step(func, var_params, i, *func_args, phi=0.0, m_phi=None):
     """Gradient free optimization step using specific points to
-    characterize objective function w.r.t to parameter values. Based on
-    formulas in arXiv:1905.09692, Mateusz Ostaszewski
+    characterize objective function w.r.t to parameter values.
+    Based on formulas in arXiv:1905.09692, Mateusz Ostaszewski
 
         Args:
             func (function handle): The function that performs energy
@@ -27,35 +52,44 @@ def rotosolve_step(func, var_params, i, *func_args):
             var_params (list of float): The variational parameters.
             i (int): Index of the variational parameter to update.
             *func_args (tuple): Optional arguments to pass to func.
+            phi (float): Optional angle phi for extrapolation (default is 0.0).
+            m_phi (float): Optional estimated value of m_phi
         Returns:
-            list of floats: Optimal parameters.
+            list of floats: Optimal parameters
+            float: Estimated optimal value of func
     """
 
     # Charaterize sinusoid of objective function using specific parameters
-    var_params[i] = 0
-    m_1 = func(var_params, *func_args)
+    var_params[i] = phi
+    m_0 = func(var_params, *func_args) if m_phi is None else m_phi
 
-    var_params[i] = 0.5 * np.pi
-    m_2 = func(var_params, *func_args)
+    var_params[i] = phi + 0.5 * np.pi
+    m_plus = func(var_params, *func_args)
 
-    var_params[i] = -0.5 * np.pi
-    m_3 = func(var_params, *func_args)
+    var_params[i] = phi - 0.5 * np.pi
+    m_minus = func(var_params, *func_args)
 
     # Calculate theta_min based on measured values
-    theta_min = -0.5 * np.pi - np.arctan2(2. * m_1 - m_2 - m_3, m_2 - m_3)
+    theta_min = phi - 0.5 * np.pi - \
+        np.arctan2(2. * m_0 - m_plus - m_minus, m_plus - m_minus)
 
     if theta_min < -np.pi:
         theta_min += 2 * np.pi
     elif theta_min > np.pi:
         theta_min -= 2 * np.pi
 
+    # calculate extrapolated minimum energy estimate:
+    m_min_estimate = \
+        extrapolate_expval(theta_min, m_0, m_minus, m_plus, phi=phi)
+
     # Update parameter to theta_min
     var_params[i] = theta_min
 
-    return var_params
+    return var_params, m_min_estimate
 
 
-def rotosolve(func, var_params, *func_args, ftol=1e-5, maxiter=100):
+def rotosolve(func, var_params, *func_args, ftol=1e-5, maxiter=100,
+              extrapolate=False):
     """Optimization procedure for parameterized quantum circuits whose
      objective function varies sinusoidally with the parameters. Based
      on the work by arXiv:1905.09692, Mateusz Ostaszewski.
@@ -68,6 +102,12 @@ def rotosolve(func, var_params, *func_args, ftol=1e-5, maxiter=100):
         ftol (float): Convergence threshold.
         maxiter (int): The maximum number of iterations.
         *func_args (tuple): Optional arguments to pass to func.
+        extrapolate (bool): If True, the expectation value of func
+            extrapolated from previous calls to `rotosolve_step()` will
+            be used instead of a function evaluation. This requires
+            only two function evaluations per parameter per iteration,
+            but may be less stable on noisy devices. If False, three
+            evaluations are used per parameter per iteration.
 
     Returns:
         float: The optimal energy found by the optimizer.
@@ -75,10 +115,21 @@ def rotosolve(func, var_params, *func_args, ftol=1e-5, maxiter=100):
      """
     # Get intial value, and run rotosolve for up to maxiter iterations
     energy_old = func(var_params, *func_args)
+
     for it in range(maxiter):
+
         # Update parameters one at a time using rotosolve_step
-        for i in range(len(var_params)):
-            var_params = rotosolve_step(func, var_params, i, *func_args)
+        energy_est = energy_old
+        for i, theta in enumerate(var_params):
+            # Optionally re-use the extrapolated energy as m_phi
+            if extrapolate:
+                var_params, energy_est = \
+                    rotosolve_step(func, var_params, i, *func_args,
+                                   phi=theta, m_phi=energy_est)
+            else:
+                var_params, energy_est = \
+                    rotosolve_step(func, var_params, i, *func_args)
+
         energy_new = func(var_params, *func_args)
 
         # Check if convergence tolerance is met
@@ -89,3 +140,99 @@ def rotosolve(func, var_params, *func_args, ftol=1e-5, maxiter=100):
         energy_old = energy_new
 
     return energy_new, var_params
+
+
+def rotoselect_step(func, var_params, var_rot_axes, i, *func_args):
+    """Gradient free optimization step using specific points to
+    characterize objective function w.r.t to parameterized
+    rotation axes and rotation angles. Based on formulas in
+    arXiv:1905.09692, Mateusz Ostaszewski
+
+        Args:
+            func (function handle): The function that performs energy
+                estimation. This function takes variational parameters and
+                parameter rotation axes (list of "RX", "RY", "RZ" strings)
+                as input and returns a float.
+            var_params (list of float): The variational parameters.
+            var_rot_axes (list): List of strings ("RX", "RY", or "RZ")
+                corresonding to the axis of rotation for each angle in
+                the list of variational parameters.
+            i (int): Index of the variational parameter to update.
+            *func_args (tuple): Optional arguments to pass to func.
+        Returns:
+            list of floats: Optimal parameters
+            list of strs: Optimal rotation axes
+    """
+    axes = ['RX', 'RY', 'RZ']
+    m_axes = np.zeros(3)
+    theta_min_axes = np.zeros(3)
+
+    # Evaluate func at phi = 0 (same result for all axes)
+    var_params[i] = 0
+    m_0 = func(var_params, var_rot_axes, *func_args)
+
+    # Do a rotosolve step for each axis:
+    rotosolve_func_args = (var_rot_axes,) + func_args
+    for k, axis in enumerate(axes):
+        var_rot_axes[i] = axis
+        var_params, m_axes[k] = \
+            rotosolve_step(func, var_params, i,
+                           *rotosolve_func_args)
+        theta_min_axes[k] = var_params[i]
+
+    # Select optimal axis yielding minimal value
+    k_opt = np.argmin(m_axes)
+    var_rot_axes[i] = axes[k_opt]
+    var_params[i] = theta_min_axes[k_opt]
+
+    return var_params, var_rot_axes
+
+
+def rotoselect(func, var_params, var_rot_axes, *func_args, ftol=1e-5,
+               maxiter=100):
+    """Optimization procedure for parameterized quantum circuits whose
+       objective function varies sinusoidally with the parameters. This
+       routine differs from `rotosolve` by sampling expectation values
+       using the Pauli {X,Y,Z} generators instead of shifted angles of
+       rotation. Based on the work by arXiv:1905.09692, Mateusz
+       Ostaszewski.
+
+    Args:
+        func (function handle): The function that performs energy
+            estimation. This function takes variational parameters and
+            parameter rotation axes (list of "RX", "RY", "RZ" strings)
+            as input and returns a float.
+        var_params (list): The variational parameters.
+        var_rot_axes (list): List of strings ("RX", "RY", or "RZ")
+            corresonding to the axis of rotation for each angle in
+            the list of variational parameters.
+        ftol (float): Convergence threshold.
+        maxiter (int): The maximum number of iterations.
+        *func_args (tuple): Optional arguments to pass to func.
+
+    Returns:
+        float: The optimal energy found by the optimizer.
+        list of floats: Optimal parameters.
+        list of strings: Optimal rotation axes.
+    """
+    # Check parameters and rotation axes are the same length:
+    assert len(var_params) == len(var_rot_axes)
+
+    # Get intial value, and run rotosolve for up to maxiter iterations
+    energy_old = func(var_params, var_rot_axes, *func_args)
+    for it in range(maxiter):
+
+        # Update parameters one at a time using rotosolve_step
+        for i in range(len(var_params)):
+            var_params, var_rot_axes = \
+                rotoselect_step(func, var_params, var_rot_axes, i, *func_args)
+        energy_new = func(var_params, var_rot_axes, *func_args)
+
+        # Check if convergence tolerance is met
+        if abs(energy_new - energy_old) <= ftol:
+            break
+
+        # Update energy value
+        energy_old = energy_new
+
+    return energy_new, var_params, var_rot_axes
diff --git a/tangelo/toolboxes/optimizers/tests/test_rotosolve.py b/tangelo/toolboxes/optimizers/tests/test_rotosolve.py
index e251480d..bb719746 100644
--- a/tangelo/toolboxes/optimizers/tests/test_rotosolve.py
+++ b/tangelo/toolboxes/optimizers/tests/test_rotosolve.py
@@ -17,7 +17,8 @@
 import numpy as np
 
 from tangelo.linq import get_backend, Gate, Circuit
-from tangelo.toolboxes.optimizers.rotosolve import rotosolve
+from tangelo.toolboxes.optimizers.rotosolve import rotosolve, rotoselect
+from tangelo.toolboxes.operators.operators import QubitOperator
 from tangelo.molecule_library import mol_H2_sto3g
 from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
 from tangelo.toolboxes.ansatz_generator import VariationalCircuitAnsatz
@@ -58,10 +59,134 @@ def exp(var_params, ansatz, qubit_hamiltonian):
             return energy
 
         # Run rotosolve, returning energy
-        energy, _ = rotosolve(exp, ansatz.var_params_default, ansatz, qubit_hamiltonian)
+        energy, _ = rotosolve(exp, ansatz.var_params_default, ansatz, qubit_hamiltonian, extrapolate=False)
 
         self.assertAlmostEqual(energy, -1.137270422018, delta=1e-4)
 
+    def test_rotosolve_extrapolate(self):
+        """Test rotosovle on H2 without VQE, using custom variational circuit
+        and qubit Hamiltonian with JW qubit mapping on an exact simulator.
+        """
+        sim = get_backend()
+        # Create qubit Hamiltonian compatible with UCC1 Ansatz
+        qubit_hamiltonian = fermion_to_qubit_mapping(fermion_operator=mol_H2_sto3g.fermionic_hamiltonian,
+                                                     mapping="jw",
+                                                     n_spinorbitals=mol_H2_sto3g.n_active_sos,
+                                                     up_then_down=True,)
+
+        # Manual input of UCC1 circuit with extra variational parameters
+        circuit = Circuit()
+        # Create excitation ladder circuit used to build entangler
+        excit_gates = [Gate("RX", 0, parameter=np.pi/2, is_variational=True)]
+        excit_gates += [Gate("H", i) for i in {1, 2, 3}]
+        excit_gates += [Gate("CNOT", i+1, i) for i in range(3)]
+        excit_circuit = Circuit(excit_gates)
+        # Build UCC1 circuit: mean field + entangler circuits
+        circuit = Circuit([Gate("X", i) for i in {0, 2}])
+        circuit += excit_circuit
+        circuit.add_gate(Gate("RZ", 3, parameter=0, is_variational=True))
+        circuit += excit_circuit.inverse()
+        # Translate circuit into variational ansatz
+        ansatz = VariationalCircuitAnsatz(circuit)
+
+        # Define function to calculate energy and update variational parameters
+        def exp(var_params, ansatz, qubit_hamiltonian):
+            ansatz.update_var_params(var_params)
+            energy = sim.get_expectation_value(qubit_hamiltonian, ansatz.circuit)
+            return energy
+
+        # Run rotosolve, returning energy
+        energy, _ = rotosolve(exp, ansatz.var_params_default, ansatz, qubit_hamiltonian, extrapolate=True)
+
+        self.assertAlmostEqual(energy, -1.137270422018, delta=1e-4)
+
+    def test_rotoselect(self):
+        """Test rotoselect using a single-qubit Euler rotation circuit"""
+
+        sim = get_backend()
+
+        # Build an Euler rotation circuit as an ansatz
+        euler_circuit = Circuit([
+            Gate('RZ', 0, parameter=0, is_variational=True),
+            Gate('RX', 0, parameter=0, is_variational=True),
+            Gate('RZ', 0, parameter=0, is_variational=True)
+        ])
+        ansatz = VariationalCircuitAnsatz(euler_circuit)
+
+        # Build a single-qubit Hamiltonian
+        hamiltonian = \
+            QubitOperator((0,'X'), 1.0) + \
+            QubitOperator((0,'Y'), 2.0) + \
+            QubitOperator((0,'Z'), 3.0)
+
+        # Define function to calculate energy and update parameters and rotation axes
+        def exp_rotoselect(var_params, var_rot_axes, ansatz, qubit_hamiltonian):
+            ansatz.update_var_params(var_params)
+            for i, axis in enumerate(var_rot_axes):
+                ansatz.circuit._variational_gates[i].name = axis
+            energy = sim.get_expectation_value(qubit_hamiltonian, ansatz.circuit)
+            return energy
+
+        # Run rotoselect, return energy, parameters and axes of rotation:
+        init_params = ansatz.var_params_default
+        init_axes = ['RX']*len(init_params)
+        energy, _, axes = rotoselect(exp_rotoselect,
+                                init_params, init_axes, ansatz, hamiltonian)
+
+        # compare with exact energy:
+        min_energy = -np.sqrt(1**2 + 2**2 + 3**2)
+        self.assertAlmostEqual(energy, min_energy, delta=1e-4)
+
+        # Ensure axes are all valid rotation gates:
+        self.assertTrue(set(axes).issubset({'RX', 'RY', 'RZ'}))
+
+    def test_rotoselect_heisenberg(self):
+        """Test rotoselect using the 5-qubit periodic Heisenberg model"""
+
+        sim = get_backend()
+        n_qubits = 3
+        n_layers = 2
+        J = h = 1.0
+
+        # Construct a "hardware efficient" CZ-based ansatz layer
+        heisenberg_gates = [Gate('Ry', i,parameter=0, is_variational=True) for i in range(n_qubits)]
+        heisenberg_gates += [Gate('CZ', i, (i+1) % n_qubits) for i in range(0,n_qubits-1,2)]
+        heisenberg_gates += [Gate('CZ', i, (i+1) % n_qubits) for i in range(1,n_qubits,2)]
+        heisenberg_layer = Circuit(heisenberg_gates)
+
+        heisenberg_circuit = Circuit()
+        for _ in range(n_layers):
+            heisenberg_circuit += heisenberg_layer
+        ansatz = VariationalCircuitAnsatz(heisenberg_circuit)
+
+        # Construct periodic Heisenberg Hamiltonian
+        hamiltonian = QubitOperator()
+        for i in range(n_qubits):
+            hamiltonian += QubitOperator((i,'Z'), h)
+            for S in ['X','Y','Z']:
+                hamiltonian += QubitOperator([(i,S),((i+1) % n_qubits,S)],J)
+
+        # Define function to calculate energy and update parameters and rotation axes
+        def exp_rotoselect(var_params, var_rot_axes, ansatz, qubit_hamiltonian):
+            ansatz.update_var_params(var_params)
+            for i, axis in enumerate(var_rot_axes):
+                ansatz.circuit._variational_gates[i].name = axis
+            energy = sim.get_expectation_value(qubit_hamiltonian, ansatz.circuit)
+            return energy
+
+        # Run rotoselect, return energy, parameters and axes of rotation:
+        init_params = [np.pi/3]*ansatz.n_var_params
+        init_axes = ['RX']*len(init_params)
+        energy, _, axes = rotoselect(exp_rotoselect,
+                                init_params, init_axes, ansatz, hamiltonian)
+
+        # compare with known ground state energy:
+        min_energy = -4.0
+        self.assertAlmostEqual(energy, min_energy, delta=1e-4)
+
+        # Ensure axes are all valid rotation gates:
+        self.assertTrue(set(axes).issubset({'RX', 'RY', 'RZ'}))
+
 
 if __name__ == "__main__":
     unittest.main()