Optimizing water's molecular geometry


I’m working on replicating the Pennylane quantum chemistry tutorial: “Optimization of Molecular geometry”, but for water instead of the Trihydronium molecule they use. However, when I’m running the final part of the code, it seems that the gradient for the nuclear coordinates starts to diverge, despite the starting values for water’s bond length and angle being very close to their ground state values.

I’m wondering if anyone has done tried a similar adapation, and can help me out by looking at my code and see where I’m going wrong with the gradient’s calculation! Thanks so much!

Hey @ImranNasrullah, welcome to the forum!

It’s a little hard to tell what’s going on without seeing your full code, but make sure that you’re setting the problem up in a way that makes sense for water (not like what’s in the demo with \text{H}_3). Check this demo out for how to set up the Hamiltonian: Building molecular Hamiltonians | PennyLane Demos

Let me know if this helps!

I believe I was the parameters that I put into the function qchem.molecular_hamiltonian are the appropriate ones to reflect water (i.e correct number of electrons and active orbitals, charge, multiplicity).

Here is the github for the code, where I’ve adapted the majority of it from the tutorial ‘Optimizing Molecular Geometries’. However, I also included the ‘Manual Construction’ section of the ‘Adaptive Circuits for Quantum Chemistry’ tutorial to adaptively build the appropriate quantum circuit for water.

I’ve tried varying up the the ‘step size’ parameter for the opt_x optimizer to make it to 0.01. In this case, I don’t think the gradient converges, but if I make the step-size to 0.8 (which is what the tutorial uses), the calculated gradient for the nuclear coordinates diverges. So I’m not sure where the error is.

In addition, here are the versions for the libraries I used for the code:

Hi @ImranNasrullah,

In your code you’re changing many things at the same time so it’s very hard to debug. I would recommend to stay as close as possible to the demo (eg using the molecule in the demo) and changing one item at a time in order to find the point of failure.

On the other hand, I’m not sure why you choose the number of electrons and orbitals as you do. As you can see in the demo on Building molecular Hamiltonians you probably want to use:

electrons = 10
orbitals = 7

I hope this helps you! :smile:

Sounds good, I can go ahead and use the original molecule first, the Tri-hydrogen Cation.

The reason that I chose the number of electrons and orbitals to be 8 and 6 respectively is to see if I can follow along with the paper: https://arxiv.org/pdf/2106.13840.pdf
from which I believe the tutorial is based. In page 5 of the paper, (page 5, Table I) they use 8 active electrons for water, and I believe they use 6 active orbitals to exclude the non-valence s-orbital of Oxygen, so that’s what I wanted to use as well.

Ok, I’ve adapted my code to use the same molecule in the tutorial, H3+! I’ve also uploaded this script to github as well.

After changing the relevant molecular parameters from H2O to H3+, I’m actually able to reproduce the exact same output as the ‘Optimizing molecular geometries’ tutorial (including being able to select the same excitation gates through Adapt-VQE that was hardcoded in the tutorial). So if it’s working for one molecule, I’m a bit stumped for why it would not be working for water.

And just to clarify for the water script, almost everything is near the same as the tutorial. Also, the only major addition that I added was a way to keep track of water’s bond angle while the algorithm was optimizing the nuclear coordinates.

Hi @ImranNasrullah ,

I’m glad that it works for H3+!

From Table I in page 5 of the paper you referenced you can see that for H2O the number of active molecular spin orbitals is the same as the number of qubits, which is 12 for H2O. The number of active electrons, Ne, is 8 in the paper.

If you change this does the result improve?

In the function qchem.molecular_hamiltonian I did set the ‘active_electrons’ parameter to 8, and the ‘active_orbitals’ parameter to 6. But since the ‘active_orbitals’ parameter really specifies molecular orbitals, I believe this indeed refers to 12 spin-molecular orbitals.

Despite this, the code code will not converge

Hey @ImranNasrullah,

It’s difficult to say what’s going on in your code (there’s a lot going on :sweat_smile:). I made this code that is heavily based on this demo. The water Hamiltonian has 14 qubits with the default basis set (sto-3g), so it’s taking me a bit to run on my humble laptop. But, maybe it will give you some ideas for why your code is stalling out.

import pennylane as qml
import pennylane.numpy as np

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

    return gradient

