量子化学のためのハミルトニアン

まず、VQEにおけるハミルトニアンの役割について簡単に概要を説明します。

VQEにおけるハミルトニアンの概要

Victoria Lipinska博士が、ハミルトニアンとそれを量子コンピューティングで使用するためのマッピング方法について解説します。

参考文献

以下の論文は上記の動画で参照されています。

• Quantum Algorithms for Fermionic Simulations, Ortiz, et al.
• Simulating Chemistry using Quantum Computers, Kassal et al.
• A Comparison of the Bravyi–Kitaev and Jordan–Wigner Transformations for the Quantum Simulation of Quantum Chemistry, Tranter, et al.
• Quantum Chemistry in the Age of Quantum Computing, Cao, et al.
• Quantum computational chemistry, McArdle, et al.
• The Bravyi-Kitaev transformation for quantum computation of electronic structure, Seeley, et al., McArdle, et al.

量子化学のためのハミルトニアン準備

量子コンピューティングを化学問題に適用する際の最初の良いステップは、対象システムのハミルトニアンを定義することです。ここでは、量子化学ハミルトニアンに議論を限定します。これらのハミルトニアンは、同一フェルミオン系に特有のマッピングを必要とするためです。

量子化学に携わっている方であれば、分子をモデル化するためのお気に入りのソフトウェアがすでにあるかもしれません。そのソフトウェアで、関心のあるシステムを記述するハミルトニアンを生成できるでしょう。ここでは、PySCF、numpy、およびQiskitのみを使用して構築したコードを使用します。ただし、ハミルトニアン準備のプロセスはパッケージ化されたソリューションにも同様に適用できます。このアプローチと他のソフトウェアの違いは、わずかな構文の違いだけです。既存のワークフローの統合を容易にするために、これらの一部については「サードパーティソフトウェア」のサブセクションで取り上げています。

IBM Quantum® QPU上で使用する量子化学ハミルトニアンを生成するには、以下の手順が必要です：

1. 分子を定義する（形状、スピン、活性空間など）
2. フェルミオンハミルトニアンを生成する（生成演算子と消滅演算子）
3. フェルミオンハミルトニアンをボゾン演算子にマッピングする（この文脈ではパウリ演算子を使用）
4. サードパーティソフトウェアを使用する場合：生成ソフトウェアとQiskitの構文の不一致を処理する

フェルミオンハミルトニアンはフェルミオン演算子で記述され、特に電子が区別不可能なフェルミオンであることを考慮しています。これは、電子が区別可能なボゾン量子ビットとは全く異なる統計に従うことを意味します。これがマッピングプロセスが必要な理由です。

これらのプロセスをすでによく理解している方は、このセクションをスキップしても構いません。 目標：

最終的な目標は、以下の形式のハミルトニアンを得ることです：

# Added by doQumentation — required packages for this notebook
!pip install -q numpy openfermion pyscf qiskit

H = [(1, "XX"), (1, "YY"), (1, "ZZ")]
print(H)

[(1, 'XX'), (1, 'YY'), (1, 'ZZ')]

または

from qiskit.quantum_info import SparsePauliOp

H = SparsePauliOp(["XX", "YY", "ZZ"], coeffs=[1.0 + 0.0j, 1.0 + 0.0j, 1.0 + 0.0j])
print(H)

SparsePauliOp(['XX', 'YY', 'ZZ'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j])

まず、いくつかのパッケージをインポートします：

import numpy as np
from pyscf import ao2mo, gto, mcscf, scf

1. 分子を定義する

ここでは、対象分子の属性を指定します。この例では、二原子水素（H₂）を選択しています（結果として得られるハミルトニアンが表示できる程度に短いためです）。 Pythonベースの化学シミュレーションフレームワーク（PySCF）は、量子計算に適した分子ハミルトニアンの生成など、様々な用途に使用できる豊富な電子構造モジュールを備えています。PySCFクイックスタートガイドは、すべての変数と機能の詳細な説明として優れたリソースです。多くの方にとってすでに馴染みのある内容であるため、ここでは最低限の概要のみを説明します。詳細についてはPySCFをご覧ください。 簡単に説明すると：

distance：二原子分子に使用できます。または各原子のデカルト座標を直接指定します。距離の単位はオングストロームです。

gto：ガウス型軌道を生成します。

basis：分子軌道のモデル化に使用する関数を指します。ここで使用する 'sto-6g' は一般的な最小基底で、6つの基本ガウス軌道を使ってスレーター型軌道（STO）をフィッティングすることに由来しています。

spin：不対電子数を示す整数値（$2S$ に等しい）。一部のソフトウェアでは多重度（$2S+1$）を使用する点に注意してください。

charge：分子の電荷。

symmetry：分子の点対称群。文字列で指定するか、"symmetry = True" を設定して自動検出させます。ここで "Dooh" は同種の原子2つからなる二原子分子に適した対称群です。

distance = 0.735
a = distance / 2
mol = gto.Mole()
mol.build(
    verbose=0,
    atom=[
        ["H", (0, 0, -a)],
        ["H", (0, 0, a)],
    ],
    basis="sto-6g",
    spin=0,
    charge=0,
    symmetry="Dooh",
)

<pyscf.gto.mole.Mole at 0x7fc718f07610>

全エネルギー（核反発エネルギーと電子エネルギーを含む）、全電子軌道エネルギー、または電子軌道の一部のエネルギー（残りの部分を凍結した場合）を記述できることに注意してください。特に $\text{H}_2$ の場合、以下に示す異なるエネルギーを確認してください。全エネルギーから核反発エネルギーを引くと、電子エネルギーが得られることがわかります：

mf = scf.RHF(mol)
mf.scf()

print(
    mf.energy_nuc(),
    mf.energy_elec()[0],
    mf.energy_tot(),
    mf.energy_tot() - mol.energy_nuc(),
)

0.7199689944489797 -1.8455976628764188 -1.125628668427439 -1.8455976628764188

active_space = range(mol.nelectron // 2 - 1, mol.nelectron // 2 + 1)

2. フェルミオンハミルトニアンを生成する

scf：幅広い自己無撞着場法を指します。

rhf：mf = scf.RHF(mol) における mf は、制限ハートリーフォック（Restricted Hartree Fock）計算を使用するソルバーです。このカーネル（以下の E）は、核反発と分子軌道を含む全エネルギーです。

mcscf：多配置自己無撞着場（multi-configuration self-consistent fields）パッケージです。

ao2mo：原子軌道から分子軌道への変換です。

また、以下の変数を使用します：

ncas：完全活性空間（complete active space）内の軌道数

nelecas：完全活性空間内の電子数

E1 = mf.kernel()
mx = mcscf.CASCI(mf, ncas=2, nelecas=(1, 1))
mo = mx.sort_mo(active_space, base=0)
E2 = mx.kernel(mo)[:2]

ここでハミルトニアンが必要になります。これは通常、電子コアのエネルギー（ecore、最小化には関与しない）、一電子演算子（h1e）、および二電子エネルギー（h2e）に分けられます。これらは最後の2行で明示的に抽出されます。

h1e, ecore = mx.get_h1eff()
h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)

これらのハミルトニアンは現在、（区別不可能な）フェルミオン系に適用可能なフェルミオン（生成および消滅）演算子で表されており、交換に対する反対称性に従います。これにより、区別可能またはボゾン系に適用される統計とは異なる統計が生じます。IBM Quantum QPUで計算を実行するには、エネルギーを記述するボゾン演算子が必要です。このようなマッピングの結果は、パウリ演算子がエルミートかつユニタリーであるため、慣例的にパウリ演算子で記述されます。使用できるマッピングはいくつかあります。最もシンプルなものの1つがJordan-Wigner変換です。

3. ハミルトニアンをマッピングする

