Nota

Esta página foi gerada a partir do tutorials/chemistry/06_calculating_thermodynamic_observables.ipynb.

Execute interativamente no IBM Quantum lab.

Calculando Observáveis Termodinâmicos com um computador quântico

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

Um rascunho preliminar com mais informações relacionadas a este tutorial pode ser encontrado no preprint: Stober et al, arXiv 2003.02303 (2020)

Cálculo de Born Oppenheimer Potential Energy Surface (BOPES)

Para computar observáveis termodinâmicos começamos com o cálculo de energia de um único ponto que calcula a função de onda e a densidade de carga e, portanto, a energia de um determinado arranjo de núcleos. Aqui computamos a superfície de energia potencial de Born-Oppenheimer de uma molécula de hidrogênio, como exemplo, que é simplesmente a energia eletrônica em função do comprimento da ligação.

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

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)

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

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)

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

Text(0, 0.5, 'Energy')

energy_surface = EnergySurface1DSpline()

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

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)

(-1.1576826065021262, -0.9127521770514491)

Cálculo dos níveis de Energia Vibrônica molecular

A aproximação de Born-Oppenheimer remove vibrações internucleares do Hamiltoniano molecular e a energia computada a partir de cálculos de energia do estado fundamental da mecânica quântica utilizando esta aproximação contém apenas a energia eletrônica. Uma vez que mesmo no zero absoluto as vibrações internucleares ainda ocorrem, é necessária uma correção para obter a verdadeira energia de temperatura zero de uma molécula. Essa correção é chamada de energia vibracional do ponto zero (ZPE), que é computada somando a contribuição de modos vibracionais internucleares. Por isso, o próximo passo na computação dos observáveis termodinâmicos é determinar os níveis de energia vibracional. Isso pode ser feito através da construção da matriz Hessiana baseda nas energias de ponto único computadas próximas ao comprimento da ligação em equilíbrio. Os autovalores da matriz Hessiana podem, então, ser usados para determinar os níveis de energia vibracional e a energia vibracional de ponto zero

\begin{equation} {\rm ZPE} = \frac{1}{2}\, \sum_i ^M \nu_i \, , \end{equation} com \(\nu_i\) sendo as frequências vibracionais, \(M = 3N -- 6\) ou \(M = 3N -- 5\) para moléculas não-lineares ou lineares, respectivamente e \(N\) é o número de partículas.

Aqui ajustamos uma superfície de energia “completa” usando um potencial de spline 1D e utilizando-o para avaliar os níveis de energia vibracional molecular.

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

Crie uma função de partição para o cálculo da capacidade térmica

A função de partição para uma molécula é o produto das contribuições dos graus de liberdade translacional, rotacional, vibracional, electrônico e nuclear. Tendo as frequências vibracionais, agora podemos obter a função de partição vibracional \(q_{\rm vibration}\) para calcular a função de partição molecular inteira \begin{equation} q_{\rm vibration} = \prod_{i=1} ^M \frac{\exp\,(-\Theta_{\nu_i}/2T)}{1-\exp\,(-\Theta_{\nu_i}/2T} \, . \end{equation} Aqui \(\Theta_{\nu_i}= h\nu_i/k_B\), \(T\) é a temperatura e \(k_B\) é a constante de Boltzman.

Os cálculos da energia de ponto único e a função de partição resultante podem ser usados para calcular a capacidade térmica (a volume constante ou à pressão constante) das moléculas. A capacidade térmica a volume constante, por exemplo, é dada por

\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\) é a energia interna, \(V\) é o volume e \(Q\) é a função de partição.

Aqui ilustramos o uso mais simples da função de partição, a saber, criar um objeto Thermodynamics para computar propriedades como a capacidade térmica à pressão constante definida acima.

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

Aqui demonstramos como acessar componentes particulares (a parte rotacional) da função de partição, que no caso do H2 podemos dividir ainda mais nas componentes para-hidrogênio e orto-hidrogênio.

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

Vamos agora visualizar a capacidade térmica a volume constante (da parte rotacional) demonstrando como podemos chamar diretamente as funções no módulo ‘thermodynamics’, fornecendo um objeto chamável para a função de partição (ou neste caso seu componente rotacional). Observe que no gráfico nós normalizamos-o dividindo pela constante universal dos gases R (número de Avogadro vezes a constante de Boltzman) e usamos cruzes para comparar com dados experimentais encontrados na literatura.

# 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]

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

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.6.12 |Anaconda, Inc.| (default, Sep 8 2020, 17:50:39) [GCC Clang 10.0.0 ]
OS	Darwin
CPUs	6
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.