Skip to content

Molecule Eigensolver (VQE Method)

View on GitHub

Evaluating the ground state of a molecular Hamiltonian allows you to understand the chemical properties of the molecule. This tutorial demonstrates the use of Variational Quantum Eigensolver (VQE) to find the ground states and energies of 𝐻2, 𝐻2𝑂, and 𝐿𝑖𝐻 molecules.

VQE is a leading method for finding approximate values of ground state wave functions and energies for complicated quantum systems and can give solutions for complex molecular structures. The overview of the VQE method is as follows: a problem (i.e., a molecule) is defined by a Hamiltonian whose ground state is sought. Then, a choice of a parameterized ansatz is made. A hybrid quantum-classical algorithm finds a solution for the defined parameters that minimizes the expectation value for the energy. A clever ansatz leads to an estimated ground state solution.

Within the scope of Classiq's VQE algorithm, define a molecule that is translated to a concise Hamiltonian. Then, choose among types of well studied ansatzes, which are carefully selected to fit your molecule type. In the last stage, the Hamiltonian and ansatz are sent to a classical optimizer. This tutorial demonstrates the steps and options in Classiq's VQE algorithm. It presents the optimization strength of Classiq's VQE algorithm and its state-of-the-art results in terms of efficient quantum circuit, with the ultimate combination of low depth and high accuracy while minimizing the number of CX gates.

!pip install -qq "classiq[chemistry]"

Generating a Qubit Hamiltonian

Define the molecule to simulate, declaring the MolecularData class and inserting a list of atoms and their spatial positions (the distances are received in Γ… =\(10^{-10} m\)). In addition, provide basis, multiplicity, and charge.

As mentioned above, this tutorial demonstrates how to define and find the ground state and energies for these molecules:

molecule_H2_geometry = [("H", (0.0, 0.0, 0)), ("H", (0.0, 0.0, 0.735))]
molecule_O2_geometry = [("O", (0.0, 0.0, 0)), ("O", (0.0, 0.0, 1.16))]
molecule_LiH_geometry = [("H", (0.0, 0.0, 0.0)), ("Li", (0.0, 0.0, 1.596))]
molecule_H2O_geometry = [
    ("O", (0.0, 0.0, 0.0)),
    ("H", (0, 0.586, 0.757)),
    ("H", (0, 0.586, -0.757)),
]
molecule_BeH2_geometry = [
    ("Be", (0.0, 0.0, 0.0)),
    ("H", (0, 0, 1.334)),
    ("H", (0, 0, -1.334)),
]

You can construct any valid assembly of atoms in a similar manner.

from openfermion.chem import MolecularData
from openfermionpyscf import run_pyscf

geometry = molecule_H2_geometry

basis = "sto-3g"  # Basis set
multiplicity = 1  # Singlet state S=0
charge = 0  # Neutral molecule
molecule = MolecularData(molecule_H2_geometry, basis, multiplicity, charge)

molecule = run_pyscf(
    molecule,
    run_mp2=True,
    run_cisd=True,
    run_ccsd=True,
    run_fci=True,  # relevant for small, classically solvable problems
)

Define the parameters of the Hamiltonian problem (FermionHamiltonianProblem) and the mapper (FermionToQubitMapper) between Fermionic Hamiltonian and qubit Hamiltonians (Jordan Wigner or Bravyi Kitaev). If you want to use Z2-symmteries for reducing the problem size you can use Z2SymTaperMapper (see below).

\[\langle \psi_{hf}| H|\psi_{hf}\rangle\]
from classiq.applications.chemistry.mapping import FermionToQubitMapper
from classiq.applications.chemistry.problems import FermionHamiltonianProblem

# Define a  Hamiltonian in an active space
problem = FermionHamiltonianProblem.from_molecule(molecule=molecule)
mapper = FermionToQubitMapper()


