Add notebook drafts (#6)

notebooks/fanci_tutorial.ipynb

@@ -0,0 +1,216 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "FanCI is a Python library for computing various post-Hartree-Fock methods\n",
+    "(of the Flexible Ansätze for N-electron CI, or \"FanCI\" type) using the PyCI library.\n",
+    "\n",
+    "Here we will look at the following wavefunction models using **Be** as our toy model:\n",
+    "* AP1roG\n",
+    "* pCCD+S"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Lets start showing how to run an AP1roG computation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Import modules\n",
+    "import pyci\n",
+    "from pyci.fanci import AP1roG, pCCDS\n",
+    "# optional\n",
+    "import numpy as np\n",
+    "from pyci.test import datafile"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Define our system"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# System information\n",
+    "filename = datafile(\"be_ccpvdz.fcidump\")\n",
+    "ham = pyci.hamiltonian(filename)\n",
+    "e_dict = {}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compute the Hartree-Fock energy for comparison"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Get Hartree-Fock energy\n",
+    "hf_wfn = pyci.doci_wfn(ham.nbasis, 2, 2)\n",
+    "hf_wfn.add_hartreefock_det()\n",
+    "hf_op = pyci.sparse_op(ham, hf_wfn)\n",
+    "e_dict[\"HF\"] = hf_op.solve(n=1)[0][0] - ham.ecore"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Optimize the AP1roG wave function"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Initialize AP1roG instance\n",
+    "ap1rog = AP1roG(ham, 2)\n",
+    "\n",
+    "# Make initial guess\n",
+    "ap1_params = np.zeros(ap1rog.nparam, dtype=pyci.c_double)\n",
+    "ap1_params[-1] = e_dict[\"HF\"]\n",
+    "\n",
+    "# Optimize wavefunction\n",
+    "ap1rog_results = ap1rog.optimize(ap1_params, use_jac=True)\n",
+    "e_dict[\"AP1roG\"] = ap1rog_results.x[-1]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Use the AP1roG wave function to generate an initial guess for the pCCD+S wave function; optimize it"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Initialize pCCD+S instance\n",
+    "pccds = pCCDS(ham, 2, 2)\n",
+    "\n",
+    "# Make initial guess from AP1roG params\n",
+    "pccds_params = np.zeros(pccds.nparam, dtype=pyci.c_double)\n",
+    "pccds_params[:pccds.wfn.nocc_up * pccds.wfn.nvir_up] = ap1rog_results.x[:-1]\n",
+    "pccds_params[-1] = ap1rog_results.x[-1]\n",
+    "\n",
+    "# Optimize wavefunction\n",
+    "pccds_results = pccds.optimize(pccds_params, use_jac=False)\n",
+    "e_dict[\"pCCD+S\"] = pccds_results.x[-1]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compare the energies and verify that  $E_\\text{HF}$ > $E_\\text{AP1roG}$ > $E_\\text{pCCD+S}$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Print energies from various methods\n",
+    "print(f\"{'METHOD':6s}, {'ENERGY':16s}\")\n",
+    "for name, energy in e_dict.items():\n",
+    "    print(f\"{name:6s}, {energy:16.9e}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Use a larger projection space (seniority-(0,2) + excitation-(0,1,2,3) CI) to approximately evaluate the pCCD+S RDMs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Make larger projection space over which to evaluate RDMs\n",
+    "proj_wfn = pyci.fullci_wfn(ham.nbasis, 2, 2)\n",
+    "# Add seniority 0 and 2 determinants\n",
+    "pyci.add_seniorities(proj_wfn, 0, 2)\n",
+    "# Add singly, doubly, triply excited determinants\n",
+    "pyci.add_excitations(proj_wfn, 0, 1, 2, 3)\n",
+    "# Evaluate coefficients with optimized pCCD+S\n",
+    "coeffs = pccds.compute_overlap(pccds_results.x, proj_wfn.to_occs_array())\n",
+    "\n",
+    "# Compute RDMs using larger projection space and pCCD+S coefficients\n",
+    "d1, d2 = pyci.compute_rdms(proj_wfn, coeffs)\n",
+    "rdm1, rdm2 = pyci.spinize_rdms(d1, d2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Validate the RDMs; ensure that  $\\text{tr}\\left(\\gamma\\right) = n_\\text{elec}$  and that  $\\frac{1}{2}\\sum_{pq}{\\Gamma_{pqpq}} = n_\\text{pair}$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Validate the RDMs\n",
+    "nelec = np.trace(rdm1)\n",
+    "npairs = np.einsum('pqpq', rdm2) / 2.0\n",
+    "print(f\"Number of electrons = {nelec:.1f}\")\n",
+    "print(f\"Number of pairs     = {npairs:.1f}\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.11.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}