Simulating more complex molecules in Pennylane using VQE

Hello everyone,

I am trying to simulate some molecules using the variational quantum eigensolver in Pennylane. I’ve read somewhere that you can get results in reasonable time for situations in which less than about 16 qubits are required. I saw numerous examples of simulations of H2O, LiH or H2 molecules, but say I want to use VQE for a molecule with more complex atoms, CuO for example. When I try to input the molecular geometry, I end up with an error. Is there some way to get around this?

Hey @mileva2319! Welcome to the forum :rocket:!

but say I want to use VQE for a molecule with more complex atoms, CuO for example. When I try to input the molecular geometry, I end up with an error.

I can only guess as to what the error is, but maybe it’s something like this:

import pennylane as qml
from pennylane import numpy as np

symbols = ["Cu", "O"]
coords = np.array([[0.0, 3.0, 0.0], [0.0, 0.0, 0.0]])

hamiltonian, qubits = qml.qchem.molecular_hamiltonian(symbols, coords, mult=2)

'''
ValueError: Atoms in {'Cu'} are not supported.
'''

As this error suggests, Cu isn’t a supported atom in our differentiable Hartree-Fock solver. However, you can use a different backend — the OpenFermion-PySCF backend — that supports atoms like Cu. You simply need to add a special keyword argument to specify as much:

hamiltonian, qubits = qml.qchem.molecular_hamiltonian(symbols, coords, mult=2, method='pyscf')
print(hamiltonian)

'''
<Hamiltonian: terms=172654, wires=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45]>
'''

You’ll need to pip install openfermionpyscf to get this to work :slight_smile:. Just note that the resulting Hamiltonian is non-differentiable — you can still perform VQE with it, but it will rely on numerical approaches for computing gradients, which can lead to inaccuracies and instability compared to the more robust differentiable Hartree-Fock solver.

Here’s an example with a smaller molecule:

symbols = ["H", "O"]
coords = np.array([[0.0, 1.0, 0.0], [0.0, 0.0, 0.0]])

hamiltonian, qubits = qml.qchem.molecular_hamiltonian(symbols, coords, charge=-1, method='pyscf')

dev = qml.device("default.qubit", wires=12)

@qml.qnode(dev)
def circuit(x):
    for i in range(12):
        qml.RX(x[i], wires=i)
    
    return qml.expval(hamiltonian)

x = np.random.uniform(0, 2*np.pi, size=(12,))
print(circuit(x))
print(qml.grad(circuit)(x))

'''
-54.046642484365606
[10.84198255 10.38715297  2.15646718  1.49026374  1.7440145  -1.07904077
  0.58801415 -1.1730143   1.62609713 -1.67494756  0.95318472  0.16550297]
'''

Hope this helps!