Run interactively in jupyter notebook.

Electronic structure


The molecular Hamiltonian is

\[\mathcal{H} = - \sum_I \frac{\nabla_{R_I}^2}{M_I} - \sum_i \frac{\nabla_{r_i}^2}{m_e} - \sum_I\sum_i \frac{Z_I e^2}{|R_I-r_i|} + \sum_i \sum_{j>i} \frac{e^2}{|r_i-r_j|} + \sum_I\sum_{J>I} \frac{Z_I Z_J e^2}{|R_I-R_J|}\]

Because the nuclei are much heavier than the electrons they do not move on the same time scale and therefore, the behavior of nuclei and electrons can be decoupled. This is the Born-Oppenheimer approximation.

Therefore, one can first tackle the electronic problem with nuclear coordinate entering only as parameters. The energy levels of the electrons in the molecule can be found by solving the non-relativistic time independent Schroedinger equation,

\[\mathcal{H}_{\text{el}} |\Psi_{n}\rangle = E_{n} |\Psi_{n}\rangle\]


\[\mathcal{H}_{\text{el}} = - \sum_i \frac{\nabla_{r_i}^2}{m_e} - \sum_I\sum_i \frac{Z_I e^2}{|R_I-r_i|} + \sum_i \sum_{j>i} \frac{e^2}{|r_i-r_j|}.\]

In particular the ground state energy is given by:

\[E_0 = \frac{\langle \Psi_0 | H_{\text{el}} | \Psi_0 \rangle}{\langle \Psi_0 | \Psi_0 \rangle}\]

where \(\Psi_0\) is the ground state of the system.

However, the dimensionality of this problem grows exponentially with the number of degrees of freedom. To tackle this issue we would like to prepare \(\Psi_0\) on a quantum computer and measure the Hamiltonian expectation value (or \(E_0\)) directly.

So how do we do that concretely?

The Hartree-Fock initial state

A good starting point for solving this problem is the Hartree-Fock (HF) method. This method approximates a N-body problem into N one-body problems where each electron evolves in the mean-field of the others. Classically solving the HF equations is efficient and leads to the exact exchange energy but does not include any electron correlation. Therefore, it is usually a good starting point to start adding correlation.

The Hamiltonian can then be re-expressed in the basis of the solutions of the HF method, also called Molecular Orbitals (MOs):

\[\hat{H}_{elec}=\sum_{pq} h_{pq} \hat{a}^{\dagger}_p \hat{a}_q + \frac{1}{2} \sum_{pqrs} h_{pqrs} \hat{a}^{\dagger}_p \hat{a}^{\dagger}_q \hat{a}_r \hat{a}_s\]

with the 1-body integrals

\[h_{pq} = \int \phi^*_p(r) \left( -\frac{1}{2} \nabla^2 - \sum_{I} \frac{Z_I}{R_I- r} \right) \phi_q(r)\]

and 2-body integrals

\[h_{pqrs} = \int \frac{\phi^*_p(r_1) \phi^*_q(r_2) \phi_r(r_2) \phi_s(r_1)}{|r_1-r_2|}.\]

The MOs (\(\phi_u\)) can be occupied or virtual (unoccupied). One MO can contain 2 electrons. However, in what follows we actually work with Spin Orbitals which are associated with a spin up (\(\alpha\)) of spin down (\(\beta\)) electron. Thus Spin Orbitals can contain one electron or be unoccupied.

We now show how to concretely realise these steps with Qiskit.

Qiskit is interfaced with different classical codes which are able to find the HF solutions. Interfacing between Qiskit and the following codes is already available: * Gaussian * Psi4 * PyQuante * PySCF

In the following we set up a PySCF driver, for the hydrogen molecule at equilibrium bond length (0.735 angstrom) in the singlet state and with no charge.

from qiskit_nature.drivers import PySCFDriver, UnitsType, Molecule
molecule = Molecule(geometry=[['H', [0., 0., 0.]],
                              ['H', [0., 0., 0.735]]],
                     charge=0, multiplicity=1)
driver = PySCFDriver(molecule = molecule, unit=UnitsType.ANGSTROM, basis='sto3g')
/opt/hostedtoolcache/Python/3.8.10/x64/lib/python3.8/site-packages/pyscf/lib/ H5pyDeprecationWarning: Using default_file_mode other than 'r' is deprecated. Pass the mode to h5py.File() instead.
  h5py.get_config().default_file_mode = 'a'

