Qchem.electron_integrals

PennyLane Help
1

Subject: Issue with Open-shell Systems in Resource Estimation for Quantum Chemistry

Hi everyone,

I was working on this demo from Soran that caught my attention: Resource Estimation for Quantum Chemistry.

However, I encountered a constraint: the function does not support open-shell systems, displaying the following error message:

ValueError: Open-shell systems are not supported. Change the charge or spin multiplicity of the molecule.

This limitation means that we cannot fully use it for the larger molecules I’m working on. Has anyone faced this issue before? Do we have an alternative solution?

Here’s the code to replicate the error message:

# Molecule definition  
<Molecule = C3H7, Charge: 0, Basis: STO-3G, Orbitals: 22, Electrons: 25>

# Error traceback
line 31 ValueError: Open-shell systems are not supported. Change the charge or spin multiplicity of the molecule.
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[174], line 31
     29 print("Molecule CH_molecule:")
     30 print(mol)
---> 31 core, one, two = qml.qchem.electron_integrals(mol)()
     32 algo = qml.resource.DoubleFactorization(one, two)
     33 print(f'Estimated gates : {algo.gates:.2e} \nEstimated qubits: {algo.qubits}')

File /opt/anaconda3/envs/pen/lib/python3.10/site-packages/pennylane/qchem/hamiltonian.py:120, in electron_integrals.<locals>._electron_integrals(*args)
    111 def _electron_integrals(*args):
    112     r"""Compute the one- and two-electron integrals in the molecular orbital basis.
    113 
    114     Args:
   (...)
    118         tuple[array[float]]: 1D tuple containing core constant, one- and two-electron integrals
    119     """
