Nota

Esta página foi gerada a partir do tutorials/chemistry/02_vibronic_structure.ipynb.

Execute interativamente no IBM Quantum lab.

Estrutura vibrônica¶

Introdução¶

O Hamiltoniano molecular é

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

Dentro da aproximação Born-Oppenheimer, uma função de onda molecular é fatorizada como um produto de uma parte eletrônica, que é a solução da equação de Schroedinger eletrônica, e de uma vibro-rotacional, que é a solução da equação de Schroedinger nuclear na superfície de energia potencial (PES) gerada pela amostragem dos autovalores da equação de Schroedinger eletrônica para diferentes geometrias.

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.

Para enfrentar o problema de escala gostaríamos de usar algoritmos quânticos.

A equação de Schroedinger nuclear é

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

O Hamiltoniano de Watson (negligenciando os termos de acoplamento vibro-rotacional) é

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

onde \(Q_l\) são as coordenadas normais ponderadas pela massa harmônica.

\(\mathcal{H}_\text{vib}\) deve ser mapeado para um operador que atua sobre os estados de um determinado conjunto de \(N_q\) qubits para calcular suas autofunções no hardware quântico. Em cálculos da estrutura eletrônica, o mapeamento é alcançado através da expressão do Hamiltoniano não relativístico eletrônico na segunda quantização, \textit{i.e.} projetando-o no conjunto completo de vetores de número de ocupação antisimétricos (ONV) gerados por um dado (finito) conjunto de orbitais. Para codificar o Hamiltoniano vibracional em operadores de segunda quantização análogos, nós expandimos o potencial \(V(Q_1, \ldots, Q_L)\) com a expansão dos \(n\) corpos da seguinte forma:

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

Uma forma flexível de segunda quantização do Hamiltoniano de Watson é obtida na chamada representação de modo n. Vamos supor que cada modo \(l\) é descrito por um conjunto base \(N_l\)-dimensional \(S_l\) definido da seguinte forma:

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

A função de onda de modo \(n\) pode ser expandida na base de produto \(\mathcal{S} = \otimes_{i=1}^L \mathcal{S}_i\) como a seguinte expansão semelhante a CI:

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

A função base de muitos corpos \(\phi_{k_1}^{(1)}(Q_1) \cdots \phi_{k_L}^{(L)}(Q_L)\) é codificada dentro da chamada segunda quantização de modo \(n\) como vetores de números de ocupação (ONVs) da seguinte forma:

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

O ONV definido acima é, portanto, o produto de \(L\) ONVs de modo específico, cada um descrevendo um modo individual. Como cada modo é descrito por uma e apenas uma função base, a ocupação de cada ONV de modo específico é um. De uma perspectiva teórica, cada modo pode ser interpretado como uma quase-partícula distinguível (definidas como fônons na física do estado sólido). A distinguibilidade resulta do fato do PES não ser invariante à permutação de dois modos, também neste caso, ao contrário da interação de Coulomb entre duas partículas iguais. Desta perspectiva, uma molécula pode ser interpretada como uma coleção de \(L\) partículas indistinguíveis que interagem através do operador PES.

Com base nesta representação da segunda quantização, introduzimos um par de operadores de criação e aniquilação por modo \(l\) \textit{and} por função base :math:`k_l definidos como:

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

com

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

A forma da segunda quantização é obtida através da expressão do potencial como

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

Salientamos aqui a diferença entre os operadores definidos aqui acima e o da estrutura electrônica. Em primeiro lugar, como já mencionamos, o potencial contém (em princípio) termos de acoplamento de três e mais corpos que resultam em strings com seis (ou mais) operadores de segunda quantização. Além disso, os Hamiltonianos conservam o número de partículas para cada modo, como pode ser visto pelo fato do número de operadores de criação e de aniquilação para um determinado modo ser igual em cada termo. No entanto, diferentes modos são acoplados por termos de dois (e mais) corpos contendo operadores SQ pertencentes a diferentes modos \(l\) e \(m\).

Referência: Ollitrault, Pauline J., et al., arXiv:2003.12578 (2020).

Calcule o potencial eletrônico

Resolvendo o ESE para diferentes configurações nucleares para obter a função PES \(V(Q_1, \ldots, Q_L)\). Até agora, o Qiskit dá a possibilidade de aproximar o PES com um campo de força quártico.

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

A vantagem de tal forma para o PES é que os campos de força anarmônicos (\(k_{ij}\), \(k_{ijk}\), \(k_{ijkl}\)) podem ser calculados com abordagens de diferenças finitas. Para métodos para os quais a energia nuclear Hessiana pode ser calculada analiticamente com métodos baseados na teoria de resposta (como HF e DFT), o campo de força quártica pode ser calculado pela diferenciação semi-numérica como:

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

e

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

Tal procedimento numérico é implementado, por exemplo, na suíte Gaussiana de programas.

Na prática isso pode ser feito com Qiskit usando o 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',
#                  '',
#                  ''

Mapear para um Hamiltoniano de qubit¶

Agora que temos uma aproximação para o potencial, precisamos escrever o Hamiltoniano na segunda quantização. Para esse fim, precisamos selecionar uma base modal para calcular as integrais de um corpo \(\langle\phi_{k_i}| V(Q_i) | \phi_{h_i} \rangle\), integrais de dois corpos \(\langle\phi_{k_i} \phi_{k_j}| V(Q_i,Q_j) | \phi_{h_i} \phi_{h_j} \rangle\)…

No caso mais simples, as funções \(\phi\) são as autofunções do oscilador harmônico para cada modo. A principal vantagem desta escolha reside no fato das integrais de um PES expresso como uma expansão de Taylor serem fáceis de calcular com essa base. Uma rotina para calcular essas integrais está implementada no Qiskit.

O operador bosônico, \(\mathcal{H}_\text{vib}^{SQ}\), é então criado e deve ser mapeado para um operador de qubit. O mapeamento direto introduzido na primeira seção deste tutorial pode ser usado no Qiskit da seguinte forma:

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

Na célula anterior definimos uma transformação bosónica para expressar o Hamiltoniano na base modal harmônica, com 2 modais por modo com o potencial truncado na ordem 2 e o mapeamento ‘direto’ bóson para qubit. O cálculo é então executado como:

[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

Para ter um número diferente de modais por modo:

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

Agora que o Hamiltoniano está pronto, ele pode ser usado em um algoritmo quântico para encontrar informações sobre a estrutura vibrônica da molécula correspondente. Confira nossos tutoriais em Cálculo do Estado Fundamental e Cálculo dos Estados Excitados para saber mais sobre como fazer isso no 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.

[ ]: