Geometry optimization on hydrogen with distinct spatial configurations and non-integer charge

[graciaira]

Hello, Pennylane community!

I am Ira from Indonesia, I attempting to do research on Quantum Chemistry with Variational Quantum Eigensolver using Pennylane. I focus on finding the equilibrium geometry of H2-H2 and H2-H4 and then adding non-integer charges to study error delocalization using Local Density Approximation based on Kohn-Sham Density Functional Theory. However, I have difficulty turning this theory into a valid code.

I follow Pennylane’s tutorial on geometry optimization and given rotations.
I try to get the equilibrium geometry of the H6 molecule on several configurations of H2-H4 as follows.

The results on linear configuration, each atom separate far as follows.

At first, I predicted the stable configuration is on pairs of H2 molecules or chains that have the same bond length and symmetry. I analyzed these spreading results might reflect many-body phenomena on hydrogen chains. As I am doubting and can’t validate this optimized geom of H6, I am trying to work on H4 with H2-H2 configurations, which have more references.

First, I calculate PES, then I use the equilibrium bond for the initial coordinate on geom optimization. In initial configurations, each atom has the same bond and is symmetric.

The following pictures are my results. The bond length of the H-H atom that forms as H2 became longer than the initial equilibrium bond, and the H2-H2 molecule separate far.

This does not seem to be a problem, looking at the final results of linear and rectangular configurations according to the reference, where the distance is differentiated between H2/H-H bonds as d and between H2-H2 molecules as r. Also, considering the tendency of H atoms to form more stable H2 molecules.

For the linear PES H4 curve, the results are in good agreement with the reference “Calculating Electronic Correlation Energy using Linear Depth Quantum Circuits”. For the T-shaped, the final configuration resembles a variation of the paper trapezoidal (T) configuration on paper “Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz”

Seeing this prototyping code, I decided to proceed with the research by going to the deep geometry optimization procedure based on the paper “Variational quantum algorithm for molecular geometry optimization” using Pennylane.

There are many complications behind this procedure that I miss to get a good optimization plot on my use case. So that I reach out to this community and open discussion forum to get some help and follow-up guidance. Thank you in advance.

Warm Regards,
Ira

# PES H4
import pennylane as qml
from pennylane import qchem

# Hartree-Fock state
hf = qml.qchem.hf_state(electrons=4, orbitals=8)

print(hf)

from pennylane import numpy as np

# atomic symbols defining the molecule
symbols = ['H', 'H', 'H', 'H']

# list to store energies
energies = []

# set up a loop to change bond length
r_range = np.arange(0.5, 5.0, 0.25)

# keeps track of points in the potential energy surface
pes_point = 0

for r in r_range:
    # Change only the z coordinate of one atom
    coordinates = np.array([0.0, 0.0, 0.0, 0.0, 0.0, r, 0.0, 0.0, 2*r, 0.0, 0.0, 3*r])

    # Obtain the qubit Hamiltonian
    H, qubits = qchem.molecular_hamiltonian(symbols, coordinates, method='pyscf')

    # define the device, optimizer and circuit
    dev = qml.device("default.qubit", wires=qubits)
    opt = qml.GradientDescentOptimizer(stepsize=0.4)

    @qml.qnode(dev, interface='autograd')
    def circuit(parameters):
        # Prepare the HF state: |11110000>
        qml.BasisState(hf, wires=range(qubits))
        qml.DoubleExcitation(parameters[0], wires=[0, 1, 4, 5])
        qml.DoubleExcitation(parameters[1], wires=[2, 3, 6, 7])
        qml.SingleExcitation(parameters[2], wires=[0, 4])
        qml.SingleExcitation(parameters[3], wires=[1, 5])
        qml.SingleExcitation(parameters[4], wires=[2, 6])
        qml.SingleExcitation(parameters[5], wires=[3, 7])

        return qml.expval(H)  # we are interested in minimizing this expectation value

    # initialize the gate parameters
    params = np.zeros(6, requires_grad=True)

    # initialize with converged parameters from previous point
    if pes_point > 0:
        params = params_old

    prev_energy = 0.0
    for n in range(50):
        # perform optimization step
        params, energy = opt.step_and_cost(circuit, params)

        if np.abs(energy - prev_energy) < 1e-6:
            break
        prev_energy = energy

    # store the converged parameters
    params_old = params
    pes_point = pes_point + 1

    energies.append(energy)

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.plot(r_range, energies)

ax.set(
    xlabel="Bond length (Bohr)",
    ylabel="Total energy (Hartree)",
    title="Potential energy surface for H$_4$ dissociation",
)
ax.grid()
plt.show()

#Let’s use our results to compute the equilibrium bond length and the bond dissociation energy:

# equilibrium energy
e_eq = min(energies)
# energy when atoms are far apart
e_dis = energies[-1]

# Bond dissociation energy
bond_energy = e_dis - e_eq

# Equilibrium bond length
idx = energies.index(e_eq)
bond_length = r_range[idx]

print(f"The equilibrium bond length is {bond_length:.1f} Bohrs")
print(f"The bond dissociation energy is {bond_energy:.6f} Hartrees")

# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903
bond_angs = 1.8 * bohr_angs
print(f"The equilibrium bond length is {bond_angs} Angstrom")

# PES H4 T-shaped
import pennylane as qml
from pennylane import qchem

# Hartree-Fock state
hf = qml.qchem.hf_state(electrons=4, orbitals=8)

print(hf)

from pennylane import numpy as np

# atomic symbols defining the molecule
symbols = ['H', 'H', 'H', 'H']

# list to store energies
energies = []

# set up a loop to change bond length
r_range = np.arange(0.5, 5.0, 0.25)

