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 
!
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 
. 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!
Hi, I am gettting similar problems to simulate FeH bonding. I am quite new in these topics so I am wondering if I am getting the number of electrons and orbitals correctly. To this I am mainly addresing my code from Optimization of molecular geometries | PennyLane Demos. For charge, I am assuming that this is neutral but with a multiplicity of 3 because the number of unpaired electrons is 2 (following the n +1 rule). I am awared that my error it might not related with these parameters, but I think that they could be important. Also I already include the method as “openfermion”
Thank you in advance for your help
The error is,
packages/pennylane/qchem/molecule.py:104) self.coordinates = coordinates ValueError: Atoms in {'Fe'} are not supported.
The code is,
from pennylane import numpy as np
import pennylane as qml
symbols = ["Fe", "H"]
x = np.array([0.028, 0.054, 0.0, 0.986, 1.610, 0.0], requires_grad=True)
elec = 26 + 1
orbs = (15 + 1)*2
delta_sz = 0
q = 1
m = 1
def H(x):
molecule = qml.qchem.Molecule(symbols, x, charge=q,mult=m)
global qubits
hh, qubits = qml.qchem.molecular_hamiltonian(molecule,method='openfermion' )
# print(qubits)
return hh
hf = qml.qchem.hf_state(electrons=elec, orbitals=orbs)
print(hf)
singlet, doublet = qml.qchem.excitations(electrons = elec , orbitals = orbs, delta_sz = delta_sz)
dev = qml.device("lightning.qubit", wires=qubits)
@qml.qnode(dev, interface="autograd")
def circuit(params, obs, wires):
qml.BasisState(hf, wires=wires)
# [qml.SingleExcitation(phi = params[i] , wires = s) for i, s in enumerate(singles)];
[qml.DoubleExcitation(phi = params[j] , wires = d) for j, d in enumerate(doublet)];
return qml.expval(obs)
def cost(params, x):
hamiltonian = H(x)
# print(hamiltonian)
return circuit(params, obs=hamiltonian, wires=range(qubits))
def finite_diff(f, x, delta=0.01):
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(qubits)) for obs in grad_h]
return np.array(grad)
print(len(doublet))
opt_theta = qml.GradientDescentOptimizer(stepsize=0.4)
opt_x = qml.GradientDescentOptimizer(stepsize=0.8)
theta = np.array(np.zeros(len(doublet)), requires_grad=True)
from functools import partial
energy = [ ]
bond_length = [ ]
bohr_angs = 0.529177210903
for n in range(100):
theta.requires_grad = True
x.requires_grad = False
theta, _ = opt_theta.step(cost, theta, x)
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")
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}")
'''Hi @Raul_Guerrero ,
Currently qml.qchem.Molecule only supports atoms with atomic numbers between 1-10.
The molecule you’re trying to simulate is very large and your laptop will most likely crash before you can simulate it.
However, you can explore the question of how many resources you would need to simulate this. The last section of this demo on resource estimation addresses the question of estimating the cost (in terms of gates and qubits) of implementing the QPE algorithm for simulating dilithium iron silicate (Li2FeSiO4) using Hamiltonians represented in first quantization and in a plane wave basis.
I recommend that you read the whole demo to see how even a water molecule can be very complex.
I hope this helps!
Hi @CatalinaAlbornoz,
thank you for your quick answer. I was expecting something like that. I thougth that I can simulate “Fe” because in PennyLane’s Q&A I read about cases like “Cu” in using openfermion as method. Simulating more complex molecules in Pennylane using VQE - PennyLane Help - Discussion Forum — PennyLane. So I understand that “Fe” is imposible to simulate with PennyLane, even with “openfermion” method?
Thank you in advance,
Raúl
Hi @Raul_Guerrero ,
I think there may be other options of simulating it, just not using qml.qchem.Molecule but something else instead.
Let me check what could work here and get back to you.
I have some great news! We have this Pull Request under review, which updates Molecule to accept all atom types. Once it’s merged it will be possible to simulate all elements, including Iron. Our next stable release is planned for early September, so once this pull request is merged you can use the development version of PennyLane for a few weeks and then use version v0.38.0 of PennyLane when it’s released in September.
# install the latest development version of PennyLane
pip install git+https://github.com/PennyLaneAI/pennylane.git#egg=pennylane
Let me know if you have any questions about this!
The Pull Request has been merged! Please let me know if you face any issues.
Wow that is really good news I will try to implement it as asap and let you know through this channel.
Thank you again for your help!
No problem @Raul_Guerrero
Just remember that you’re trying to simulate a very large molecule so I recommend restricting your active space (a lot) to avoid running out of memory.
Let me know how it goes!
Yes, already consider it jejej. 10^{eternity} thank yous
My code is working by now for Fe (and very reduced active space) 
. Let’s see how converge 
That’s great to hear @Raul_Guerrero ! 