How does setting active electrons and orbitals work exactly?

If I’ve got

hamiltonian, qubits = qchem.molecular_hamiltonian( symbols, coordinates, active_electrons=4, active_orbitals=4, )

How are the active orbitals chosen? Are they automatically the 2 highest occupied orbitals and the 2 lowest unoccupied orbitals? In the case of water, wouldn’t this be an issue since the highest occupied orbital would be occupied by the lone pair instead of the bonding electrons? Is there a way to select exactly which orbitals are chosen as active? What if I wanted to consider only the orbital occupied by the lone pair and its excitations?

Hi @somearthling ,

The first figure in the molecular_hamiltonian documentation can help you visualize this. If you don’t set active_orbitals= then all of the orbitals will be considered active.

If you do set active_orbitals=something then the documentation for active_space can help you understand how it works. First it calculates the number of core orbitals with ncore_orbs = (electrons - active_electrons) // 2. These are the lowest orbitals and they are always occupied by two electrons. If you specified active_orbitals=4 then the next 4 orbitals will be considered active orbitals. The remaining orbitals (if any) will be considered external orbitals, hence non-active.

We might have a function to select the orbitals arbitrarily. I’ll look into it and get back to you on this.

Hi @somearthling ,

My colleague Alain has this great answer for you!

The total number of active orbitals is the number of orbitals selected below and above the Highest-Occupied Molecular Orbital (HOMO) to be included in the simulation (see Fig. here). Note that the number of orbitals below and above the HOMO does not need to be same.
Example:

  • water molecule, minimal basis set => 10 electrons, 7 orbitals
  • For spin multiplicity (2S+1)=1 => 5 doubly-occupied and 2 unoccupied orbitals
  • The code snippet below shows how to build an active space with 3 orbitals (the fifth doubly-occupied orbital + 2 orbitals above the HOMO)
>>> from pennylane import qchem

>>> electrons = 10
>>> orbitals = 7
>>> core, active = qchem.active_space(electrons, orbitals, active_electrons=2, active_orbitals=3)
>>> print("Core orbitals:", core)
>>> print("Active orbitals:", active)
Core orbitals: [0, 1, 2, 3]
Active orbitals: [4, 5, 6]  

You can build the corresponding Hamiltonian as:

>>> import numpy as np

>>> symbols = ["H", "O", "H"]
>>> coordinates = np.array([-0.0399, -0.0038, 0.0, 1.5780, 0.8540, 0.0, 2.7909, -0.5159, 0.0])

>>> H, qubits = qchem.molecular_hamiltonian(symbols, coordinates, active_electrons=2, active_orbitals=3)

Hope this helps!

Hey @CatalinaAlbornoz,

Thanks for the answer, but I’m trying to understand how I can specify active and core orbitals - if you take a look at this paper, they use HOMO -2, HOMO -1, LUMO and LUMO+1 as the active space, since the HOMO mainly corresponds to lone pair effects rather than the bonding electrons. I’m trying to do something similar so was wondering if there’s a way for me to specify - this would be

core = [0,1,4] active = [2, 3, 5, 6]

but there’s no way to feed this into molecular_hamiltonian which only accepts integers for the active orbitals argument

That’s right @somearthling.

We recognize that this functionality can be useful for cases like yours but unfortunately we currently don’t have capacity to implement it. Is this a functionality that you would like to implement yourself?