# keeps track of points in the potential energy surface
pes_point = 0

for r in r_range:
    # Change only the z coordinate of one atom
    coordinates = np.array([0.0, r, 0.0, 0.0, 0.0, 0.0, 0.0, 1/2*r, r, 0.0, 1/2*r, 2*r])

    # Obtain the qubit Hamiltonian
    H, qubits = qchem.molecular_hamiltonian(symbols, coordinates, method='pyscf')

    # define the device, optimizer and circuit
    dev = qml.device("default.qubit", wires=qubits)
    opt = qml.GradientDescentOptimizer(stepsize=0.4)

    @qml.qnode(dev, interface='autograd')
    def circuit(parameters):
        # Prepare the HF state: |11110000>
        qml.BasisState(hf, wires=range(qubits))
        qml.DoubleExcitation(parameters[0], wires=[0, 1, 4, 5])
        qml.DoubleExcitation(parameters[1], wires=[2, 3, 6, 7])
        qml.SingleExcitation(parameters[2], wires=[0, 4])
        qml.SingleExcitation(parameters[3], wires=[1, 5])
        qml.SingleExcitation(parameters[4], wires=[2, 6])
        qml.SingleExcitation(parameters[5], wires=[3, 7])

        return qml.expval(H)  # we are interested in minimizing this expectation value

    # initialize the gate parameters
    params = np.zeros(6, requires_grad=True)

    # initialize with converged parameters from previous point
    if pes_point > 0:
        params = params_old

    prev_energy = 0.0
    for n in range(50):
        # perform optimization step
        params, energy = opt.step_and_cost(circuit, params)

        if np.abs(energy - prev_energy) < 1e-6:
            break
        prev_energy = energy

    # store the converged parameters
    params_old = params
    pes_point = pes_point + 1

    energies.append(energy)

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.plot(r_range, energies)

ax.set(
    xlabel="Bond length (Bohr)",
    ylabel="Total energy (Hartree)",
    title="Potential energy surface for T-shaped H$_4$ dissociation",
)
ax.grid()
plt.show()

#Let’s use our results to compute the equilibrium bond length and the bond dissociation energy:

# equilibrium energy
e_eq = min(energies)
# energy when atoms are far apart
e_dis = energies[-1]

# Bond dissociation energy
bond_energy = e_dis - e_eq

# Equilibrium bond length
idx = energies.index(e_eq)
bond_length = r_range[idx]

print(f"The equilibrium bond length is {bond_length:.1f} Bohrs")
print(f"The bond dissociation energy is {bond_energy:.6f} Hartrees")

# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903
bond_angs = 1.8 * bohr_angs
print(f"The equilibrium bond length is {bond_angs} Angstrom")

#PES H6 linier
import pennylane as qml
from pennylane import qchem

# Hartree-Fock state
hf = qml.qchem.hf_state(electrons=6, orbitals=12)

print(hf)

from pennylane import numpy as np

# atomic symbols defining the molecule
symbols = ['H', 'H', 'H', 'H', 'H', 'H']

# list to store energies
energies = []

# set up a loop to change bond length
r_range = np.arange(0.5, 5.0, 0.25)

# keeps track of points in the potential energy surface
pes_point = 0

for r in r_range:
    # Change only the z coordinate of one atom
    coordinates = np.array([0.0, 0.0, 0.0, 0.0, 0.0, r, 0.0, 0.0, 2*r, 0.0, 0.0, 3*r, 0.0, 0.0, 4*r, 0.0, 0.0, 5*r])

    # Obtain the qubit Hamiltonian
    H, qubits = qchem.molecular_hamiltonian(symbols, coordinates, method='pyscf')

    # define the device, optimizer and circuit
    dev = qml.device("default.qubit", wires=qubits)
    opt = qml.GradientDescentOptimizer(stepsize=0.4)

    @qml.qnode(dev, interface='autograd')
    def circuit(parameters):
        # Prepare the HF state: |111111000000>
        qml.BasisState(hf, wires=range(qubits))
        qml.DoubleExcitation(parameters[0], wires=[0, 1, 6, 7])
        qml.DoubleExcitation(parameters[1], wires=[2, 3, 8, 9])
        qml.DoubleExcitation(parameters[2], wires=[4, 5, 10, 11])
        qml.SingleExcitation(parameters[3], wires=[0, 6])
        qml.SingleExcitation(parameters[4], wires=[1, 7])
        qml.SingleExcitation(parameters[5], wires=[2, 8])
        qml.SingleExcitation(parameters[6], wires=[3, 9])
        qml.SingleExcitation(parameters[7], wires=[4, 10])
        qml.SingleExcitation(parameters[8], wires=[5, 11])

        return qml.expval(H)  # we are interested in minimizing this expectation value

    # initialize the gate parameters
    params = np.zeros(9, requires_grad=True)

    # initialize with converged parameters from previous point
    if pes_point > 0:
        params = params_old

    prev_energy = 0.0
    for n in range(50):
        # perform optimization step
        params, energy = opt.step_and_cost(circuit, params)

        if np.abs(energy - prev_energy) < 1e-6:
            break
        prev_energy = energy

    # store the converged parameters
    params_old = params
    pes_point = pes_point + 1

    energies.append(energy)

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.plot(r_range, energies)

ax.set(
    xlabel="Bond length (Bohr)",
    ylabel="Total energy (Hartree)",
    title="Potential energy surface for linear H$_6$ dissociation",
)
ax.grid()
plt.show()

#Let’s use our results to compute the equilibrium bond length and the bond dissociation energy:

# equilibrium energy
e_eq = min(energies)
# energy when atoms are far apart
e_dis = energies[-1]

# Bond dissociation energy
bond_energy = e_dis - e_eq

