Note

This page was generated from docs/tutorials/08_property_framework.ipynb.

# The Property Framework¶

```
[1]:
```

```
import numpy as np
```

Qiskit Nature 0.2.0 introduces the *Property* framework. This framework replaces the legacy driver return types like `QMolecule`

and `WatsonHamiltonian`

with a more modular and extensible approach.

In this tutorial, we will walk through the framework, explain its most important components and show you, how you can extend it with a custom *property* yourself.

## What is a `Property`

?¶

At its core, a `Property`

is an object complementing some raw data with functions that allow you to transform/manipulate/interpret this raw data. This “definition” is kept rather abstract on purpose, but what it means is essentially the following. A `Property`

:

represents a physical observable (that’s the raw data)

can be expressed as an operator

can be evaluated with a wavefunction

provides an

`interpret`

method which gives meaning to the eigenvalue of the evaluated qubit operator

```
[2]:
```

```
from qiskit_nature.properties import Property, GroupedProperty
```

The `qiskit_nature.properties`

module provides two classes:

`Property`

: this is the basic interface. It requires only a`name`

and an`interpret`

method to be implemented.`GroupedProperty`

: this class is an implementation of the Composite pattern which allows you to*group*multiple properties into one.

**Note:** Grouped properties must have unique `name`

attributes!

## Second Quantization Properties¶

At the time of writing, Qiskit Nature ships with a single variant of properties: the `SecondQuantizedProperty`

objects.

This sub-type adds one additional requirement because any `SecondQuantizedProperty`

**must**implement a`second_q_ops`

method which constructs a (list of)`SecondQuantizedOp`

s.

The `qiskit_nature.properties.second_quantization`

module is further divided into `electronic`

and `vibrational`

modules (similar to other modules of Qiskit Nature). Let us dive into the `electronic`

sub-module first.

### Electronic Second Quantization Properties¶

Out-of-the-box Qiskit Nature ships the following electronic properties:

```
[3]:
```

```
from qiskit_nature.properties.second_quantization.electronic import (
ElectronicEnergy,
ElectronicDipoleMoment,
ParticleNumber,
AngularMomentum,
Magnetization,
)
```

```
/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/pyscf/lib/misc.py:46: H5pyDeprecationWarning: Using default_file_mode other than 'r' is deprecated. Pass the mode to h5py.File() instead.
h5py.get_config().default_file_mode = 'a'
```

`ElectronicEnergy`

is arguably the most important property because it contains the *Hamiltonian* representing the electronic structure problem whose eigenvalues we are interested in when solving an `ElectronicStructureProblem`

. The way in which it stores this Hamiltonian is, just like the rest of the framework, highly modular. The initializer of the `ElectronicEnergy`

has a single required argument, `electronic_integrals`

, which is a `List[ElectronicIntegrals]`

.

#### ElectronicIntegrals¶

```
[4]:
```

```
from qiskit_nature.properties.second_quantization.electronic.integrals import (
ElectronicIntegrals,
OneBodyElectronicIntegrals,
TwoBodyElectronicIntegrals,
IntegralProperty,
)
```

The `ElectronicIntegrals`

are a container class for *n-body* interactions in a given basis, which can be any of the following:

```
[5]:
```

```
from qiskit_nature.properties.second_quantization.electronic.bases import ElectronicBasis
```

```
[6]:
```

```
list(ElectronicBasis)
```

```
[6]:
```

```
[<ElectronicBasis.AO: 'atomic'>,
<ElectronicBasis.MO: 'molecular'>,
<ElectronicBasis.SO: 'spin'>]
```

Let us take the `OneBodyElectronicIntegrals`

as an example. As the name suggests, this container is for 1-body interaction integrals. You can construct an instance of it like so:

```
[7]:
```

```
one_body_ints = OneBodyElectronicIntegrals(
ElectronicBasis.MO,
(
np.eye(2),
2 * np.eye(2),
),
)
print(one_body_ints)
```