qubit_hamiltonian = mapper.map(problem.fermion_hamiltonian)
print("Your Hamiltonian is", qubit_hamiltonian, sep="\n")
num_qubits = mapper.get_num_qubits(problem)
print(f"number of qubits {num_qubits}")

Your Hamiltonian is
(-0.09057898608834769+0j) [] +
(0.04523279994605784+0j) [X0 X1 X2 X3] +
(0.04523279994605784+0j) [X0 X1 Y2 Y3] +
(0.04523279994605784+0j) [Y0 Y1 X2 X3] +
(0.04523279994605784+0j) [Y0 Y1 Y2 Y3] +
(0.17218393261915538+0j) [Z0] +
(0.12091263261776627+0j) [Z0 Z1] +
(0.16892753870087907+0j) [Z0 Z2] +
(0.1661454325638241+0j) [Z0 Z3] +
(-0.2257534922240238+0j) [Z1] +
(0.1661454325638241+0j) [Z1 Z2] +
(0.17464343068300453+0j) [Z1 Z3] +
(0.1721839326191554+0j) [Z2] +
(0.12091263261776627+0j) [Z2 Z3] +
(-0.22575349222402386+0j) [Z3]
number of qubits 4

Constructing and Synthesizing a Ground State Solver

A ground state solver model consists of a parameterized eigenfunction ("the ansatz"), on which to run a VQE.

Start with a Hardware (HW) efficient ansatz:

HW Efficient Ansatz

The suggested HW efficient ansatz solution is generated to fit a specific hardware [1]. The ansatz creates a state with a given number of parameters according to your choice of the number of qubits that fits the Hamiltonian, and creates entanglement between the qubits using the inputed connectivity map. This example uses a four qubit map, which is specifically made for \(H_2\) without using qubit tapering.

After constructing the model, synthesize it and view the output circuit.

For groundstate solvers, it is typical to initialize the ansatz with the Hartree-Fock state. Use the get_hf_state and the prepare_basis_state qfunc.

from classiq import *
from classiq.applications.chemistry.hartree_fock import get_hf_state
from classiq.applications.chemistry.op_utils import qubit_op_to_pauli_terms

reps = 3
num_params = reps * num_qubits
hf_state = get_hf_state(problem, mapper)
vqe_hamiltonian = qubit_op_to_pauli_terms(mapper.map(problem.fermion_hamiltonian))


@qfunc
def main(params: CArray[CReal, num_params], state: Output[QArray]):
    prepare_basis_state(hf_state, state)
    full_hea(
        num_qubits=num_qubits,
        operands_1qubit=[lambda _, q: X(q), lambda theta, q: RY(theta, q)],
        operands_2qubit=[lambda _, q1, q2: CX(q1, q2)],
        is_parametrized=[0, 1, 0],
        angle_params=params,
        connectivity_map=[(0, 1), (1, 2), (2, 3)],
        reps=reps,
        x=state,
    )


qmod_hwea = create_model(
    main, execution_preferences=ExecutionPreferences(num_shots=1000)
)
write_qmod(qmod_hwea, "molecule_eigensolver_hwea", symbolic_only=False)
qprog_hwea = synthesize(qmod_hwea)
show(qprog_hwea)

Quantum program link: https://platform.classiq.io/circuit/30Y1nbgeb5vBuSlqGFMxdkZqCYh

Unitary Coupled Cluster (UCC) Ansatz

Create the commonly used chemistry-inspired UCC ansatz, which is a unitary version of the classical coupled cluster (CC) method [2].

The parameter that defines the UCC ansatz: excitations (List[int] or List[str]): list of desired excitations, e.g.,

  • 1 for singles

  • 2 for doubles

  • 3 for triples

  • 4 for quadruples

Once again, after running the code lines below, you can view the output circuit that creates the state with an interactive interface and print the depth of the circuit.