# Equilibrium bond length
idx = energies.index(e_eq)
bond_length = r_range[idx]

print(f"The equilibrium bond length is {bond_length:.1f} Bohrs")
print(f"The bond dissociation energy is {bond_energy:.6f} Hartrees")

# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903
bond_angs = 1.8 * bohr_angs
print(f"The equilibrium bond length is {bond_angs} Angstrom")

# optimize geometry H4 linear

from pennylane import numpy as np

symbols = ["H", "H", "H", "H"]
x = np.array([0.000, 0.000, 0.000, 0.000, 0.000, 0.953, 0.000, 0.000, 2*0.953, 0.000, 0.000, 3*0.953], requires_grad=True)

import pennylane as qml


def H(x):
    return qml.qchem.molecular_hamiltonian(symbols, x, charge=0)[0]

hf = qml.qchem.hf_state(electrons=4, orbitals=8)
print(hf)

num_wires = 8
dev = qml.device("lightning.qubit", wires=num_wires)


@qml.qnode(dev, interface="autograd")
def circuit(params, obs, wires):
    qml.BasisState(hf, wires=wires)
    qml.DoubleExcitation(params[0], wires=[0, 1, 4, 5])
    qml.DoubleExcitation(params[1], wires=[2, 3, 6, 7])

    return qml.expval(obs)

def cost(params, x):
    hamiltonian = H(x)
    return circuit(params, obs=hamiltonian, wires=range(num_wires))

def finite_diff(f, x, delta=0.01):
    """Compute the central-difference finite difference of a function"""
    gradient = []

    for i in range(len(x)):
        shift = np.zeros_like(x)
        shift[i] += 0.5 * delta
        res = (f(x + shift) - f(x - shift)) * delta**-1
        gradient.append(res)

    return gradient


def grad_x(params, x):
    grad_h = finite_diff(H, x)
    grad = [circuit(params, obs=obs, wires=range(num_wires)) for obs in grad_h]
    return np.array(grad)

opt_theta = qml.GradientDescentOptimizer(stepsize=0.4)
opt_x = qml.GradientDescentOptimizer(stepsize=0.8)

theta = np.array([0.0, 0.0], requires_grad=True)

from functools import partial

# store the values of the cost function
energy = []

# store the values of the bond length
bond_length = []

# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903

for n in range(200):

    # Optimize the circuit parameters
    theta.requires_grad = True
    x.requires_grad = False
    theta, _ = opt_theta.step(cost, theta, x)

    # Optimize the nuclear coordinates
    x.requires_grad = True
    theta.requires_grad = False
    _, x = opt_x.step(cost, theta, x, grad_fn=grad_x)

    energy.append(cost(theta, x))
    bond_length.append(np.linalg.norm(x[0:4] - x[4:8]) * bohr_angs)

    if n % 4 == 0:
        print(f"Step = {n},  E = {energy[-1]:.8f} Ha,  bond length = {bond_length[-1]:.5f} A")

    # Check maximum component of the nuclear gradient
    if np.max(grad_x(theta, x)) <= 1e-05:
        break

print("\n" f"Final value of the ground-state energy = {energy[-1]:.8f} Ha")

print("\n" "Ground-state equilibrium geometry")
print("%s %4s %8s %8s" % ("symbol", "x", "y", "z"))
for i, atom in enumerate(symbols):
    print(f"  {atom}    {x[3 * i]:.4f}   {x[3 * i + 1]:.4f}   {x[3 * i + 2]:.4f}")

import matplotlib.pyplot as plt

fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(12)

# Add energy plot on column 1
E_fci = -2.175339244358
E_vqe = np.array(energy)
ax1 = fig.add_subplot(121)
ax1.plot(range(n + 1), E_vqe - E_fci, "go", ls="dashed")
ax1.plot(range(n + 1), np.full(n + 1, 0.001), color="red")
ax1.set_xlabel("Optimization step", fontsize=13)
ax1.set_ylabel("$E_{VQE} - E_{FCI}$ (Hartree)", fontsize=13)
ax1.text(5, 0.0013, r"Chemical accuracy", fontsize=13)
plt.yscale("log")
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)

# Add bond length plot on column 2
d_fci = 0.953
ax2 = fig.add_subplot(122)
ax2.plot(range(n + 1), bond_length, "go", ls="dashed")
ax2.plot(range(n + 1), np.full(n + 1, d_fci), color="red")
ax2.set_ylim([0.885, 0.97])
ax2.set_xlabel("Optimization step", fontsize=13)
ax2.set_ylabel("bond length ($\AA$)", fontsize=13)
ax2.text(5, 0.9865, r"Equilibrium bond length", fontsize=13)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)

plt.subplots_adjust(wspace=0.3)
plt.show()

#Plot belum disesuaikan, ini untuk test dan lihat polanya
import matplotlib.pyplot as plt
E_vqe = np.array(energy)
plt.plot(bond_length, E_vqe)

plt.plot(range(n + 1), bond_length)

E_vqe = np.array(energy)
plt.plot(range(n + 1), E_vqe - E_fci)

np.savetxt('./optimize_geom_h4_linear_bondlength.txt', bond_length)
np.savetxt('./optimize_geom_h4_linear_evqe.txt', E_vqe)

# optimize geometry H4 T-shaped

from pennylane import numpy as np

symbols = ["H", "H", "H", "H"]
x = np.array([0.0, 0.953, 0.0, 0.0, 0.0, 0.0, 0.0, 1/2*0.953, 0.953, 0.0, 1/2*0.953, 2*0.953], requires_grad=True)

import pennylane as qml


def H(x):
    return qml.qchem.molecular_hamiltonian(symbols, x, charge=0)[0]