```
(MO) 1-Body Terms:
Alpha
<(2, 2) matrix with 2 non-zero entries>
[0, 0] = 1.0
[1, 1] = 1.0
Beta
<(2, 2) matrix with 2 non-zero entries>
[0, 0] = 2.0
[1, 1] = 2.0
```

As you can see, the first argument simply specifies the basis of the integrals. The second argument requires a little further explanation:

In the case of the

`AO`

or`MO`

basis, this argument**must**be a pair of numpy arrays, where the first and second one are the alpha- and beta-spin integrals, respectively.

**Note:** The second argument may be `None`

, for the case where the beta-spin integrals are the same as the alpha-spin integrals (so there is no need to provide the same values twice).

In the case of the

`SO`

basis, this argument**must**be a single numpy array, storing the alpha- and beta-spin integrals in blocked order, i.e. like so:

```
spin_basis = np.block([[alpha_spin, zeros], [zeros, beta_spin]])
```

The `TwoBodyElectronicIntegrals`

work pretty much the same except that they contain all possible spin-combinations in the tuple of numpy arrays. For example:

```
[8]:
```

```
two_body_ints = TwoBodyElectronicIntegrals(
ElectronicBasis.MO,
(
np.arange(1, 17).reshape((2, 2, 2, 2)),
np.arange(16, 32).reshape((2, 2, 2, 2)),
np.arange(-16, 0).reshape((2, 2, 2, 2)),
None,
),
)
print(two_body_ints)
```

```
(MO) 2-Body Terms:
Alpha-Alpha
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 1.0
[0, 0, 0, 1] = 2.0
[0, 0, 1, 0] = 3.0
[0, 0, 1, 1] = 4.0
[0, 1, 0, 0] = 5.0
... skipping 11 entries
Beta-Alpha
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 16.0
[0, 0, 0, 1] = 17.0
[0, 0, 1, 0] = 18.0
[0, 0, 1, 1] = 19.0
[0, 1, 0, 0] = 20.0
... skipping 11 entries
Beta-Beta
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = -16.0
[0, 0, 0, 1] = -15.0
[0, 0, 1, 0] = -14.0
[0, 0, 1, 1] = -13.0
[0, 1, 0, 0] = -12.0
... skipping 11 entries
Alpha-Beta
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 16.0
[0, 0, 0, 1] = 24.0
[0, 0, 1, 0] = 20.0
[0, 0, 1, 1] = 28.0
[0, 1, 0, 0] = 18.0
... skipping 11 entries
```

We should take note of a few observations:

the numpy arrays shall be ordered as

`("alpha-alpha", "beta-alpha", "beta-beta", "alpha-beta")`

the

`alpha-alpha`

matrix may**not**be`None`

if the

`alpha-beta`

matrix is`None`

, but`beta-alpha`

is not, its transpose will be used (like above)in any other case, matrices which are

`None`

will be filled with the`alpha-alpha`

matrixin the

`SO`

basis case, a single numpy array must be specified. Refer to`TwoBodyElectronicIntegrals.to_spin()`

for its exact formatting.

#### ElectronicEnergy¶

Now we are ready to construct an `ElectronicEnergy`

instance:

```
[9]:
```

```
electronic_energy = ElectronicEnergy(
[one_body_ints, two_body_ints],
)
print(electronic_energy)
```