For the current example, use the Z2SymTaperMapper that exploits Z2-symmetries of the molecule Hamiltonian to reduce the problem size. You can confirm that using Z2SymTaperMapper.from_problem compared to FermionToQubitMapper, the number of qubits is reduced as:

  • for \(H_2\) - from 4 to 1

  • for \(LiH\) from 12 to 8 (together with freezing the core orbital first_active_index=1)

  • for \(H_{2}O\) from 14 to 10 (together with freezing the core orbital first_active_index=1)

from classiq.applications.chemistry.ucc import get_ucc_hamiltonians
from classiq.applications.chemistry.z2_symmetries import Z2SymTaperMapper

problem = FermionHamiltonianProblem.from_molecule(molecule=molecule)
mapper = Z2SymTaperMapper.from_problem(problem)


qubit_hamiltonian = mapper.map(problem.fermion_hamiltonian)
print("Your Hamiltonian is", qubit_hamiltonian, sep="\n")
num_qubits = mapper.get_num_qubits(problem)
print(f"number of qubits {num_qubits}")


hf_state = get_hf_state(problem, mapper)
uccsd_hamiltonians = get_ucc_hamiltonians(problem, mapper, excitations=[1, 2])
num_params = len(uccsd_hamiltonians)
vqe_hamiltonian = qubit_op_to_pauli_terms(mapper.map(problem.fermion_hamiltonian))


@qfunc
def main(params: CArray[CReal, num_params], state: Output[QArray]):
    prepare_basis_state(hf_state, state)
    multi_suzuki_trotter(uccsd_hamiltonians, params, 1, 1, state)


qmod_ucc = create_model(main, execution_preferences=ExecutionPreferences(num_shots=1e6))
write_qmod(qmod_ucc, "molecule_eigensolver_ucc", symbolic_only=False)
qprog_ucc = synthesize(qmod_ucc)

show(qprog_ucc)

print(f"circuit depth: {qprog_ucc.transpiled_circuit.depth}")

Your Hamiltonian is
-0.32112414706764497 [] +
0.18093119978423144 [X0] +
0.7958748496863588 [Z0]
number of qubits 1
Quantum program link: https://platform.classiq.io/circuit/30Y1nuIOUd8XbNKwQIIfhumJqFp
circuit depth: 3

The Classiq UCC algorithm provides a highly efficient solution in terms of circuit depth and number of CX gates. These ultimately reduce the gate's time and amount of resources needed for operation.

Executing to Find the Ground State

After synthesizing the model you can execute it:

After you specified a Hamiltonian and an ansatz, send the resulting quantum program to the VQE algorithm to find the Hamiltonian's ground state. In the process, the algorithm sends requests to a classical server, whose 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. You can control the max_iteration value so the solution reaches a stable convergence. In addition, the num_shots value sets the number of measurements performed after each iteration, thus influencing the accuracy of the solutions.

with ExecutionSession(qprog_ucc) as es:
    result_ucc = es.minimize(
        cost_function=vqe_hamiltonian,
        initial_params={"params": [0.0] * num_params},
        max_iteration=200,
    )

optimizer_res = result_ucc[-1][0]
optimal_params = result_ucc[-1][1]
print(f"optimizer result: {optimizer_res}")
print(f"optimal parameter: {optimal_params}")

optimizer result: -1.1379656848485658
optimal parameter: {'params': [-0.21640000000000004]}

Note that energy is presented in units of Hartree.

Finally, compare the VQE solution to the classical solution:

expected_energy = molecule.fci_energy
print("exact result:", expected_energy)
print("vqe result:", optimizer_res)

exact result: -1.1373060357533995
vqe result: -1.1379656848485658

[1][Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M. Chow, Jay M. Gambetta Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets. Nature 549, 242 (2017).](https://arxiv.org/abs/1704.05018)

[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: Particle-hole Hamiltonian and optimized wave-function expansions. Phys. Rev. A 98, 022322 (2018).](https://arxiv.org/abs/1805.04340)