Thanks @Raul_Guerrero, I created a simple working example that you can use to debug the original code. Please also have a look at this discussion which is a bout a similar topic. Please let me know if you have any further questions. Thanks!
from autograd import grad
import pennylane as qml
from pennylane import numpy as np
symbols = ["O", "H", "H"]
geometry = np.array([[0.028, 0.054, 0.10],
[0.986, 1.610, -0.10],
[1.855, 0.002, 0.20]], requires_grad=True)
mol = qml.qchem.Molecule(symbols, geometry)
dev = qml.device("default.qubit")
def energy(mol):
@qml.qnode(dev, interface="autograd")
def circuit(*args):
# note that active_electrons=2, active_orbitals=2 in this example
qml.BasisState(np.array([1, 1, 0, 0]), wires=range(4))
# note that this simple circuit has only one gate
qml.DoubleExcitation(*args[0][0], wires=[0, 1, 2, 3])
# note that active_electrons=2, active_orbitals=2 in this example
H = qml.qchem.molecular_hamiltonian(mol.symbols, mol.coordinates,
active_electrons=2, active_orbitals=2, args=args[1:])[0]
return qml.expval(H)
return circuit
# number of zeros should match the number of circuit gates
circuit_param = [np.array([0.0], requires_grad=True)]
for n in range(36):
args = [circuit_param, geometry]
mol = qml.qchem.Molecule(symbols, geometry)
g_param = grad(energy(mol), argnum = 0)(*args)
circuit_param = circuit_param - 0.25 * g_param[0]
forces = -grad(energy(mol), argnum = 1)(*args)
geometry = geometry + 0.5 * forces
print(f'n: {n}, E: {energy(mol)(*args):.8f}, Force-max: {abs(forces).max():.8f}')