```
ElectronicEnergy
(MO) 1-Body Terms:
Alpha
<(2, 2) matrix with 2 non-zero entries>
[0, 0] = 1.0
[1, 1] = 1.0
Beta
<(2, 2) matrix with 2 non-zero entries>
[0, 0] = 2.0
[1, 1] = 2.0
(MO) 2-Body Terms:
Alpha-Alpha
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 1.0
[0, 0, 0, 1] = 2.0
[0, 0, 1, 0] = 3.0
[0, 0, 1, 1] = 4.0
[0, 1, 0, 0] = 5.0
... skipping 11 entries
Beta-Alpha
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 16.0
[0, 0, 0, 1] = 17.0
[0, 0, 1, 0] = 18.0
[0, 0, 1, 1] = 19.0
[0, 1, 0, 0] = 20.0
... skipping 11 entries
Beta-Beta
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = -16.0
[0, 0, 0, 1] = -15.0
[0, 0, 1, 0] = -14.0
[0, 0, 1, 1] = -13.0
[0, 1, 0, 0] = -12.0
... skipping 11 entries
Alpha-Beta
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 16.0
[0, 0, 0, 1] = 24.0
[0, 0, 1, 0] = 20.0
[0, 0, 1, 1] = 28.0
[0, 1, 0, 0] = 18.0
... skipping 11 entries
```

This property can now be used to construct a `SecondQuantizedOp`

(which can then be mapped to a `QubitOperator`

):

```
[10]:
```

```
hamiltonian = electronic_energy.second_q_ops()[0] # here, output length is always 1
print(hamiltonian)
```

```
Fermionic Operator
register length=4, number terms=20
(22+0j) * ( +_0 -_1 +_2 -_3 )
+ (-26+0j) * ( +_0 -_1 -_2 +_3 )
+ (30+0j) * ( +_0 -_1 +_3 -_3 )
+ (18+0j) * ( +_0 -_1 +_2 -_2 )
+ (-21+0j) * ( -_0 +_1 +_2 -_3 )
+ (25+0j) * ( -_0 +_1 -_2 +_3 )
+ ...
```

#### Result interpretation¶

An additional benefit which we gain from the `Property`

framework, is that the result interpretation of a computed eigenvalue can be handled by each property itself. This results in nice and logically consistent classes because the result gets interpreted in the same context where the raw data is available.

```
[11]:
```

```
from qiskit_nature.results import ElectronicStructureResult
# some dummy result
result = ElectronicStructureResult()
result.eigenenergies = np.asarray([-1])
result.computed_energies = np.asarray([-1])
# now, let's interpret it
electronic_energy.interpret(result)
print(result)
```

```
=== GROUND STATE ENERGY ===
* Electronic ground state energy (Hartree): -1.
- computed part: -1.
```

While this particular example is not yet very impressive, wait until we use more properties at once.

#### ParticleNumber¶

The `ParticleNumber`

property also takes a special place among the builtin properties because it serves a dual purpose of:

storing the number of particles in the electronic system

providing the operators to evaluate the number of particles for a given eigensolution of the problem

Therefore, this property is required if you want to use additional functionality like the `ActiveSpaceTransformer`

or the `ElectronicStructureProblem.default_filter_criterion()`

.

```
[12]:
```

```
particle_number = ParticleNumber(
num_spin_orbitals = 4,
num_particles = (1, 1),
)
print(particle_number)
```

```
ParticleNumber:
4 SOs
1 alpha electrons
orbital occupation: [1. 0.]
1 beta electrons
orbital occupation: [1. 0.]
```

#### GroupedProperty¶

Rather than iterating all of the other properties one by one, let us simply consider a group of properties as provided by any `ElectronicStructureDriver`

from the `qiskit_nature.drivers.second_quantization`

module.

```
[13]:
```

```
from qiskit_nature.drivers.second_quantization.pyscfd import PySCFDriver
```

```
[14]:
```

```
electronic_driver = PySCFDriver(atom="H 0 0 0; H 0 0 0.735", basis="sto3g")
electronic_driver_result = electronic_driver.run()
```

```
[15]:
```

```
print(electronic_driver_result)
```