def optimize_coordinates():

    electrons = 10
    orbitals = 7
    active_orbitals = 4
    active_electrons = 4
    qubits = 14 # from previous calculation

    def water_hamiltonian(coordinates):
        H = qml.qchem.molecular_hamiltonian(
            ["H", "O", "H"],

        return H

    init_coords = np.array([-1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0], requires_grad=True)

    hf_state = qml.qchem.hf_state(active_electrons, active_orbitals)
    singles, doubles = qml.qchem.excitations(electrons, qubits)

    init_params = np.random.uniform(0, 1, size=(len(singles) + len(doubles),), requires_grad=True)

    coords = init_coords
    params = init_params

    hf_state = qml.qchem.hf_state(electrons, qubits)
    singles, doubles = qml.qchem.excitations(electrons, qubits)

    dev = qml.device("default.qubit", wires=qubits)

    def ansatz(params, obs):
        qml.templates.AllSinglesDoubles(params, range(qubits), hf_state, singles, doubles)
        return qml.expval(obs)

    def cost(params, coords):
        H = water_hamiltonian(coords)
        return ansatz(params, H)
    def grad_coords(params, coords):
        grad_h = finite_diff(water_hamiltonian, coords)
        grad = [ansatz(params, obs) for obs in grad_h]
        return np.array(grad)

    opt_params = qml.AdamOptimizer(0.05)
    opt_coords = qml.AdamOptimizer(0.1)

    for i in range(500):
        coords.requires_grad = False
        params.requires_grad = True
        params, _ = opt_params.step(cost, params, coords)

        coords.requires_grad = True
        params.requires_grad = False
        _, coords = opt_coords.step(cost, params, coords, grad_fn=grad_coords)

        print("cost", cost(params, coords))
    return coords


Hope this helps!

Hello Isaac,
Thanks for the response! I tried out the code, and I believe it wasn’t able to optimize (the bond length and angle ended up diverging). However, I notice this code is the same as the code I wrote myself (except I also had the Adapt-VQE algorithm to select the different excitation gates), you included all the gates. Despite this, I still can’t get my code to optimize for water

Any ideas for where to go from here?

Furthermore, I have another script where I optimizing just the bond length r, and the bond angle for H3+. When i do this, the results are very consistent:

and this as well:

As you can see here, it seems to be that both bond length and angle for the H3+ are converging correctly.

However, when I repurpose this exact same script for H2O (the only thing I change is the “symbols” argument to be “O,H,H” and change the molecular_hamiltonian(). When I fix the angle, here is the output for the bond length of water:

It diverges from the accepted value, which is around 1

Any thoughts on how I can approach this?? The bond length is diverging from the accepted value for water!

Hey @ImranNasrullah,

I got in touch with one of our qchem specialists and they’ll write back here by early next week :slight_smile:

Hi @ImranNasrullah and thanks for using PennyLane!

Could you please run the code in this demo but just change the molecular symbols and Hamiltonian as provided in the following? Depending on your results, we go further and find the problem.

symbols = ["O", "H", "H"]
x = np.array([0.028, 0.054, 0.0,
              0.986, 1.610, 0.0,
              1.855, 0.002, 0.0], requires_grad=True)

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

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

num_wires = 6

# the rest of the code is the same as "Optimization of molecular geometries" demo.

Hello Soran! Thanks so much for the suggestion. I did as you said, and the code in the tutorial, this time using the smaller-active space that you suggested. The code was able to converge with a OH1 bond-radius of 0.98958 Angstroms, and an OH1-OH2 bond angle of 100 degrees.

However, I’m curious for why this convergence can’t be replicated for the active-space size used in this paper: https://arxiv.org/pdf/2106.13840.pdf ? I’m wondering if you could try it on your end, or show me the results of the code that was used for the paper and let me know how it goes. Thanks so much

Hi @ImranNasrullah. Great that the geometry converges for the small active space.

Could you please increase the size of the activate space step-by-step, e.g., 2 electron and 1 orbital at a time, and see when the optimisation goes wrong?

Also, please change the initial coordinates, e.g., by rotating the molecule or slightly changing the bond lengths / angle.

x = np.array([0.028,  0.054,  0.10,
              0.986,  1.610, -0.10,
              1.855,  0.002,  0.20], requires_grad=True)

The initial geometry and possible molecular symmetries might affect the accuracy of the gradient calculations with the finite difference method. All of the calculations in the paper were done with PennyLane using a code similar to that in the demo. So, following the two hints provided above should fix the issue or at least help us to find where the problem is. Please try them and let me know the results.


Sounds good, I’ll give it a shot

1 Like

Hello Soran, apologies for the delay. As you suggested, I ran my code starting from 4 electrons and 4-orbitals, increasing the # of electrons by 2 and the orbitals by 1 every time. 4 e- and 4 m.o., as well as 6 e- and 5 m.o. seem to converge just fine. But when I get to 8-electrons and 6-m.o, here are the results for:

And finally, here is the associated energy when the algorithm takes these steps:

As you can see, both the bond-length and the bond-angle start to diverge pretty quickly. Furthermore, I’m confused why the energy is still around -73- (-74) Ha even with a drastically small angle. I remember this happening for the initial choice of water-coordinates x i chose a while back, and now I believe these are consistent results for your suggestion of the rotated-coordinates x you provided me above. So even with 2 initializations, the algorithm consistently diverges.

If you have any suggestions from where to go from here that’d be amazing. And would it be possible for anyone to look at my code?
Since I basically used the tutorial, it should be pretty readable

hello, actually would it have to do with the difference in calculation of the dhf method? Soran, I saw ur answer about how using the psi4 method would work, linked here:

Hi @ImranNasrullah, could you please use the ‘pyscf’ backend and let me know the results? It will help us to find the actual reason for the diverging coordinates.

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

Regarding the energy fluctuations, they are actually pretty large because 1 Ha = 627.5096 kcal/mol.

Please also attach your code here, so I can have a look. Thanks!

Hello Soran,
I tried using the ‘pyscf’ backend like you suggested, and the bond length becomes very large,
while the angle becomes super small. Not sure why this is the case. I’ve attached this code for you to take a look over.

Both files are pretty much the same thing, one has results with the default ‘dhf’ method, one is using the ‘pyscf’ method, but both are optimizing water with the same size of active space.

I remember you replied to another person having the same issue as me a while back ago, and their issue was fixed by using the ‘psi4’ library, which I believe is now deprecated. Is there any sort of backend like that as well, or can we only use the ‘dhf’ and ‘pyscf’ backends.

And thank you so much by the way!