For further information about the drivers see

The mapping from fermions to qubits


The Hamiltonian given in the previous section is expressed in terms of fermionic operators. To encode the problem into the state of a quantum computer, these operators must be mapped to spin operators (indeed the qubits follow spin statistics).

There exist different mapping types with different properties. Qiskit already supports the following mappings: * The Jordan-Wigner ‘jordan_wigner’ mapping (über das paulische äquivalenzverbot. In The Collected Works of Eugene Paul Wigner (pp. 109-129). Springer, Berlin, Heidelberg (1993)). * The Parity ‘parity’ (The Journal of chemical physics, 137(22), 224109 (2012)) * The Bravyi-Kitaev ‘bravyi_kitaev’ (Annals of Physics, 298(1), 210-226 (2002))

The Jordan-Wigner mapping is particularly interesting as it maps each Spin Orbital to a qubit (as shown on the Figure above).

Here we set up the Electronic Structure Problem to generate the Second quantized operator and a qubit converter that will map it to a qubit operator.

from qiskit_nature.problems.second_quantization import ElectronicStructureProblem
from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit_nature.mappers.second_quantization import JordanWignerMapper, ParityMapper
es_problem = ElectronicStructureProblem(driver)
second_q_op = es_problem.second_q_ops()
  +-+- * (0.18093119978423158+0j)
+ +--+ * (-0.18093119978423156+0j)
+ -++- * (-0.18093119978423156+0j)
+ -+-+ * (0.18093119978423156+0j)
+ IIIN * (-0.47189600728114073+0j)
+ IINI * (-1.2563390730032507+0j)
+ IINN * (0.48365053047106565+0j)
+ INII * (-0.47189600728114073+0j)
+ ININ * (0.6985737227320187+0j)
+ INNI * (0.6645817302552972+0j)
+ NIII * (-1.2563390730032507+0j)
+ NIIN * (0.6645817302552972+0j)
+ NINI * (0.675710154803517+0j)
+ NNII * (0.48365053047106565+0j)

If we now transform this Hamiltonian for the given driver defined above we get our qubit operator:

qubit_converter = QubitConverter(mapper=JordanWignerMapper())
qubit_op = qubit_converter.convert(second_q_op[0])
-0.8105479805373268 * IIII
- 0.22575349222402502 * ZIII
+ 0.17218393261915543 * IZII
+ 0.12091263261776641 * ZZII
- 0.22575349222402505 * IIZI
+ 0.17464343068300467 * ZIZI
+ 0.1661454325638243 * IZZI
+ 0.17218393261915546 * IIIZ
+ 0.1661454325638243 * ZIIZ
+ 0.16892753870087926 * IZIZ
+ 0.12091263261776641 * IIZZ
+ 0.04523279994605789 * XXXX
+ 0.04523279994605789 * YYXX
+ 0.04523279994605789 * XXYY
+ 0.04523279994605789 * YYYY

In the minimal (STO-3G) basis set 4 qubits are required. We can reduce the number of qubits by using the Parity mapping, which allows for the removal of 2 qubits by exploiting known symmetries arising from the mapping.

qubit_converter = QubitConverter(mapper = ParityMapper(), two_qubit_reduction=True)
qubit_op = qubit_converter.convert(second_q_op[0], num_particles=es_problem.num_particles)
-1.0523732457728592 * II
+ (-0.3979374248431804+1.3877787807814457e-17j) * ZI
+ 0.3979374248431804 * IZ
- 0.011280104256235296 * ZZ
+ (0.18093119978423153+3.469446951953614e-18j) * XX

This time only 2 qubits are needed.

Now that the Hamiltonian is ready, it can be used in a quantum algorithm to find information about the electronic structure of the corresponding molecule. Check out our tutorials on Ground State Calculation and Excited States Calculation to learn more about how to do that in Qiskit!


Version Information

Qiskit SoftwareVersion
IBM Q ProviderNone
System information
Python3.8.10 (default, May 4 2021, 07:16:51) [GCC 9.3.0]
Memory (Gb)6.791343688964844
Tue May 11 14:50:36 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

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.