```
ElectronicStructureDriverResult:
DriverMetadata:
Program: PYSCF
Version: 1.7.6
Config:
atom=H 0 0 0; H 0 0 0.735
unit=Angstrom
charge=0
spin=0
basis=sto3g
method=rhf
conv_tol=1e-09
max_cycle=50
init_guess=minao
max_memory=4000
ElectronicBasisTransform:
Initial basis: atomic
Final basis: molecular
Alpha coefficients:
[0, 0] = 0.5483020229014736
[0, 1] = -1.2183273138546826
[1, 0] = 0.548302022901473
[1, 1] = 1.2183273138546828
Beta coefficients:
[0, 0] = 0.5483020229014736
[0, 1] = -1.2183273138546826
[1, 0] = 0.548302022901473
[1, 1] = 1.2183273138546828
ParticleNumber:
4 SOs
1 alpha electrons
orbital occupation: [1. 0.]
1 beta electrons
orbital occupation: [1. 0.]
ElectronicEnergy
(AO) 1-Body Terms:
Alpha
<(2, 2) matrix with 4 non-zero entries>
[0, 0] = -1.1242175791954514
[0, 1] = -0.9652573993472758
[1, 0] = -0.9652573993472758
[1, 1] = -1.1242175791954512
Beta
<(2, 2) matrix with 4 non-zero entries>
[0, 0] = -1.1242175791954514
[0, 1] = -0.9652573993472758
[1, 0] = -0.9652573993472758
[1, 1] = -1.1242175791954512
(AO) 2-Body Terms:
Alpha-Alpha
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 0.7746059439198978
[0, 0, 0, 1] = 0.4474457245330949
[0, 0, 1, 0] = 0.447445724533095
[0, 0, 1, 1] = 0.5718769760004512
[0, 1, 0, 0] = 0.4474457245330951
... skipping 11 entries
Beta-Alpha
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 0.7746059439198978
[0, 0, 0, 1] = 0.4474457245330949
[0, 0, 1, 0] = 0.447445724533095
[0, 0, 1, 1] = 0.5718769760004512
[0, 1, 0, 0] = 0.4474457245330951
... skipping 11 entries
Beta-Beta
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 0.7746059439198978
[0, 0, 0, 1] = 0.4474457245330949
[0, 0, 1, 0] = 0.447445724533095
[0, 0, 1, 1] = 0.5718769760004512
[0, 1, 0, 0] = 0.4474457245330951
... skipping 11 entries
Alpha-Beta
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 0.7746059439198978
[0, 0, 0, 1] = 0.4474457245330949
[0, 0, 1, 0] = 0.447445724533095
[0, 0, 1, 1] = 0.5718769760004512
[0, 1, 0, 0] = 0.4474457245330951
... skipping 11 entries
(MO) 1-Body Terms:
Alpha
<(2, 2) matrix with 2 non-zero entries>
[0, 0] = -1.2563390730032507
[1, 1] = -0.47189600728114073
Beta
<(2, 2) matrix with 2 non-zero entries>
[0, 0] = -1.2563390730032507
[1, 1] = -0.47189600728114073
(MO) 2-Body Terms:
Alpha-Alpha
<(2, 2, 2, 2) matrix with 8 non-zero entries>
[0, 0, 0, 0] = 0.675710154803517
[0, 0, 1, 1] = 0.6645817302552974
[0, 1, 0, 1] = 0.18093119978423153
[0, 1, 1, 0] = 0.18093119978423136
[1, 0, 0, 1] = 0.18093119978423164
... skipping 3 entries
Beta-Alpha
<(2, 2, 2, 2) matrix with 8 non-zero entries>
[0, 0, 0, 0] = 0.675710154803517
[0, 0, 1, 1] = 0.6645817302552974
[0, 1, 0, 1] = 0.18093119978423153
[0, 1, 1, 0] = 0.18093119978423136
[1, 0, 0, 1] = 0.18093119978423164
... skipping 3 entries
Beta-Beta
<(2, 2, 2, 2) matrix with 8 non-zero entries>
[0, 0, 0, 0] = 0.675710154803517
[0, 0, 1, 1] = 0.6645817302552974
[0, 1, 0, 1] = 0.18093119978423153
[0, 1, 1, 0] = 0.18093119978423136
[1, 0, 0, 1] = 0.18093119978423164
... skipping 3 entries
Alpha-Beta
<(2, 2, 2, 2) matrix with 8 non-zero entries>
[0, 0, 0, 0] = 0.675710154803517
[0, 0, 1, 1] = 0.6645817302552974
[0, 1, 0, 1] = 0.18093119978423153
[0, 1, 1, 0] = 0.18093119978423136
[1, 0, 0, 1] = 0.18093119978423164
... skipping 3 entries
ElectronicDipoleMoment:
DipoleMomentX
(AO) 1-Body Terms:
Alpha
<(2, 2) matrix with 0 non-zero entries>
Beta
<(2, 2) matrix with 0 non-zero entries>
(MO) 1-Body Terms:
Alpha
<(2, 2) matrix with 0 non-zero entries>
Beta
<(2, 2) matrix with 0 non-zero entries>
DipoleMomentY
(AO) 1-Body Terms:
Alpha
<(2, 2) matrix with 0 non-zero entries>
Beta
<(2, 2) matrix with 0 non-zero entries>
(MO) 1-Body Terms:
Alpha
<(2, 2) matrix with 0 non-zero entries>
Beta
<(2, 2) matrix with 0 non-zero entries>
DipoleMomentZ
(AO) 1-Body Terms:
Alpha
<(2, 2) matrix with 3 non-zero entries>
[0, 1] = 0.4605377079660319
[1, 0] = 0.4605377079660319
[1, 1] = 1.3889487015553204
Beta
<(2, 2) matrix with 3 non-zero entries>
[0, 1] = 0.4605377079660319
[1, 0] = 0.4605377079660319
[1, 1] = 1.3889487015553204
(MO) 1-Body Terms:
Alpha
<(2, 2) matrix with 4 non-zero entries>
[0, 0] = 0.69447435077766
[0, 1] = 0.9278334704592323
[1, 0] = 0.9278334704592324
[1, 1] = 0.6944743507776605
Beta
<(2, 2) matrix with 4 non-zero entries>
[0, 0] = 0.69447435077766
[0, 1] = 0.9278334704592323
[1, 0] = 0.9278334704592324
[1, 1] = 0.6944743507776605
AngularMomentum:
4 SOs
Magnetization:
4 SOs
Molecule:
Multiplicity: 1
Charge: 0
Geometry:
H [0.0, 0.0, 0.0]
H [0.0, 0.0, 1.3889487015553204]
Masses:
H 1
H 1
```

