Skip to content

Creating a Molecule's Potential Energy Curve

View on GitHub Experiment in the IDE

A potential energy curve gives the ground energy of an assembly of atoms as a function of the distances between them. The global minima of the curve indicates the binding energy and internuclear distance for the stable molecule. Therefore, such a curve can be powerful tool in computational chemistry for predicting the molecular structure and spectrum.

This tutorial demonstrates how to use the Classiq chemistry package to create a molecule's potential energy curve. It compares the result with the Hartree-Fock approximation method and the FCI (Full Configuration Interaction) energy.

Definitions and Initialization

Define the range of internuclear distances for the model to simulate, and choose the number of sampling points in this range, which determines the graph's resolution.

import numpy as np

# define the sampling params
num1 = 5  # how many sampling points - determines your resolution
start1 = 0.20  # what is your sampling start distance
stop1 = 1  # what is your sampling end distance
num2 = 7  # how many sampling points - determines your resolution
start2 = 1.4  # what is your sampling start distance
stop2 = 3.5  # what is your sampling end distance


# prepare x,y vectors
distance = np.append(np.linspace(start1, stop1, num1), np.linspace(start2, stop2, num2))
VQE_energy = []
HF_energy = []
exact_energy = []

print(distance)

[0.2  0.4  0.6  0.8  1.   1.4  1.75 2.1  2.45 2.8  3.15 3.5 ]

Energy Estimations

Define a for-loop, which takes these steps:

  1. Creating a molecule at changing distances between the atoms.

  2. Constructing a chemistry model for the corresponding Hamiltonian, using the Hartree-Fock initial state and UCC ansatz.

  3. Synthesizing the model to get a quantum program.

  4. Executing the quantum program to extract a solution for the ground energy.

  5. Obtaining the exact solution and Hartree-Fock solution.

import time

from classiq import *

from openfermion.chem import MolecularData
from openfermionpyscf import run_pyscf

from classiq.applications.chemistry.hartree_fock import get_hf_state
from classiq.applications.chemistry.op_utils import qubit_op_to_pauli_terms
from classiq.applications.chemistry.problems import FermionHamiltonianProblem
from classiq.applications.chemistry.ucc import get_ucc_hamiltonians
from classiq.applications.chemistry.z2_symmetries import Z2SymTaperMapper

# create the molecule, insert the distance, prepare H, create UCC anzats and solve in energy
qmods = []
qprogs = []
results = []
durations = []

for x in distance:
    time1 = time.time()

    # Define a molecule
    geometry = [("H", (0.0, 0.0, 0)), ("H", (0.0, 0.0, float(x)))]
    molecule = MolecularData(geometry, basis="sto-3g", multiplicity=1, charge=0)
    molecule = run_pyscf(
        molecule,
        run_mp2=True,
        run_cisd=True,
        run_ccsd=True,
        run_fci=True,  # relevant for small, classically solvable problems
    )

    # Define a problem and a mapper
    problem = FermionHamiltonianProblem.from_molecule(molecule)
    mapper = Z2SymTaperMapper.from_problem(problem)

    # Construct a model
    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)

    # Synthesize
    qprog = synthesize(main)
    qprog = set_quantum_program_execution_preferences(
        qprog,
        preferences=ExecutionPreferences(
            num_shots=1000,
            backend_preferences=ClassiqBackendPreferences(
                backend_name="simulator_statevector"
            ),
        ),
    )
    qprogs.append(qprog)

    # Execute
    with ExecutionSession(qprog) as es:
        result = es.minimize(
            cost_function=vqe_hamiltonian,
            initial_params={"params": [0] * num_params},
            max_iteration=500,
        )

    VQE_energy.append(result[-1][0])
    HF_energy.append(molecule.hf_energy)
    exact_energy.append(molecule.fci_energy)

    time2 = time.time()
    duration = time2 - time1
    durations.append(duration)
    print(duration)

5.403509855270386
3.093651056289673
3.4077019691467285
3.023664712905884
3.3719558715820312
3.364731788635254
3.5703423023223877
3.436528205871582
3.3499598503112793
3.333562135696411
3.357571840286255
3.066948175430298


9.940804243087769


10.691434860229492


11.07170820236206


9.649993658065796


11.017099618911743


10.790261507034302


9.763712167739868


9.770148515701294


9.85509705543518


15.547704696655273

# # save the last model to a qmod file
write_qmod(main, "molecular_energy_curve", symbolic_only=False)

Graph Creation

import matplotlib.pyplot as plt

plt.plot(
    distance, VQE_energy, "r--", distance, HF_energy, "bs", distance, exact_energy, "g^"
)

plt.xlabel("distance [Å]")
plt.ylabel("energy [Ha]")
plt.legend(["Classiq VQE", "Hartree-Fock", "Exact solution"])
plt.title("Binding Curve H_{2}")

plt.show()

png

This graph presents the ground state for the \(H_{2}\) molecule as a function of the distance between the two hydrogen atoms. Note that both the HF solution and Classiq VQE present decent results around the global minima. For further distances, Classiq VQE stays close to the exact solution while the HF solution gradually deviates. The source of this lack of correspondence is with the lack of flexible correlations within the HF model, which is enabled within the VQE scope. You can explore more curves, creating graphs for different molecules (even n-dimensional or larger atom assemblies) in a similar fashion.