Note

Run interactively in the IBM Quantum lab.

# Vibronic 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.

Within the Born-Oppenheimer approximation, a molecular wave function is factorized as a product of an electronic part, which is the solution of the electronic Schroedinger equation, and a vibro-rotational one, which is the solution of the nuclear Schroedinger equation in the potential energy surface (PES) generated by sampling the eigenvalues of the electronic Schroedinger equation for different geometries.

The nuclear Schroedinger equation is usually solved in two steps, in analogy with its electronic counterpart. A single-particle basis (the basis functions are called, in this case, modals) is obtained either by the harmonic approximation applied to the PES or from a vibrational self-consistent field (VSCF) calculation. Vibrational anharmonic correlations are added a-posteriori with perturbative or variational approaches. The latter include Vibrational Configuration Interaction (VCI) and Vibrational Coupled Cluster (VCC) for highly-accurate anharmonic energies. The main advantage of VCI and VCC over alternative approaches (such as perturbation theories) is that their accuracy can be systematically improved towards the complete basis set limit for a given PES. However, their applicability is limited to small molecules with up to about 10 atoms due to their unfavorable scaling with system size.

To tackle the scaling problem we would like to use quantum algorithms.

The nuclear Schroedinger equation is

$\mathcal{H}_{\text{vib}} |\Psi_{n}\rangle = E_{n} |\Psi_{n}\rangle$

The so-called Watson Hamiltonian (neglecting vibro-rotational coupling terms) is

$\mathcal{H}_\text{vib}(Q_1, \ldots, Q_L) = - \frac{1}{2} \sum_{l=1}^{L} \frac{\partial^2}{\partial Q_l^2} + V(Q_1, \ldots, Q_L)$

where $$Q_l$$ are the harmonic mass-weighted normal coordinates.

$$\mathcal{H}_\text{vib}$$ must be mapped to an operator that acts on the states of a given set of $$N_q$$ qubits in order to calculate its eigenfunctions on quantum hardware. In electronic structure calculations, the mapping is achieved by expressing the non-relativistic electronic Hamiltonian in second quantization, \textit{i.e.} by projecting it onto the complete set of antisymmetrized occupation number vectors (ONV) generated by a given (finite) set of orbitals. To encode the vibrational Hamiltonian in an analogous second quantization operators, we expand the potential $$V(Q_1, \ldots, Q_L)$$ with the $$n$$-body expansion as follows:

$V(Q_1, \ldots, Q_L) = V_0 + \sum_{l=1}^L V^{[l]}(Q_l) + \sum_{l<m}^L V^{[l,m]}(Q_l, Q_m) + \sum_{l<m<n}^L V^{[l,m,n]}(Q_l, Q_m, Q_n) + \ldots$

where $$V_0$$ is the electronic energy of the reference geometry, the one-mode term $$V^{[l]}(Q_l)$$ represents the variation of the PES upon change of the $$l$$-th normal coordinate from the equilibrium position. Similarly, the two-body potential $$V^{[l,m]}(Q_l, Q_m)$$ represents the change in the exact PES upon a simultaneous displacement along the $$l$$-th and $$m$$-th coordinates. Often, including terms up to three-body in the $$L$$-body expansion is sufficient to obtain an accuracy of about 1~cm$$^{-1}$$. We highlight that the many-body expansion of the potential operator defining the Watson Hamiltonian contains arbitrarily high coupling terms. This is a crucial difference compared to the non-relativistic electronic-structure Hamiltonian that contains only pairwise interactions.

A flexible second quantization form of the Watson Hamiltonian is obtained within the so-called n-mode representation. Let us assume that each mode $$l$$ is described by a $$N_l$$-dimensional basis set $$S_l$$ defined as follows:

$\mathcal{S}_l = \{ \phi_1^{(l)} (Q_l) , \ldots , \phi_{N_l}^{(l)} (Q_l) \} \, .$