There is a lot going on but with the explanations above you should not have any problems with understanding this output.

#### Constructing a Hamiltonian from raw integrals¶

If you have obtained some raw one- and two-body integrals by means other than through a driver provided by Qiskit Nature, you can still easily construct an `ElectronicEnergy`

property to serve as your access point into the stack:

```
[16]:
```

```
one_body_ints = np.arange(1, 5).reshape((2, 2))
two_body_ints = np.arange(1, 17).reshape((2, 2, 2, 2))
electronic_energy_from_ints = ElectronicEnergy.from_raw_integrals(ElectronicBasis.MO, one_body_ints, two_body_ints)
print(electronic_energy_from_ints)
```

```
ElectronicEnergy
(MO) 1-Body Terms:
Alpha
<(2, 2) matrix with 4 non-zero entries>
[0, 0] = 1.0
[0, 1] = 2.0
[1, 0] = 3.0
[1, 1] = 4.0
Beta
<(2, 2) matrix with 4 non-zero entries>
[0, 0] = 1.0
[0, 1] = 2.0
[1, 0] = 3.0
[1, 1] = 4.0
(MO) 2-Body Terms:
Alpha-Alpha
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 1.0
[0, 0, 0, 1] = 2.0
[0, 0, 1, 0] = 3.0
[0, 0, 1, 1] = 4.0
[0, 1, 0, 0] = 5.0
... skipping 11 entries
Beta-Alpha
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 1.0
[0, 0, 0, 1] = 2.0
[0, 0, 1, 0] = 3.0
[0, 0, 1, 1] = 4.0
[0, 1, 0, 0] = 5.0
... skipping 11 entries
Beta-Beta
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 1.0
[0, 0, 0, 1] = 2.0
[0, 0, 1, 0] = 3.0
[0, 0, 1, 1] = 4.0
[0, 1, 0, 0] = 5.0
... skipping 11 entries
Alpha-Beta
<(2, 2, 2, 2) matrix with 16 non-zero entries>
[0, 0, 0, 0] = 1.0
[0, 0, 0, 1] = 2.0
[0, 0, 1, 0] = 3.0
[0, 0, 1, 1] = 4.0
[0, 1, 0, 0] = 5.0
... skipping 11 entries
```

