Calculating number of active electrons and active orbitals for a given molecule

PennyLane Help
1 
symbols, coordinates = (['O','N','C','H','H','H','C'], 
                    np.array([1.1280, 0.2091, 0.0000, 
                              -1.1878, 0.1791, 0.0000, 
                              0.0598, -0.3882, 0.0000, 
                              -1.3085, 1.1864, 0.0001, 
                              -2.0305, -0.3861, -0.0001, 
                              -0.0014, -1.4883, -0.0001, 
                              -0.1805, 1.3955, 0.0000]))
h, qubits = qchem.molecular_hamiltonian(
    symbols,
    coordinates,
    name="hiv",
    charge=0,
    mult=1,
    basis='sto-3g',
    active_electrons=2,
    active_orbitals=2
    )
print(h,qubits)

electrons = 2
S2 = qchem.spin2(electrons, qubits)
print(S2)

I have the above molecule with 7 atoms (1O, 1N, 3H, 2C). I am attempting to find its ground state energy by implementing the code provided as one of Amazon Braket’s examples for Pennylane (https://github.com/aws/amazon-braket-examples/blob/main/examples/pennylane/3_Quantum_chemistry_with_VQE.ipynb).

Is there a right way to calculate the number of active_electrons, active_orbitals, and electrons? If yes, in case of the molecule I am using, what values should I use in order to get the ground state energy that is closer to the experimental value?

2

Hi @akb,

Thank you very much for sharing your question and welcome to the Xanadu forum!

Is there a right way to calculate the number of active_electrons, active_orbitals, and electrons?

Not really. Let me start by clarifying that defining an active space is an approximation to keep the simulation within the computational resources at hand. In principle, one has to include all the electrons in the molecule. On the other hand, the total number of molecular orbitals M determines the number of qubits N_\mathrm{qubits} = 2 \times M (see this tutorial for more details) .

For the “hiv” molecule you have:

  • Total number of electrons N_e=30

  • For a minimal basis set (e.g., “sto-3g”) M=23

  • This requires 46 qubits to simulate the molecular Hamiltonian. This is hard due to the size and depth of the variational circuit, and the number of terms in the decomposed Hamiltonian.


A natural approximation is to neglect the core orbitals (orbitals with deep energies that are typically not important for chemical bonding). In your case there are 4 core orbitals populated by 8 core electrons. So you would end up with a system of 22 electrons and 19 valence orbitals which still requires significant resources (38 qubits). In order to reduce further the size of your simulation you can decrease the number of active electrons and orbitals. Currently, PennyLane simulators should be able to handle up to ~ 16 qubits. This corresponds to a symmetric active space with active_electrons = 8 and active_orbitals = 8. Note that if you have a chemical intuition about which orbitals you want to include in the active space you can use the decompose function (see here) that allows you to choose the active orbitals.

what values should I use in order to get the ground state energy that is closer to the experimental value?

Increasing the active space is crucial to improve the accuracy of the ground-state energy. Having said that, I would like to stress that a minimal basis set provides a poor representation of the molecular orbitals. Note that you can improve your basis set while keeping the size of your active space fixed.

You can also optimize your simulation by defining a quantum circuit that prepares a molecular state as close as possible to the true ground state of the system. For example, the most important contributions to the ground-state correlation energy can be captured by applying DoubleExcitation operations to the M-qubit system initialized to encode the Hartree-Fock state of the molecule. I encourage you to take a look at the tutorial Adaptive circuits for quantum chemistry showing how you can reduce the number of quantum gates in your circuit and improve the performance of your simulation.

Hope this helps.

Do not hesitate to get back to us if you have further questions.

Thank you!

PS: Please, make sure that the coordinates you are using to build the molecular Hamiltonian are in atomic units (Bohr radii). Based on the interatomic distances, they seem to be in Angstroms.

3 Likes
3

Thank You Alain for the detailed explanation,
also please advise where can one read/learn in more detail about how it works the atoms, orbitals, i.e more in depth intricate details…

4

Hi @pratjz,

For these topics there are many books/papers. I would recommend the book by Prof. Frank Jensen “Introduction to Computational Chemistry”.

If you have questions, feel free to get back to us.

Best regards,

Alain.

1 Like
5

Thank You Alain for the suggestion

6

HI @Alain_Delgado_Gran

One more Question , is there a way I can find the ground state energy classically in pennylane ? for comparison b/w quantum & classical

Thanks
Prateek

7

Hi @pratjz, you can use sparse_hamiltonian to build a sparse matrix which represents the hamiltonian, and then the smaller eigenvalue will be the ground state energy.

This is an example of how you could do it:

import pennylane as qml
from pennylane import numpy as np
from pennylane import qchem

symbols = ["Li", "H"]
geometry = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 2.969280527])

H, qubits = qchem.molecular_hamiltonian(
    symbols,
    geometry,
    active_electrons=2,
    active_orbitals=5
)

H_sparse = qml.utils.sparse_hamiltonian(H)
eigenvalues, eigenvectors = np.linalg.eig(sparse_matrix)
ground_state_energy=min(np.real(eigenvalues))
print('ground_state_energy: ',ground_state_energy)

Let me know if this helps!

2 Likes
8

Thank You very much @CatalinaAlbornoz this looks like what I was looking for …will give it a try

1 Like
9

HI @CatalinaAlbornoz I was trying the suggested approach but it always gives the following error , Please help what am I missing

LinAlgError: 0-dimensional array given. Array must be at least two-dimensional

10

Please Ignore found the solution , Need to explicitly convert the sparse_matrix to an array like

eigenvalues, eigenvectors = np.linalg.eig(sparse_matrix.toarray())

Thank You :slight_smile:

11

I’m glad you found a solution @pratjz!