Pennylane takes too long to build a Hamiltonian

Dear Pennylane team,

I’ve buid a number of molecular Hamiltonians and quantum circuits by using the calls
qchem.molecular_hamiltonian() and ApproxTimeEvolution().

Most molecules I’ve tried are toy examples and featured in the Pennylane tutorials or documentation pages. I’ve tried: LiH, H2O and H2.

For these small molecules Pennylane computes Hamiltonians and respective quantum circuits in few minutes.

However, recently I’ve tried to get a Hamiltonian and a quantum circuit for
the C2HF3 molecule. The geometry is as follows:

    symbols = ["C", "C", "F", "F", "F", "H"]
    coordinates = np.array(
        [
            0.0, 0.430631, 0.000000,
            -0.704353, -0.686847, 0.0,
            1.314791, 0.500587, 0.0,
            -0.554822, 1.630861, 0.0,
            -0.092517, -1.882156, 0.0,
            -1.780948, -0.706341, 0.0
        ]
    )
hamiltonian, num_qubits = qchem.molecular_hamiltonian(symbols, coordinates)
print(f"Number of qubits: {num_qubits}")
print("Qubit Hamiltonian")
print(hamiltonian)

This is a larger molecule and Pennylane is not able to compute a Hamiltonian let alone the circuit. That is, calling qchem.molecular_hamiltonian() runs for hours and never terminates.
Any ideas on how to get a Hamiltonian and a circuit in reasonable time (say few minutes)? What do I miss here?

P.S. I tried the same molecule in a different quantum chem python package and it was able to produce a Hamiltonian and a respective circuit in under 2 minutes.

Hi @Einar_Gabbassov,

Thank you for sharing this information with us. The molecule you’re looking at is quite large so my guess is that any other packages you’re using are doing approximations or using qubit-reduction techniques under the hood.

Something you could try is doing qubit tapering. This can reduce your computational needs and speed up the simulation. You can learn more about this in this demo.

Please let me know if you have any questions about this and let me know if this gives you a speedup!

Hi @Einar_Gabbassov,

My colleague Soran suggested that you will get much better speed if you use method="pyscf" when you build the molecular Hamiltonian. So you would do something like:

H, qubits = qchem.molecular_hamiltonian(symbols, coordinates, method="pyscf")

You can learn more about this in this demo. I hope this is helpful!

Hi @CatalinaAlbornoz
Thanks for the suggestion. I have tried method=pyscf. Unfortunattely, it crashes my jupyter notebook and on top of that it produces an error.

On a positive side… it takes about 5 minutes to crash my jupyter notebook. So if the bug is fixed, it might take about 5 minutes to compute the Hamiltonian?? :man_shrugging:

Below is a partial traceback. Note that even though all coefficients of a Hamiltonian are real, the datatype seems to be a complex number.

File ~/.pyenv/versions/circuit_comp_env/lib/python3.8/site-packages/pennylane/qchem/openfermion_obs.py:956, in molecular_hamiltonian(symbols, coordinates, name, charge, mult, basis, method, active_electrons, active_orbitals, mapping, outpath, wires, alpha, coeff, args, grouping_type, grouping_method)
    950 core, active = qml.qchem.active_space(
    951     molecule.n_electrons, molecule.n_orbitals, mult, active_electrons, active_orbitals
    952 )
    954 h_of, qubits = (decompose(hf_file, mapping, core, active), 2 * len(active))
--> 956 h_pl = qml.qchem.convert.import_operator(h_of, wires=wires)
    958 return (
    959     qml.Hamiltonian(h_pl.coeffs, h_pl.ops, grouping_type=grouping_type, method=grouping_method),
    960     qubits,
    961 )

File ~/.pyenv/versions/circuit_comp_env/lib/python3.8/site-packages/pennylane/qchem/convert.py:325, in import_operator(qubit_observable, format, wires, tol)
    320     raise TypeError(f"Converter does not exist for {format} format.")
    322 if any(
    323     np.iscomplex(np.real_if_close(coef, tol=tol)) for coef in qubit_observable.terms.values()
    324 ):
--> 325     raise TypeError(
    326         f"The coefficients entering the QubitOperator must be real;"
    327         f" got complex coefficients in the operator {qubit_observable}"
    328     )
    330 return qml.Hamiltonian(*_openfermion_to_pennylane(qubit_observable, wires=wires))

TypeError: The coefficients entering the QubitOperator must be real; got complex coefficients in the operator (-169.0923877964394+0j) [] +
(-3.41902239433203e-05+0j) [X0 X1 Y2 Y3] +
(-2.996469373464669e-06+0j) [X0 X1 Y2 Z3 Z4 Y5] +  
...  
...  

Hi @Einar_Gabbassov, thank you for posting this here. Is this the full error traceback?

Also, can you please post the output of qml.about()?

Hi @Einar_Gabbassov, I managed to reproduce your problem. I will look into it and come back with an answer.

Hi @Einar_Gabbassov and thanks for using PennyLane to simulate this interesting molecule. The molecule is relatively large and it is more efficient to use the non-differentiable workflow of PennyLane, which uses OpenFermion-PySCF, to construct the molecular Hamiltonian. Please see here and here for more details.

The error you get is because of small imaginary components in the Hamiltonian coefficients. The default threshold for raising this error is defined here. We are working on making this default threshold slightly larger. In the meantime, you can directly construct your Hamiltonian using the following code. Please note that the atomic coordinates should be in atomic units and a conversion factor must be applied to the original values.

import pennylane as qml
from pennylane import numpy as np
import openfermion
from pennylane.qchem.openfermion_obs import meanfield, decompose


symbols = ["C", "C", "F", "F", "F", "H"]
coordinates = np.array(
        [
            0.0, 0.430631, 0.000000,
            -0.704353, -0.686847, 0.0,
            1.314791, 0.500587, 0.0,
            -0.554822, 1.630861, 0.0,
            -0.092517, -1.882156, 0.0,
            -1.780948, -0.706341, 0.0
        ]
    ) / 0.529177210903  # convert to Bohr


hf_file = meanfield(symbols, coordinates)
molecule = openfermion.MolecularData(filename=hf_file)
h_of = decompose(hf_file)
h_pl = qml.qchem.convert.import_operator(h_of, tol=1e09)  # use a larger threshold for the imaginary part

print(h_pl)

Please feel free to let us know if you have any further questions.

1 Like