hf = qml.qchem.hf_state(electrons=4, orbitals=8)
print(hf)

num_wires = 8
dev = qml.device("lightning.qubit", wires=num_wires)


@qml.qnode(dev, interface="autograd")
def circuit(params, obs, wires):
    qml.BasisState(hf, wires=wires)
    qml.DoubleExcitation(params[0], wires=[0, 1, 4, 5])
    qml.DoubleExcitation(params[1], wires=[2, 3, 6, 7])

    return qml.expval(obs)

def cost(params, x):
    hamiltonian = H(x)
    return circuit(params, obs=hamiltonian, wires=range(num_wires))

def finite_diff(f, x, delta=0.01):
    """Compute the central-difference finite difference of a function"""
    gradient = []

    for i in range(len(x)):
        shift = np.zeros_like(x)
        shift[i] += 0.5 * delta
        res = (f(x + shift) - f(x - shift)) * delta**-1
        gradient.append(res)

    return gradient


def grad_x(params, x):
    grad_h = finite_diff(H, x)
    grad = [circuit(params, obs=obs, wires=range(num_wires)) for obs in grad_h]
    return np.array(grad)

opt_theta = qml.GradientDescentOptimizer(stepsize=0.4)
opt_x = qml.GradientDescentOptimizer(stepsize=0.8)

theta = np.array([0.0, 0.0], requires_grad=True)

from functools import partial

# store the values of the cost function
energy = []

# store the values of the bond length
bond_length = []

# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903

for n in range(300):

    # Optimize the circuit parameters
    theta.requires_grad = True
    x.requires_grad = False
    theta, _ = opt_theta.step(cost, theta, x)

    # Optimize the nuclear coordinates
    x.requires_grad = True
    theta.requires_grad = False
    _, x = opt_x.step(cost, theta, x, grad_fn=grad_x)

    energy.append(cost(theta, x))
    bond_length.append(np.linalg.norm(x[0:4] - x[4:8]) * bohr_angs)

    if n % 4 == 0:
        print(f"Step = {n},  E = {energy[-1]:.8f} Ha,  bond length = {bond_length[-1]:.5f} A")

    # Check maximum component of the nuclear gradient
    if np.max(grad_x(theta, x)) <= 1e-05:
        break

print("\n" f"Final value of the ground-state energy = {energy[-1]:.8f} Ha")

print("\n" "Ground-state equilibrium geometry")
print("%s %4s %8s %8s" % ("symbol", "x", "y", "z"))
for i, atom in enumerate(symbols):
    print(f"  {atom}    {x[3 * i]:.4f}   {x[3 * i + 1]:.4f}   {x[3 * i + 2]:.4f}")

#Plot belum disesuaikan, ini untuk test dan lihat polanya
import matplotlib.pyplot as plt
E_vqe = np.array(energy)
plt.plot(bond_length, E_vqe)

plt.plot(range(n + 1), bond_length)

E_fci = -2.095776493960
E_vqe = np.array(energy)
plt.plot(range(n + 1), E_vqe - E_fci)

import matplotlib.pyplot as plt

fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(12)

# Add energy plot on column 1
E_fci = -2.095776493960
E_vqe = np.array(energy)
ax1 = fig.add_subplot(121)
ax1.plot(range(n + 1), E_vqe - E_fci, "go", ls="dashed")
ax1.plot(range(n + 1), np.full(n + 1, 0.001), color="red")
ax1.set_xlabel("Optimization step", fontsize=13)
ax1.set_ylabel("$E_{VQE} - E_{FCI}$ (Hartree)", fontsize=13)
ax1.text(5, 0.0013, r"Chemical accuracy", fontsize=13)
plt.yscale("log")
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)

# Add bond length plot on column 2
d_fci = 0.953
ax2 = fig.add_subplot(122)
ax2.plot(range(n + 1), bond_length, "go", ls="dashed")
ax2.plot(range(n + 1), np.full(n + 1, d_fci), color="red")
ax2.set_ylim([0.885, 2.97])
ax2.set_xlabel("Optimization step", fontsize=13)
ax2.set_ylabel("bond length ($\AA$)", fontsize=13)
ax2.text(5, 0.9865, r"Equilibrium bond length", fontsize=13)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)

plt.subplots_adjust(wspace=0.3)
plt.show()

np.savetxt('./optimize_geom_h4_tshaped_bondlength.txt', bond_length)
np.savetxt('./optimize_geom_h4_tshaped_evqe.txt', E_vqe)

# Energy FCI H4

import pyscf

mol = pyscf.M(
    atom = 'H 0.000 0.000 0.000; H 0.000 0.000 0.953; H 0.000 0.000 2*0.953; H 0.000 0.000 3*0.953',  # in Angstrom
    basis = 'sto3g',
    symmetry = True,
)
myhf = mol.RHF().run()

#
# create an FCI solver based on the SCF object
#
cisolver = pyscf.fci.FCI(myhf)
print('E(FCI) = %.12f' % cisolver.kernel()[0])

#
# create an FCI solver based on the SCF object
#
myuhf = mol.UHF().run()
cisolver = pyscf.fci.FCI(myuhf)
print('E(UHF-FCI) = %.12f' % cisolver.kernel()[0])

#
# create an FCI solver based on the given orbitals and the num. electrons and
# spin of the mol object
#
cisolver = pyscf.fci.FCI(mol, myhf.mo_coeff)
print('E(FCI) linear = %.12f' % cisolver.kernel()[0])

import pyscf

mol = pyscf.M(
    atom = 'H 0.000 0.953 0.000; H 0.000 0.000 0.000; H 0.000 1/2*0.953 0.953; H 0.000 1/2*0.953 2*0.953',  # in Angstrom
    basis = 'sto3g',
    symmetry = True,
)
myhf = mol.RHF().run()