化学ハミルトニアンを量子コンピューター上で実行するのに適した形式にマッピングするためのツールは多数あります。ここでは、PySCF、numpy、およびQiskitのみを使用してJordan-Wignerマッピングを直接実装します。他のソリューションの構文上の考慮事項については以下でコメントします。 Cholesky関数は、ハミルトニアン中の二電子項の低ランク分解を得るのに役立ちます。

def cholesky(V, eps):
    # see https://arxiv.org/pdf/1711.02242.pdf section B2
    # see https://arxiv.org/abs/1808.02625
    # see https://arxiv.org/abs/2104.08957
    no = V.shape[0]
    chmax, ng = 20 * no, 0
    W = V.reshape(no**2, no**2)
    L = np.zeros((no**2, chmax))
    Dmax = np.diagonal(W).copy()
    nu_max = np.argmax(Dmax)
    vmax = Dmax[nu_max]
    while vmax > eps:
        L[:, ng] = W[:, nu_max]
        if ng > 0:
            L[:, ng] -= np.dot(L[:, 0:ng], (L.T)[0:ng, nu_max])
        L[:, ng] /= np.sqrt(vmax)
        Dmax[: no**2] -= L[: no**2, ng] ** 2
        ng += 1
        nu_max = np.argmax(Dmax)
        vmax = Dmax[nu_max]
    L = L[:, :ng].reshape((no, no, ng))
    print(
        "accuracy of Cholesky decomposition ",
        np.abs(np.einsum("prg,qsg->prqs", L, L) - V).max(),
    )
    return L, ng

identity 関数と creators_destructors 関数は、フェルミオンハミルトニアンの生成演算子と消滅演算子をパウリ演算子に置き換えます。creators_destructors はJordan-Wignerマッピングを使用します。

def identity(n):
    return SparsePauliOp.from_list([("I" * n, 1)])

def creators_destructors(n, mapping="jordan_wigner"):
    c_list = []
    if mapping == "jordan_wigner":
        for p in range(n):
            if p == 0:
                ell, r = "I" * (n - 1), ""
            elif p == n - 1:
                ell, r = "", "Z" * (n - 1)
            else:
                ell, r = "I" * (n - p - 1), "Z" * p
            cp = SparsePauliOp.from_list([(ell + "X" + r, 0.5), (ell + "Y" + r, -0.5j)])
            c_list.append(cp)
    else:
        raise ValueError("Unsupported mapping.")
    d_list = [cp.adjoint() for cp in c_list]
    return c_list, d_list

最後に、build_hamiltonian は cholesky、identity、および creators_destructors 関数を使用して、量子コンピューター上で実行するのに適した最終的なハミルトニアンを作成します。

def build_hamiltonian(ecore: float, h1e: np.ndarray, h2e: np.ndarray) -> SparsePauliOp:
    ncas, _ = h1e.shape

    C, D = creators_destructors(2 * ncas, mapping="jordan_wigner")
    Exc = []
    for p in range(ncas):
        Excp = [C[p] @ D[p] + C[ncas + p] @ D[ncas + p]]
        for r in range(p + 1, ncas):
            Excp.append(
                C[p] @ D[r]
                + C[ncas + p] @ D[ncas + r]
                + C[r] @ D[p]
                + C[ncas + r] @ D[ncas + p]
            )
        Exc.append(Excp)

    # low-rank decomposition of the Hamiltonian
    Lop, ng = cholesky(h2e, 1e-6)
    t1e = h1e - 0.5 * np.einsum("pxxr->pr", h2e)

    H = ecore * identity(2 * ncas)
    # one-body term
    for p in range(ncas):
        for r in range(p, ncas):
            H += t1e[p, r] * Exc[p][r - p]
    # two-body term
    for g in range(ng):
        Lg = 0 * identity(2 * ncas)
        for p in range(ncas):
            for r in range(p, ncas):
                Lg += Lop[p, r, g] * Exc[p][r - p]
        H += 0.5 * Lg @ Lg

    return H.chop().simplify()

