{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "frontmatter",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Scale SQD chemistry workflows with Dice solver\"\n",
        "description: \"Scale SQD chemistry workflows with Dice solver for the latest version of Sample-based quantum diagonalization (SQD)\"\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "65adf860-4443-494f-ad90-5d8ea444ac4e",
      "metadata": {},
      "source": [
        "# Scale SQD chemistry workflows with Dice solver\n",
        "\n",
        "{/* cspell:ignore hcore */}\n",
        "\n",
        "This guide shows how to use the Dice solver as an alternative SCI solver to diagonalize larger Fermionic problems beyond what the default PySCF solver can support. For information on how to install and use `qiskit-addon-dice-solver`, [visit the documentation](https://qiskit.github.io/qiskit-addon-dice-solver/). For more details on the SQD code used in this example, see the [chemistry Hamiltonian tutorial](/docs/tutorials/sample-based-quantum-diagonalization).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "baf7d074-4443-4695-875e-65f7fb8cef25",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "converged SCF energy = -108.835236570774\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_dice_solver import solve_sci_batch\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=300,\n",
        "    norb=num_orbitals,\n",
        "    nelec=(num_elec_a, num_elec_b),\n",
        "    num_batches=5,\n",
        "    max_iterations=5,\n",
        "    sci_solver=solve_sci_batch,\n",
        "    symmetrize_spin=True,\n",
        "    seed=rng,\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "99709a88-f9ec-4dbc-a561-e682f90aa4c5",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Exact energy: -109.04667177808028\n",
            "Estimated energy: -109.03402667558743\n"
          ]
        }
      ],
      "source": [
        "print(f\"Exact energy: {exact_energy}\")\n",
        "print(f\"Estimated energy: {result.energy + nuclear_repulsion_energy}\")"
      ]
    },
    {
      "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
}