--> 120     _, coeffs, _, h_core, repulsion_tensor = scf(mol)(*args)
    121     one = qml.math.einsum("qr,rs,st->qt", coeffs.T, h_core, coeffs)
    122     two = qml.math.swapaxes(
    123         qml.math.einsum(
    124             "ab,cd,bdeg,ef,gh->acfh", coeffs.T, coeffs.T, repulsion_tensor, coeffs, coeffs
   (...)
    127         3,
    128     )

File /opt/anaconda3/envs/pen/lib/python3.10/site-packages/pennylane/qchem/hartree_fock.py:120, in scf.<locals>._scf(*args)
    110 r"""Perform the self-consistent-field iterations.
    111 
    112 Args:
   (...)
    117     Fock matrix, core matrix
    118 """
    119 if mol.n_electrons % 2 == 1 or mol.mult != 1:
--> 120     raise ValueError(
    121         "Open-shell systems are not supported. Change the charge or spin multiplicity of the molecule."
    122     )
    124 basis_functions = mol.basis_set
    125 charges = mol.nuclear_charges

Any insights or suggestions would be greatly appreciated.

Many thanks,
Parfait

2

Hi @Parfait_Atchade ,

Are you able to share the code you’re using to generate the molecule? The open-shell error is probably not related to the size of the molecule. Although this is indeed a huge molecule so I want to try to replicate your error.

1 Like
3

Hi @Catalina,

I thought I had shared this before—sorry for that! Here is the code.

Many thanks,
Parfait.

import numpy as np
import pennylane as qml

# Atomic coordinates (in bohr)
coordinates_list = np.array([
    [1.3122, -0.3072, -0.0711],  # Carbon (C)
    [0.0728,  0.5625,  0.0627],  # Carbon (C)
    [-1.2351, -0.2499, -0.0399], # Carbon (C)
    [-2.0976,  0.4019,  0.0633], # Hydrogen (H)
    [-1.3009, -0.7534, -0.9997], # Hydrogen (H)
    [-1.2810, -1.0023,  0.7420], # Hydrogen (H)
    [0.0901,   1.0865,  1.0210], # Hydrogen (H)
    [0.0843,   1.3280, -0.7122], # Hydrogen (H)
    [2.2655,   0.1678, -0.2663], # Hydrogen (H)
    [1.3400,   1.2607,  0.4413]  # Hydrogen (H)
])

# Element symbols in the same order as the coordinates
symbols = ["C", "C", "C", "H", "H", "H", "H", "H", "H", "H"]

# Flatten coordinates into a 1D array
coordinates = coordinates_list.flatten()

# Create the molecule
mol = qml.qchem.Molecule(
    symbols=symbols,
    coordinates=coordinates,
    charge=0,
    mult=1,
    basis_name='sto-3g',
    name='CH_molecule',
    load_data=False,
    normalize=True,
    unit='bohr'
)

# Print molecule information
print("CH_molecule structure:")
print(mol)

# Compute electronic integrals
core, one, two = qml.qchem.electron_integrals(mol)()

# Apply double factorization
algo = qml.resource.DoubleFactorization(one, two)

# Display estimated quantum resources
print(f'Estimated gates: {algo.gates:.2e}')
print(f'Estimated qubits: {algo.qubits}')

# Plot factorization progress (handling missing function error)
try:
    plot_factorization_progressive(one, two, num_points_coarse=10, num_points_fine=20)
except NameError:
    print("Error: The function 'plot_factorization_progressive' is not defined.")
4

Thanks @Parfait_Atchade !

So the issue here is that C3H7 has an extra electron lying around, meaning that it’s an open-shell system, which isn’t supported. The way to fix this is by computing C3H7+ instead. All you need to do is set charge=1 when you define your molecule with qml.qchem.Molecule. Alternatively, you could try computing C3H8 in case you prefer to compute a molecule with neutral charge.

That being said, C3H7+ and C3H8 are very large molecules. I suggest you set several orbitals as core when you calculate the electron_integrals because otherwise it may take a long time to calculate them and you may run into memory issues. Below I share one example of how you can do this:

# Set active space
electrons = mol.n_electrons
orbitals = mol.n_orbitals
core, active = qml.qchem.active_space(electrons, orbitals, active_electrons=4, active_orbitals=4)

# Compute electronic integrals
core, one, two = qml.qchem.electron_integrals(mol, core=core, active=active)()

I’m not sure why but it’s still taking a long time to run even with the much reduced active space, so let me know if you manage tu run this successfully.

At least you won’t be getting the error you were seeing before.

I hope this helps!

1 Like
5

Hi @Catalina,

Thank you very much for your suggestions. I followed your advice by setting charge=1 in qml.qchem.Molecule and also tried using C₃H₈ to avoid dealing with an open-shell system. The calculations are still taking quite a bit of time, which I understand is normal given the combinatorial cost associated with the Slater determinant.

Thank you again for your help and time.

Best regards,

6

No problem @Parfait_Atchade !

I ran it with the reduced active space and it managed to complete the computation. I didn’t time it though, so I don’t know how long it took.

I hope your computation works out too!

7

Hi @Parfait_Atchade ,

I realized that there might be a way to reduce the computation time for the resource estimation that you’re doing.

I’m not sure that it will work with your workflow, I haven’t tested it yet, but we have a device called NullQubit. This doesn’t do computations (so the answer from any computation would be wrong) but if you just want to calculate resources needed this can be the way to go. You can use it by defining your device as dev = qml.device("null.qubit").

You can then use qml.Tracker or qml.specs.

Here’s how you can use qml.Tracker

with qml.Tracker(dev) as tracker:
    circuit(params)

I hope this helps!

8

Hi Catalina,

Thank you for reaching out and for the suggestion! I wasn’t aware of NullQubit, but it sounds like a great way to speed up the resource estimation process. I’ll definitely look into integrating it into my workflow and see if it aligns with my current approach.

I appreciate the pointers on using qml.Tracker and qml.specs—this could be a useful way to assess the computational requirements without running full calculations. I’ll test it out and let you know if I have any questions or run into any issues.

Thanks again for the insight!

Best,
Parfait

1 Like