最後に、build_hamiltonian を使用して、Jordan-Wigner変換によりパウリ演算子からクビットハミルトニアンを構築します。また、使用したCholesky分解の精度も確認できます。

H = build_hamiltonian(ecore, h1e, h2e)
print(H)

accuracy of Cholesky decomposition  2.220446049250313e-16
SparsePauliOp(['IIII', 'IIIZ', 'IZII', 'IIZI', 'ZIII', 'IZIZ', 'IIZZ', 'ZIIZ', 'IZZI', 'ZZII', 'ZIZI', 'YYYY', 'XXYY', 'YYXX', 'XXXX'],
              coeffs=[-0.09820182+0.j, -0.1740751 +0.j, -0.1740751 +0.j,  0.2242933 +0.j,
  0.2242933 +0.j,  0.16891402+0.j,  0.1210099 +0.j,  0.16631441+0.j,
  0.16631441+0.j,  0.1210099 +0.j,  0.17504456+0.j,  0.04530451+0.j,
  0.04530451+0.j,  0.04530451+0.j,  0.04530451+0.j])

この分子サンプルノートブックには、様々な複雑さの分子のセットアップとハミルトニアンが示されています。少し修正するだけで、ほとんどの小分子を調べることができます。

分子のフェルミオン演算子を構築する際に考慮すべき2つの重要な点を簡単に確認しておきましょう。分子の種類が変わると、対称性も変わります。同様に、円筒対称な "A1" のように、様々な対称性を持つ軌道の数も変わります。LiHへの単純な拡張でもこれらの変化が明らかになります：

distance = 1.56
mol = gto.Mole()
mol.build(
    verbose=0,
    atom=[["Li", (0, 0, 0)], ["H", (0, 0, distance)]],
    basis="sto-6g",
    spin=0,
    charge=0,
    symmetry="Coov",
)
mf = scf.RHF(mol)
E1 = mf.kernel()

# %% ----------------------------------------------------------------------------------------------

mx = mcscf.CASCI(mf, ncas=5, nelecas=(1, 1))
cas_space_symmetry = {"A1": 3, "E1x": 1, "E1y": 1}
mo = mcscf.sort_mo_by_irrep(mx, mf.mo_coeff, cas_space_symmetry)
E2 = mx.kernel(mo)[:2]
h1e, ecore = mx.get_h1eff()
h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)

Jordan-Wignerマッパーを使用したLiHのハミルトニアンがすでに276項から成ることにも注目する価値があります。最終的なハミルトニアンに対する直感は、すぐに失われてしまいます。

len(build_hamiltonian(ecore, h1e, h2e))

accuracy of Cholesky decomposition  1.1102230246251565e-16

276

対称性に迷った場合は、symmetry = True と verbose = 4 を設定することで、分子の対称性情報を生成することもできます：

distance = 1.56
mol = gto.Mole()
mol.build(
    verbose=4,
    atom=[["Li", (0, 0, 0)], ["H", (0, 0, distance)]],
    basis="sto-6g",
    spin=0,
    charge=0,
    symmetry=True,
)

System: uname_result(system='Linux', node='IBM-R912JTRT', release='5.10.102.1-microsoft-standard-WSL2', version='#1 SMP Wed Mar 2 00:30:59 UTC 2022', machine='x86_64')  Threads 16
Python 3.11.12 (main, May 16 2025, 02:33:32) [GCC 11.4.0]
numpy 2.3.1  scipy 1.16.0  h5py 3.14.0
Date: Mon Jun 30 12:56:55 2025
PySCF version 2.9.0
PySCF path  /home/porter284/.pyenv/versions/3.11.12/lib/python3.11/site-packages/pyscf

