{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "frontmatter",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"軌道最適化によるSQD推定値の改善\"\n",
        "description: \"最新バージョンのサンプルベース量子対角化（SQD）における軌道最適化を用いたSQD推定値の改善\"\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "bb5a576d",
      "metadata": {},
      "source": [
        "<span id=\"improve-an-sqd-estimate-with-orbital-optimization\" />\n",
        "\n",
        "# 軌道最適化によるSQD推定値の改善\n",
        "\n",
        "サンプルベース量子対角化（SQD）は、\n",
        "電子配置の固定された部分空間においてハミルトニアンを対角化することで、基底状態のエネルギーを近似する。 その\n",
        "推定値は、ハミルトニアンが表現される軌道基底に依存しており、\n",
        "*軌道最適化* （OO）はこの自由度を利用して、\n",
        "部分空間を拡大することなくエネルギーを低減させる。\n",
        "\n",
        "このガイドでは、 $N_2$ 分子に対してSQDを実行し、その後、軌道\n",
        "最適化を用いて結果を改善します。ここでは、\n",
        "ハミルトニアンを表現し[`ffsim`](https://qiskit-community.github.io/ffsim/)、エネルギーを最小化する軌道回転を求めるためにを使用します。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "77ab953b",
      "metadata": {},
      "source": [
        "<span id=\"run-sqd\" />\n",
        "\n",
        "## SQDを実行する\n",
        "\n",
        "分子軌道（MO）基底において、 $N_2$ の分子積分を構築し、\n",
        "一様乱数サンプルを生成した上で、SQDを実行して基底状態の近似を求める。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "b8d5618e",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-07-16T01:36:48.334261Z",
          "iopub.status.busy": "2026-07-16T01:36:48.334063Z",
          "iopub.status.idle": "2026-07-16T01:37:42.585853Z",
          "shell.execute_reply": "2026-07-16T01:37:42.584303Z"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "converged SCF energy = -108.835236570775\n",
            "CASCI E = -109.046671778080  E(CI) = -32.8155692383187  S^2 = 0.0000000\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "import pyscf\n",
        "import pyscf.cc\n",
        "import pyscf.mcscf\n",
        "from qiskit_addon_sqd.counts import generate_bit_array_uniform\n",
        "from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian\n",
        "\n",
        "# Specify molecule properties\n",
        "num_orbitals = 16\n",
        "num_elec_a = num_elec_b = 5\n",
        "spin_sq = 0\n",
        "\n",
        "# Build N2 molecule\n",
        "mol = pyscf.gto.Mole()\n",
        "mol.build(\n",
        "    atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]],\n",
        "    basis=\"6-31g\",\n",
        "    symmetry=\"Dooh\",\n",
        ")\n",
        "\n",
        "# Define active space\n",
        "n_frozen = 2\n",
        "active_space = range(n_frozen, mol.nao_nr())\n",
        "\n",
        "# Get molecular integrals\n",
        "scf = pyscf.scf.RHF(mol).run()\n",
        "num_orbitals = len(active_space)\n",
        "n_electrons = int(sum(scf.mo_occ[active_space]))\n",
        "num_elec_a = (n_electrons + mol.spin) // 2\n",
        "num_elec_b = (n_electrons - mol.spin) // 2\n",
        "cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))\n",
        "mo = cas.sort_mo(active_space, base=0)\n",
        "hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n",
        "eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)\n",
        "\n",
        "# Compute exact energy\n",
        "exact_energy = cas.run().e_tot\n",
        "\n",
        "# Create a seed to control randomness throughout this workflow\n",
        "rng = np.random.default_rng(24)\n",
        "\n",
        "\n",
        "# Generate random samples\n",
        "bit_array = generate_bit_array_uniform(\n",
        "    10_000, num_orbitals * 2, rand_seed=rng\n",
        ")\n",
        "\n",
        "# Run SQD\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    hcore,\n",
        "    eri,\n",
        "    bit_array,\n",
        "    samples_per_batch=100,\n",
        "    norb=num_orbitals,\n",
        "    nelec=(num_elec_a, num_elec_b),\n",
        "    num_batches=1,\n",
        "    max_iterations=5,\n",
        "    symmetrize_spin=True,\n",
        "    seed=rng,\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "0ca70028",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-07-16T01:37:42.590790Z",
          "iopub.status.busy": "2026-07-16T01:37:42.589477Z",
          "iopub.status.idle": "2026-07-16T01:37:42.595563Z",
          "shell.execute_reply": "2026-07-16T01:37:42.595133Z"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Exact energy:  -109.04667178\n",
            "SQD energy:    -108.98469255\n"
          ]
        }
      ],
      "source": [
        "sqd_energy = result.energy + nuclear_repulsion_energy\n",
        "print(f\"Exact energy:  {exact_energy:.8f}\")\n",
        "print(f\"SQD energy:    {sqd_energy:.8f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9160d09c",
      "metadata": {},
      "source": [
        "<span id=\"optimize-the-orbitals\" />\n",
        "\n",
        "## 軌道を最適化する\n",
        "\n",
        "軌道最適化とは、変分エネルギーを低減させる軌道回転を探すことである\n",
        "\n",
        "$$\n",
        "E = \\langle \\psi | \\mathcal{U}^\\dagger\\, H\\, \\mathcal{U} | \\psi \\rangle\n",
        "$$\n",
        "\n",
        "SQD基底状態近似 $|\\psi\\rangle$ に基づく。軌道回転は、\n",
        "$N \\times N$ ユニタリ行列 $\\mathbf{U}$ （ $N$ は空間軌道の数）によって指定され、\n",
        "これは演算子を通じて多体状態に作用する。\n",
        "\n",
        "$$\n",
        "\\mathcal{U} = \\exp\\left[\\sum_{pq, \\sigma} \\log(\\mathbf{U})_{pq}\\, a^\\dagger_{p\\sigma} a_{q\\sigma}\\right].\n",
        "$$\n",
        "\n",
        "`ffsim.optimize_orbitals` **行列**$\\mathbf{U}$ を返し、これを\n",
        "軌道基底に（ via を用いて `hamiltonian.rotated`）適用することは、 $\\mathcal{U}$ を\n",
        "その状態に適用することと同等である。 詳細については、\n",
        "[ffsimの軌道回転に関する説明](https://qiskit-community.github.io/ffsim/explanations/orbital-rotation.html)\n",
        "を参照してください。\n",
        "\n",
        "軌道回転を行うと、その部分空間から見てハミルトニアンが変化するため、\n",
        "エネルギーがこれ以上改善しなくなるまで、以下の2つの手順を交互に繰り返します：\n",
        "\n",
        "1. 固定された一連の\n",
        "   構成について、現在の基底におけるハミルトニアンを**対角化する**。\n",
        "2. 結果として得られる状態のエネルギーを最小化する回転を見つけ、 **軌道関数を最適化し**、\n",
        "   その後、積分を新しい基底に変換する。\n",
        "\n",
        "軌道回転のステップは、\n",
        "状態の1体および2体の縮約密度行列（RDM）から、\n",
        "エネルギーを最小化する回転を求める処理に\n",
        "[`ffsim.optimize_orbitals`](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.optimize_orbitals)委ねます。 「第○条」を参照\n",
        "[。 詳細はII A 4](https://arxiv.org/pdf/2405.05068) を参照のこと。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "006745b5",
      "metadata": {},
      "source": [
        "<span id=\"why-orbital-optimization-helps-here\" />\n",
        "\n",
        "### なぜこの場面で軌道最適化が役立つのか\n",
        "\n",
        "SCF分子軌道（MO）基底は、*完全な*-CI問題において軌道回転に対して静止している。 しかし、SQDは小さな切り詰め部分空間（ここでは、およそ1,900万個の完全CI行列式のうち、\n",
        "数百個のCI文字列）で動作するため、その部分空間ではMO\n",
        "基底は一般的に最適ではない。そのため、軌道を回転させると、その部分空間が\n",
        "表現できるエネルギーが低下してしまう。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "46eb2f5a",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-07-16T01:37:42.597573Z",
          "iopub.status.busy": "2026-07-16T01:37:42.597399Z",
          "iopub.status.idle": "2026-07-16T01:37:42.654809Z",
          "shell.execute_reply": "2026-07-16T01:37:42.654267Z"
        }
      },
      "outputs": [],
      "source": [
        "import ffsim\n",
        "from pyscf import fci\n",
        "\n",
        "# ffsim's ``MolecularHamiltonian`` uses the same \"chemist\" ordering for the two-body\n",
        "# tensor as PySCF's ``eri``, and stores the nuclear repulsion energy as the constant\n",
        "# term so that expectation values come out as total energies."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "93179dc4",
      "metadata": {},
      "source": [
        "<span id=\"alternate-diagonalization-and-orbital-optimization\" />\n",
        "\n",
        "### 交互対角化と軌道最適化\n",
        "\n",
        "SQDによって上記で発見された構成に対して、対角化部分空間を**一定に**保つことで、\n",
        "各反復計算において軌道回転の影響のみを分離できるようにしています。 各\n",
        "反復ごとに：\n",
        "\n",
        "1. PySCF'sの選択CIソルバーを用いて、現在の基底における固定部分空間上でハミルトニアンを**対角化します**。\n",
        "\n",
        "2. 結果として得られる状態の**RDMを構築します**\n",
        "   。これだけで十分`ffsim.optimize_orbitals`\n",
        "   です。\n",
        "\n",
        "3. **軌道を最適化します** ：エネルギーを最小化する\n",
        "   回転を返し`ffsim.optimize_orbitals`、これを積分に適用して、改良された基底に移行させます。\n",
        "\n",
        "各最適化ステップ*の前に*、エネルギー値を記録します。 基底は反復ごとに\n",
        "改善されるため、この数列は、固定された部分空間において達成可能な最良のエネルギーに向かって\n",
        "単調に減少する。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "a0783e88",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-07-16T01:37:42.657394Z",
          "iopub.status.busy": "2026-07-16T01:37:42.657213Z",
          "iopub.status.idle": "2026-07-16T01:38:48.374706Z",
          "shell.execute_reply": "2026-07-16T01:38:48.373697Z"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Iteration 0: energy = -108.98452447\n",
            "Iteration 1: energy = -108.99981993\n",
            "Iteration 2: energy = -109.00585329\n",
            "Iteration 3: energy = -109.00816569\n",
            "Iteration 4: energy = -109.00936616\n",
            "Iteration 5: energy = -109.01014322\n",
            "Iteration 6: energy = -109.01069439\n",
            "Iteration 7: energy = -109.01109308\n",
            "Iteration 8: energy = -109.01138928\n",
            "Iteration 9: energy = -109.01161411\n"
          ]
        }
      ],
      "source": [
        "# Fix the diagonalization subspace to the configurations found by SQD.\n",
        "ci_strings = (result.sci_state.ci_strs_a, result.sci_state.ci_strs_b)\n",
        "nelec = (num_elec_a, num_elec_b)\n",
        "\n",
        "# Start from the MO basis in which we ran SQD.\n",
        "hamiltonian_opt = ffsim.MolecularHamiltonian(\n",
        "    hcore, eri, constant=nuclear_repulsion_energy\n",
        ")\n",
        "\n",
        "num_iters = 10\n",
        "for i in range(num_iters):\n",
        "    # Diagonalize over the fixed subspace in the current basis.\n",
        "    myci = fci.selected_ci.SelectedCI()\n",
        "    myci = fci.addons.fix_spin_(myci, ss=spin_sq)\n",
        "    _, amplitudes = fci.selected_ci.kernel_fixed_space(\n",
        "        myci,\n",
        "        hamiltonian_opt.one_body_tensor,\n",
        "        hamiltonian_opt.two_body_tensor,\n",
        "        num_orbitals,\n",
        "        nelec,\n",
        "        ci_strs=ci_strings,\n",
        "    )\n",
        "\n",
        "    # Build the RDMs and record the energy before re-optimizing the orbitals.\n",
        "    dm1, dm2 = myci.make_rdm12(amplitudes, num_orbitals, nelec)\n",
        "    rdm = ffsim.ReducedDensityMatrix(dm1, dm2)\n",
        "    energy = rdm.expectation(hamiltonian_opt).real\n",
        "    print(f\"Iteration {i}: energy = {energy:.8f}\")\n",
        "\n",
        "    # Rotate the Hamiltonian into the energy-minimizing basis for the next iteration.\n",
        "    # optimize_orbitals returns the unitary matrix U minimizing\n",
        "    # rdm.rotated(U).expectation(hamiltonian), equivalently\n",
        "    # rdm.expectation(hamiltonian.rotated(U.conj().T)), so we rotate by U^dagger.\n",
        "    orbital_rotation = ffsim.optimize_orbitals(rdm, hamiltonian_opt)\n",
        "    hamiltonian_opt = hamiltonian_opt.rotated(orbital_rotation.T.conj())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a61dbda8",
      "metadata": {},
      "source": [
        "<span id=\"compare-the-results\" />\n",
        "\n",
        "### 結果の比較\n",
        "\n",
        "軌道最適化により、固定部分空間による推定値が改善され、\n",
        "正確なエネルギーとの差の大部分が縮まり、かつその値を上回る状態が維持される。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "762b0903",
      "metadata": {
        "execution": {
          "iopub.execute_input": "2026-07-16T01:38:48.377620Z",
          "iopub.status.busy": "2026-07-16T01:38:48.377404Z",
          "iopub.status.idle": "2026-07-16T01:38:53.111168Z",
          "shell.execute_reply": "2026-07-16T01:38:53.110312Z"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Exact energy:      -109.04667178\n",
            "SQD energy (MO):   -108.98469255\n",
            "Energy after OO:   -109.01178727\n"
          ]
        }
      ],
      "source": [
        "# Diagonalize once more in the final optimized basis to report the improved energy.\n",
        "myci = fci.selected_ci.SelectedCI()\n",
        "myci = fci.addons.fix_spin_(myci, ss=spin_sq)\n",
        "_, amplitudes = fci.selected_ci.kernel_fixed_space(\n",
        "    myci,\n",
        "    hamiltonian_opt.one_body_tensor,\n",
        "    hamiltonian_opt.two_body_tensor,\n",
        "    num_orbitals,\n",
        "    nelec,\n",
        "    ci_strs=ci_strings,\n",
        ")\n",
        "dm1, dm2 = myci.make_rdm12(amplitudes, num_orbitals, nelec)\n",
        "energy_after_oo = (\n",
        "    ffsim.ReducedDensityMatrix(dm1, dm2).expectation(hamiltonian_opt).real\n",
        ")\n",
        "\n",
        "print(f\"Exact energy:      {exact_energy:.8f}\")\n",
        "print(f\"SQD energy (MO):   {sqd_energy:.8f}\")\n",
        "print(f\"Energy after OO:   {energy_after_oo:.8f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}