Hamiltonian generation with active space

Hi !

I am working on creating a Hamiltonian using molecular integrals. I checked out this demo ( Using PennyLane with PySCF and OpenFermion | PennyLane Demos ). I am working on larger molecule, I need to introduce active space component while calculating the 1, 2 body integrals. Can you please help me in this regard ?

My next question is how to load the basis set (eg. ANO-RCC basis set) for FeS. I checked out the qchem.load_basisset(basis, element). Say for FeS, there are multiple orbitals, coefficients and exponents. How to load them into coefficients and alpha to create the Hamiltonian ?

import pennylane as qml
from pennylane import numpy as np

symbols = ["Fe", "S"]
geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], requires_grad=True)

# Load basis sets for both atoms
basis_Fe = qml.qchem.load_basisset('ANO-RCC-MB', 'Fe')
basis_S  = qml.qchem.load_basisset('ANO-RCC-MB', 'S')

# Concatenate all basis function exponents and coefficients for all atoms
exponents    = basis_Fe["exponents"]    + basis_S["exponents"]    # list concatenation, not array!
coefficients = basis_Fe["coefficients"] + basis_S["coefficients"]

# Now convert to object arrays (optional for `molecular_hamiltonian`; lists-of-lists is OK)
alpha = exponents
coeff = coefficients

molecule= qml.qchem.Molecule(symbols, geometry, alpha=alpha, coeff=coeff, load_data = True)


My error message:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[14], line 19
     16 alpha = exponents
     17 coeff = coefficients
---> 19 molecule= qml.qchem.Molecule(
     20     symbols, geometry, alpha=alpha, coeff=coeff, load_data = True
     21 )
     22 print(hamiltonian)

File ~/miniconda3/envs/na/lib/python3.11/site-packages/pennylane/qchem/molecule.py:181, in Molecule.__init__(self, symbols, coordinates, charge, mult, basis_name, name, load_data, l, alpha, coeff, normalize, unit)
    178 self.coeff = coeff
    179 self.r = r
--> 181 self.basis_set = [
    182     BasisFunction(self.l[i], self.alpha[i], self.coeff[i], self.r[i]) for i in range(len(l))
    183 ]
    184 self.n_orbitals = len(self.l)
    186 self.mo_coefficients = None

File ~/miniconda3/envs/na/lib/python3.11/site-packages/pennylane/qchem/molecule.py:182, in <listcomp>(.0)
    178 self.coeff = coeff
    179 self.r = r
    181 self.basis_set = [
--> 182     BasisFunction(self.l[i], self.alpha[i], self.coeff[i], self.r[i]) for i in range(len(l))
    183 ]
    184 self.n_orbitals = len(self.l)
    186 self.mo_coefficients = None

IndexError: list index out of range
Name: pennylane
Version: 0.42.1
Summary: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Home-page: 
Author: 
Author-email: 
License-Expression: Apache-2.0
Location: /Users/s.poyyapakkam/miniconda3/envs/na/lib/python3.11/site-packages
Requires: appdirs, autograd, autoray, cachetools, diastatic-malt, networkx, numpy, packaging, pennylane-lightning, requests, rustworkx, scipy, tomlkit, typing_extensions
Required-by: PennyLane-qiskit, pennylane_lightning

Platform info:           macOS-15.4-arm64-arm-64bit
Python version:          3.11.8
Numpy version:           2.3.2
Scipy version:           1.16.1
Installed devices:
- qiskit.aer (PennyLane-qiskit-0.42.0)
- qiskit.basicaer (PennyLane-qiskit-0.42.0)
- qiskit.basicsim (PennyLane-qiskit-0.42.0)
- qiskit.remote (PennyLane-qiskit-0.42.0)
- lightning.qubit (pennylane_lightning-0.42.0)
- default.clifford (pennylane-0.42.1)
- default.gaussian (pennylane-0.42.1)
- default.mixed (pennylane-0.42.1)
- default.qubit (pennylane-0.42.1)
- default.qutrit (pennylane-0.42.1)
- default.qutrit.mixed (pennylane-0.42.1)
- default.tensor (pennylane-0.42.1)
- null.qubit (pennylane-0.42.1)
- reference.qubit (pennylane-0.42.1)

Hi @Srivathsan_Sundar ,

My colleague Diksha had a look at your problem and made some suggestions! Overall it looks like you might be overcomplicating things a bit. I’m adding Diksha’s thoughts below:

It seems like you’re not using the basis set function in the intended way. Your molecule is quite big but you should be able to simulate it by using the active space functionality (it may take a long time to run though).

When using a certain basis set to generate the molecular Hamiltonian there’s no need to separate the exponents and coefficients nor use the load_basisset function. We can do something like this:

symbols = ["Fe", "S"]
geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], requires_grad=True)

mol = qml.qchem.Molecule(
    symbols,
    geometry,
    basis_name="ANO-RCC-MB",  # Specify the basis set name                                                                                           
    load_data=True        # Crucially, set this flag to True                                                                                          
)

# 3. Pass the Molecule object to molecular_hamiltonian                                                                                                
H, n_qubits = qml.qchem.molecular_hamiltonian(mol)
print(H, n_qubits)

For defining the active space, we need to define the active electrons and orbitals:

H, n_qubits = qml.qchem.molecular_hamiltonian(mol,active_electrons=2,active_orbitals=2)

Note:
You can also use the load_basisset functionality but that can be more complicated. For simple cases like H, when there is only one type of orbital, it’s easier, but in this example you would need to separate the exponents and coefficients based on the orbital type. Currently, it’s done under the hood in qchem.basis_set, with some post-processing done in the Molecule class.

Note 2:
We support only the basis sets mentioned in basissetexchange.org, so anything not mentioned there won’t work.

Let us know if you have any follow-up questions.

I hope this helps!

Thanks for your message @CatalinaAlbornoz. Yeah I am going to use the active space component and make it to 12 qubits system. I din’t know such an easy way of including the basis set name existed. Thanks for highlighting it.

My small suggestion is that in the qml.qchem.Molecule function, the way to add new basis set isn’t clear. If you put this above case as an example. It would be really helpful.

Thanks
Sri

Hi @Srivathsan_Sundar , I’m glad this helped!

I agree that adding this example to the docs could be very helpful. We’ll note it down for improving the docs. Thanks for the suggestion and helping us improve :heart: .

Enjoy using PennyLane!