Creating a Molecule's Potential Energy Curve

View on GitHub

Open this notebook in GitHub to run it yourself

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.

!pip install -qq "classiq[chemistry]"

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)

Output:

[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:

Creating a molecule at changing distances between the atoms.
Constructing a chemistry model for the corresponding Hamiltonian, using the Hartree-Fock initial state and UCC ansatz.
Synthesizing the model to get a quantum program.
Executing the quantum program to extract a solution for the ground energy.
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_qmod
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_qmod(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)
    qprogs.append(qprog)

    # Execute
    result = variational_minimize(
        qprog,
        cost_function=vqe_hamiltonian,
        initial_params={"params": [0] * num_params},
        max_iteration=500,
        verbose=False,
    )

    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)

Output:

9.569384336471558

Output:

8.400216341018677

Output:

6.803268194198608

Output:

8.204976081848145

Output:

7.538207054138184

Output:

8.526277303695679

Output:

7.766880035400391

Output:

7.56188702583313

Output:

7.304831504821777

Output:

8.095475912094116

Output:

8.021538734436035

Output:

7.48616099357605

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()

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.

Algorithms

Applications

Tutorials

Creating a Molecule's Potential Energy Curve

View on GitHub

Definitions and Initialization

Energy Estimations

Graph Creation

View on GitHub

​Definitions and Initialization

​Energy Estimations

​Graph Creation

Definitions and Initialization

Energy Estimations

Graph Creation