Hi @ImranNasrullah, Great that the code works. The issue could be with the large step in the optimisation. Could you please use a more reasonable geometry for water and print the geometry, the bond length and the bond angle at each step? That will help to see why the geometry diverges. Thanks!
Ok actually, I think the issue was this line:
#z_matrix = z_matrix + (0.5 * forces)
z_matrix = z_matrix + (1 * forces)
The values were fine, but the gradients were very slow:
Starting bond length = 0.9727221432914391, bond angle = 60.64705337449555
n: 0, E: -74.87727582
r gradient: 0.30221041, Bond-length: 1.12383
theta gradient: 0.00487345, Bond-angle: 60.64949
n: 1, E: -74.87611641
r gradient: 0.24895165, Bond-length: 0.99935
theta gradient: 0.00318285, Bond-angle: 60.65108
n: 2, E: -74.88441809
r gradient: 0.16882385, Bond-length: 1.08376
theta gradient: 0.00455928, Bond-angle: 60.65336
n: 3, E: -74.88448594
r gradient: 0.14527409, Bond-length: 1.01113
theta gradient: 0.00361623, Bond-angle: 60.65517
n: 4, E: -74.88625368
r gradient: 0.11589304, Bond-length: 1.06907
theta gradient: 0.00442514, Bond-angle: 60.65738
n: 5, E: -74.88638265
r gradient: 0.10065774, Bond-length: 1.01874
theta gradient: 0.00377800, Bond-angle: 60.65927
n: 6, E: -74.88705213
r gradient: 0.08337407, Bond-length: 1.06043
theta gradient: 0.00433874, Bond-angle: 60.66144
n: 7, E: -74.88716148
r gradient: 0.07262312, Bond-length: 1.02412
theta gradient: 0.00387317, Bond-angle: 60.66338
n: 8, E: -74.88746120
r gradient: 0.06118253, Bond-length: 1.05471
theta gradient: 0.00427770, Bond-angle: 60.66552
n: 9, E: -74.88754066
r gradient: 0.05332279, Bond-length: 1.02805
theta gradient: 0.00393601, Bond-angle: 60.66748
n: 10, E: -74.88768821
r gradient: 0.04533859, Bond-length: 1.05072
theta gradient: 0.00423296, Bond-angle: 60.66960
n: 11, E: -74.88774280
r gradient: 0.03950059, Bond-length: 1.03097
theta gradient: 0.00397973, Bond-angle: 60.67159
In an attempt to make the gradients be faster (as shown above, in 11 steps, the gradient is very small), I removed the factor of 0.5, and I think this causes the algorithm to diverge.
Would there be any better to way to make the gradients be larger? Or is it an issue in itself that the gradient of ‘r’ and ‘theta’ are smaller despite the values being not close to water’s ground state?
Hey @ImranNasrullah,
I think at this stage you’ve reached a point of needing to test different hyperparameters (e.g., learning rates, number of parameters in your circuit), optimization methods, parameter initialization strategies, etc. If you have any questions about PennyLane features, be they errors you encounter or features you’d like to understand better, please let us know!
Sure! In that case, I’m trying out the QNG Optimizer for VQE according to this tutorial. Can the QNGOptimizer be used to update the nuclear coordinates of the Hamiltonian (in my case the ‘r’ and ‘theta’) or can they only be used to update the parameters of the circuit to create the ansatz? I know thatthe QNGOptimizer has to do with searching along the wavefunction’s parameter space, which suggests the circuit parameters, so not sure if it inherently only works for just the circuit params.
I’m able to get the QNGOptimizer to work for the circuit parameters, but when I get it to update the nuclear coordinate parameters, I get this issue:
File ~/anaconda3/lib/python3.11/site-packages/pennylane/transforms/core/transform_program.py:324, in TransformProgram._set_all_classical_jacobians.<locals>.jacobian(classical_function, program, argnums, *args, **kwargs)
321 argnums = 0 if argnums is None else argnums
323 if not indices and argnums is None:
--> 324 raise qml.QuantumFunctionError("No trainable parameters.")
326 classical_function = partial(classical_function, program)
328 if qnode.interface == "autograd":
QuantumFunctionError: No trainable parameters.
Also another question. Is ther any way to try using the ‘PSI4’ package, or is Pennylane long since deprecated its usage?
Because using PSI4 fixed a lot of errors for this dicussion post (second to last post) so I figure it is worth a shot. If not PS4, are there recommendations on using any other methods/packages?
Could you attach your code example that tries to use QNGOptimizer?
Regarding PSI4, we have plans for showing people how to use it with pennylane, similar to what we did for pyscf in this demo. But nothing right now 
Sure, here is the code (and correct me if I’m wrong, but it would be pretty powerful to implement the QNGOptimizer for the Optimization of Molecular Geometries tutorial because its independent of geometry initializations.
Here is the self-contained code:
# Working and Self-Contained Example:
import pennylane as qml
from pennylane import qchem
from pennylane import numpy as np
import time
symbols = ["O", "H", "H"]
r = np.array([0.95769], requires_grad=True) # Z-matrix parameters:
angle = np.array([90.5739], requires_grad=True)
def H(x):
return qml.qchem.molecular_hamiltonian(symbols, x, charge=0, active_electrons = 4, active_orbitals = 4)[0]
def build_xyz(r, angle): # prepares an equilateral triagnle's coordinates of H3+!
r = r/bohr_angs
x = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
# O at (0,0,0)
x[0] = 0.0
x[1] = 0.0
x[2] = 0.0
# H1 at (r*cos(angle),r*sin(angle),0)
x[3] = r * np.cos(angle * np.pi/180)
x[4] = r * np.sin(angle * np.pi/180)
x[5] = 0.0
# H2 at (r, 0, 0)
x[6] = r
x[7] = 0.0
x[8] = 0.0
return x
hf = qml.qchem.hf_state(4, 8)
num_wires = 8
dev = qml.device("lightning.qubit", wires=num_wires)
@qml.qnode(dev, interface="autograd")
def circuit(params, obs, wires):
# prepares reference state
qml.BasisState(hf, wires=wires)
# One gate for simplicity
qml.DoubleExcitation(params[0], wires=[0,1,4,5])
# returns expectation value of the ansatz prepared from this quantum circuit:
return qml.expval(obs)
########### Costs: ###########
def cost(params, r, angle):
hamiltonian = H(build_xyz(r[0], angle[0]))
return circuit(params, obs=hamiltonian, wires=range(num_wires))
# Cost for QNG Optimizer:
def circuit1(params, obs, wires):
# prepares reference state
qml.BasisState(hf, wires=wires)
# One gate for simplicity
qml.DoubleExcitation(params[0], wires=[0,1,4,5])
# returns expectation value of the ansatz prepared from this quantum circuit:
return qml.expval(obs)
dev = qml.device("default.qubit", wires=num_wires)
@qml.qnode(dev, interface="autograd")
def cost_QNG(params): # Is exact same as "cost", but is decorated with QNode -> for use of the QNG Optimizer
hamiltonian = H(build_xyz(r[0], angle[0]))
return circuit1(params, obs=hamiltonian, wires=range(num_wires))
########### Bond length Finite Diffs ###########
def finite_diff_r(f, r, angle, delta=1e-07):
'''Compute the central-difference finite difference of a function'''
gradient = []
shift = 0.5 * delta
temp1 = build_xyz(r+shift, angle)
temp2 = build_xyz(r-shift, angle)
res = (f(temp1) - f(temp2)) * delta**-1 # dH/dx
gradient.append(res)
return gradient
def grad_r(params):
grad_h = finite_diff_r(H, r,angle)
grad = [circuit(params, obs=obs, wires=range(num_wires)) for obs in grad_h]
return np.array(grad)
############## Angles Finite Diffs ###################
def finite_diff_angle(f, r, angle, delta=1e-07):
'''Compute the central-difference finite difference of a function'''
gradient = []
shift = 0.5 * delta
temp1 = build_xyz(r, angle+shift)
temp2 = build_xyz(r, angle-shift)
res = (f(temp1) - f(temp2)) * delta**-1 # dH/dx
gradient.append(res)
return gradient
def grad_angle(params):
grad_h = finite_diff_angle(H, r, angle)
grad = [circuit(params, obs=obs, wires=range(num_wires)) for obs in grad_h]
return np.array(grad)
#############################################################################
# Optimizers
opt_theta = qml.QNGOptimizer(0.5, lam=0.001, approx="block-diag")
opt_r = qml.QNGOptimizer(0.5, lam=0.001, approx="block-diag")
opt_angle = qml.QNGOptimizer(0.5, lam=0.001, approx="block-diag")
# circuit parameters:
theta = np.array([0.0], requires_grad=True)
################# Optimization Loop: #################
from functools import partial
# store the values of the cost function
energy = []
# store the values of the bond length
bond_length = []
bond_angle = []
# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903
for n in range(100):
# Optimize the circuit parameters
theta.requires_grad = True
r.requires_grad = False
angle.requires_grad = False
theta = opt_theta.step(cost_QNG, theta)
print("Done with theta optimization!")
# Optimize the nuclear coordinate bond length
theta.requires_grad = False
r.requires_grad = True
angle.requires_grad = False
r = opt_r.step(cost_QNG, theta, grad_fn=grad_r)
# Optimize the nuclear coordinate bond angle
theta.requires_grad = False
r.requires_grad = False
angle.requires_grad = True
angle = opt_angle.step(cost_QNG, theta, grad_fn=grad_angle)
energy.append(cost_QNG(theta))
bond_length.append(r[0])
bond_angle.append(angle[0])
if n % 1 == 0:
print(f"Step = {n}, E = {energy[-1]:.8f} Ha")
print(f"Bond length = {bond_length[-1]:.5f} A, Bond angle = {bond_angle[-1]:.5f}" + '\u00b0')
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}")
Hey @ImranNasrullah,
The issue is because of this block:
theta.requires_grad = False
r.requires_grad = True
angle.requires_grad = False
r = opt_r.step(cost_QNG, theta, grad_fn=grad_r)
The QNode cost_QNG strictly depends on params, which is theta in the block above, and you wrote theta.requires_grad = False. This is why you see QuantumFunctionError: No trainable parameters.
If you put theta.requires_grad = True, everything runs without any errors.
Hello Isaac,
Can I be sure that setting theta.requires_grad = True would also cause the value that is returned in this context to be the new value of r? It just seems like it wouldn’t be the case because we set theta being optimized to True
Thanks
If you don’t want to differentiate w.r.t. theta, then your cost function needs to depend on something else that is differentiable. Would that be r?
I believe so, because the hamiltonian is a function of r. It would be nice to apply the principle of the QNGOptimizer to both r and theta so that they would be able to be optimized independent of the geometry initialization. That’s why I was looking to apply the QNGOptimizer to both r and theta
Your cost_QNG function doesn’t depend on r explicitly, which is why it’s not being differentiated with respect to it. I’m not sure if just adding the dependency fixes it, but worth a try!