分子上的变分量子本征求解器 (VQE)#

概述#

VQE 是一种变分算法，用于计算满足 \(H \left|\psi_g\right> =E_g\left|\psi_g\right>\) 的给定哈密顿 H 的基态，我们称之为 \(\psi_g\)。对于任意归一化波函数 \(\psi_f\)，期望值 \(\left<\psi_f|H|\psi_f \right>\) 总是不低于基态能量，除非 \(\psi_f = \psi_g\) (这里我们假设基态没有简并）。基于这个事实，如果我们使用参数化波函数 \(\psi_\theta\)，例如由具有参数 \(\theta\) 的参数化量子电路 (PQC) 给出，我们可以通过最小化 \(H\) 的期望值来给出基态能量和波函数的近似值。在实际的量子硬件中，该算法可以在量子-经典混合范式中实现，在量子硬件中使用有限差分或参数移位计算梯度，在经典计算机中使用梯度下降法进行优化。在数值模拟中，我们可以使用自动微分计算梯度。

计算分子的基态能量对于量子化学任务通常很重要，因为它可以用来找出分子的原子结构。在分子的模拟中，我们不考虑原子核的运动，这意味着我们固定了其组成原子的原子核坐标。我们只考虑分子中的电子，因为原子核比电子重得多，因此声子携带的能量可以忽略不计，或者可以使用 Born-Oppenheimer 近似重新考虑。严格来说，电子存在于连续空间中，因此希尔伯特空间是无限维的。为了进行实际计算，我们只保留一些重要的单粒子基，例如低能原子轨道。在二次量子化形式中，我们可以将这些原子轨道表示为 \(c_i^\dagger|0>\)。通过考虑原子核和电子的相互作用作为背景以及电子 - 电子相互作用，分子哈密顿量通常可以表示为 \(H = \sum_{i, j} h_{i,j} c_i^\dagger c_j + \sum_{i, j, k, l} \alpha_{i, j, k, l} c_i^\dagger c_j^\dagger c_k c_l\)。请注意，自旋指数也被吸收到轨道指数中。有很多软件可以在 H 中给出这些参数，例如我们将在本教程后面使用的 pyscf。现在我们有了对分子的费米子描述。通过使用从费米子算子到自旋算子的映射，例如 Jordan-Wigner 变换或 Bravyi-Kitaev 变换，我们可以将费米子哈密顿算子映射到与量子计算机更兼容的自旋哈密顿算子。对于自旋哈密顿算子，我们可以很容易地使用 PQC 来构造轨迹波函数并进行 VQE 算法。在本教程的以下部分，我们将演示如何使用 TensorCircuit 在分子上模拟 VQE 算法的完整示例。

设置#

我们应该首先 pip install openfermion openfermionpyscf 根据 openfermion 和 pyscf 提供的量子化学计算生成H2O分子的费米子和量子比特哈密顿量。

[1]:

import numpy as np
from openfermion.chem import MolecularData
from openfermion.transforms import (
    get_fermion_operator,
    jordan_wigner,
    binary_code_transform,
    checksum_code,
    reorder,
)
from openfermion.chem import geometry_from_pubchem
from openfermion.utils import up_then_down
from openfermion.linalg import LinearQubitOperator
from openfermionpyscf import run_pyscf
import tensorflow as tf

import tensorcircuit as tc

K = tc.set_backend("tensorflow")

生成哈密顿量#

  • 获取分子能量信息和分子轨道

[2]:

multiplicity = 1
basis = "sto-3g"
# H2O 的 14 个自旋轨道
geometry = geometry_from_pubchem("h2o")
description = "h2o"
molecule = MolecularData(geometry, basis, multiplicity, description=description)
molecule = run_pyscf(molecule, run_mp2=True, run_cisd=True, run_ccsd=True, run_fci=True)
print(molecule.fci_energy, molecule.ccsd_energy, molecule.hf_energy)

-75.0155301894916 -75.01540899923558 -74.96444758276998

  • 获取费米子哈密顿量

[3]:

mh = molecule.get_molecular_hamiltonian()

[4]:

fh = get_fermion_operator(mh)

[5]:

print(fh.terms[((0, 1), (0, 0))])  # 获取费米子哈密顿量

-32.68991541360029

  • 转换为量子比特哈密顿量

[6]:

# 对于 H2O 的 14 个轨道，诸如 JW 或 BK 之类的正常变换需要 14 个量子位

a = jordan_wigner(fh)
LinearQubitOperator(a).n_qubits

[6]:

14

我们可以使用二进制编码来保存另外两个量子位，因为自旋向上和自旋向下填充的数量都是 5（5/7 个轨道中的奇数电子）

[7]:

b = binary_code_transform(reorder(fh, up_then_down), 2 * checksum_code(7, 1))
# 7 是 7 个自旋极化轨道，1 是奇数占用
LinearQubitOperator(b).n_qubits

[7]:

12

[8]:

print(b.terms[((0, "Z"),)])  # Z_0 泡利字符串的系数

12.412562749393349

  • 将 openfermion 中的量子比特哈密顿量转换为 TensorCircuit 中的格式

[9]:

lsb, wb = tc.templates.chems.get_ps(b, 12)
lsa, wa = tc.templates.chems.get_ps(a, 14)

  • 以矩阵形式检查哈密顿量

[10]:

ma = tc.quantum.PauliStringSum2COO_numpy(lsa, wa)

[11]:

mb = tc.quantum.PauliStringSum2COO_numpy(lsb, wb)

[12]:

mad, mbd = ma.todense(), mb.todense()

这两种哈密顿量对应的 Hartree Fock 乘积状态

[13]:

bin(np.argmin(np.diag(mad)))

[13]:

'0b11111111110000'

[14]:

bin(np.argmin(np.diag(mbd)))

[14]:

'0b111110111110'

VQE 设置#

原则上，我们可以将哈密顿量的每个泡利串评估为期望测量，但它会花费大量模拟时间，相反，我们将它们融合为如上所示的哈密顿矩阵来运行 VQE。

  • 使用密集矩阵期望

[15]:

n = 12
depth = 4
mbd_tf = tc.array_to_tensor(mbd)


def vqe(param):
    c = tc.Circuit(n)
    for i in [0, 1, 2, 3, 4, 6, 7, 8, 9, 10]:
        c.X(i)
    for j in range(depth):
        for i in range(n - 1):
            c.exp1(i, i + 1, unitary=tc.gates._xx_matrix, theta=param[j, i, 0])
        for i in range(n):
            c.rx(i, theta=param[j, i, 1])
        for i in range(n):
            c.ry(i, theta=param[j, i, 2])
        for i in range(n):
            c.rx(i, theta=param[j, i, 3])
    return tc.templates.measurements.operator_expectation(c, mbd_tf)

[16]:

vags = tc.backend.jit(tc.backend.value_and_grad(vqe))
lr = tf.keras.optimizers.schedules.ExponentialDecay(
    decay_rate=0.5, decay_steps=300, initial_learning_rate=0.5e-2
)
opt = tc.backend.optimizer(tf.keras.optimizers.Adam(lr))

param = tc.backend.implicit_randn(shape=[depth, n, 4], stddev=0.02, dtype="float32")
for i in range(600):
    e, g = vags(param)
    param = opt.update(g, param)
    if i % 100 == 0:
        print(e)

tf.Tensor(-74.76671, shape=(), dtype=float32)
tf.Tensor(-74.95493, shape=(), dtype=float32)
tf.Tensor(-74.95319, shape=(), dtype=float32)
tf.Tensor(-74.954315, shape=(), dtype=float32)
tf.Tensor(-74.956116, shape=(), dtype=float32)
tf.Tensor(-74.95809, shape=(), dtype=float32)

  • 使用稀疏矩阵期望

我们还可以使用稀疏哈密顿矩阵进行电路期望评估，唯一的区别是将 mbd_tf 替换为 mb_tf

[17]:

mb_tf = tc.backend.coo_sparse_matrix(
    np.transpose(np.stack([mb.row, mb.col])), mb.data, shape=(2**n, 2**n)
)

稀疏矩阵评估和密集矩阵评估之间的一个微基准，用于比较计算期望的时间，当然，稀疏总是在空间方面获胜。

[18]:

def dense_expt(param):
    c = tc.Circuit(n)
    for i in range(n):
        c.H(i)
        c.rx(i, theta=param[i])
    return tc.templates.measurements.operator_expectation(c, mbd_tf)


def sparse_expt(param):
    c = tc.Circuit(n)
    for i in range(n):
        c.H(i)
        c.rx(i, theta=param[i])
    return tc.templates.measurements.operator_expectation(c, mb_tf)

[19]:

dense_vag = tc.backend.jit(tc.backend.value_and_grad(dense_expt))
sparse_vag = tc.backend.jit(tc.backend.value_and_grad(sparse_expt))

v0, g0 = dense_vag(tc.backend.ones([n]))
v1, g1 = sparse_vag(tc.backend.ones([n]))

# 一致性检查

np.testing.assert_allclose(v0, v1, atol=1e-5)
np.testing.assert_allclose(g0, g1, atol=1e-5)

[20]:

%timeit dense_vag(tc.backend.ones([n]))

30.7 ms ± 1.45 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

[21]:

%timeit sparse_vag(tc.backend.ones([n]))

3.6 ms ± 63 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

因此，稀疏矩阵求值除了节省空间外，还可以节省时间，这总是被推荐的。