#
# create an FCI solver based on the SCF object
#
cisolver = pyscf.fci.FCI(myhf)
print('E(FCI) = %.12f' % cisolver.kernel()[0])

#
# create an FCI solver based on the SCF object
#
myuhf = mol.UHF().run()
cisolver = pyscf.fci.FCI(myuhf)
print('E(UHF-FCI) = %.12f' % cisolver.kernel()[0])

#
# create an FCI solver based on the given orbitals and the num. electrons and
# spin of the mol object
#
cisolver = pyscf.fci.FCI(mol, myhf.mo_coeff)
print('E(FCI) T shaped = %.12f' % cisolver.kernel()[0])

# optimize geometry H6 linear

from pennylane import numpy as np

symbols = ["H", "H", "H", "H", "H", "H"]
x = np.array([0.000, 0.000, 0.000, 0.000, 0.000, 0.953, 0.000, 0.000, 2*0.953, 0.000, 0.000, 3*0.953, 0.000, 0.000, 4*0.953, 0.000, 0.000, 5*0.953], requires_grad=True)

import pennylane as qml


def H(x):
    return qml.qchem.molecular_hamiltonian(symbols, x, charge=0)[0]

hf = qml.qchem.hf_state(electrons=6, orbitals=12)
print(hf)

num_wires = 12
dev = qml.device("lightning.qubit", wires=num_wires)


@qml.qnode(dev, interface="autograd")
def circuit(params, obs, wires):
    qml.BasisState(hf, wires=wires)
    qml.DoubleExcitation(params[0], wires=[0, 1, 6, 7])
    qml.DoubleExcitation(params[1], wires=[2, 3, 8, 9])
    qml.DoubleExcitation(params[2], wires=[4, 5, 10, 11])

    return qml.expval(obs)

def cost(params, x):
    hamiltonian = H(x)
    return circuit(params, obs=hamiltonian, wires=range(num_wires))

def finite_diff(f, x, delta=0.01):
    """Compute the central-difference finite difference of a function"""
    gradient = []

    for i in range(len(x)):
        shift = np.zeros_like(x)
        shift[i] += 0.5 * delta
        res = (f(x + shift) - f(x - shift)) * delta**-1
        gradient.append(res)

    return gradient


def grad_x(params, x):
    grad_h = finite_diff(H, x)
    grad = [circuit(params, obs=obs, wires=range(num_wires)) for obs in grad_h]
    return np.array(grad)

opt_theta = qml.GradientDescentOptimizer(stepsize=0.4)
opt_x = qml.GradientDescentOptimizer(stepsize=0.8)

theta = np.array([0.000, 0.000, 0.000], requires_grad=True)

from functools import partial

# store the values of the cost function
energy = []

# store the values of the bond length
bond_length = []

# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903

for n in range(300):

    # Optimize the circuit parameters
    theta.requires_grad = True
    x.requires_grad = False
    theta, _ = opt_theta.step(cost, theta, x)

    # Optimize the nuclear coordinates
    x.requires_grad = True
    theta.requires_grad = False
    _, x = opt_x.step(cost, theta, x, grad_fn=grad_x)

    energy.append(cost(theta, x))
    bond_length.append(np.linalg.norm(x[0:6] - x[6:12]) * bohr_angs)

    if n % 4 == 0:
        print(f"Step = {n},  E = {energy[-1]:.8f} Ha,  bond length = {bond_length[-1]:.5f} A")

    # Check maximum component of the nuclear gradient
    if np.max(grad_x(theta, x)) <= 1e-05:
        break

print("\n" f"Final value of the ground-state energy = {energy[-1]:.8f} Ha")

# optimize geometry H6 T-shaped

from pennylane import numpy as np

symbols = ["H", "H", "H", "H", "H", "H"]
x = np.array([0.0, 0.953, 0.0, 0.0, 0.0, 0.0, 0.0, 1/2*0.953, 0.953, 0.0, 1/2*0.953, 2*0.953, 0.0, 1/2*0.953, 3*0.953, 0.0, 1/2*0.953, 4*0.953], requires_grad=True)

import pennylane as qml


def H(x):
    return qml.qchem.molecular_hamiltonian(symbols, x, charge=0)[0]

hf = qml.qchem.hf_state(electrons=6, orbitals=12)
print(hf)

num_wires = 12
dev = qml.device("lightning.qubit", wires=num_wires)


@qml.qnode(dev, interface="autograd")
def circuit(params, obs, wires):
    qml.BasisState(hf, wires=wires)
    qml.DoubleExcitation(params[0], wires=[0, 1, 6, 7])
    qml.DoubleExcitation(params[1], wires=[2, 3, 8, 9])
    qml.DoubleExcitation(params[2], wires=[4, 5, 10, 11])

    return qml.expval(obs)

def cost(params, x):
    hamiltonian = H(x)
    return circuit(params, obs=hamiltonian, wires=range(num_wires))

def finite_diff(f, x, delta=0.01):
    """Compute the central-difference finite difference of a function"""
    gradient = []

    for i in range(len(x)):
        shift = np.zeros_like(x)
        shift[i] += 0.5 * delta
        res = (f(x + shift) - f(x - shift)) * delta**-1
        gradient.append(res)

    return gradient


def grad_x(params, x):
    grad_h = finite_diff(H, x)
    grad = [circuit(params, obs=obs, wires=range(num_wires)) for obs in grad_h]
    return np.array(grad)

opt_theta = qml.GradientDescentOptimizer(stepsize=0.4)
opt_x = qml.GradientDescentOptimizer(stepsize=0.8)

theta = np.array([0.0, 0.0, 0.0], requires_grad=True)

