Hey @ImranNasrullah,

It’s difficult to say what’s going on in your code (there’s a lot going on ). 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
gradient.append(res)
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"],
coordinates,
charge=0,
method="dhf",
)[0]
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)
@qml.qnode(dev)
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
optimize_coordinates()
```

Hope this helps!