{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "461ae8c6-4564-43e1-8e6a-6dde924e0c74",
   "metadata": {},
   "source": [
    "# Search with fallback mode\n",
    "\n",
    "Active space selection with the ASF package is largely based on pair information criteria derived from the two-electron cumulant, augmented with single-orbital entropies. We designed procedures to select active spaces that are likely to converge in CASSCF calculations. To fulfill the latter criterion, it helps to omit weakly entangled orbitals from the active space. For some systems, however, correlation of relevant orbitals may be relatively weak; in particular, if those systems lack any pronounced static correlation. The following examples show how to obtain a reasonable active space by controlling and relaxing the search criteria via fallback modes.\n",
    "\n",
    "The example covers:\n",
    "- searching an active space using `find_one_entropy` with varying fallback modes\n",
    "- searching active spaces using `find_many` with a fallback option\n",
    "- visualizing results with Jmol\n",
    "\n",
    "Before turning to the example calculations, we briefly review the selection process of `find_many`, `find_one_sized`, and `find_one_entropy`.\n",
    "The ASF classes are able to generate a large set of candidate spaces (\"unfiltered spaces\"); only a subset of those will converge at all in a CASSCF calculation, and only a tiny subset will converge to a solution that is consistent with the guess (which is the desirable outcome). Procedures in the ASF aim to identify such a tiny subset. Search functions, such as `find_one_entropy` or `find_many`, first narrow down the large unfiltered set to a collection of reasonable spaces by applying the following three criteria:\n",
    "\n",
    "  1) Removing any orbital from the active space lowers the amount of correlation in the active space more strongly than adding any orbital to the active space increases it. This is quantified by the smallest decrement `min_decrement` exceeding the largest increment `max_increment` (to within a comparison tolerance).\n",
    "  2) It is not possible to increase the amount of correlation contained in an active space by swapping out an active orbital for an inactive orbital (while keeping the number of active orbitals and electrons constant). Numerically, for each group of active spaces (N, M) with the same number of active orbitals and active electrons, a space from that group is kept if it includes the orbitals with the largest node weight to within a comparison tolerance. The considered node weight may be the minimal single-orbital entropy, the minimal edge sum (minimal sum of all edge values associated with an orbital node), or no node weight (`\"none\"`) which is equivalent of disabling this criterion.\n",
    "  3) Too weakly correlated orbitals are omitted from the active space. The associated requirement is that the minimal entropy of any orbital must be above the given minimal entropy threshold (default: 0.07).\n",
    "\n",
    "We make a further distinction between orbitals that are members of a minimal space always needed to accommodate unpaired electrons, and further (\"exchangeable\") MOs. The reason for this is that orbitals needed to host unpaired electrons may be weakly correlated in some systems, but they always have to be included in the active space; and despite being weakly correlated, they are unable to rotate out of the active space in CASSCF calculations as a consequence of hosting the spin density. Generally, the selection criteria set out above are applied to \"exchangeable\" MOs that are not needed to accomodate unpaired electrons. Minimal space MOs that are identified as hosting the unpaired electrons are exempt. An exception from this exception is the minimal space itself, which has a dedicated flag and contains the actual values, not the values determined for the exchangeable MOs only.\n",
    "\n",
    "The `find_one` variants then apply additional criteria to select a single space, while `find_many` returns the collection of reasonable spaces directly.\n",
    "\n",
    "\n",
    "## Example 1: Hexaaqua-Iron(III) cation\n",
    "\n",
    "The first example is the hexaaqua-iron(III) complex in its sextet ground state using a molecular structure taken from *Phys. Chem. Chem. Phys. 21*, 4854-4870, (**2019**) [[doi:10.1039/C9CP00105K]](https://doi.org/10.1039/C9CP00105K), which was optimized in the gas phase for the sextet ground state at the BP86-D3/def2-TZVP level of theory. The complex has ideal octahedral symmetry with all Fe-O distances measuring 2.06 Å."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11b148c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Suppress PySCF warning...\n",
    "import pyscf\n",
    "pyscf.__config__.B3LYP_WITH_VWN5 = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b2d20c2-7148-4bf9-a6e5-d4760d4ddb91",
   "metadata": {},
   "outputs": [],
   "source": [
    "from asf.asfbase import print_mo_table\n",
    "from asf.filters import ActiveSpaceSelectionError\n",
    "from asf.scf import stable_scf\n",
    "from asf.preselection import MP2PairinfoPreselection, MP2NatorbPreselection\n",
    "from asf.utility import pictures_Jmol, show_mos_grid\n",
    "from asf import ASFDMRG, ASFCI\n",
    "from pathlib import Path\n",
    "from pprint import pprint\n",
    "from pyscf.gto import M"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "671feca4-da48-498a-a95a-7eee180f64bf",
   "metadata": {},
   "source": [
    "### 1.1. Setting up the molecular system\n",
    "\n",
    "The ASF package assumes that molecular systems are represented as PySCF Mole objects. The first step is therefore to initialize the hexaaqua-iron(III) complex as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf51d818",
   "metadata": {},
   "outputs": [],
   "source": [
    "mol = M(\n",
    "    atom=\"\"\"\n",
    "    Fe   -0.0000000    0.0000000    0.0000000\n",
    "    O     0.0000000    2.0622910    0.0000000\n",
    "    H     0.7919274    2.6471973    0.0000000\n",
    "    H    -0.7919274    2.6471973    0.0000000\n",
    "    O    -0.0000000    0.0000000    2.0622910\n",
    "    H    -0.0000000    0.7919274    2.6471973\n",
    "    H     0.0000000   -0.7919274    2.6471973\n",
    "    O     2.0622910   -0.0000000   -0.0000000\n",
    "    H     2.6471973   -0.0000000    0.7919274\n",
    "    H     2.6471973   -0.0000000   -0.7919274\n",
    "    O    -0.0000000   -2.0622910    0.0000000\n",
    "    H    -0.7919274   -2.6471973   -0.0000000\n",
    "    H     0.7919274   -2.6471973   -0.0000000\n",
    "    O     0.0000000    0.0000000   -2.0622910\n",
    "    H     0.0000000   -0.7919274   -2.6471973\n",
    "    H    -0.0000000    0.7919274   -2.6471973\n",
    "    O    -2.0622910    0.0000000    0.0000000\n",
    "    H    -2.6471973    0.0000000   -0.7919274\n",
    "    H    -2.6471973    0.0000000    0.7919274\n",
    "    \"\"\",\n",
    "    charge=3,\n",
    "    spin=5,\n",
    "    basis=\"def2-svp\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "482d2cc6-3e6e-4224-af1f-57232dc539b4",
   "metadata": {},
   "source": [
    "### 1.2. Initializing the ASF class\n",
    "\n",
    "Next, we prepare an ASF class instance that will later compute the single-orbital entropies and pair information based on an approximate DMRG-CI calculation. We start with an unrestricted Hartree-Fock calculation using the convenient `stable_scf` function from the `asf.scf` module. Using the `from_preselection` method, the selection of the initial space and instantiation of the `ASFDMRG` is carried out with one function call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6afa97dc-4e1b-4e79-9ade-f68ccfc05776",
   "metadata": {},
   "outputs": [],
   "source": [
    "mf = stable_scf(mol)\n",
    "sf = ASFDMRG.from_preselection(MP2NatorbPreselection(mf), maxM=150)\n",
    "print(f\"-> Selected an initial space of {sf.nel} electrons and {sf.norb} orbitals.\")\n",
    "sf.calculate()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "841175fa-ad30-418a-8f96-8066456d7959",
   "metadata": {},
   "source": [
    "Inspecting the single-orbital densities (w) and single-orbital entropy (S) for the initial space, we notice that *all* orbitals display very low entropies. For comparison, the default entropy threshold applied by `find_one_entropy` is ca. 0.139 and the minimal threshold applied to select the reasonable set of spaces is 0.07. We can therefore expect that at least one of the three criteria mentioned above is not met and the collection of reasonable spaces is empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dee9bd10-8baf-4595-bad7-75b32fc6908f",
   "metadata": {},
   "outputs": [],
   "source": [
    "print_mo_table(orbital_densities=sf.one_orbital_density(), mo_list=sf.mo_list.tolist(), selections={\"initial\": sf.mo_list})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fd57b51-ab33-448d-b8b8-61577de9d30d",
   "metadata": {},
   "source": [
    "### 1.3. Searching for multiple spaces\n",
    "\n",
    "Following the explanations outlined above, it is unsurprising if `find_many` does not suggest any non-trivial active space. Indeed, this is the case if we strictly apply the three criteria mentioned above without exemption for the minimal space, i.e. with `keep_minimal=False`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c18319f-7feb-46a5-914a-47ae697aa8c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "spaces_strict = sf.find_many(keep_minimal=False)\n",
    "print(f\"Active spaces matching criteria 1) - 3): {spaces_strict}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b936e75-e164-4f94-b401-ea07711866e0",
   "metadata": {},
   "source": [
    "Running `find_many` with the default settings, i.e. `keep_minimal=True`, the minimal space is exempt from the criteria 1) to 3) outlined previously. For open-shell systems with *N* unpaired electrons it is always possible to construct an *(N, N)* minimal space, containing *N* electrons in *N* orbitals: using it to perform a CASSCF calculation is, effectively, equivalent with restricted open-shell Hartree-Fock (ROHF). Therefore, the default behavior of `find_many` is to preserve this *(N, N)* space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae319bc3-0d5a-41f5-8262-81f464e2c8b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "spaces_minimal = sf.find_many()\n",
    "print(f\"Active spaces matching criteria 1) - 3) including minimal space:\")\n",
    "pprint(spaces_minimal)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd0932ca-4967-4ba2-bc2f-16fd48e9cb08",
   "metadata": {},
   "source": [
    "### 1.4. Searching for a single active space\n",
    "\n",
    "The method `find_one_entropy` invokes the same procedure to extract sensible spaces from the initial unfiltered collection as `find_many`, so no active space is suggested if the search criteria are applied in a strict manner. To obtain a reasonable alternative in such cases, `find_one_entropy` provides several fallback modes with slightly altered search parameters. These fallback criteria follow a hierarchy outlined below:\n",
    "\n",
    "1) `\"strict\"` applies the original search criteria: among the set of reasonable spaces containing only MOs with entropies above a threshold (0.139 by default), it selects the one that maximizes the pair information sum. \n",
    "2) `\"below_threshold\"` chooses among spaces that were discarded with the strict entropy threshold. Among the collection of discarded but otherwise sensible spaces, it identifies the one where the lowest entropy of any MO is as large as possible. If multiple spaces have identical values of the minimal entropy, it selects the one among them that maximizes the pair information sum.\n",
    "3) `\"minimal\"` always includes the minimal space in set of reasonable spaces; criterion `\"below_threshold\"` excludes the minimal space, unless it satisfies the criteria for filtering of sensible spaces.\n",
    "4) `\"unfiltered\"` extends the search to the large set of unfiltered spaces if the selection attempt with the reasonable spaces was unsuccessful. Similarly to options 1)-3) in this enumeration, it will first try to select an unfiltered space with a minimal entropy above the given threshold that maximizes the pair information sum and if it fails, then attempts to find the unfiltered space with the highest minimal entropy below the entropy threshold.\n",
    "\n",
    "The `fallback_mode` argument determines to which point this sequence should be executed if the \"strict\" search does not yield a result. By default it is set to \"minimal\", i.e. at most the \"minimal\" fallback search is run before an exception is thrown. `\"unfiltered\"` is an expert option that is disabled by default, as the active spaces returned thereby are not at all likely to lead to converging CASSCF calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "322fc3fd-b0bd-4701-99e9-d8ff4963c5fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "FALLBACK_MODE = \"strict\"\n",
    "\n",
    "try:\n",
    "    selected = sf.find_one_entropy(fallback_mode=FALLBACK_MODE)\n",
    "    print(selected)\n",
    "except ActiveSpaceSelectionError:\n",
    "    print(f\"-> Could not determine active space with {FALLBACK_MODE=}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c88d169-7c57-43fd-8d13-18e92cf17d3d",
   "metadata": {},
   "source": [
    "Since the collection of reasonable spaces is empty, also the subsequent fallback search below the entropy threshold does not yield a result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92ecd3f1-03a3-49b8-8f0d-876892bebdff",
   "metadata": {},
   "outputs": [],
   "source": [
    "FALLBACK_MODE = \"below_threshold\"\n",
    "\n",
    "try:\n",
    "    selected = sf.find_one_entropy(fallback_mode=FALLBACK_MODE)\n",
    "    print(selected)\n",
    "except ActiveSpaceSelectionError:\n",
    "    print(f\"-> Could not determine active space with {FALLBACK_MODE=}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11992bcc-129b-4f97-ac5e-579d0235ba3e",
   "metadata": {},
   "source": [
    "Once the \"minimal\" fallback mode is activated, the set of reasonable spaces has one member, the minimal space, which is returned as the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "315662dc-d087-40c2-bc05-9b6070945971",
   "metadata": {},
   "outputs": [],
   "source": [
    "FALLBACK_MODE = \"minimal\"\n",
    "\n",
    "try:\n",
    "    selected = sf.find_one_entropy(fallback_mode=FALLBACK_MODE)\n",
    "    print(selected)\n",
    "except ActiveSpaceSelectionError:\n",
    "    print(f\"-> Could not determine active space with {FALLBACK_MODE=}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47a57849-599b-4b17-9a6f-3d889fb84345",
   "metadata": {},
   "source": [
    "### 1.5. Visualizing the results\n",
    "\n",
    "To visualize the suggested (5, 5) space, we generate molecular orbital images with [Jmol](https://chemapps.stolaf.edu/jmol/) via the `pictures_Jmol` utility function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d512bc20-5812-4ace-9173-e4a5ff273731",
   "metadata": {},
   "outputs": [],
   "source": [
    "# If plots have been generated in a previous run of the Notebook, they will be reused here.\n",
    "# Set redo_plots = True, if new orbital images should be created regardless of existing files.\n",
    "redo_plots = True\n",
    "\n",
    "def mo_images(name_glob: str = \"mo_*.png\") -> list[Path]:\n",
    "    \"\"\"Return the list of all MO images in this directory.\"\"\"\n",
    "    return list(Path.cwd().glob(name_glob))\n",
    "\n",
    "if not mo_images() or redo_plots:\n",
    "    pictures_Jmol(\n",
    "        mol=mol,\n",
    "        mo_coeff=sf.mo_coeff,\n",
    "        mo_list=selected.mo_list,\n",
    "        rotate=(60, 45, 0),\n",
    "        mo_index_title=\" \",\n",
    "        zoom=100,\n",
    "    )\n",
    "show_mos_grid(images=mo_images(), mo_list=selected.mo_list, columns=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "501d45fb-b796-4e59-88b9-ce2643fbd217",
   "metadata": {},
   "source": [
    "## Example 2: Hydroxyl radical\n",
    "\n",
    "The hydroxyl radical is another case where the strict application entropy-based search criteria does not lead to an active space suggestion. It serves as an example why the ASF applies not one but several fallback modes for `find_one_entropy`. \n",
    "\n",
    "### 2.1. Setting up the molecular system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74ca66d9-79bc-4b9d-bddd-d70bb25dcd2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "mol2 = M(\n",
    "    atom=\"\"\"\n",
    "    O     0.0000000    0.0000000    0.0000000 \n",
    "    H     0.0000000    0.0000000    0.9700000\n",
    "    \"\"\",\n",
    "    charge=0,\n",
    "    spin=1,\n",
    "    basis=\"def2-svp\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b9fd813",
   "metadata": {},
   "source": [
    "### 2.2. Initializing the ASF class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ad9ad15-5612-42fe-a3f8-b2ba0c51d6e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mf2 = stable_scf(mol2)\n",
    "sf2 = ASFCI.from_preselection(MP2PairinfoPreselection(mf2))\n",
    "print(f\"Selected an initial space of {sf2.nel} electrons and {sf2.norb} orbitals.\")\n",
    "sf2.calculate()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "563bec81-31ff-41aa-977f-f003820666e7",
   "metadata": {},
   "source": [
    "From the single-orbital densities and single-orbital entropy associated with the initial space, we notice that no entropy is above the default threshold (0.139) but most entropies are above the minimal entropy threshold used to select the sensible spaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d44fe786-9ccd-42ec-a24f-504b9af1c609",
   "metadata": {},
   "outputs": [],
   "source": [
    "print_mo_table(orbital_densities=sf2.one_orbital_density(), mo_list=sf2.mo_list.tolist(), selections={\"initial\": sf2.mo_list})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d135e4e1-cbe0-49d7-b8e2-ffda8dbf2313",
   "metadata": {},
   "source": [
    "### 2.3. Searching for a single active space\n",
    "\n",
    "As expected, applying the search criteria without a fallback does not yield a suggestion for an active space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7de8b642-e41b-4a65-b2b7-c7a03470f93e",
   "metadata": {},
   "outputs": [],
   "source": [
    "FALLBACK_MODE = \"strict\"\n",
    "\n",
    "try:\n",
    "    selected2 = sf2.find_one_entropy(fallback_mode=FALLBACK_MODE)\n",
    "    print(selected2)\n",
    "except ActiveSpaceSelectionError:\n",
    "    print(f\" - Could not determine active space with {FALLBACK_MODE=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7192dbed-640e-41a6-99c5-4f386dc38a9f",
   "metadata": {},
   "source": [
    "With the subsequent fallback mode, `\"below_threshold\"`, there are two possible spaces to choose from (note that `find_many` is called here to show the sensible spaces below the entropy threshold):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fafc9b7-e54d-4930-b42f-85befce06b57",
   "metadata": {},
   "outputs": [],
   "source": [
    "pprint(sf2.find_many(keep_minimal=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbd688e3-5197-4726-a36e-aad97b8595e7",
   "metadata": {},
   "source": [
    "Since this fallback option selects a space below the threshold where the lowest entropy of any MO is largest (as opposed to maximizing the pair information above the threshold), we obtain the (3, 3) space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24b1384a-7509-455e-b439-de7b10d03366",
   "metadata": {},
   "outputs": [],
   "source": [
    "FALLBACK_MODE = \"below_threshold\"\n",
    "\n",
    "try:\n",
    "    selected2 = sf2.find_one_entropy(fallback_mode=FALLBACK_MODE)\n",
    "    print(selected2)\n",
    "except ActiveSpaceSelectionError:\n",
    "    print(f\" - Could not determine active space with {FALLBACK_MODE=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b102856f-4dd1-4c67-a00a-f84d1bbc8daa",
   "metadata": {},
   "source": [
    "### 2.4. Visualizing the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ca8c414-fa4a-4995-ac7e-404c74c72505",
   "metadata": {},
   "outputs": [],
   "source": [
    "redo_plots = True\n",
    "\n",
    "if not mo_images(\"oh_mo_*.png\") or redo_plots:\n",
    "    pictures_Jmol(\n",
    "        mol=mol2,\n",
    "        mo_coeff=sf2.mo_coeff,\n",
    "        mo_list=selected2.mo_list,\n",
    "        rotate=(90, 45, 90),\n",
    "        name=\"oh_mo\",\n",
    "        mo_index_title=\" \",\n",
    "        zoom=100,\n",
    "    )\n",
    "show_mos_grid(images=mo_images(\"oh_mo_*.png\"), mo_list=selected2.mo_list, columns=4, image_regex=r\"oh_mo_(\\d+)\")"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}