The $$n$$-mode wave function can be expanded in the product basis $$\mathcal{S} = \otimes_{i=1}^L \mathcal{S}_i$$ as the following CI-like expansion:

$|\Psi\rangle = \sum_{k_1=1}^{N_1} \cdots \sum_{k_L=1}^{N_L} C_{k_1,\ldots,k_L} \phi_{k_1}^{(1)}(Q_1) \cdots \phi_{k_L}^{(L)}(Q_L) \, ,$

The many-body basis function $$\phi_{k_1}^{(1)}(Q_1) \cdots \phi_{k_L}^{(L)}(Q_L)$$ are encoded within the so-called $$n$$-mode second quantization as occupation-number vectors (ONVs) as follows:

$\phi_{k_1}(Q_1) \cdots \phi_{k_L}(Q_L) \equiv |0_1 \cdots 1_{k_1} \cdots 0_{N_1}, 0_1 \cdots 1_{k_2} \cdots 0_{N_2}, \cdots , 0_1 \cdots 1_{k_L} \cdots 0_{N_L}\rangle \, .$

The ONV defined above is, therefore, the product of $$L$$ mode-specific ONVs, each one describing an individual mode. Since each mode is described by one and only one basis function, the occupation of each mode-specific ONV is one. From a theoretical perspective, each mode can be interpreted as a distinguishable quasi-particle (defined as phonons in solid-state physics). Distinguishability arises from the fact that the PES is not invariant by permutation of two modes, also in this case unlike the Coulomb interaction between two equal particles. From this perspective, a molecule can be interpreted as a collection of $$L$$ indistinguishable particles that interact through the PES operator.

Based on this second-quantization representation we introduce a pair of creation and annihilation operators per mode $$l$$ \textit{and} per basis function $$k_l$$ defined as:

\begin{split}\begin{aligned} a_{k_l}^\dagger |\cdots, 0_1 \cdots 0_{k_l} \cdots 0_{N_l}, \cdots\rangle &= | \cdots, 0_1 \cdots 1_{k_l} \cdots 0_{N_l}, \cdots\rangle \\ a_{k_l}^\dagger | \cdots, 0_1 \cdots 1_{k_l} \cdots 0_{N_l}, \cdots\rangle &= 0 \\ a_{k_l} | \cdots, 0_1 \cdots 1_{k_l} \cdots 0_{N_l}, \cdots\rangle &= | \cdots, 0_1 \cdots 0_{k_l} \cdots 0_{N_l}, \cdots\rangle \\ a_{k_l} | \cdots, 0_1 \cdots 0_{k_l} \cdots 0_{N_l}, \cdots\rangle &= 0 \\ \end{aligned}\end{split}

with

\begin{split}\begin{aligned} \left[ a_{k_l}^\dagger, a_{h_m}^\dagger \right] &= 0 \\ \left[ a_{k_l}, a_{h_m} \right] &= 0 \\ \left[ a_{k_l}^\dagger, a_{h_m} \right] &= \delta_{l,m} \, , \delta_{k_l,h_m} \end{aligned}\end{split}

The second quantization form is obtained by expressing the potential as

\begin{split}\begin{aligned} \mathcal{H}_\text{vib}^{SQ} =& \sum_{l=1}^L \sum_{k_l,h_l}^{N_l} \langle \phi_{k_l} | T(Q_l) + V^{[l]}(Q_l) | \phi_{h_l} \rangle a_{k_l}^+ a_{h_l} \\ +& \sum_{l<m}^L \sum_{k_l,h_l}^{N_l} \sum_{k_m,h_m}^{N_m} \langle \phi_{k_l} \phi_{k_m} | V^{[l,m]}(Q_l, Q_m) | \phi_{h_l} \phi_{h_m} \rangle a_{k_l}^+ a_{k_m}^+ a_{h_l} a_{h_m} + \cdots \end{aligned}\end{split}

We highlight here the difference between the operators defined here above and the electronic structure one. First, as we already mentioned, the potential contains (in principle) three- and higher-body coupling terms that lead to strings with six (or more) second-quantization operators. Moreover, the Hamiltonian conserves the number of particles for each mode, as can be seen from the fact that the number of creation and annihilation operators for a given mode is the same in each term. Nevertheless, different modes are coupled by two- (and higher) body terms containing SQ operators belonging to different modes $$l$$ and $$m$$.

Reference: Ollitrault, Pauline J., et al., arXiv:2003.12578 (2020).

Compute the electronic potential

Solving the ESE for different nuclear configurations to obtain the PES function $$V(Q_1, \ldots, Q_L)$$. So far Qiskit gives the possibility to approximate the PES with a quartic force field.

$V(Q_1, \ldots, Q_L) = \frac{1}{2} \sum_{ij} k_{ij} Q_i Q_j + \frac{1}{6} \sum_{ijk} k_{ijk} Q_i Q_j Q_k + \frac{1}{16} \sum_{ijkl} k_{ijkl} Q_i Q_j Q_k Q_l$

The advantage of such form for the PES is that the anharmonic force fields ($$k_{ij}$$, $$k_{ijk}$$, $$k_{ijkl}$$) can be calculated by finite-difference approaches. For methods for which the nuclear energy Hessian can be calculated analytically with response theory-based methods (such as HF and DFT), the quartic force field can be calculated by semi-numerical differentiation as:

$k_{ijk} = \frac{H_{ij}(+\delta Q_k) - H_{ij}(-\delta Q_k)}{2\delta Q_k}$

and

$k_{ijkl} = \frac{H_{ij}(+\delta Q_k+\delta Q_l) - H_{ij}(+\delta Q_k-\delta Q_l) -H_{ij}(-\delta Q_k+\delta Q_l) + H_{ij}(-\delta Q_k+\delta Q_l)} {4\delta Q_k \delta Q_l}$

Such numerical procedure is implemented, for instance, in the Gaussian suite of programs.

In practice this can be done with Qiskit using the GaussianForceDriver.

[1]:

from qiskit.chemistry.drivers import GaussianForcesDriver
# if you ran Gaussian elsewhere and already have the output file
driver = GaussianForcesDriver(logfile='aux_files/CO2_freq_B3LYP_ccpVDZ.log')

# if you want to run the Gaussian job from Qiskit
# driver = GaussianForcesDriver(
#                 ['#p B3LYP/6-31g Freq=(Anharm) Int=Ultrafine SCF=VeryTight',
#                  '',
#                  'CO2 geometry optimization B3LYP/6-31g',
#                  '',
#                  '0 1',
#                  'C  -0.848629  2.067624  0.160992',
#                  'O   0.098816  2.655801 -0.159738',
#                  'O  -1.796073  1.479446  0.481721',
#                  '',
#                  ''

/home/runner/work/qiskit/qiskit/.tox/docs/lib/python3.8/site-packages/qiskit/chemistry/__init__.py:170: DeprecationWarning: The package qiskit.chemistry is deprecated. It was moved/refactored to qiskit_nature (pip install qiskit-nature). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
warn_package('chemistry', 'qiskit_nature', 'qiskit-nature')


## Map to a qubit Hamiltonian¶

Now that we have an approximation for the potential, we need to write the Hamiltonian in second quantization. To this end we need to select a modal basis to calculate the one-body integrals $$\langle\phi_{k_i}| V(Q_i) | \phi_{h_i} \rangle$$, two-body integrals $$\langle\phi_{k_i} \phi_{k_j}| V(Q_i,Q_j) | \phi_{h_i} \phi_{h_j} \rangle$$

In the simplest case, the $$\phi$$ functions are the harmonic-oscillator eigenfunctions for each mode. The main advantage of this choice is that the integrals of a PES expressed as a Taylor expansion are easy to calculate with such basis. A routine for computing these integrals is implemented in Qiskit.

The bosonic operator, $$\mathcal{H}_\text{vib}^{SQ}$$, is then created and must be mapped to a qubit operator. The direct mapping introduced in the first section of this tutorial can be used is Qiskit as follows:

