探究

实验室

计算

公开

量子计算｜高斯玻色采样模拟分子振动光谱

量子化学

量子计算

量子变分算法

量子算法

量子化学量子计算量子变分算法量子算法

JiaweiMiao

更新于 2025-03-09

推荐镜像 :DeePMD-kit:2.2.8-cuda12.0

推荐机型 :c4_m15_1 * NVIDIA T4

数据集

DeepQuamtum_dataset(v2)

量子计算｜高斯玻色采样模拟分子振动光谱

©️ Copyright 2023 @ Authors
作者： 图灵算法
日期：2024-05-23
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
Note: 本 notebook 搬运自 DeepQuantum 网站

此研究涉及两个主要领域。其中之一是关于量子信息和量子计算机的研究。我们选择玻色采样进行讨论。目前，有多种物理实现方式可以进行玻色子取样，但最常见的是光子解决方案。在该方案中，单光子态被用作输入到干涉仪，在输出端测量光子数。干涉装置相对简单；而制备单光子和测量光子数对仪器的要求较高。因此，如果我们考虑使用相干激光，可以让实验简便一些。由于相干态可以用高斯分布描述，这种类型的玻色子采样被称为高斯玻色采样。

另一个主要领域是光谱学。我们以分子振动光谱为例，模拟不同带电态之间转变的概率，从而近似确定分子振动谱。这种方法可以使用高斯玻色采样来实现。参考 Huh 的研究工作，我们在基于玻色采样的模拟器上实现该模拟过程。

代码

文本

高斯玻色采样

量子比特是一个抽象的二维量子物理系统。其希尔伯特空间通常可由两个基态 $∣ 0 ⟩$ 和 $∣ 1 ⟩$ 来描述。不同量子比特可以处于叠加态和纠缠状态。我们可将其推广到3维、4维以至有限的 $d$ -维空间中。

以上内容可以进一步推广到无限维度。玻色子粒子是这种类型的量子物理系统的一个例子，其中在给定状态中可以有无限多个粒子。这样的系统仍然是量子系统，仍具有叠加和纠缠性。在可数无限维度情况下，量子态被称为 qumodes。

在真实物理系统中，任何物理量都需要通过多次测量来确定其状态分布。因此量子计算机都是对抽样问题的物理实现。在可数无限维度的情况下，这被称为玻色子采样。

在多种玻色采样的物理实现方式中，最简单的仍然是光子方案。量子光学设备在 qumodes 上为幺正算符。常用高斯算符如下：

$P_{i} (ϕ) = e^{- i ϕ a_{i}^{†} a_{i}}$

$D_{i} (α) = e^{α a_{i}^{†} - α^{*} a_{i}}$

$S_{i} (ζ) = e^{\frac{1}{2} (ζ^{*} a_{i}^{2} - ζ a_{i}^{† 2})}$

$B_{ij} (θ) = e^{θ (a_{i}^{†} a_{j} - a_{i} a_{j}^{†})}$

事实上，一个纯的光子数态，即 $∣ n ⟩$ , 在技术上很难做到。 $∣ 1 ⟩$ 意味着特定量子模式中只有一个光子。精确地产生一个光子是一项巨大挑战。而相干光是每个量子模式的可数无限状态之和，高斯玻色采样以相干光作为输入，是一种更简单的解决方案。

代码

文本

分子的振动光谱简要描述

Franck–Condon原理能级图

代码

文本

光谱学的一个主要问题是检查原子和分子之间电荷状态的跃迁。分子吸收的光频率便是取决于不同电子态之间的允许跃迁。这些电子跃迁可能伴随着分子振动能量的变化。在这种情况下，代表更强烈吸收光的频率的吸收线被称为振动光谱。然而，除了氢原子可以被精确解决外，随着电子数量增加，量子力学方程变得指数级复杂，使得对它们的计算变得不可能。因此必须使用近似方法。

如果假设核和电子状态是独立的，对于极其复杂的分子系统，一个很好的近似称为玻恩-奥本海默近似。由于原子核的质量较大，其对电子壳层变化反应更慢，因此在 Franck-Condon 原理中可以将原子核视为恒定。

通过用简单的量子谐振子来近似电子的振动，可以将电子的哈密顿算符写成反应坐标的函数。然后，电子势能表面已经可以用抛物线描述。

$H = p^{2} + q^{2}$

Duschinsky 提出一种近似方法，不同电子态的简正坐标之间存在线性关系

$q^{'} = U_{D} q + d$

Doktorov 证明，在这些条件下，对于这样的量子谐振子，不同势能面上的量子态 $∣ ϕ ⟩$ 和 $∣ ϕ^{'} ⟩$ 之间的关系可以用以下方式给出:

$∣ ϕ^{'} ⟩ = U_{Do k t} ∣ ϕ ⟩$

其中，Doktorov 算子定义为

$U_{Do k} = D (d) S^{†} (Ω^{'}) R (U_{D u sc h}) S_{Ω} (Ω)$

旋转算符的输入参数是Duschinsky混合矩阵本身
位移算符D的参数是Duschinsky位移
压缩算符是从分子的物理特性推导出来的。 $Ω$ 和 $Ω^{'}$ 是分子内原子在发生跃迁前后谐波的角频率

上面哈密顿量的本征值问题的解是 Fock 空间中的相干态。高斯算符作用于它们之后，将一个相干态转换为另一个相干态。不同电子态之间的转移概率，即所谓的 Franck-Condon 因子，可以近似表示为

$FC F_{n^{'}, n} = ∣ ⟨ ϕ^{'} ∣ ϕ ⟩ ∣^{2} = ∣ ⟨ n^{'} ∣ U_{Do k t} ∣ n ⟩ ∣^{2}$

玻色子采样的计算在于对 qumodes 进行旋转矩阵操作。在我们计算振动光谱时，使用上面的 Doktorov 算符。其过程为

制备相干态
将相干态传输到一个玻色采样设备，以获得转移概率
进行第二次压缩，测量各模式分布，进一步计算FCF

代码

文本

Hessian matrix

对于一双原子分子，假设谐振子势能

$V = \frac{k}{2} (r - r_{e})^{2}$

其中 $r$ 是两个核之间的距离， $r_{e}$ 是它们的平衡距离。简单起见，可只考虑X轴。我们可以用两个原子的笛卡尔坐标 $x_{1}$ （对应原子1）和 $x_{2}$ （对应原子2）来表示这一点。则， $r = x_{2} - x_{1}$ 。

势能的二阶导数是：

$\frac{d ^{2} V}{d x _{1}^{2}} = k$

$\frac{d ^{2} V}{d x _{2}^{2}} = k$

$\frac{d ^{2} V}{d x _{1} d x _{2}} = \frac{d ^{2} V}{d x _{2} d x _{1}} = - k$

Hessian矩阵（仅针对x方向）为：

$H = (k - k - k k)$

质量加权的Hessian矩阵是：

$F = (\frac{k}{m _{1}} - \frac{k}{m _{1} m _{2}} - \frac{k}{m _{1} m _{2}} \frac{k}{m _{2}})$

如果我们有了系统的Hessian矩阵和质量加权Hessian矩阵，就可以由特征向量定义质量加权简正模式。接着，利用分子的平衡结构坐标数据、质量加权简正坐标等数据可以求解所需高斯算符参数。

Jankowiak的补充材料提供了一些分子的相关数据。我们以甲酸（ $1^{1} A^{'} \to 1^{2} A^{'}$ ）为例计算分子振动光谱：

代码

文本

[1]

%%capture

import sys

# !{sys.executable} -m pip install torchvision

!{sys.executable} -m pip install deepquantum -i https://pypi.tuna.tsinghua.edu.cn/simple

代码

文本

[2]

import numpy as np

import deepquantum as dq

import torch

import matplotlib.pyplot as plt

代码

文本

$1^{1} A^{'}$ 态及 $1^{2} A^{'}$ 态的平衡结构坐标如下

代码

文本

[3]

ri = np.genfromtxt("/bohr/imagenet-w3dd/v2/data_vibronic_excita/formic_ri.csv", delimiter=",", skip_header=0)[:, np.newaxis]

rf = np.genfromtxt("/bohr/imagenet-w3dd/v2/data_vibronic_excita/formic_rf.csv", delimiter=",", skip_header=0)[:, np.newaxis]

代码

文本

$1^{1} A^{'}$ 态及 $1^{2} A^{'}$ 态的质量加权简正坐标如下

代码

文本

[4]

li = np.genfromtxt("/bohr/imagenet-w3dd/v2/data_vibronic_excita/formic_li.csv", delimiter=",", skip_header=0)

lf = np.genfromtxt("/bohr/imagenet-w3dd/v2/data_vibronic_excita/formic_lf.csv", delimiter=",", skip_header=0)

代码

文本

$1^{1} A^{'}$ 态及 $1^{2} A^{'}$ 态简正模式频率如下

代码

文本

[5]

omega = np.genfromtxt("/bohr/imagenet-w3dd/v2/data_vibronic_excita/formic_omega.csv", delimiter=",", skip_header=0)

omegap = np.genfromtxt("/bohr/imagenet-w3dd/v2/data_vibronic_excita/formic_omegap.csv", delimiter=",", skip_header=0)

代码

文本

计算中用到的物理学常量

代码

文本

[6]

c = 299792458.0 # 光速

mu = 1.6605390666 * 10**-27 # 原子质量单位

h = 6.62607015 * 10**-34 # 普朗克常数

m_c = 12 # 碳原子相对原子质量

m_h = 1.007825037 # 氢原子相对原子质量

m_o = 15.994914640 # 氧原子相对原子质量

代码

文本

Duschinsky 矩阵 $U$ 及位移矢量 $δ$ 计算如下

代码

文本

[7]

u = []

for li_ele in li:

for lf_elf in lf:

u.append(np.sum(li_ele * lf_elf))

u = np.array(u[-1::-1]).reshape(7, 7).T

array([[ 0.99343181,  0.01440011,  0.01532633,  0.02861045,  0.06378083,
         0.07513988, -0.04280796],
       [-0.01485231,  0.99314303,  0.07419536,  0.0769291 , -0.03610189,
        -0.00248855,  0.01727882],
       [-0.01185718, -0.09164895,  0.84227531,  0.17994129, -0.38567948,
         0.30738928,  0.08008576],
       [ 0.03813492,  0.04087165, -0.34031972, -0.52311845, -0.66785609,
         0.38477092,  0.11420873],
       [-0.04133251, -0.03419326, -0.40036477,  0.76362651, -0.10356638,
         0.48381836,  0.09406589],
       [ 0.0907918 , -0.04184572, -0.09066643,  0.31511355, -0.59003401,
        -0.719268  ,  0.13037051],
       [-0.03245558,  0.00500992, -0.02063661,  0.06935002, -0.20181377,
         0.01732396, -0.97588364]])

代码

文本

[8]

delta = []

m = np.diag([m_c, m_c, m_c, m_o, m_o, m_o, m_o, m_o, m_o, m_h, m_h, m_h, m_h, m_h, m_h])

for i in range(len(omegap)):

d = lf[i].T @ np.sqrt(m) @ (ri - rf)

l = np.sqrt(h / (4 * np.pi**2 * 100 * omegap[i] * c * mu)) / (10**-10)

delta.append(d / l)

delta = np.array(delta[-1::-1])

delta

array([[ 0.22536617],
       [ 0.14689208],
       [ 1.55989779],
       [-0.37838396],
       [ 0.45525871],
       [-0.34391138],
       [ 0.06184607]])

代码

文本

在用于计算振动光谱的GBS算法中，以上这些化学参数足以确定GBS设备的配置。利用这些数据，我们即可计算甲酸分子振动光谱。

代码

文本

实际上，在可能仅涉及有限数量的光子的情况下需要非线性相互作用，第二次挤压操作在光学装置中难以实现。通常，我们需要将两次挤压操作压缩为一次：

$U_{Do k} = R_{C_{L}} S_{Σ}^{†} R_{C_{R}} D_{2^{- 1/2} J^{- 1} δ}$

其中 $J = Ω^{'} U Ω^{- 1}$ $δ = ℏ^{- 1/2} Ω^{'} d$ $J = C_{L} Σ C_{R}^{t}$

各 GBS 参数计算如下：

代码

文本

[9]

pre_transition_squeezing = np.sqrt(omega[-1::-1])

post_transition_squeezing = np.sqrt(omegap[-1::-1])

j_mat = (

np.diag(post_transition_squeezing)

@ u

@ np.linalg.inv(np.diag(pre_transition_squeezing))

)

cl, lambda_1, cr = np.linalg.svd(j_mat)

delta_2 = np.linalg.inv(j_mat) @ delta / np.sqrt(2)

delta_2 = delta_2.flatten()

lambda_2 = np.log(lambda_1)

代码

文本

甲酸分子振动光谱计算如下所示

代码

文本

[10]

modes = 7 # 简正模式数量

cutoff = 3

shots = 500000

cir = dq.photonic.QumodeCircuit(

nmode=modes,

init_state="vac",

# init_state=init_state,

cutoff=cutoff,

backend="gaussian",

)

for i in range(modes):

cir.d(wires=[i], r=delta_2[i])

cir.any(cr, wires=list(range(modes)))

for i in range(modes):

cir.s(wires=[i], r=-lambda_2[i])

cir.any(cl, wires=list(range(modes)))

state = cir()

# 线路可视化

cir.draw()

代码

文本

[11]

sample = cir.measure(shots=shots)

Automatically using MCMC to sample the final states!
chain 1: 100%|██████████████████████████| 99999/99999 [00:34<00:00, 2900.03it/s]
chain 2: 100%|█████████████████████████| 99999/99999 [00:02<00:00, 46761.71it/s]
chain 3: 100%|█████████████████████████| 99999/99999 [00:02<00:00, 41915.68it/s]
chain 4: 100%|█████████████████████████| 99999/99999 [00:02<00:00, 41636.58it/s]
chain 5: 100%|█████████████████████████| 99999/99999 [00:02<00:00, 42822.48it/s]

代码

文本

[12]

wave_number = []

counts = []

for ele in sample.items():

# print(ele[0].state)

wave_number.append(torch.sum(ele[0].state * omegap))

counts.append(ele[1])

plt.bar(wave_number, counts, width=100, color="g")

plt.xlabel(r"Energy $(cm^{-1})$")

plt.ylabel(r"Counts")

plt.xlim(-1000, 8000)

plt.show()

代码

文本

参考结果

甲酸分子振动光谱图

代码

文本

量子化学

量子计算

量子变分算法

量子算法

量子化学量子计算量子变分算法量子算法

点个赞吧

量子计算 ｜ 高斯玻色采样模拟分子振动光谱

高斯玻色采样

分子的振动光谱简要描述

Hessian matrix

甲酸分子振动光谱计算如下所示

参考结果

量子计算｜高斯玻色采样模拟分子振动光谱