Hello @isaacdevlugt ,
Thanks for the response!
I was able to get a non-zero output, and here is the self-contained code:
# Minimally working and self-contained example:
# imports and functions:
import pennylane as qml
from pennylane import qchem
from pennylane import numpy as np
from autograd import grad
import time
import math
bohr_angs = 0.529177210903
#################### Functions: ####################
def unit_vector(vector):
""" Returns the unit vector of the vector. """
return vector / np.linalg.norm(vector)
def angle_between(v1, v2):
""" Returns the angle in radians between vectors 'v1' and 'v2'::
"""
v1_u = unit_vector(v1)
v2_u = unit_vector(v2)
return (np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0)) * 180/math.pi) #convert to degrees
def build_geometry(r, angle): # prepares water's xyz coordinates given a bond length and angle
r = r/bohr_angs # converts from units of angstroms back to units of Bohrs
x = np.array([[0.0, 0.0, 0.0],
[0.0, 0.0, 0.0],
[0.0, 0.0, 0.0]], requires_grad=True)
# O1's nuclear coordinates, is positioned at x,y,z: (0,0,0)
# H1's nuclear coordinates, is positioned at x,y,z: (r,0,0)
# x:
zeros = np.zeros((x.shape)) # Doing it this way works with 'autograd':
zeros[1][0] = 1
x = x + (zeros * r)
# H2's nuclear coordinates, is position at x,y,z: (r*cos(angle), r*sin(angle), 0)
# x:
zeros = np.zeros((x.shape)) # Doing it this way works with 'autograd':
zeros[2][0] = 1
x = x + (zeros * r * np.cos(angle * np.pi/180))
# y:
zeros = np.zeros((x.shape)) # Doing it this way works with 'autograd':
zeros[2][1] = 1
x = x + (zeros * r * np.sin(angle * np.pi/180))
return x
dev = qml.device("default.qubit")
def energy(mol):
@qml.qnode(dev, interface="autograd")
def circuit(*args):
# re-define geometry according via our parametrizations, r and theta:
geometry = build_geometry(args[1][0][0], args[1][1][0]) # re-initialize geometry:
# note that active_electrons=8, active_orbitals=6 in this example
qml.BasisState([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0], wires=range(12))
# Just apply one gate for simplicity:
qml.DoubleExcitation(*args[0][0], wires=[0, 1, 10, 11])
H = qml.qchem.molecular_hamiltonian(mol.symbols, mol.coordinates,
active_electrons=8, active_orbitals=6, args=[geometry])[0]
return qml.expval(H)
return circuit
############################################################
##################### Initializations: #####################
symbols = ["O", "H", "H"]
# construct the z_matrix vector (which will store 'r' and 'theta') according to the geometry below:
geometry = np.array([[0.10000000, 0.10000000, 0.28377432],
[0.10000000, 1.45278171, -1.00662237],
[0.10000000, -1.45278171, -1.00662237]], requires_grad=True)
# Initialize starting r and starting theta
OH1 = geometry[0] - geometry[1]
OH2 = geometry[0] - geometry[2]
r = np.linalg.norm(geometry[0] - geometry[1]) * bohr_angs
theta = angle_between(OH1, OH2)
# Finally, initialize z_matrix, which we will use to construct the "args" vector for the first time:
z_matrix = [np.array([r], requires_grad=True)] + [np.array([theta], requires_grad=True)]
##################### Begin optimization: #######################################
circuit_param = [np.array([0.0], requires_grad=True)] * 1
for n in range(10): # Just run 1 iteration
args = [circuit_param, z_matrix]
geometry = build_geometry(args[1][0][0], args[1][1][0]) # re-initialize geometry:
mol = qml.qchem.Molecule(symbols, geometry)
# Update circuit parameters:
g_param = grad(energy(mol), argnum = 0)(*args)
circuit_param = circuit_param - 0.25 * g_param[0]
# Update z_matrix:
forces = -np.array(grad(energy(mol), argnum = 1)(*args))
z_matrix = z_matrix + (0.5 * forces)
print(f'n: {n}, E: {energy(mol)(*args):.8f}')
print(f'r gradient: {abs(forces[0]).max():.8f}, Bond-length: {z_matrix[0][0]:.5f}')
print(f'theta gradient: {abs(forces[1]).max():.8f}, Bond-angle: {z_matrix[1][0]:.5f}')
print("")
I was able to get the optimization in terms of r and theta. However, I think something is slightly off with the code as the energy doesn’t consistently decrease across the steps (sometimes it goes up slightly before heading back down)
. Implementing this into PennyLane as a proper differentiation method is something that would take a good chunk of time. But, you can most certainly define the function like you have outside of PennyLane and use it as a differentiation method.
.
) and I got: