Calculating ground energy using VQE for benzene

Hi,

I’ve tried the code for calculating ground energy for water molecule and it worked fine but when i’m trying more complex molecule such as benzene c6h6, it get stuck on building the hamiltonian . i have tried changing the parameters such as number of orbitals, basis set etc but no progress.

Hey @Dassen_Sathan! It’s probably that the Hamiltonian is so big that it just takes a while to generate the qubit representation of it. Indeed, if I try generating the qubit Hamiltonian of benzene on my laptop, it takes a while!

Do you need to run this calculation on your laptop? It might be better if you have access to a remote server that you can run this on.

Thanks Isaac, for the reply as i have docking molecules that is bigger than benzene. Can you please give me the steps to follow to build hamiltonian on a remote server.

Is there a way to speed up the process on laptop?
cheers

You might actually be able to get away with generating the molecular Hamiltonian on your laptop with the pyscf backend:

import pennylane as qml
from pennylane import qchem 

labels, coords = qchem.mol_data("C6H6")
benzene_hamiltonian = qchem.molecular_hamiltonian(labels, coords, method='pyscf')

As far as finding its ground state energy, though, you will definitely need to define an active space — the qubit Hamiltonian is 72 qubits :astonished:! In that case, I’d still recommend using the pyscf backend, but specify a reasonable number of active electrons:

benzene_hamiltonian = qchem.molecular_hamiltonian(labels, coords, method='pyscf', active_electrons=4)

That line above ran for me in under 2 minutes on my laptop :grin:. Hope this helps!

Hi,

Thanks a lot for your suggestions. What an improvement under 2 mins. Will definitely try that.

have a good day.

2 Likes

Hi isaac,

I’ve tried the pyscf, it is much quicker but getting the following error when trying with the docking molecule