from functools import partial

# store the values of the cost function
energy = []

# store the values of the bond length
bond_length = []

# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903

for n in range(300):

    # Optimize the circuit parameters
    theta.requires_grad = True
    x.requires_grad = False
    theta, _ = opt_theta.step(cost, theta, x)

    # Optimize the nuclear coordinates
    x.requires_grad = True
    theta.requires_grad = False
    _, x = opt_x.step(cost, theta, x, grad_fn=grad_x)

    energy.append(cost(theta, x))
    bond_length.append(np.linalg.norm(x[0:6] - x[6:12]) * bohr_angs)

    if n % 4 == 0:
        print(f"Step = {n},  E = {energy[-1]:.8f} Ha,  bond length = {bond_length[-1]:.5f} A")

    # Check maximum component of the nuclear gradient
    if np.max(grad_x(theta, x)) <= 1e-05:
        break

print("\n" f"Final value of the ground-state energy = {energy[-1]:.8f} Ha")

And, finally, make sure to include the versions of your packages. Specifically, show us the output of qml.about().

Name: PennyLane
Version: 0.36.0
Summary: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Home-page: GitHub - PennyLaneAI/pennylane: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Author:
Author-email:
License: Apache License 2.0
Location: /home/ira/anaconda3/envs/penny/lib/python3.11/site-packages
Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, pennylane-lightning, requests, rustworkx, scipy, semantic-version, toml, typing-extensions
Required-by: PennyLane_Lightning

Platform info: Linux-6.5.0-41-generic-x86_64-with-glibc2.35
Python version: 3.11.9
Numpy version: 1.26.4
Scipy version: 1.13.1
Installed devices:

• lightning.qubit (PennyLane_Lightning-0.36.0)
• default.clifford (PennyLane-0.36.0)
• default.gaussian (PennyLane-0.36.0)
• default.mixed (PennyLane-0.36.0)
• default.qubit (PennyLane-0.36.0)
• default.qubit.autograd (PennyLane-0.36.0)
• default.qubit.jax (PennyLane-0.36.0)
• default.qubit.legacy (PennyLane-0.36.0)
• default.qubit.tf (PennyLane-0.36.0)
• default.qubit.torch (PennyLane-0.36.0)
• default.qutrit (PennyLane-0.36.0)
• default.qutrit.mixed (PennyLane-0.36.0)
• null.qubit (PennyLane-0.36.0)

I would like to tag Ms. Ivana so my thread could be found easily @Ivana_at_Xanadu
Thank you again, Ivana

#geometryoptimization

[CatalinaAlbornoz]

Hi @graciaira, welcome to the Forum!

It looks like you’re doing some awesome research. I love your geometry plots.

However I’m having trouble understanding what your issue is. My understanding is that you’re getting no errors but the results aren’t what you expect. Can you please let me know what the result was vs what you expected? For example: “When simulating the PES for H4 I expected to get a minimum energy of -2Ha at a bond length ob 1.8 Bohr but instead I got a minimal energy of X at a bond length of Y.”

Something like this (maybe starting with a single configuration) can help me understand the problem and look into it.

[graciaira]

Hello, Catalina, thankyou for welcoming me. I am glad to have you here. How do you!

Yes, correct i am having a research on molecular simulation.

Actually we plan to using LDA-DFT approach on VQE to exposing delocalization error. Getting equilibrium geometry is basic for us to understand how spatial-dependent configurations effecting it’s energy level pior stepping to more complex analysis.

I plot the geom on xCrysden by it’s coordinate, i also thought they look cute together in one frame, hoho ^.^

That’s true i am not having any error, Catalina. It also possible for me to start writing the results based on the pattern it occurs. But, it’s important for me to go deeper than the skin and reveals subtle error when it comes to quantify it’s optimized bond length and study the stability phenomena here. Considering i would adding fractional charge as well.

Let’s go with H2-H2 linear configuration first, Catalina.

import pennylane as qml

from pennylane import qchem

# Hartree-Fock state

hf = qml.qchem.hf_state(electrons=4, orbitals=8)

print(hf)

from pennylane import numpy as np

# atomic symbols defining the molecule

symbols = ['H', 'H', 'H', 'H']

# list to store energies

energies = []

# set up a loop to change bond length

r_range = np.arange(0.5, 5.0, 0.25)

# keeps track of points in the potential energy surface

pes_point = 0
for r in r_range:
   # Change only the z coordinate of one atom
   coordinates = np.array([0.0, 0.0, 0.0, 0.0, 0.0, r, 0.0, 0.0, 2*r, 0.0, 0.0, 3*r])

   # Obtain the qubit Hamiltonian
   H, qubits = qchem.molecular_hamiltonian(symbols, coordinates, method='pyscf')

   # define the device, optimizer and circuit
   dev = qml.device("default.qubit", wires=qubits)
   opt = qml.GradientDescentOptimizer(stepsize=0.4)

   @qml.qnode(dev, interface='autograd')
   def circuit(parameters):
       # Prepare the HF state: |11110000>
       qml.BasisState(hf, wires=range(qubits))
       qml.DoubleExcitation(parameters[0], wires=[0, 1, 4, 5])
       qml.DoubleExcitation(parameters[1], wires=[2, 3, 6, 7])
       qml.SingleExcitation(parameters[2], wires=[0, 4])
       qml.SingleExcitation(parameters[3], wires=[1, 5])
       qml.SingleExcitation(parameters[4], wires=[2, 6])
       qml.SingleExcitation(parameters[5], wires=[3, 7])

       return qml.expval(H)  # we are interested in minimizing this expectation value

   # initialize the gate parameters
   params = np.zeros(6, requires_grad=True)

   # initialize with converged parameters from previous point
   if pes_point > 0:
       params = params_old

   prev_energy = 0.0
   for n in range(50):
       # perform optimization step
       params, energy = opt.step_and_cost(circuit, params)

       if np.abs(energy - prev_energy) < 1e-6:
           break
       prev_energy = energy

   # store the converged parameters
   params_old = params
   pes_point = pes_point + 1

   energies.append(energy)
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.plot(r_range, energies)