### Vibrational Second Quantization Properties¶

The `vibrational`

stack for vibrational structure calculations also integrates with the Property framework. After the above introduction you should be able to understand the following directly:

```
[17]:
```

```
from qiskit_nature.drivers.second_quantization.gaussiand import GaussianForcesDriver
```

```
[18]:
```

```
# if you ran Gaussian elsewhere and already have the output file
vibrational_driver = GaussianForcesDriver(logfile='aux_files/CO2_freq_B3LYP_ccpVDZ.log')
vibrational_driver_result = vibrational_driver.run()
```

For vibrational structure calculations we always need to define the basis which we want to work in, separately:

```
[19]:
```

```
from qiskit_nature.properties.second_quantization.vibrational.bases import HarmonicBasis
```

```
[20]:
```

```
harmonic_basis = HarmonicBasis([2] * 4)
```

```
[21]:
```

```
vibrational_driver_result.basis = harmonic_basis
print(vibrational_driver_result)
```

```
VibrationalStructureDriverResult:
HarmonicBasis:
Modals: [2, 2, 2, 2]:
VibrationalEnergy:
HarmonicBasis:
Modals: [2, 2, 2, 2]
1-Body Terms:
<sparse integral list with 13 entries>
(1, 1) = 605.3643675
(-1, -1) = -605.3643675
(2, 2) = 340.5950575
(-2, -2) = -340.5950575
(3, 3) = 163.7595125
... skipping 8 entries
2-Body Terms:
<sparse integral list with 15 entries>
(2, 1, 1) = -89.09086530649508
(3, 3, 2) = 21.644966371722838
(4, 4, 2) = 6.412754934114705
(2, 2, 1, 1) = 5.03965375
(3, 1, 1, 1) = -2.4473854166666666
... skipping 10 entries
3-Body Terms:
<sparse integral list with 9 entries>
(3, 2, 1) = 44.01468537435673
(4, 2, 1) = -78.71701132125833
(4, 3, 2) = 17.15529085952822
(3, 2, 2, 1) = -3.73513125
(4, 2, 2, 1) = 6.68000625
... skipping 4 entries
OccupiedModals:
HarmonicBasis:
Modals: [2, 2, 2, 2]
```

## Writing custom Properties¶

You can extend the Property framework with your own implementations. Here, we will walk through the simple example of constructing a Property for the *electronic density*. Since this observable is essentially based on matrices, we will be leveraging the `OneBodyElectronicIntegrals`

container to store the actual density matrix.

```
[22]:
```

```
from itertools import product
from typing import List
from qiskit_nature.drivers import QMolecule
from qiskit_nature.operators.second_quantization import FermionicOp
from qiskit_nature.properties.second_quantization.electronic.bases import ElectronicBasis
from qiskit_nature.properties.second_quantization.electronic.types import ElectronicProperty
from qiskit_nature.properties.second_quantization.electronic.integrals import OneBodyElectronicIntegrals
```

```
[23]:
```

