6-1. OpenFermion basics
In this chapter, we introduce a method to convert the Hamiltonian of an interacting electron system into a form that is easy for quantum computers to handle using OpenFermion[1], a python library for quantum chemistry calculation. OpenFermion provides connectivity to the open source libraries for quantum chemical calculations Psi4 and PySCF. Even if you do not understand the details of how to use these libraries, you can obtain the Hamiltonian of the electron system that appears in quantum chemical calculations simply by entering the structure of the molecule. Here, we use PySCF.
[ ]:
## If various libraries are not installed, please execute.
## If you run on Google Colaboratory, ignore the following message waring, 'You must restart the runtime in order to use newly installed versions.'
## Resuming runtime will lead to a crash.
!pip install qulacs pyscf openfermion openfermionpyscf
## Run only if you are in Google Colaboratory or in a local environment where Qulacs is not installed
## Qulacs errors will now be output normally.
!pip3 install wurlitzer
%load_ext wurlitzer
[1]:
#Import required libraries
#If you get an error, install openfermion later than v.1.0.0
import numpy as np
import matplotlib.pyplot as plt
from openfermion.chem import MolecularData
from openfermion.transforms import get_fermion_operator, jordan_wigner, bravyi_kitaev
from openfermion.linalg import get_sparse_operator
from openfermion.ops import FermionOperator
from openfermionpyscf import run_pyscf
from pyscf import fci
Calculation of hydrogen molecule
In openfermion, we enter data describing molecules into MolecularData class.
[2]:
#define constants
basis = "sto-3g" #basis set
multiplicity = 1 #spin multiplicity
charge = 0 #total charge for the molecule
distance = 0.65
geometry = [("H",(0,0,0)),("H", (0,0,distance))] #xyz coordinates for atoms
description = str(distance) #description for the psi4 output file
molecule = MolecularData(geometry, basis, multiplicity, charge, description)
Variable description
The meaning of the variables appearing in the above code is explained below.
basis: basis function
Set basis functions for representing molecular orbitals. There are various basis sets such as sto-3g, 6-31G.
sto-3g (Slater Type Orbital - 3 gaussian) used here is a basis function that approximates Slater type orbital with 3 gaussians.
Slater type orbital is an orbital that resembles the solution of a hydrogen atom. As a radial function, we use
In angular direction, we use spherical harmonic function \(Y_{lm}(\theta,\phi)\).
sto-3g uses a function which approximates this radial wavefunction \(R_{nl}(r)\) with three gaussians.
multiplicity: spin multiplicity
Since electrons have a spin of 1/2, the spin multiplicity is 2 when one electron exists in isolation. However, in the case of a hydrogen molecule, electrons form a singlet in the ground state, and the total spin is 0. Since spin 0 has only one state, the spin multiplicity is set to 1 in this case.
charge: total charge
Input the total charge. When considering an ion, it becomes + or −.
geometry: nuclear configuration
Specify atomic species and their coordinates in x,y,z.
description
The output results calculated by pyscf are saved in the directory where the openfermion library is saved. This is a variable that determines the name of the file.
Calculatin with PySCF
Input the MolecularData set above to the function run_pyscf and perform quantum chemical calculations with PySCF. It should finish in a few seconds.
[3]:
molecule = run_pyscf(molecule,run_scf=1,run_fci=1)
HF & Full-CI energy
Let’s take a look at the Hartree-Fock energy and Full-CI energy (=exact basis energy) calculated by PySCF. (1 Hartree = 27.2116 eV)
[4]:
print("HF energy: {} (Hartree)".format(molecule.hf_energy))
print("FCI energy: {} (Hartree)".format(molecule.fci_energy))
HF energy: -1.1129965456691684 (Hartree)
FCI energy: -1.1299047843229135 (Hartree)
One-Electron integral \(h_{ij}\)・Two-Electron integral \(h_{ijkl}\)
Quantities such as the one-electron and two-electron integrals are also stored in the MolecularData class.
[5]:
print(molecule.one_body_integrals)
[[-1.30950987e+00 1.98461056e-17]
[ 8.14840166e-17 -4.10026381e-01]]
[6]:
print(molecule.two_body_integrals)
[[[[ 6.91904405e-01 -4.16333634e-17]
[-2.77555756e-17 1.76318452e-01]]
[[-2.77555756e-17 1.76318452e-01]
[ 6.79683914e-01 0.00000000e+00]]]
[[[-4.16333634e-17 6.79683914e-01]
[ 1.76318452e-01 8.32667268e-17]]
[[ 1.76318452e-01 8.32667268e-17]
[ 0.00000000e+00 7.14671111e-01]]]]
Hamiltonian in second quantized form
openfermion computes the Hamiltonian of the second quantized form from these integrals.
(For the second quantization, see e.g. reference [2]). The Hamiltonian can be calculated by calling the get_molecular_hamiltonian method.
For example, (3,1)represents \(c_3^\dagger\), and (1,0) represents \(c_1\) .
[7]:
print(molecule.get_molecular_hamiltonian())
() 0.8141187860307693
((0, 1), (0, 0)) -1.309509868464871
((1, 1), (1, 0)) -1.309509868464871
((2, 1), (2, 0)) -0.4100263808117837
((3, 1), (3, 0)) -0.4100263808117837
((0, 1), (0, 1), (0, 0), (0, 0)) 0.3459522026149022
((0, 1), (0, 1), (2, 0), (2, 0)) 0.0881592258051036
((0, 1), (1, 1), (1, 0), (0, 0)) 0.3459522026149022
((0, 1), (1, 1), (3, 0), (2, 0)) 0.0881592258051036
((0, 1), (2, 1), (0, 0), (2, 0)) 0.0881592258051036
((0, 1), (2, 1), (2, 0), (0, 0)) 0.3398419569652305
((0, 1), (3, 1), (1, 0), (2, 0)) 0.0881592258051036
((0, 1), (3, 1), (3, 0), (0, 0)) 0.3398419569652305
((1, 1), (0, 1), (0, 0), (1, 0)) 0.3459522026149022
((1, 1), (0, 1), (2, 0), (3, 0)) 0.0881592258051036
((1, 1), (1, 1), (1, 0), (1, 0)) 0.3459522026149022
((1, 1), (1, 1), (3, 0), (3, 0)) 0.0881592258051036
((1, 1), (2, 1), (0, 0), (3, 0)) 0.0881592258051036
((1, 1), (2, 1), (2, 0), (1, 0)) 0.3398419569652305
((1, 1), (3, 1), (1, 0), (3, 0)) 0.0881592258051036
((1, 1), (3, 1), (3, 0), (1, 0)) 0.3398419569652305
((2, 1), (0, 1), (0, 0), (2, 0)) 0.33984195696523034
((2, 1), (0, 1), (2, 0), (0, 0)) 0.0881592258051036
((2, 1), (1, 1), (1, 0), (2, 0)) 0.33984195696523034
((2, 1), (1, 1), (3, 0), (0, 0)) 0.0881592258051036
((2, 1), (2, 1), (0, 0), (0, 0)) 0.0881592258051036
((2, 1), (2, 1), (2, 0), (2, 0)) 0.35733555551906837
((2, 1), (3, 1), (1, 0), (0, 0)) 0.0881592258051036
((2, 1), (3, 1), (3, 0), (2, 0)) 0.35733555551906837
((3, 1), (0, 1), (0, 0), (3, 0)) 0.33984195696523034
((3, 1), (0, 1), (2, 0), (1, 0)) 0.0881592258051036
((3, 1), (1, 1), (1, 0), (3, 0)) 0.33984195696523034
((3, 1), (1, 1), (3, 0), (1, 0)) 0.0881592258051036
((3, 1), (2, 1), (0, 0), (1, 0)) 0.0881592258051036
((3, 1), (2, 1), (2, 0), (3, 0)) 0.35733555551906837
((3, 1), (3, 1), (1, 0), (1, 0)) 0.0881592258051036
((3, 1), (3, 1), (3, 0), (3, 0)) 0.35733555551906837
Convert to tractable operators for quantum computers
The easiest to handle on a quantum computer are the Pauli operators \(I, X, Y, Z\) and their tensor products. Therefore, in order to treat the Hamiltonian of ordinary electrons with a quantum computer, convert the Hamiltonian of the second quantization form
to
Various transformation methods have been proposed, but here we use the simplest one, called the Jordan-Wigner transformation. In the Jordan-Wigner transformation, the molecular orbital \(i\) corresponds to the \(i\)-th qubit. The state that the molecular orbital is occupied by an electron is represented with \(|1\rangle\), otherwise \(|0\rangle\).
Under such a convention, constructing the Pauli operator so as to satisfy the anticommutation relation of the following creation-annihilation operator of fermion,
we obtain
For transformation methods other than the Jordan-Wigner transformation, see [2][3].
Jordan-Wigner transformation is implemented in openfermion. By passing FermionOperator to jordan_wigner function, it will return QubitOperator corresponding to the Jordan-Wigner transformation of that operator. Below, the FermionOperator is created from the MolecularData of the hydrogen molecule created above, and the Jordan-Wigner transform is used to convert the Hamiltonian of the hydrogen molecule into a form that can be easily handled by a quantum computer.
[8]:
jw_hamiltonian = jordan_wigner(get_fermion_operator(molecule.get_molecular_hamiltonian()))
print(jw_hamiltonian)
(0.0377511039464572+0j) [] +
(-0.0440796129025518+0j) [X0 X1 Y2 Y3] +
(0.0440796129025518+0j) [X0 Y1 Y2 X3] +
(0.0440796129025518+0j) [Y0 X1 X2 Y3] +
(-0.0440796129025518+0j) [Y0 Y1 X2 X3] +
(0.18601648886230573+0j) [Z0] +
(0.1729761013074511+0j) [Z0 Z1] +
(0.12584136558006342+0j) [Z0 Z2] +
(0.16992097848261523+0j) [Z0 Z3] +
(0.18601648886230576+0j) [Z1] +
(0.16992097848261523+0j) [Z1 Z2] +
(0.12584136558006342+0j) [Z1 Z3] +
(-0.26941693141632095+0j) [Z2] +
(0.17866777775953419+0j) [Z2 Z3] +
(-0.26941693141632095+0j) [Z3]
Let’s calculate the Hartree-Fock (HF) energy from this Hamiltonian. In the Jordan-Wigner transform, the qubit \(\left|0\right\rangle, \left|1\right\rangle\) and the orbital occupancy numbers correspond one-to-one. Therefore, in order to calculate the HF energy, calculate the expected value for \(\left|1100\right\rangle\). Because theare 2 electrons, the state is filled with 1 twice from the bottom.
[9]:
# function for computing the tensor product
def kron_N(*ops):
tmp = ops[0]
for op in ops[1:]:
tmp = np.kron(tmp,op)
return tmp
bra0 = np.array([[1,0]])
bra1 = np.array([[0,1]])
HFbra = kron_N(bra1, bra1, bra0, bra0)
HFket = HFbra.T
print(HFbra)
jw_matrix = get_sparse_operator(jw_hamiltonian)
print(np.real(HFbra.dot(jw_matrix.dot(HFket))), molecule.hf_energy)
[[0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0]]
[[-1.11299655]] -1.1129965456691684
It can be confirmed that it almost matches the calculation of pyscf.
Now let’s diagonalize the Hamiltonian and see that the result agrees with the Full-CI (exact solution) energy.
[10]:
eigenenergies, eigenvecs = np.linalg.eigh(jw_matrix.toarray())
print(eigenenergies[0], molecule.fci_energy)
-1.129904784322913 -1.1299047843229135
It can be confirmed that this also almost matches. The ground state wavefunction \(\left|\psi_g\right\rangle\) is
[11]:
print(eigenvecs[:,0])
[ 0. +0.j 0. +0.j 0. +0.j 0.09545811+0.j
0. +0.j 0. +0.j 0. +0.j 0. +0.j
0. +0.j 0. +0.j 0. +0.j 0. +0.j
-0.99543345+0.j 0. +0.j 0. +0.j 0. +0.j]
\(\left|\psi_g \right\rangle \approx - 0.995\left|1100\right\rangle + 0.0955 \left|0011\right\rangle\). It turns out that the Hartree-Fock solution \(\left|1100\right\rangle\) is slightly combined with the two-electron excited state \(\left|0011\right\rangle\).