[2]:

from qiskit.chemistry.transformations import (BosonicTransformation,
BosonicTransformationType, # type of modal basis
BosonicQubitMappingType) # type of boson to qubit mapping

bosonic_transformation = BosonicTransformation(qubit_mapping=BosonicQubitMappingType.DIRECT,
transformation_type=BosonicTransformationType.HARMONIC,
basis_size=2,
truncation=2)


In the previous cell we defined a bosonic transformation to express the Hamiltonian in the harmonic modal basis, with 2 modals per mode with the potential truncated at order 2 and the ‘direct’ boson to qubit mapping. The calculation is then ran as:

[3]:

qubit_op, _ = bosonic_transformation.transform(driver)
print(qubit_op.print_details())

print('\nThe total number of modes in this system is {}'.format(bosonic_transformation.num_modes))

/home/runner/work/qiskit/qiskit/.tox/docs/lib/python3.8/site-packages/qiskit/chemistry/bosonic_operator.py:180: DeprecationWarning: The package qiskit.aqua.operators is deprecated. It was moved/refactored to qiskit.opflow (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
qubit_op = WeightedPauliOperator([])

IIIIIIII        (5077.2365606250005+0j)
IIIIIIIZ        (-601.9000340625+0j)
IIIIIIZI        (-1810.6538965625004+0j)
IIIIIZII        (-342.04261812500005+0j)
IIIIXXII        (-84.41897986727395+0j)
IIIIYYII        (-84.41897986727395+0j)
IIIIZIII        (-1027.4430731250004+0j)
IIIZIIII        (-155.20426031250003+0j)
IIZIIIII        (-467.50585156250014+0j)
IZIIIIII        (-157.22939093750003+0j)
ZIIIIIII        (-472.28961593750006+0j)
IIIIIZIZ        (1.2599134375000005+0j)
IIIIIZZI        (3.779740312500002+0j)
IIIIXXIZ        (22.272716326623776+0j)
IIIIYYIZ        (22.272716326623776+0j)
IIIIXXZI        (66.81814897987134+0j)
IIIIYYZI        (66.81814897987134+0j)
IIIIZIIZ        (3.779740312500002+0j)
IIIIZIZI        (11.339220937500007+0j)
IIIZIIIZ        (-1.5962606250000009+0j)
IIIZIIZI        (-4.788781875000002+0j)
IIIZIZII        (-0.6414307812500002+0j)
IIIZXXII        (-5.411241592930711+0j)
IIIZYYII        (-5.411241592930711+0j)
IIIZZIII        (-1.9242923437500008+0j)
IIXXIIXX        (-0.5021409375000001+0j)
IIYYIIXX        (-0.5021409375000001+0j)
IIXXIIYY        (-0.5021409375000001+0j)
IIYYIIYY        (-0.5021409375000001+0j)
IIZIIIIZ        (-4.788781875000002+0j)
IIZIIIZI        (-14.366345625000008+0j)
IIZIIZII        (-1.924292343750001+0j)
IIZIXXII        (-16.233724778792137+0j)
IIZIYYII        (-16.233724778792137+0j)
IIZIZIII        (-5.772877031250004+0j)
IZIIIIIZ        (-1.1489604687500004+0j)
IZIIIIZI        (-3.4468814062500015+0j)
IZIIIZII        (-0.42099484375000024+0j)
IZIIXXII        (-1.6031887335286772+0j)
IZIIYYII        (-1.6031887335286772+0j)
IZIIZIII        (-1.2629845312500008+0j)
IZIZIIII        (-0.13775546875000005+0j)
IZZIIIII        (-0.4132664062500002+0j)
XXIIIIXX        (0.8980418750000001+0j)
YYIIIIXX        (0.8980418750000001+0j)
XXIIIIYY        (0.8980418750000001+0j)
YYIIIIYY        (0.8980418750000001+0j)
XXXXIIII        (1.986637812500001+0j)
YYXXIIII        (1.986637812500001+0j)
XXYYIIII        (1.986637812500001+0j)
YYYYIIII        (1.986637812500001+0j)
ZIIIIIIZ        (-3.4468814062500015+0j)
ZIIIIIZI        (-10.340644218750004+0j)
ZIIIIZII        (-1.2629845312500005+0j)
ZIIIXXII        (-4.809566200586032+0j)
ZIIIYYII        (-4.809566200586032+0j)
ZIIIZIII        (-3.7889535937500023+0j)
ZIIZIIII        (-0.4132664062500002+0j)
ZIZIIIII        (-1.2397992187500007+0j)

The total number of modes in this system is 4


To have a different number of modals per mode:

[4]:

bosonic_transformation_2 = BosonicTransformation(qubit_mapping=BosonicQubitMappingType.DIRECT,
transformation_type=BosonicTransformationType.HARMONIC,
basis_size=[0,0,4,2],
truncation=2)
qubit_op_2, _ = bosonic_transformation_2.transform(driver)
print(qubit_op_2.print_details())

IIIIII  (3309.8257831250007+0j)
IIIIIZ  (-164.15502593750006+0j)
IIIXIX  (1.8979435436158592+0j)
IIIYIY  (1.8979435436158592+0j)
IIIIZI  (-494.3581484375001+0j)
IIXIXI  (6.378706032974414+0j)
IIYIYI  (6.378706032974414+0j)
IIIZII  (-828.3474121875001+0j)
IIZIII  (-1166.1228171875002+0j)
IZIIII  (-161.85614656250007+0j)
ZIIIII  (-486.16988281250013+0j)
IZIIIZ  (-0.13775546875000005+0j)
IZIXIX  (0.19481565219731314+0j)
IZIYIY  (0.19481565219731314+0j)
IZIIZI  (-0.4132664062500002+0j)
IZXIXI  (0.3374306077154137+0j)
IZYIYI  (0.3374306077154137+0j)
IZIZII  (-0.6887773437500002+0j)
IZZIII  (-0.9642882812500004+0j)
XXIIXX  (1.986637812500001+0j)
YYIIXX  (1.986637812500001+0j)
XXIIYY  (1.986637812500001+0j)
YYIIYY  (1.986637812500001+0j)
XXXIIX  (0.8288402231450827+0j)
YYXIIX  (0.8288402231450827+0j)
XXYIIY  (0.8288402231450827+0j)
YYYIIY  (0.8288402231450827+0j)
XXIXXI  (4.245123515804727+0j)
YYIXXI  (4.245123515804727+0j)
XXIYYI  (4.245123515804727+0j)
YYIYYI  (4.245123515804727+0j)
XXXXII  (6.957428881323849+0j)
YYXXII  (6.957428881323849+0j)
XXYYII  (6.957428881323849+0j)
YYYYII  (6.957428881323849+0j)
ZIIIIZ  (-0.4132664062500002+0j)
ZIIXIX  (0.5844469565919395+0j)
ZIIYIY  (0.5844469565919395+0j)
ZIIIZI  (-1.2397992187500007+0j)
ZIXIXI  (1.0122918231462412+0j)
ZIYIYI  (1.0122918231462412+0j)
ZIIZII  (-2.066332031250001+0j)
ZIZIII  (-2.8928648437500013+0j)



Now that the Hamiltonian is ready, it can be used in a quantum algorithm to find information about the vibronic 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!

[5]:

import qiskit.tools.jupyter
%qiskit_version_table


### Version Information

Qiskit SoftwareVersion
Qiskit0.26.0
Terra0.17.3
Aer0.8.2
Ignis0.6.0
Aqua0.9.1
IBM Q Provider0.13.1
System information
Python3.8.10 (default, May 4 2021, 07:16:51) [GCC 9.3.0]
OSLinux
CPUs2
Memory (Gb)6.791339874267578
Mon May 17 18:57:13 2021 UTC

### This code is a part of Qiskit

[ ]: