注釈

当ページは tutorials/chemistry/02_vibronic_structure.ipynb から生成されました。

IBM Quantum lab でインタラクティブに実行します。

振動構造¶

はじめに¶

ある分子のハミルトニアンは、以下の数式で表されます。

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

ボルン・オッペンハイマー近似では、分子の波動関数は各電子の波動関数と振動回転スペクトルの積として表現されています。電子の波動関数は、電子のシュレーディンガー方程式の解に相当します。また、振動回転スペクトルは、様々な構造に対する電子のシュレディンガー方程式の固有値サンプリングによって生成されるポテンシャルエネルギー曲面(PES) での原子核のシュレーディンガー方程式の解に相当します。

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.

こうした系が大きくなると不利になる問題に対して、量子アルゴリズムの適用が望まれています。

原子核のシュレーディンガー方程式は以下で表されます。

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

いわゆるWatsonハミルトニアン（振動回転結合を無視）は、以下で表されます。

\[\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)\]

ここで \(Q_l\) は調和質量加重正規座標です。

\(\mathcal{H}_\text{vib}\) は、量子ハードウェア上で固有関数を計算するために、与えられた:math:N_q 量子ビットの状態に作用する演算子にマッピングする必要があります。電子構造計算では、非相対論的電子ハミルトニアンを第二量子化で表現することによってマッピングを実現することができます。つまり、与えられた(有限の)分子軌道によって生成された反対称化された粒子数表示ベクトル(ONV)の完全なセットとして投影します。同様に第二量子化演算子に調和振動子ハミルトニアンをエンコードするためには、ポテンシャル \(V(Q_1, \ldots, Q_L)\) を :math:`n`体に拡張して以下のように展開します。

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

Watson ハミルトニアンの柔軟な第二量子化形式は、いわゆる n モード表現から得られます。各モード \(l\) が、次のように定義された \(N_l\) 次元基底の \(S_l\) によって記述されているとしましょう。

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

\(n\) モード波動関数は、積基底 \(\mathcal{S} = \otimes_{i=1}^L \mathcal{S}_i\) によって、 Cl のように展開することができます。

\[|\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) \, ,\]

多体基底関数 \(\phi_{k_1}^{(1)}(Q_1) \cdots \phi_{k_L}^{(L)}(Q_L)\) は、いわゆる \(n\) モード第二量子化において、以下のように占有数ベクトル (occupation-number vectors, ONVs) として符号化されています。

\[\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 \, .\]

したがって上のように定義された ONV は、それぞれが独立モードを記述する \(L\) モードに固有なONV の積となります。各モードはただ一つの基底関数によって記述されているため、各モードに固有な ONV の占有数は 1 です。理論的に各モードは、区別可能な準粒子 (固体物理学ではフォノンとして定義されます) として解釈することができます。区別可能性は、 PESが二つのモードの並び替えに関して不変でないという事実から導かれます。今回の場合、これは二つの等しい粒子間にはたらくクーロン相互作用とは異なっています。この観点から、分子は PES 演算子を介して相互作用する、 \(L\) 個の区別できない粒子の集合として解釈することができます。

第二量子化表現に基づいて、モード \(l\) と基底関数 \(k_l\) ごとに、以下のように生成演算子と消滅演算子の対を導入します。

\[\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}\]

ここで

\[\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}\]

第二量子化形式は、ポテンシャルを次のように表現することで得られます。

\[\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}\]

ここで、上に定義した演算子と電子構造の演算子の違いを強調しておきます。まず既に述べたように、ポテンシャルには (原理的には) 三体やそれより高次の多体結合項が含まれており、これにより 6 つ (あるいはそれ以上の) 第二量子化演算子をもった演算子となっています。さらに、与えられたモードの生成演算子と消滅演算子の数が各項で同じであることからもわかるように、ハミルトニアンは各モードの粒子数を保存しています。一方で異なるモードは、相異なるモード \(l\) と \(m\) に属する第二量子化演算子を含んだ、二体 (あるいはそれより高次) の項によって結合されています。

参考文献: Ollitrault, Pauline J., et al., arXiv:2003.12578 (2020).

電子ポテンシャルの計算

異なる原子核配置に対して ESE を解くことで、 PES 関数 \(V(Q_1, \ldots, Q_L)\) を求めることができます。現在では Qiskit によって、 PES を四次力場で近似することが可能です。

\[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\]

このような形式の PES の利点は、非調和力場 ( \(k_{ij}\), \(k_{ijk}\), \(k_{ijkl}\) ) が有限差分法によって求められることです。応答理論に基づく手法 (HF や DFT) によって原子核エネルギーのヘシアンが解析的に求まる手法については、準数値微分によって四次力場を次のように計算することができます。

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

および

\[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}\]

例えばこのような数値的手法は、Gaussian のプログラム・セットに実装されています。

実際には、これは GaussianForceDriver を使用して、Qiskit で行うことができます。

[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',
#                  '',
#                  ''

量子ビット・ハミルトニアンへのマッピング¶

ポテンシャルの近似ができたので、第二量子化でハミルトニアンを書き下す必要があります。そのためには、一体積分 \(\langle\phi_{k_i}| V(Q_i) | \phi_{h_i} \rangle\) や二体積分 \(\langle\phi_{k_i} \phi_{k_j}| V(Q_i,Q_j) | \phi_{h_i} \phi_{h_j} \rangle\)などを計算するためのモーダル基底を選ぶ必要があります。

最も単純な場合、 \(\phi\) 関数は各モードの調和振動子固有関数です。この選択の主な利点は、テイラー展開によって表現された PES の積分が、このような基底で簡単に計算できることです。これらの積分を計算するためのルーチンは Qiskit で実装されています。

その後、ボソン演算子 \(\mathcal{H}_\text{vib}^{SQ}\) を生成し、量子ビット演算子にマッピングしなければなりません。このチュートリアルの最初のセクションで紹介した直接的なマッピングは Qiskitを用いて次のようにできます。

[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)

一つ前のセルでは、調和モーダル基底でハミルトニアンを表現するためにボソン変換を定義しました。ここでは、一つのモードごとに二つのモーダルが対応付けられ、ポテンシャルは二次で打ち切られており、ボソンの量子ビットに対する「直接的な」マッピングとなっています。計算は次のように実行されます。

[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))

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

モードごとに異なるモーダルをもたせるためには、以下のようにします。

[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)

これでハミルトニアンの準備ができたので、対応する分子の振動構造に関する情報を得るために、量子アルゴリズムでこのハミルトニアンを利用することができます。 Qiskit での計算方法については、基底状態計算と励起状態計算のチュートリアルをご覧ください。

[5]:

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
Qiskit	0.23.0
Terra	0.16.0
Aer	0.7.0
Ignis	0.5.0
Aqua	0.8.0
IBM Q Provider	0.11.0
System information
Python	3.7.4 (default, Aug 13 2019, 15:17:50) [Clang 4.0.1 (tags/RELEASE_401/final)]
OS	Darwin
CPUs	6
Memory (Gb)	16.0
Wed Oct 21 12:22:51 2020 CEST

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.

[ ]: