Molecule Eigensolver (VQE method)
Evaluating the ground state of a molecular Hamiltonian allows us to understand the chemical properties of the molecule. In this demo, we demonstrate the usage of Variational Quantum Eigensolver (VQE) for finding the ground states and energies of several molecules: π»2 , π»2π and πΏππ» .
VQE is a leading method for finding approximate values of ground state wavefuncions and energies for complicated quantum systems, and as such can give solutions for complex molecular structures. The overview of the VQE method is as following: a problem (i.e. a molecule) is defined by a Hamiltonian which ground state is sought. Then, a choice of a parameterized ansatz is made. Using a hybrid quantumclassical algorithm, a solution for the defined parameters which minimize the expectation value for the energy is found. A clever ansatz will lead to an estimated ground state solution.
Within the scope of Classiq's VQE algorithm, the user defines a Molecule which is being translated to a concise Hamiltonian. Then, a choice between several types of well studied ansatz is given, which can be carefully picked to fit your Molecule type. In the last stage, the Hamiltonian and ansatz are sent to a classical optimizer. During this tutorial we will demonstrate the steps and user options in Classiq's VQE algorithm. Furthermore, the demonstration will present the optimization strength Classiq's VQE algorithm holds, and it's state of the art results in term of efficient quantum circuit  with ultimate combination of low depth and high accuracy, while minimizing the number of CX gates.
0. Prerequirments
The model is using several Classiq's libraries.
import numpy as np
from classiq import QuantumProgram, construct_chemistry_model, execute, show, synthesize
from classiq.applications.chemistry import (
ChemistryExecutionParameters,
HEAParameters,
Molecule,
MoleculeProblem,
UCCParameters,
)
from classiq.execution import (
ClassiqBackendPreferences,
ClassiqSimulatorBackendNames,
ExecutionPreferences,
OptimizerType,
)
from classiq.synthesis import set_execution_preferences
1. Generate Qubit Hamiltonian
The first step is to define the molecule we wish to simulate. We hereby declare the class Molecule and insert a list of atoms and their spacial positions. The algorithm will automatically regard relevant attributes such as the atom's mass, charge and spin.
As mentioned above, during this tutorial, we demonstrate how to define and find the ground state and energies for 3 molecules:
molecule_H2 = Molecule(atoms=[("H", (0.0, 0.0, 0)), ("H", (0.0, 0.0, 0.735))])
molecule_O2 = Molecule(atoms=[("O", (0.0, 0.0, 0)), ("O", (0.0, 0.0, 1.16))])
molecule_LiH = Molecule(atoms=[("H", (0.0, 0.0, 0.0)), ("Li", (0.0, 0.0, 1.596))])
molecule_H2O = Molecule(
atoms=[("O", (0.0, 0.0, 0.0)), ("H", (0, 0.586, 0.757)), ("H", (0, 0.586, 0.757))]
)
molecule_BeH2 = Molecule(
atoms=[("Be", (0.0, 0.0, 0.0)), ("H", (0, 0, 1.334)), ("H", (0, 0, 1.334))]
)
Similarly, the user is able to construct any valid essambly of atoms. The distances are recived in Γ
(\(10^{10} m\)). We will continue this demonstration with a specific molecule. The user can change the molecule
below to study other cases.
molecule = molecule_H2
Next, we define the parameters of the Hamiltonian generation program. The user has a choice over the following options:

mapping (str): the mapping between the fermionic Hamiltonian and an qubits Hamiltonian. Supported types:
 "jordan_wigner"  "parity"  "bravyi_kitaev"  "fast_bravyi_kitaev"

freeze_core (bool): remove the "core" orbitals of the atoms defining the molecule.

z2_symmetries (bool): whether to perform z2 symmetries reduction. If symmetries in the molecules exist, this option will decrease the number of qubits used and will efficient the Hamiltonian and thus the calculations.
Finally, the Hamiltonian is generated from MoleculeProblem
.
chemistry_problem = MoleculeProblem(
molecule=molecule,
mapping="jordan_wigner", #'bravyi_kitaev'
z2_symmetries=True,
freeze_core=True,
)
operator = chemistry_problem.generate_hamiltonian()
gs_problem = chemistry_problem.update_problem(operator.num_qubits)
print("Your Hamiltonian is", operator.show(), sep="\n")
Your Hamiltonian is
1.041 * I
0.796 * Z
+0.181 * X
The output of the above code lines is the Hamiltonian presented as a superposition of Pauli matrices multiplication. One can easily confirm that using z2*symmetries=True, the number of qubits are reduced (compered to z2_symmetries=False): for \(H_2\)  from 4 to 1, for \(LiH\) from 12 to 8, and for \(H*{2}O\) from 14 to 10.
2. Constructing and Synthesizing a Ground State Solver
A ground state solver model consists of a parameterized eigenfunction ("the ansatz"), on which we run a VQE. In addition, a postprocess of the result allows to return the total energy (combining the ground state energy of the Hamiltonian, the nuclear repulsion and the static nuclear energy).
Once we've specified an Hamiltonian and a desired Ansatz, we send them to the VQE algorithm in order to find the Hamiltonian's ground state. In the process, the algorithm will send requests to a classical server, which task is to minimize the energy expectation value and return the optimized parameters. The simulator and optimizing parameters are defined as part of the VQE part of the model. The user should control the max_iteration
value in a manner so the solution has reached a stable convergence. In addition, the value num_shots
sets the number of measurements performed after each iteration, thus influence the accuracy of the solutions.
We demonstrate two different proposal for the wavefunction solution ansatz: (1) Hardware (HW) efficient, and (2) Unitary Coupled Cluster (UCC). For groundstate solvers it is typical to initialize the Ansatz with the HartreeFock state.
2.1 HWEfficient Ansatz
Hardwareefficient ansatz is a suggested solution that is generated to fit a specific hardware [1]. The ansatz creates a state with given number of parameters by user choice (number of qubits, that should fit the Hamiltonian), and creates entanglement between the qubits by the inputed connectivity map. In this example, a 4 qubit map is given, which is specifically made of \(H_2\) with z2_symmetries=False.
After constructing the model, we can synthesize it and view the outputted circuit, in charged on creating the state with an interactive interface.
chemistry_problem = MoleculeProblem(
molecule=molecule,
mapping="jordan_wigner", #'bravyi_kitaev'
z2_symmetries=False,
freeze_core=True,
)
hwea_params = HEAParameters(
num_qubits=4,
connectivity_map=[(0, 1), (1, 2), (2, 3)],
reps=3,
one_qubit_gates=["x", "ry"],
two_qubit_gates=["cx"],
)
qmod = construct_chemistry_model(
chemistry_problem=chemistry_problem,
use_hartree_fock=True,
ansatz_parameters=hwea_params,
execution_parameters=ChemistryExecutionParameters(
optimizer=OptimizerType.COBYLA,
max_iteration=30,
initial_point=None,
),
)
backend_preferences = ClassiqBackendPreferences(
backend_name=ClassiqSimulatorBackendNames.SIMULATOR
)
qmod = set_execution_preferences(
qmod,
execution_preferences=ExecutionPreferences(
num_shots=1000, backend_preferences=backend_preferences
),
)
from classiq import write_qmod
write_qmod(qmod, name="molecule_eigensolver")
qprog = synthesize(qmod)
show(qprog)
Opening: https://platform.classiq.io/circuit/dd765d79bcdb49e085be74af8c1b3f86?version=0.41.0.dev39%2B79c8fd0855
2.2. UCC Ansatz
Next, we show how to create the commonly used chemistryinspired UCC ansatz, which is a unitary version of the classical coupled cluster (CC) method [2] .
The parameter that defines the UCC ansatz is:
 excitations (List[int] or List[str]): list of desired excitations. Allowed excitations:
 1 for singles  2 for doubles  3 for triples  4 for quadruples
Once again, after the code lines bellow run, the user is able to view the outputted circuit, in charged on creating the state with an interactive interface. In addition, the depth of the circuit is printed.
chemistry_problem = MoleculeProblem(
molecule=molecule,
mapping="jordan_wigner", #'bravyi_kitaev'
z2_symmetries=True,
freeze_core=True,
)
serialized_chemistry_model = construct_chemistry_model(
chemistry_problem=chemistry_problem,
use_hartree_fock=True,
ansatz_parameters=UCCParameters(excitations=[1, 2]),
execution_parameters=ChemistryExecutionParameters(
optimizer=OptimizerType.COBYLA,
max_iteration=30,
initial_point=None,
),
)
backend_preferences = ClassiqBackendPreferences(
backend_name=ClassiqSimulatorBackendNames.SIMULATOR
)
serialized_chemistry_model = set_execution_preferences(
serialized_chemistry_model,
execution_preferences=ExecutionPreferences(
num_shots=1000, backend_preferences=backend_preferences
),
)
qprog = synthesize(serialized_chemistry_model)
show(qprog)
circuit = QuantumProgram.from_qprog(qprog)
print(f"circuit depth: {circuit.transpiled_circuit.depth}")
Opening: https://platform.classiq.io/circuit/39cf374d1965408185bf75c71cd73296?version=0.41.0.dev39%2B79c8fd0855
circuit depth: 3
Classiq's UCC algorithm provides an highly efficient solution in aspects of circuit depth and number of CX gates. Those ultimately reduce the gate's time and amount of resources needed for its operation.
3. Execute to Find Ground State
Once we've synthesized the model we can execute it.
result = execute(qprog).result()
chemistry_result_dict = result[1].value
Execution of the quantum program returns several useful outputs:

energy : the output of the VQE algorithm  the electronic energy simulated.

nuclear_repulsion : the electrostatic energy generated by the atom's nuclei.

hartree_fock_energy : the Hartree Fock energy.

total_energy : this is the ground state energy of the Hamiltonian (combining the energy, the nuclear repulsion and the static nuclear energy).
It also contains the full VQE result from which we can get, for example:

optimal_parameters : gives the results for the anzats parameters minimizing that expectation value.

eigenstate : gives the ground state wave function.
Note the all energy are presented in units of Hartree.
chemistry_result_dict["total_energy"]
1.1382090268147285
chemistry_result_dict["vqe_result"]["optimal_parameters"]
{'param_0': 3.2624618497328566}
Finally, we can compare the VQE solution to the classical solution by employing exact diagonalization:
mat = operator.to_matrix()
w, v = np.linalg.eig(mat)
print("exact result:", np.real(min(w)))
print("vqe result:", chemistry_result_dict["energy"])
exact result: 1.8572750302023786
vqe result: 1.8581780212637082
[1] Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M. Chow, Jay M. Gambetta Hardwareefficient variational quantum eigensolver for small molecules and quantum magnets. Nature 549, 242 (2017)
[2] Panagiotis Kl. Barkoutsos, Jerome F. Gonthier, Igor Sokolov, Nikolaj Moll, Gian Salis, Andreas Fuhrer, Marc Ganzhorn, Daniel J. Egger, Matthias Troyer, Antonio Mezzacapo, Stefan Filipp, and Ivano Tavernelli Quantum algorithms for electronic structure calculations: Particlehole Hamiltonian and optimized wavefunction expansions Phys. Rev. A 98, 022322 (2018)