Note

This page was generated from tutorials/chemistry/01_electronic_structure.ipynb.

Run interactively in the IBM Quantum lab.

Electronic structure¶

Introduction¶

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 independant Schroedinger equation,

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

where

\[\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 Hamtiltonian 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 concreatly 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.

[1]:

from qiskit.chemistry.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')

/home/computertreker/git/qiskit/qiskit-dup/.tox/docs/lib/python3.7/site-packages/pyscf/lib/misc.py:47: 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 https://qiskit.org/documentation/apidoc/qiskit.chemistry.drivers.html

The mapping from fermions to qubits¶

d17b42bd0c16441e83504a3baedf694c

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 Bravyi-Kitaev Super Fast ‘bksf’ (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 an object which contains all the information about any transformation of the fermionic Hamiltonian to the qubits Hamiltonian. In this example we simply ask for the Jordan-Wigner mapping.

[2]:

from qiskit.chemistry.transformations import (FermionicTransformation,
                                              FermionicTransformationType,
                                              FermionicQubitMappingType)

fermionic_transformation = FermionicTransformation(
            transformation=FermionicTransformationType.FULL,
            qubit_mapping=FermionicQubitMappingType.JORDAN_WIGNER,
            two_qubit_reduction=False,
            freeze_core=False)

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

[3]:

qubit_op, _ = fermionic_transformation.transform(driver)
print(qubit_op)
print(fermionic_transformation.molecule_info)

SummedOp([
  -0.8105479805373264 * IIII,
  0.17218393261915527 * IIIZ,
  -0.22575349222402502 * IIZI,
  0.17218393261915532 * IZII,
  -0.22575349222402497 * ZIII,
  0.12091263261776641 * IIZZ,
  0.16892753870087926 * IZIZ,
  0.04523279994605789 * XXYY,
  0.04523279994605789 * YYYY,
  0.04523279994605789 * XXXX,
  0.04523279994605789 * YYXX,
  0.1661454325638243 * ZIIZ,
  0.1661454325638243 * IZZI,
  0.17464343068300467 * ZIZI,
  0.12091263261776641 * ZZII
])
{'num_particles': [1, 1], 'num_orbitals': 4, 'two_qubit_reduction': False, 'z2_symmetries': <qiskit.aqua.operators.legacy.weighted_pauli_operator.Z2Symmetries object at 0x7fedcda35790>}

In the minimal (STO-3G) basis set 4 qubits are required. We could even lower the qubit count by using the Parity mapping which allows to get rid of to qubits by symmetry considerations.

[4]:

fermionic_transformation_2 = FermionicTransformation(
            transformation=FermionicTransformationType.FULL,
            qubit_mapping=FermionicQubitMappingType.PARITY,
            two_qubit_reduction=True,
            freeze_core=False)
qubit_op_2, _ = fermionic_transformation_2.transform(driver)
print(qubit_op_2)

SummedOp([
  -1.0523732457728587 * II,
  0.3979374248431802 * IZ,
  -0.3979374248431802 * ZI,
  -0.011280104256235324 * ZZ,
  0.18093119978423147 * XX
])

This time only 2 qubits are needed.

Another possibility is to use the Particle-Hole tranformation (Physical Review A, 98(2), 022322 (2018)). This shifts the vacuum state to a state lying in the N-particle Fock space. In this representation the HF (reference) state has a null energy and the optimization procedure is more faster.

[5]:

fermionic_transformation_3 = FermionicTransformation(
            transformation=FermionicTransformationType.PARTICLE_HOLE,
            qubit_mapping=FermionicQubitMappingType.JORDAN_WIGNER,
            two_qubit_reduction=False,
            freeze_core=False)
qubit_op_3, _ = fermionic_transformation_3.transform(driver)
print(qubit_op_3)

SummedOp([
  1.026420010665659 * IIII,
  0.17218393261915532 * IIIZ,
  -0.22575349222402522 * IIZI,
  0.17218393261915532 * IZII,
  -0.22575349222402522 * ZIII,
  0.16892753870087926 * IZIZ,
  0.04523279994605789 * YYYY,
  0.04523279994605789 * XXYY,
  0.04523279994605789 * YYXX,
  0.04523279994605789 * XXXX,
  0.12091263261776641 * IIZZ,
  0.1661454325638243 * IZZI,
  0.1661454325638243 * ZIIZ,
  0.17464343068300467 * ZIZI,
  0.12091263261776641 * ZZII
])

The list of available mappings and transformations are

[6]:

print('*Transformations')
for fer_transform in FermionicTransformationType:
    print(fer_transform)

print('\n*Mappings')
for fer_mapping in FermionicQubitMappingType:
    print(fer_mapping)

*Transformations
FermionicTransformationType.FULL
FermionicTransformationType.PARTICLE_HOLE

*Mappings
FermionicQubitMappingType.JORDAN_WIGNER
FermionicQubitMappingType.PARITY
FermionicQubitMappingType.BRAVYI_KITAEV

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!

[7]:

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
Qiskit	0.24.1
Terra	0.16.4
Aer	0.7.6
Ignis	0.5.2
Aqua	0.8.2
IBM Q Provider	0.12.2
System information
Python	3.7.7 (default, Apr 22 2020, 19:15:10) [GCC 9.3.0]
OS	Linux
CPUs	32
Memory (Gb)	125.71903228759766
Tue May 25 15:05:11 2021 EDT

This code is a part of Qiskit

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.

[ ]: