Japanese
言語
English
Japanese
German
Korean
Portuguese, Brazilian
French
Table of Contents
Shortcuts

注釈

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

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

量子コンピューターによる熱力学的観測量の計算

[1]:

# imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from functools import partial

from qiskit.aqua import QuantumInstance
from qiskit import BasicAer
from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE, IQPE
from qiskit.aqua.components.optimizers import SLSQP
from qiskit.chemistry.components.initial_states import HartreeFock
from qiskit.circuit.library import ExcitationPreserving

from qiskit.circuit.library import ExcitationPreserving

from qiskit import BasicAer
from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE
from qiskit.aqua.components.optimizers import SLSQP
from qiskit.chemistry.components.initial_states import HartreeFock
from qiskit.chemistry.components.variational_forms import UCCSD
from qiskit.chemistry.drivers import PySCFDriver, UnitsType
from qiskit.chemistry.core import TransformationType, QubitMappingType
from qiskit.aqua.algorithms import VQAlgorithm, VQE, MinimumEigensolver
from qiskit.chemistry.transformations import FermionicTransformation
from qiskit.chemistry.drivers import PySCFDriver
from qiskit.chemistry.algorithms.ground_state_solvers import GroundStateEigensolver
from qiskit.chemistry.algorithms.pes_samplers.extrapolator import *
from qiskit.chemistry.algorithms.pes_samplers.bopes_sampler import BOPESSampler
from qiskit.chemistry.drivers.molecule import Molecule
import qiskit.chemistry.constants as const
from qiskit.chemistry.algorithms.pes_samplers.potentials.energy_surface_spline import EnergySurface1DSpline
from thermodynamics_utils.thermodynamics import constant_volume_heat_capacity
from thermodynamics_utils.vibronic_structure_fd import VibronicStructure1DFD
from thermodynamics_utils.partition_function import DiatomicPartitionFunction
from thermodynamics_utils.thermodynamics import Thermodynamics

import warnings
warnings.simplefilter('ignore', np.RankWarning)

A preliminary draft with more information related to this tutorial can be found in preprint: Stober et al, arXiv 2003.02303 (2020)

Calculation of the Born Oppenheimer Potential Energy Surface (BOPES)

熱力学的な観測量を計算するために、波動関数と電荷密度、つまり特定の原子核配置におけるエネルギーを求める、一点計算から始めます。ここでは例として、単に結合長の関数としての電子のエネルギーである、水素分子の Born-Oppenheimer ポテンシャル・エネルギー面を計算します。

[2]:

ft = FermionicTransformation()
driver = PySCFDriver()
solver = VQE(quantum_instance=
             QuantumInstance(backend=BasicAer.get_backend('statevector_simulator')))
me_gss = GroundStateEigensolver(ft, solver)

[3]:

stretch1 = partial(Molecule.absolute_distance, atom_pair=(1, 0))
mol = Molecule(geometry=[('H', [0., 0., 0.]),
                        ('H', [0., 0., 0.2])],
                       degrees_of_freedom=[stretch1],
                       masses=[1.6735328E-27, 1.6735328E-27]
                       )


# pass molecule to PSYCF driver
driver = PySCFDriver(molecule=mol)

[4]:

# BOPES sampler testing
bs = BOPESSampler(gss=me_gss,
                  bootstrap=True)
points = np.linspace(0.45, 5, 50)
res = bs.sample(driver,points)

[5]:

energies = []
bs_res_full = res.raw_results
for point in points:
    energy = bs_res_full[point]['computed_energies'][0] + bs_res_full[point]['nuclear_repulsion_energy']
    energies.append(energy)

[6]:

fig = plt.figure()
plt.plot(points,energies)
plt.title('Dissociation profile')
plt.xlabel('Interatomic distance')
plt.ylabel('Energy')

[6]:

Text(0, 0.5, 'Energy')
../../_images/tutorials_chemistry_06_calculating_thermodynamic_observables_9_1.png 
[7]:

energy_surface = EnergySurface1DSpline()

xdata = res.points
ydata = res.energies
energy_surface.fit(xdata=xdata, ydata=ydata)

[8]:

plt.plot(xdata, ydata, 'kx')
x = np.arange(min(xdata) - 0.25, max(xdata) + 0.25, 0.05)
plt.plot(x, energy_surface.eval(x), 'r-')
plt.xlabel(r'distance, $\AA$')
plt.ylabel('energy, Hartree')
dist = max(ydata) - min(ydata)
plt.ylim(min(ydata) - 0.1 * dist, max(ydata) + 0.1 * dist)

[8]:

(-1.1576826065021262, -0.9127521770514491)
../../_images/tutorials_chemistry_06_calculating_thermodynamic_observables_11_1.png

分子の振動エネルギー・レベルの計算

Born-Oppeheimer 近似では、分子のハミルトニアンで原子核の振動を無視するため、このように近似された量子力学的な基底エネルギー計算によって求められたエネルギーは、電子的なエネルギーのみが含まれています。絶対零度においても原子核は振動しているため、分子の真の零点エネルギーを得るためには補正が必要です。この補正は、ゼロ点振動エネルギー (zero-point vibrational energy, ZPE) と呼ばれ、原子核の振動モードの寄与をあわせて計算されます。したがって、熱力学的な観測量を計算するための次のステップは、振動エネルギー・レベルを決定することです。これは、平衡結合長に近い単一点エネルギーに基づいてヘッセ行列を作成することにより行うことができます。ヘッセ行列の固有値は、振動エネルギー・レベルと零点振動エネルギーを決定するために用いることができます。
\begin{equation} {\rm ZPE} = \frac{1}{2}\, \sum_i ^M \nu_i \, , \end{equation} ここで \(\nu_i\) は振動周波数、\(N\) は粒子数を表しており、非線形な分子に対しては \(M = 3N − 6\) 、線形な分子に対しては \(M = 3N − 5\) です。

ここで、一次元スプライン・ポテンシャルを用いてエネルギー面「全体」のフィッティングを行い、それを分子の振動エネルギー・レベルの評価に使用します。

[9]:

vibronic_structure = VibronicStructure1DFD(mol, energy_surface)

plt.plot(xdata, ydata, 'kx')
x = np.arange(min(xdata) - 0.25, max(xdata) + 0.25, 0.05)
plt.plot(x, energy_surface.eval(x), 'r-')
plt.xlabel(r'distance, $\AA$')
plt.ylabel('energy, Hartree')
dist = max(ydata) - min(ydata)
plt.ylim(min(ydata) - 0.1 * dist, max(ydata) + 0.1 * dist)
for N in range(15):
    on = np.ones(x.shape)
    on *= (energy_surface.eval(energy_surface.get_equilibrium_geometry()) +
           vibronic_structure.vibrational_energy_level(N))
    plt.plot(x, on, 'g:')
on = np.ones(x.shape)

plt.show()
../../_images/tutorials_chemistry_06_calculating_thermodynamic_observables_15_0.png

熱容量を計算するための分配関数の作成

分子の分配関数は、並進・回転・振動・電子・原子核の自由度からの寄与の積です。振動の周波数がわかっているので、分子全体の分配関数 \begin{equation} q_{\rm vibration} = \prod_{i=1} ^M \frac{\exp\,(-\Theta_{\nu_i}/2T)}{1-\exp\,(-\Theta_{\nu_i}/2T} \, . \end{equation} を計算するための振動分配関数 \(q_{\rm vibration}\) を求めることができます。ここで \(\Theta_{\nu_i}= h\nu_i/k_B\) は温度、\(k_B\) はボルツマン定数です。

一点エネルギー計算とその結果得られる分配関数は、分子の (定積あるいは定圧における) 熱容量を計算することができます。定積熱容量は、例えば次のようになります。

\begin{equation} C_v = \left.\frac{\partial U}{\partial T}\right|_{N,V}\, , \qquad {\rm with} \quad U=k_B T^2 \left.\frac{\partial {\rm ln} Q}{\partial T}\right|_{N,V} . \end{equation}

\(U\) は内部エネルギー、 \(V\) は体積、\(Q\) は分配関数です。

ここでは、上で定義した定圧熱容量のような性質を計算するための熱力学オブジェクトを作成するという、分配関数の最も簡単な使い方を説明します。

[10]:

Q = DiatomicPartitionFunction(mol, energy_surface, vibronic_structure)

P = 101350  # Pa
temps = np.arange(10, 1050, 5)  # K

mol.spins = [1/2,1/2]

td = Thermodynamics(Q, pressure = 101350)
td.set_pressure(101350)
temps = np.arange(10, 1500, 5)
ymin = 5
ymax = 11

plt.plot(temps,
         td.constant_pressure_heat_capacity(temps) / const.CAL_TO_J)
plt.xlim(0, 1025)
plt.ylim(ymin, ymax)
plt.xlabel('Temperature, K')
plt.ylabel('Cp, cal mol$^{-1}$ K$^{-1}$')

plt.show()
../../_images/tutorials_chemistry_06_calculating_thermodynamic_observables_19_0.png

ここでは、分配関数の特定の成分 (回転部分) にアクセスする方法を示します。これは H2 の場合ではさらに、パラ水素の成分とオルト水素の成分に分けることができます。

[11]:

eq = Q.get_partition(part="rot", split="eq")
para = Q.get_partition(part="rot", split="para")
ortho = Q.get_partition(part="rot", split="ortho")

どのように ‘thermodynamics’ モジュールの関数を直接呼び出し、分配関数 (今回の場合はその回転部分) に呼び出し可能オブジェクトを与えるかを示すために、(回転部分の) 定積熱容量をプロットします。このプロットでは、普遍的な気体定数 R (アボガドロ数にボルツマン定数を掛けたもの) で割ることで正規化を行い、文献の実験データと比較するために、十字印を用いていることに注意してください。

[15]:

# REFERENCE DATA from literature
df_brink_T = [80.913535,135.240157,176.633783,219.808499,246.226899]
df_brink_Cv = [0.118605,0.469925,0.711510,0.833597,0.895701]

df_eucken_T = [25.120525, 30.162485, 36.048121, 41.920364, 56.195875, 62.484934, 72.148692, 73.805910, 73.804236, 92.214423,180.031917,230.300866]
df_eucken_Cv  = [0.012287,0.012354,0.008448,0.020478,0.032620,0.048640,0.048768,0.076678,0.078670,0.170548,0.667731,0.847681]

df_gia_T = [190.919338,195.951254,202.652107,204.292585,209.322828,225.300754,234.514217,243.747768]
df_gia_Cv = [0.711700,0.723719,0.749704,0.797535,0.811546,0.797814,0.833793,0.845868]

df_parting_T = [80.101665, 86.358919,185.914204,239.927797]
df_parting_Cv = [0.084730,0.138598,0.667809,0.891634]

df_ce_T = [80.669344,135.550569,145.464190,165.301153,182.144856,203.372528,237.993108,268.696642,294.095771,308.872014]
df_ce_Cv = [0.103048,0.467344,0.541364,0.647315,0.714078,0.798258,0.891147,0.944848,0.966618,0.985486]

[13]:

HeatCapacity = constant_volume_heat_capacity

R = const.N_A * const.KB_J_PER_K
plt.plot(temps,
         HeatCapacity(eq, temps) / R, '-k',
         label='Cv_rot Equilibrium')
plt.plot(temps, HeatCapacity(para, temps) / R, '-b',
         label='Cv_rot Para')
plt.plot(temps, HeatCapacity(ortho, temps) / R, '-r',
         label='Cv_rot Ortho')
plt.plot(temps, 0.25 * HeatCapacity(para, temps) / R
         + 0.75 * HeatCapacity(ortho, temps) / R, '-g',
         label='Cv_rot 1:3 para:ortho')
plt.plot(df_brink_T, df_brink_Cv, '+g')
plt.plot(df_eucken_T, df_eucken_Cv, '+g')
plt.plot(df_gia_T, df_gia_Cv, '+g')
plt.plot(df_parting_T, df_parting_Cv, '+g')
plt.plot(df_ce_T, df_ce_Cv, '+g', label = 'experimental data')
plt.legend(loc='upper right', framealpha=100, frameon=False)
plt.xlim(10, 400)
plt.ylim(-0.1, 2.8)
plt.xlabel('Temperature, K')
plt.ylabel('Cv (rotational)/R')
plt.tight_layout()
plt.show()
../../_images/tutorials_chemistry_06_calculating_thermodynamic_observables_24_0.png 
[14]:

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit SoftwareVersion
Qiskit0.23.0
Terra0.16.0
Aer0.7.0
Ignis0.5.0
Aqua0.8.0
IBM Q Provider0.11.0
System information
Python3.6.12 |Anaconda, Inc.| (default, Sep 8 2020, 17:50:39) [GCC Clang 10.0.0 ]
OSDarwin
CPUs6
Memory (Gb)16.0
Wed Oct 21 18:45:59 2020 CEST

This code is a part of Qiskit

© Copyright IBM 2017, 2020.

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.