[CONFIG] conf_file None
[INPUT] verbose = 4
[INPUT] num. atoms = 2
[INPUT] num. electrons = 4
[INPUT] charge = 0
[INPUT] spin (= nelec alpha-beta = 2S) = 0
[INPUT] symmetry True subgroup None
[INPUT] Mole.unit = angstrom
[INPUT] Symbol           X                Y                Z      unit          X                Y                Z       unit  Magmom
[INPUT]  1 Li     0.000000000000   0.000000000000   0.000000000000 AA    0.000000000000   0.000000000000   0.000000000000 Bohr   0.0
[INPUT]  2 H      0.000000000000   0.000000000000   1.560000000000 AA    0.000000000000   0.000000000000   2.947972754321 Bohr   0.0

nuclear repulsion = 1.01764848253846
point group symmetry = Coov
symmetry origin: [0.         0.         0.73699319]
symmetry axis x: [1. 0. 0.]
symmetry axis y: [0. 1. 0.]
symmetry axis z: [0. 0. 1.]
num. orbitals of irrep A1 = 4
num. orbitals of irrep E1x = 1
num. orbitals of irrep E1y = 1
number of shells = 4
number of NR pGTOs = 36
number of NR cGTOs = 6
basis = sto-6g
ecp = {}
CPU time:         9.85

<pyscf.gto.mole.Mole at 0x7fc719f94850>

この出力には多くの有用な情報のほか、point group symmetry = Coov と各既約表現における軌道数の両方が含まれています。

point group symmetry = Coov
num. orbitals of irrep A1 = 4
num. orbitals of irrep E1x = 1
num. orbitals of irrep E1y = 1
number of shells = 4

これは必ずしも活性空間に含める軌道数を教えてくれるわけではありませんが、どのような軌道が存在し、それぞれの対称性を確認するのに役立ちます。

対称性と軌道を指定することはしばしば有用ですが、含めたい軌道数を直接指定することもできます。以下のエチレンの例を考えてみましょう。verbose = 4 を使用することで、様々な軌道の対称性を出力できます：

# Replace these variables with correct distances:
a = 1
b = 1
c = 1

# Build
mol = gto.Mole()
mol.build(
    verbose=4,
    atom=[
        ["C", (0, 0, a)],
        ["C", (0, 0, -a)],
        ["H", (0, c, b)],
        ["H", (0, -c, b)],
        ["H", (0, c, -b)],
        ["H", (0, -c, -b)],
    ],
    basis="sto-6g",
    spin=0,
    charge=0,
    symmetry=True,
)

System: uname_result(system='Linux', node='IBM-R912JTRT', release='5.10.102.1-microsoft-standard-WSL2', version='#1 SMP Wed Mar 2 00:30:59 UTC 2022', machine='x86_64')  Threads 16
Python 3.11.12 (main, May 16 2025, 02:33:32) [GCC 11.4.0]
numpy 2.3.1  scipy 1.16.0  h5py 3.14.0
Date: Mon Jun 30 12:57:07 2025
PySCF version 2.9.0
PySCF path  /home/porter284/.pyenv/versions/3.11.12/lib/python3.11/site-packages/pyscf

[CONFIG] conf_file None
[INPUT] verbose = 4
[INPUT] num. atoms = 6
[INPUT] num. electrons = 16
[INPUT] charge = 0
[INPUT] spin (= nelec alpha-beta = 2S) = 0
[INPUT] symmetry True subgroup None
[INPUT] Mole.unit = angstrom
[INPUT] Symbol           X                Y                Z      unit          X                Y                Z       unit  Magmom
[INPUT]  1 C      0.000000000000   0.000000000000   1.000000000000 AA    0.000000000000   0.000000000000   1.889726124565 Bohr   0.0
[INPUT]  2 C      0.000000000000   0.000000000000  -1.000000000000 AA    0.000000000000   0.000000000000  -1.889726124565 Bohr   0.0
[INPUT]  3 H      0.000000000000   1.000000000000   1.000000000000 AA    0.000000000000   1.889726124565   1.889726124565 Bohr   0.0
[INPUT]  4 H      0.000000000000  -1.000000000000   1.000000000000 AA    0.000000000000  -1.889726124565   1.889726124565 Bohr   0.0
[INPUT]  5 H      0.000000000000   1.000000000000  -1.000000000000 AA    0.000000000000   1.889726124565  -1.889726124565 Bohr   0.0
[INPUT]  6 H      0.000000000000  -1.000000000000  -1.000000000000 AA    0.000000000000  -1.889726124565  -1.889726124565 Bohr   0.0

nuclear repulsion = 29.3377079104231
point group symmetry = D2h
symmetry origin: [0. 0. 0.]
symmetry axis x: [0. 1. 0.]
symmetry axis y: [1. 0. 0.]
symmetry axis z: [-0. -0. -1.]
num. orbitals of irrep Ag = 4
num. orbitals of irrep B2g = 2
num. orbitals of irrep B3g = 1
num. orbitals of irrep B1u = 4
num. orbitals of irrep B2u = 1
num. orbitals of irrep B3u = 2
number of shells = 10
number of NR pGTOs = 84
number of NR cGTOs = 14
basis = sto-6g
ecp = {}
CPU time:         9.92

<pyscf.gto.mole.Mole at 0x7fc719fa9290>

次の結果が得られます：

num. orbitals of irrep Ag = 4

num. orbitals of irrep B2g = 2

num. orbitals of irrep B3g = 1

num. orbitals of irrep B1u = 4

num. orbitals of irrep B2u = 1

num. orbitals of irrep B3u = 2

しかし、すべての軌道を対称性で指定する代わりに、単純に次のように記述することもできます：

active_space = range(mol.nelectron // 2 - 2, mol.nelectron // 2 + 2)

このアプローチでは、充填レベル付近（価電子軌道と非占有軌道）のいくつかの軌道を選択します。ここでは、活性空間に5つの軌道が選択されています（6番目から10番目）。

print(
    mol.nelectron // 2 - 2,
    mol.nelectron // 2 + 2,
)

6 10

4. サードパーティソフトウェア

量子化学向けに開発されたソフトウェアパッケージは複数あり、複数のマッパーや活性空間を制限するためのツールを提供しているものもあります。上記の手順は一般的なものであり、サードパーティソフトウェアにも同様に適用されます。ただし、他のソフトウェアがQiskitで受け入れられない形式でハミルトニアンを返す場合があります。例えば、一部のソフトウェアは以下の形式でハミルトニアンを返します：

H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3] 特に、ゲートに番号が付いており、恒等演算子が表示されていないことに注目してください。これは、項 [Z2 Z3] を ZZII（クビット0と1に恒等演算子、クビット2と3にZ演算子が作用し、クビット0が右端に来るよう並べられている）と記述するQiskitのハミルトニアン形式とは対照的です。

既存のワークフローに対応するために、以下のコードブロックで一方の構文からもう一方への変換を行います。convert_openfermion_to_qiskit 関数は、OpenFermionまたはTangeloで生成されたハミルトニアン（任意の利用可能なマッパーを使用してすでにパウリ演算子にマッピングされているもの）と、分子に必要なクビット数を引数として受け取ります。

from openfermion import QubitOperator
from qiskit.quantum_info import SparsePauliOp

def convert_openfermion_to_qiskit(
    openfermion_operator: QubitOperator, num_qubits: int
) -> SparsePauliOp:
    terms = openfermion_operator.terms

    labels = []
    coefficients = []

    for term, constant in terms.items():
        # Default set to identity
        operator = list("I" * num_qubits)

        # Iterate through PauliSum and replace I with Pauli
        for index, pauli in term:
            operator[index] = pauli
        label = "".join(operator)
        labels.append(label)
        coefficients.append(constant)

    return SparsePauliOp(labels, coefficients)

さらに、このPythonノートブックには、上記の変換を含む、他のソフトウェアワークフローからQiskitへのハミルトニアン移行に関する完全なサンプルコードが含まれています。

これで、IBM®量子コンピューターで量子化学計算を実行するために必要なハミルトニアンを取得するための豊富なツールが揃いました。