```
class ElectronicDensity(ElectronicProperty):
"""A simple electronic density property.
This basic example works only in the MO basis!
"""
def __init__(self, num_molecular_orbitals: int) -> None:
super().__init__(self.__class__.__name__)
self._num_molecular_orbitals = num_molecular_orbitals
def __str__(self) -> str:
string = [super().__str__() + ":"]
string += [f"\t{self._num_molecular_orbitals} MOs"]
return "\n".join(string)
@classmethod
def from_legacy_driver_result(cls, result) -> "ElectronicDensity":
cls._validate_input_type(result, QMolecule)
qmol = cast(QMolecule, result)
return cls(qmol.num_molecular_orbitals)
def second_q_ops(self) -> List[FermionicOp]:
ops = []
# iterate all pairs of molecular orbitals
for mo_i, mo_j in product(range(self._num_molecular_orbitals), repeat=2):
# construct an auxiliary matrix where the only non-zero entry is at the current pair of MOs
number_op_matrix = np.zeros((self._num_molecular_orbitals, self._num_molecular_orbitals))
number_op_matrix[mo_i, mo_j] = 1
# leverage the OneBodyElectronicIntegrals to construct the corresponding FermionicOp
one_body_ints = OneBodyElectronicIntegrals(ElectronicBasis.MO, (number_op_matrix, number_op_matrix))
ops.append(one_body_ints.to_second_q_op())
return ops
def interpret(self, result) -> None:
# here goes the code which interprets the eigenvalues returned for the auxiliary operators
pass
```

```
[24]:
```

```
density = ElectronicDensity(2)
```

```
[25]:
```

```
print(density)
```

```
ElectronicDensity:
2 MOs
```

```
[26]:
```

```
for idx, op in enumerate(density.second_q_ops()):
print(idx, ":", op)
```

```
0 : Fermionic Operator
register length=4, number terms=2
(1+0j) * ( +_0 -_0 )
+ (1+0j) * ( +_2 -_2 )
1 : Fermionic Operator
register length=4, number terms=2
(1+0j) * ( +_0 -_1 )
+ (1+0j) * ( +_2 -_3 )
2 : Fermionic Operator
register length=4, number terms=2
(1+0j) * ( +_1 -_0 )
+ (1+0j) * ( +_3 -_2 )
3 : Fermionic Operator
register length=4, number terms=2
(1+0j) * ( +_1 -_1 )
+ (1+0j) * ( +_3 -_3 )
```

Of course, the above example is very minimal and can be extended at will.

**Note:** as of Qiskit Nature version 0.2.0, the direct integration of custom Property objects into the stack is not implemented yet, due to limitations of the auxiliary operator parsing. See https://github.com/Qiskit/qiskit-nature/issues/312 for more details.

For the time being, you can still evaluate a custom Property, by passing it’s generated operators directly to the `Eigensolver.solve`

method by means of constructing the `aux_operators`

. Note, however, that you will have to deal with transformations applied to your properties manually, until the above issue is resolved.

```
[27]:
```

```
# set up some problem
problem = ...
# set up a solver
solver = ...
# when solving the problem, pass additional operators in like so:
aux_ops = density.second_q_ops()
# solver.solve(problem, aux_ops)
```

```
[28]:
```

```
import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright
```

### Version Information

Qiskit Software | Version |
---|---|

`qiskit-terra` | 0.18.3 |

`qiskit-aer` | 0.9.1 |

`qiskit-nature` | 0.2.2 |

System information | |

Python | 3.8.12 (default, Sep 13 2021, 08:28:12) [GCC 9.3.0] |

OS | Linux |

CPUs | 2 |

Memory (Gb) | 6.790920257568359 |

Wed Oct 13 20:53:44 2021 UTC |

### This code is a part of Qiskit

© Copyright IBM 2017, 2021.

This code is licensed under the Apache License, Version 2.0. You may

obtain a copy of this license in the LICENSE.txt file in the root directory

of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this

copyright notice, and modified files need to carry a notice indicating

that they have been altered from the originals.