raise RuntimeError(‘Electron number %d and spin %d are not consistent\n’
RuntimeError: Electron number 63 and spin 0 are not consistent
Note mol.spin = 2S = Nalpha - Nbeta, not 2S+1

Any help from you part will be greatly appreciated, have a nice week-end

cheers

Hey @Dassen_Sathan,

It sounds like you’re trying to simulate an unphysical electronic structure. The number of electrons and the corresponding spin you set aren’t physically possible. Can you share your code?

Hi isaac,

Please find the code and the structure i was testing it( 3htb and jz4 docking). Actually i’ve used siremol to calculate the binding energy(Free energy perturbation) and it worked fine though it took me 4 days to get the results which is normal according to siremol documents. I was thinking of speeding the process and thought about using VQE ( some papers mentioned about it but did not give the steps to do it). i don’t know if it is the right way to do it but I don’t have any result yet to be able to compare.

thanks for your help

(Attachment 3htbjz4_dock.xyz is missing)

FEP_VQE_pyscf.py (5.8 KB)

Hi isaac,

Please find the code and the structure i was testing it( 3htb and jz4 docking). Actually i’ve used siremol to calculate the binding energy(Free energy perturbation) and it worked fine though it took me 4 days to get the results which is normal according to siremol documents. I was thinking of speeding the process and thought about using VQE ( some papers mentioned about it but did not give the steps to do it). i don’t know if it is the right way to do it but I don’t have any results yet to be able to compare.

thanks for your help

3htbjz4_dock.txt (567 Bytes)

FEP_VQE_pyscf.py (5.8 KB)

The error here is stemming from this method. The number of alpha/beta electrons here means number of spin up/down electrons, and the error is saying that n_\alpha + n_\beta \neq 63.

Looking at your code, I can indeed replicate the error you’re getting. But, there are a few things that I’m confused about:

  • benzene_hamiltonian isn’t used beyond grabbing the number of qubits, but later on you define a different variable H with more electronic structure details :thinking:. Not sure which one you want to use, but I assume it’s H because that’s what your VQE is optimizing against.

  • symbols isn’t defined, but I assume you meant labels.

  • If I get rid of benzene_hamiltonian and move the electronic structure variable definitions above H (e.g., charge, multiplicity, etc.) so that it can be defined, I get a new error: ValueError: Openshell systems are not supported. Change the charge or spin multiplicity of the molecule. This means that there are unpaired electrons, which, as the error suggests, isn’t supported. I tried several multiplicity values, assuming that you want charge to remain 0, and none of them worked.

In any case, I’m going to get in touch with one of our qchem experts to see if they can help here :slight_smile:. Hang tight!

Hi,

Thanks a lot for your feedback, i was testing it on benzene, when it worked i tried it on the 3htb_jz4 without modifying the code, that’s why the code is a bit messy. Hope your qchem colleagues can help out, i’ve been stuck for some time on these issues.

cheers

Hi @Dassen_Sathan and thanks for using PennyLane!

The error you are getting is because of an inconsistency in the number of electrons and spin multiplicity. The molecule defined in the 3htbjz4_dock.txt file has 63 electrons (9 * 6 + 8 + 1, for C, O, and H, respectively) but your spin multiplicity is the default value of 1 which corresponds to a system with no spin-unpaird electron. To fix the error, please make sure that you construct your molecular Hamiltonian with the correct multiplicity of 2 to account for an odd number of electrons in your molecule:

H, qubits = qchem.molecular_hamiltonian(
    labels,
    coords,
    charge=0,
    mult=2,
    basis='sto-3g',
    method='pyscf',
    active_electrons=5,
    active_orbitals=4,
)

Please see here for more details on multiplicity. Please also note that I have changed the number of active_electrons from 4 to 5 to make sure that the core electrons are not unpaired.

Also, your molecule seems to be missing some hydrogen atoms as I think what you aim to simulate is 2-Propylphenol. If this is the case, please make sure that you are using the correct charge and multiplicity after adding the missing hydrogen atoms. Finlay, if you need to construct your Hamiltonian few times in the code, please make sure that you use the correct input parameters everywhere. Feel free to let us know if you have any other questions!

Hi,

thanks a lot for your explanations, i will try what you have suggested and hopefully it will work .

Can you please let me know how to add the missing hydrogen atoms, is it using pymol??

cheers

cheers

Hi @Dassen_Sathan. You can use the mol_data function to obtain the symbols and geometry of a compound from the PubChem database. For instance, the symbols and geometry of 2-Propylphenol can be obtained with

qml.qchem.mol_data("2-Propylphenol")

Please inspect the symbols and the geometry properly to make sure that the molecule is the correct one and let us know if you have any other questions.

Hi,

thanks a lot for your help Will definitely contact if i have any other question.

1 Like

Hi,

I’ve tried what you have suggested using the 2- Propylphenol, launched it on friday and still no result but still in running state, i don’t know if there is an issue in building the hamiltonian.

cheers

Hi @Dassen_Sathan. The following code should work and give you the Hamiltonian in a few seconds, assuming that the desired molecule is 2-Propylphenol.

import pennylane as qml
from pennylane import numpy as np

labels, coords = qml.qchem.mol_data("2-Propylphenol")

H, qubits = qml.qchem.molecular_hamiltonian(
    labels,
    coords,
    charge=0,
    mult=1,
    basis='sto-3g',
    method='pyscf',
    active_electrons=4,
    active_orbitals=4,
)

Please make sure that you build the Hamiltonian correctly everywhere. You might also change the number of active_electrons and active_orbitals depending on the problem.

1 Like

hi,

Ok it seems i made a mistake somewhere, will check my code gain.Thanks a lot! have a good day

1 Like

Let us know if you still need help!

Hi,

I’ve made the changes as suggested and managed to get the energy for a docking of 3 htb and 2 propylphenol but the result i’m getting is -400 for the energy using VQE compared to classical methods with is -29 (siremol tool). Any suggestion what might have gone wrong…