ax.set(
   xlabel="Bond length (Bohr)",
   ylabel="Total energy (Hartree)",
   title="Potential energy surface for H$_4$ dissociation",
)
ax.grid()
plt.show()
#Let’s use our results to compute the equilibrium bond length and the bond dissociation energy:

# equilibrium energy

e_eq = min(energies)

# energy when atoms are far apart

e_dis = energies[-1]

# Bond dissociation energy

bond_energy = e_dis - e_eq

# Equilibrium bond length

idx = energies.index(e_eq)

bond_length = r_range[idx]

print(f"The equilibrium bond length is {bond_length:.1f} Bohrs")

print(f"The bond dissociation energy is {bond_energy:.6f} Hartrees")

# Factor to convert from Bohrs to Angstroms

bohr_angs = 0.529177210903

bond_angs = 1.8 * bohr_angs

print(f"The equilibrium bond length is {bond_angs} Angstrom")

I didn’t make much change from tutorial, just the array coordinates, and single-double excitation gates.

First of all, i would like to validate, is my code given above correct? I arrange the gates based on my understanding that we need to move each atom to unoccupied orbitals.

The equilibrium bond gives me 1.8 Bohr and I convert it to Angstrom as the input geom to optimize the H2-H2 linear.

from pennylane import numpy as np

symbols = ["H", "H", "H", "H"]

x = np.array([0.000, 0.000, 0.000, 0.000, 0.000, 0.953, 0.000, 0.000, 2*0.953, 0.000, 0.000, 3*0.953], requires_grad=True)

import pennylane as qml

def H(x):
    return qml.qchem.molecular_hamiltonian(symbols, x, charge=0)[0]

hf = qml.qchem.hf_state(electrons=4, orbitals=8)

print(hf)

num_wires = 8
dev = qml.device("lightning.qubit", wires=num_wires)


@qml.qnode(dev, interface="autograd")
def circuit(params, obs, wires):
    qml.BasisState(hf, wires=wires)
    qml.DoubleExcitation(params[0], wires=[0, 1, 4, 5])
    qml.DoubleExcitation(params[1], wires=[2, 3, 6, 7])

    return qml.expval(obs)

def cost(params, x):
    hamiltonian = H(x)
    return circuit(params, obs=hamiltonian, wires=range(num_wires))

def finite_diff(f, x, delta=0.01):
    """Compute the central-difference finite difference of a function"""
    gradient = []

    for i in range(len(x)):
        shift = np.zeros_like(x)
        shift[i] += 0.5 * delta
        res = (f(x + shift) - f(x - shift)) * delta**-1
        gradient.append(res)

    return gradient

def grad_x(params, x):
    grad_h = finite_diff(H, x)
    grad = [circuit(params, obs=obs, wires=range(num_wires)) for obs in grad_h]
    return np.array(grad)

opt_theta = qml.GradientDescentOptimizer(stepsize=0.4)

opt_x = qml.GradientDescentOptimizer(stepsize=0.8)

theta = np.array([0.0, 0.0], requires_grad=True)

from functools import partial

# store the values of the cost function
energy = []

# store the values of the bond length
bond_length = []

# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903

for n in range(200):

    # Optimize the circuit parameters
    theta.requires_grad = True
    x.requires_grad = False
    theta, _ = opt_theta.step(cost, theta, x)

    # Optimize the nuclear coordinates
    x.requires_grad = True
    theta.requires_grad = False
    _, x = opt_x.step(cost, theta, x, grad_fn=grad_x)

    energy.append(cost(theta, x))
    bond_length.append(np.linalg.norm(x[0:4] - x[4:8]) * bohr_angs)

    if n % 4 == 0:
        print(f"Step = {n},  E = {energy[-1]:.8f} Ha,  bond length = {bond_length[-1]:.5f} A")

    # Check maximum component of the nuclear gradient
    if np.max(grad_x(theta, x)) <= 1e-05:
        break

print("\n" f"Final value of the ground-state energy = {energy[-1]:.8f} Ha")

print("\n" "Ground-state equilibrium geometry")
print("%s %4s %8s %8s" % ("symbol", "x", "y", "z"))
for i, atom in enumerate(symbols):
    print(f"  {atom}    {x[3 * i]:.4f}   {x[3 * i + 1]:.4f}   {x[3 * i + 2]:.4f}")

#Plot belum disesuaikan, ini untuk test dan lihat polanya

import matplotlib.pyplot as plt

E_vqe = np.array(energy)

plt.plot(bond_length, E_vqe)

plt.plot(range(n + 1), bond_length)

E_vqe = np.array(energy)

plt.plot(range(n + 1), E_vqe - E_fci)

import matplotlib.pyplot as plt

fig = plt.figure()

fig.set_figheight(5)

fig.set_figwidth(12)

# Add energy plot on column 1

E_fci = -2.175339244358

E_vqe = np.array(energy)

ax1 = fig.add_subplot(121)

ax1.plot(range(n + 1), E_vqe - E_fci, "go", ls="dashed")

ax1.plot(range(n + 1), np.full(n + 1, 0.001), color="red")

ax1.set_xlabel("Optimization step", fontsize=13)

ax1.set_ylabel("$E_{VQE} - E_{FCI}$ (Hartree)", fontsize=13)

ax1.text(5, 0.0013, r"Chemical accuracy", fontsize=13)

plt.yscale("log")

plt.xticks(fontsize=12)

plt.yticks(fontsize=12)

# Add bond length plot on column 2

d_fci = 0.953

ax2 = fig.add_subplot(122)

ax2.plot(range(n + 1), bond_length, "go", ls="dashed")

ax2.plot(range(n + 1), np.full(n + 1, d_fci), color="red")

ax2.set_ylim([0.885, 0.97])

ax2.set_xlabel("Optimization step", fontsize=13)

ax2.set_ylabel("bond length ($\AA$)", fontsize=13)

ax2.text(5, 0.9865, r"Equilibrium bond length", fontsize=13)

plt.xticks(fontsize=12)

plt.yticks(fontsize=12)

plt.subplots_adjust(wspace=0.3)

plt.show()

Now, the challenge begin, the error occurs when i gives single excitation, so i simply remove all singles. And the optimization doesn’t actually converges after reach maximum iteration with 200 steps.
This is the optimized geom given:

Ground-state equilibrium geometry
symbol x y z
H 0.0000 0.0000 -1.6867
H 0.0000 0.0000 -0.3834
H 0.0000 0.0000 3.1861
H 0.0000 0.0000 4.6019

The visualization looks okay, if we consider H2 tend to comes together and stabilize itself as molecule and separate from other H2 pairs. But if we pay attention to the distances here, all exceed the equilibrium bonds, that is around 1.8 bohr or 0.953 A.
I am also start confusing which one on bohr which one on angstroms. Is it correct if i put the initial geom on Angstorms? And the printed is on Angstorms right?

Then next challenge is:
I couldn’t analyze the correlation of equilibrium geom with energy level because the quantities differ as well from initial PES curve. Here, the minimum seems on -2.24 Ha at a bond length 0.9 A, while PES is less than -2Ha at the bond length 0.95A
(ah just realize i didn’t print the energy but based on the curve is below -2Ha)

Another challenge is,
Is it okay the optimization zig-zag at the tails? I worried there’s something i miss.

And regarding accuracy, I am not sure how to calculate FCI, i use this code

import pyscf

mol = pyscf.M(

atom = 'H 0.000 0.000 0.000; H 0.000 0.000 0.953; H 0.000 0.000 2*0.953; H 0.000 0.000 3*0.953', # in Angstrom

basis = 'sto3g',

symmetry = True,

)

myhf = mol.RHF().run()

#

# create an FCI solver based on the SCF object

#

cisolver = pyscf.fci.FCI(myhf)

print('E(FCI) = %.12f' % cisolver.kernel()[0])

#

# create an FCI solver based on the SCF object

#

myuhf = mol.UHF().run()

cisolver = pyscf.fci.FCI(myuhf)

print('E(UHF-FCI) = %.12f' % cisolver.kernel()[0])

#

# create an FCI solver based on the given orbitals and the num. electrons and

# spin of the mol object

#

cisolver = pyscf.fci.FCI(mol, myhf.mo_coeff)

print('E(FCI) = %.12f' % cisolver.kernel()[0])

gives me results

converged SCF energy = -2.11330013628225
E(FCI) = -2.175339244358
converged SCF energy = -2.11330013628175 <S^2> = 3.5080827e-12 2S+1 = 1
E(UHF-FCI) = -2.175339244359
E(FCI) = -2.175339244358

And i plot the curve like this

Ah, about equilibrium bond length plot on red straight line on optimization step curve, i confuse which one to insert. Is it by FCI calculation or by PES i obtained at the very beginning? Here i insert based on PES.
(Appologies for the scaling, i didn’t modifying those yet, but full curve already shown above).

Catalina, we see a lot issue already from the parameters set up, i need your help for sure to proceed as i generate many data for distinct spatial configuration.

Thank you so much for being there.

Warm Regards,
Ira

[CatalinaAlbornoz]

Hi @graciaira,

Thank you for clarifying your questions!

1. qchem.molecular_hamiltonian takes coordinates in atomic units. The conversion is: 1 Angstrom = 1.8897259885789 a.u. This seems to be the same as what you’re using in the first section of your code so I see no issue there. However, you seem to be using Angstrom instead of atomic units in the second part of your code, so this would result in a different simulation from the one you want to do.
2. I noticed you’re manually calculating finite differences in the second part of your code. While you can indeed do this, it can lead to numerical error. You can use other methods such as grad as shown in the last section of this demo. In that demo we perform a full optimization where the circuit parameters, the nuclear coordinates and the basis set parameters are all differentiable parameters that can be optimized simultaneously.
3. You’re using a GradientDescentOptimizer with a stepsize that seems a bit large and can lead to instabilities like the ones you’re seeing. I’d recommend trying qml.AdamOptimizer, which has an adaptive learning rate and can help you in these cases.
4. You can download the FCI energy (and a lot of other data) from this dataset. You just have to set the configuration that you want and the information that you want to download and then just copy and run the code snippet that is shown in the download button. The code below shows the result of downloading the fci energies for the full range of bondlengths available in the dataset.

import pennylane as qml
import numpy as np
import matplotlib.pyplot as plt

H4datasets = qml.data.load("qchem", molname="H4", bondlength="full", basis="STO-3G", attributes=["fci_energy"])

data = np.array([[float(elem.bondlength),elem.fci_energy] for elem in H4datasets])
bondlengths = data[:,0]
fci_energies = data[:,1]

fig, ax = plt.subplots()
ax.scatter(bondlengths, fci_energies)

ax.set(
   xlabel="Bond length (Bohr)",
   ylabel="FCI energy (Hartree)"
)
ax.grid()
plt.show()

Let me know if this answers your questions!