Help in tutorial optimization of molecules

PennyLane Help
1 
from pennylane import numpy as np

symbols = ["Be", "Be"]
x = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.5], requires_grad=True)
import pennylane as qml


def H(x):
    return qml.qchem.molecular_hamiltonian(symbols, x, charge=0)[0]
hf = qml.qchem.hf_state(electrons=2, orbitals=8)
print(hf)
num_wires = 8
dev = qml.device("lightning.qubit", wires=num_wires)


@qml.qnode(dev, interface="autograd")
def circuit(params, obs, wires):
    qml.BasisState(hf, wires=wires)
    qml.DoubleExcitation(params[0], wires=[0, 1, 2, 3])
    qml.DoubleExcitation(params[1], wires=[0, 1, 4, 5])
    qml.DoubleExcitation(params[2], wires=[2, 3, 6, 7])

    return qml.expval(obs)
def cost(params, x):
    hamiltonian = H(x)
    return circuit(params, obs=hamiltonian, wires=range(num_wires))
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 grad_x(params, x):
    grad_h = finite_diff(H, x)
    grad = [circuit(params, obs=obs, wires=range(num_wires)) for obs in grad_h]
    return np.array(grad)
opt_theta = qml.GradientDescentOptimizer(stepsize=0.4)
opt_x = qml.GradientDescentOptimizer(stepsize=0.8)
theta = np.array([0.0, 0.0], requires_grad=True)
from functools import partial

# store the values of the cost function
energy = []

# store the values of the bond length
bond_length = []

# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903

for n in range(100):

    # Optimize the circuit parameters
    theta.requires_grad = True
    x.requires_grad = False
    theta, _ = opt_theta.step(cost, theta, x)

    # Optimize the nuclear coordinates
    x.requires_grad = True
    theta.requires_grad = False
    _, x = opt_x.step(cost, theta, x, grad_fn=grad_x)

    energy.append(cost(theta, x))
    bond_length.append(np.linalg.norm(x[0:3] - x[3:6]) * bohr_angs)

    if n % 4 == 0:
        print(f"Step = {n},  E = {energy[-1]:.8f} Ha,  bond length = {bond_length[-1]:.5f} A")

    # Check maximum component of the nuclear gradient
    if np.max(grad_x(theta, x)) <= 1e-05:
        break

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}")

this code is showing error

IndexError: index 2 is out of bounds for axis 0 with size 2

Can you explain if I use this tutorial code for other molecules , how can I change the values ?

2

Hi @Kaushik,

What the error is saying is that you’re trying to access params[2] but params only has two values (params[0] and params[1]) so it can’t access params[2]. In your code, params is given by theta, which you will notice has only two values. You need to add a third value to theta in order to run your code.

I’m noticing that you have a very big molecule. If you look at the qubits needed to run the Hamiltonian you’ll notice that you actually need 20 qubits instead of 8. This is in the limit of what you can do with a laptop. You can try running it but if you run out of memory you may need to choose a smaller molecule or try different techniques to lower the number of qubits needed.

ham, qub = qml.qchem.molecular_hamiltonian(symbols, x, charge=0)
print(qub)

I hope this helps.

3

How can I use active space approximation and freeze core approximation to get optimized geometry for big molecules like C2,C3,C4,C5 etc. Actually I want ground state energy value and excied state value for various molecules using ssvqe algorithm . I am using molecules like Be2, Be3, Be4 etc for various basis sets. For this I need optimized geometry and then energy values

4

Hi @Kaushik ,

I recommend that you go to the docs for qchem.molecular_hamiltonian.

In the SSVQE code which I shared here you will notice that we call qchem.molecular_hamiltonian and we set a number of active electrons and active orbitals.

Let’s imagine you have a molecule with 10 electrons. Most likely the first 6 electrons will stay in the lowest-energy orbitals, so you can set your active electrons as 4. You can then decide how many active orbitals you want to have. It makes sense to have a maximum of 4 active orbitals because most likely the electrons will need a lot of extra energy to occupy higher-energy orbitals, so it’s unlikely to be a natural state for that molecule.

Conclusion, instead of simulating all configurations of all 10 electrons in a ton of orbitals you’re restricting your problem to just a set of combinations that will most likely give you the answer you want.

For some very big molecules however you may still need too many qubits to simulate active electrons in the active orbitals, and a further reduction in the size of the simulation may compromise accuracy. There’s a reason why we’re excited about building actual quantum hardware! We hope this will allow us to tackle these complicated simulations in the future.

1 Like