dev_guide.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction: Welcome to the PySCF tutorial\n",
    "\n",
    "This component of the tutorial is primarily concerned with computational chemistry calculations and will not cover the theoretical details of the methods utilized [i.e., Hartree–Fock (HF), Kohn–Sham Density Functional Theory (KS-DFT), Møller–Plesset Peturbation Theory (MP), Coupled Cluster Theory (CC), etc.].\n",
    "\n",
    "We will go through the steps necessary to carry out a variety of relevant calculations. These include:\n",
    "\n",
    "* Single Point Energies\n",
    "  * Ground State\n",
    "    1. atomization energies and bond dissociation energies\n",
    "    2. reaction energies\n",
    "    3. electron affinities, ionization potentials, and proton affinities (vertical and adiabatic)\n",
    "    4. relative energies of isomers\n",
    "    5. binding energies of non-covalent complexes (weak/hydrogen/rare-gas)\n",
    "    6. barrier heights\n",
    "    7. potential energy curves\n",
    "  * Excited State (Time-Dependent)\n",
    "    1. TD-HF/TD-DFT\n",
    "    2. wavefunction stability analysis\n",
    "    3. orbital visualization\n",
    "\n",
    "* Forces/Geometry Optimizations\n",
    "  * Ground State\n",
    "    1. force calculation of equilibrium structure\n",
    "    2. force calculation of non-equilibrium structure\n",
    "    3. transition state search\n",
    "    4. geometry optimization\n",
    "  * Excited State (Time-Dependent)\n",
    "    1. TODO\n",
    "\n",
    "* Frequencies\n",
    "  * Ground State\n",
    "    1. vibrational modes\n",
    "    2. free-energy calculation\n",
    "    3. is the result of a geometry optimization a local minimum?\n",
    "  * Excited State (Time-Dependent)\n",
    "    1. TODO\n",
    "\n",
    "Furthermore, we will demonstrate the effect of the most important computational settings on the final calculated results using a handful of the above interactions:\n",
    "  * Relevant Settings\n",
    "    1. basis set (minimal vs. double vs. triple. vs. quadruple vs. basis set limit) BSIE\n",
    "    2. DFT integration grid (coarse vs. fine, etc.)\n",
    "    3. method (HF/DFT/MP2/CC/etc.)\n",
    "    4. counterpoise corrections (which is an extension of basis set category) talk about BSIE/BSSE\n",
    "    5. integral threshes and convergence criteria\n",
    "    6. cart vs spherical?\n",
    "    7. restricted vs unrestricted\n",
    "    8. density fitting\n",
    "    9. opt algo, like diis gdm Newton etc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Getting started\n",
    "\n",
    "In order to get started, we need to import the pyscf module, which is accomplished by the command in the Code cell below. If this command fails, please see the [quick setup guide](http://sunqm.github.io/pyscf/tutorial.html#quick-setup) or the [detailed installation instructions](http://sunqm.github.io/pyscf/install.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyscf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to define the molecule, which requires importing the gto (gaussian type orbital) submodule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscf import gto"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Creating a molecule (mol) object with gto.Mole() <a name=\"molecule\"></a>\n",
    "\n",
    "We will now go through the steps involved in setting up a molecule, including a description of the eight most important attributes. The first molecule used in this demonstration is the water dimer (two hydrogen-bonded water molecules) from the widely-used [S22 dataset](http://pubs.rsc.org/en/Content/ArticleLanding/2006/CP/B600027D#!divAbstract) of non-covalent interactions, and can be found [here](http://www.begdb.com/index.php?action=oneMolecule&state=show&id=82). The molecule is pictured here:\n",
    "\n",
    "<img style=\"float: left;\" src=\"images/water_dimer.jpg\"><br clear=\"all\" />\n",
    "\n",
    "The first step is to create the molecule object, mol:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol=gto.Mole()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: The molecule object can be named nearly anything (e.g., mol1, mol2, mymol, etc.).\n",
    "\n",
    "## atom attribute <a name=\"atom\"></a>\n",
    "#### (NO DEFAULT)\n",
    "\n",
    "PySCF allows for a variety of [molecular input formats](http://sunqm.github.io/pyscf/gto.html#geometry), but the one that is most suitable for the current example is the use of triple-quotes. This allows one to simply copy and paste the geometry from the BEGDB website without further modification:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.atom=\"\"\"\n",
    "O  -1.551007  -0.114520   0.000000\n",
    "H  -1.934259   0.762503   0.000000\n",
    "H  -0.599677   0.040712   0.000000\n",
    "O   1.350625   0.111469   0.000000\n",
    "H   1.680398  -0.373741  -0.758561\n",
    "H   1.680398  -0.373741   0.758561\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The triple-quotes can be used to read four unique molecular input formats:\n",
    "\n",
    "1. XYZ (see geom/water_dimer.xyz) [TODO]\n",
    "2. Z-matrix (see geom/water_dimer.zmat) [TODO]\n",
    "3. Cartesian coordinates only (see geom/water_dimer.mol and example above)\n",
    "4. Cartesian coordinates or Z-matrix prepended by a line containing the charge and spin separated by a space (see geom/water_dimer.qc) [TODO]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to set up a molecule is to read the coordinates from a file on disk. The file containing the molecule can be formatted in any of the four ways listed above. As long as the file adheres to one of these formats, PySCF can automatically parse the file regardless of the extension. Thus, reading in a molecule file can be accomplished simply via [TODO]:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mol.atom=read_file('water_dimer.mol')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## basis attribute\n",
    "#### (NO DEFAULT)\n",
    "\n",
    "The simplest way to set up the basis is to specify the same basis set for all atoms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.basis='aug-cc-pVDZ'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the aug-cc-pVDZ basis set is used for all of the hydrogen and oxygen atoms.\n",
    "\n",
    "Alternatively, one can use a different basis set for different elements by using Python's dictionary data structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.basis={'H': 'cc-pVDZ', 'O': 'aug-cc-pVDZ'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the cc-pVDZ basis set is used for all of the hydrogen atoms, while the diffuse aug-cc-pVDZ basis set is used for all of the oxygen atoms.\n",
    "\n",
    "One can also use different basis sets for different atoms of the same element by labeling the atoms with integers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.atom=\"\"\"\n",
    "O1  -1.551007  -0.114520   0.000000\n",
    "H2  -1.934259   0.762503   0.000000\n",
    "H3  -0.599677   0.040712   0.000000\n",
    "O4   1.350625   0.111469   0.000000\n",
    "H5   1.680398  -0.373741  -0.758561\n",
    "H6   1.680398  -0.373741   0.758561\"\"\"\n",
    "mol.basis={'O1': 'cc-pVDZ', 'H2': 'cc-pVDZ', 'H3': 'aug-cc-pVDZ', 'O4': 'aug-cc-pVDZ', 'H5': 'cc-pVDZ', 'H6': 'cc-pVDZ'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the cc-pVDZ basis set is used for all of the atoms except for those involved in the hydrogen bond (H3 and O4). The diffuse aug-cc-pVDZ basis set is used for H3 and O4.\n",
    "\n",
    "Finally, one common way of reducing basis set superposition error (BSSE) is to use counterpoise corrections by employing ghost functions on ghost atoms. A useful guide to BSSE can be found [here](http://vergil.chemistry.gatech.edu/notes/cp.pdf). In PySCF, counterpoise corrections can be applied via simple modifications to the atom and basis attributes.\n",
    "\n",
    "For example, running the water dimer calculation with ghost functions on the second water monomer can be accomplished via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.atom=\"\"\"\n",
    "O  -1.551007  -0.114520   0.000000\n",
    "H  -1.934259   0.762503   0.000000\n",
    "H  -0.599677   0.040712   0.000000\n",
    "X-O   1.350625   0.111469   0.000000\n",
    "X-H   1.680398  -0.373741  -0.758561\n",
    "X-H   1.680398  -0.373741   0.758561\"\"\"\n",
    "mol.basis='aug-cc-pVDZ'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the parser understands that the latter three atoms are meant to be ghost atoms since the atomic symbol is prepended with the string, 'X-'.\n",
    "\n",
    "Advanced users who are interested in using custom basis sets or basis sets that are not pre-defined in PySCF are directed [here](http://sunqm.github.io/pyscf/gto.html#input-basis). Useful resources for the latter include the [EMSL basis set exchange](https://bse.pnl.gov/bse/portal), the [PSI4 basis set page](http://www.psicode.org/psi4manual/master/basissets_byfamily.html), Grant Hill's [ccRepo](http://www.grant-hill.group.shef.ac.uk/ccrepo/), and the [TURBOMOLE basis set library](http://www.cosmologic-services.de/basis-sets/basissets.php). PySCF can parse NWChem- and Gaussian94-formatted basis sets [TODO].\n",
    "\n",
    "For example, the 'doubly-augmented' d-aug-cc-pVDZ basis set is not pre-defined in PySCF. If one wants to use the d-aug-cc-pVDZ basis set for the oxygen atoms and the cc-pVDZ basis set for the hydrogen atoms, then the appropriate basis setup would be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.basis={'H': 'cc-pVDZ','O': gto.basis.parse('''\n",
    "O    S\n",
    "  11720.0000000              0.0007100             -0.0001600        \n",
    "   1759.0000000              0.0054700             -0.0012630        \n",
    "    400.8000000              0.0278370             -0.0062670        \n",
    "    113.7000000              0.1048000             -0.0257160        \n",
    "     37.0300000              0.2830620             -0.0709240        \n",
    "     13.2700000              0.4487190             -0.1654110        \n",
    "      5.0250000              0.2709520             -0.1169550        \n",
    "      1.0130000              0.0154580              0.5573680        \n",
    "O    S\n",
    "      0.3023000              1.0000000        \n",
    "O    S\n",
    "      0.0789600              1.0000000        \n",
    "O    S\n",
    "      0.0206000              1.0000000        \n",
    "O    P\n",
    "     17.7000000              0.0430180        \n",
    "      3.8540000              0.2289130        \n",
    "      1.0460000              0.5087280        \n",
    "O    P\n",
    "      0.2753000              1.0000000        \n",
    "O    P\n",
    "      0.0685600              1.0000000        \n",
    "O    P\n",
    "      0.0171000              1.0000000        \n",
    "O    D\n",
    "      1.1850000              1.0000000        \n",
    "O    D\n",
    "      0.3320000              1.0000000        \n",
    "O    D\n",
    "      0.0930000              1.0000000        \n",
    "''')}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additionally, PySCF can read a file containing the desired basis set [TODO]:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mol.basis={'H': gto.basis.read('basis/STO-2G.dat'),'O': 'STO-3G'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## cart attribute\n",
    "#### DEFAULT: `mol.cart=0`\n",
    "\n",
    "The cart attribute determines whether the d and/or higher basis functions are taken to be *spherical* (i.e., d=5 functions, f=7 functions, g=9 functions, etc.) or *Cartesian* (i.e., d=6 functions, f=10 functions, g=15 functions, etc.). `mol.cart=0` specifices *spherical* functions and `mol.cart=1` specifies *Cartesian* functions. If an existing basis set in PySCF is being used, the cart attribute is automatically set based on how the basis set itself was optimized. Thus, it is not necessary to define this attribute for most calculations.\n",
    "\n",
    "There are two instances where specifying the cart attribute is necessary:\n",
    "\n",
    "1.  If either a custom basis set is being used or one that is not pre-defined in PySCF, it is important to determine the appropriate value of cart and set it accordingly. For example, the [UGBS basis set](https://aip.scitation.org/doi/abs/10.1063/1.475959) is commonly used to determine absolute atomic energies. This basis set was optimized for use with spherical functions. Since UGBS is not pre-defined in PySCF, it is important to set `mol.cart=0`.\n",
    "\n",
    "2.  If the user is interested in experimenting with the effect of using spherical vs. Cartesian functions for a given basis set, it is possible to override the default for pre-defined basis sets. For example, the Dunning basis sets (cc-pVXZ and aug-cc-pVXZ) are meant to be used with spherical functions (`mol.cart=0`), yet one can set `mol.cart=1` to gauge the sensitivity of absolute and relative energies to this setting.\n",
    "\n",
    "## charge attribute\n",
    "#### DEFAULT: `mol.charge=0`\n",
    "\n",
    "The charge attribute sets the charge of the molecule. One should set `mol.charge=0` for neutral systems, `mol.charge=1` for cations, `mol.charge=2` for dications, `mol.charge=-1` for anions, `mol.charge=-2` for dianions, etc.\n",
    "\n",
    "## ecp attribute\n",
    "#### DEFAULT: `mol.ecp={}`\n",
    "\n",
    "When heavy elements are present in a molecule, standard Gaussian basis sets are oftentimes both insufficient (because they are unable to capture the relativistic nature of heavier elements) and impractical (because as one descends the rows of the periodic table, Gaussian basis sets tend to become very large and heavily-contracted due to the increased number of electrons). For this reason, a great deal of work has been devoted to developing effective core potentials (ECP) which replace core electrons around a nucleus by pseudopotentials.\n",
    "\n",
    "In PySCF, the ecp attribute is very similar to the basis attribute. To demonstrate the use of the ecp attribute, we will switch to a system with heavy elements, namely, TeH<sub>2</sub>–HI. This system is from the [HEAVY28](http://www.thch.uni-bonn.de/tc.old/downloads/GMTKN/GMTKN55/HEAVY28.html) dataset of the [GMTKN55](https://www.chemie.uni-bonn.de/pctc/mulliken-center/software/GMTKN/gmtkn55) database. Specifically, it is [data point](http://www.thch.uni-bonn.de/tc.old/downloads/GMTKN/GMTKN55/HEAVY28ref.html) \\#27, an intermolecular bond between tellerium hydride and hydrogen iodide:\n",
    "\n",
    "<img style=\"float: left;\" src=\"images/TeH2_HI.jpg\" width=\"50%\"><br clear=\"all\" />\n",
    "\n",
    "The atom and basis attributes can be specified as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.atom=\"\"\"\n",
    "Te         0.03866320      2.16215120      0.00000000\n",
    "H         -0.10453280      0.53339420      0.00000000\n",
    "H         -1.58803480      2.30795520      0.00000000\n",
    "H          1.61524120     -2.39760280      0.00000000\n",
    "I          0.03866320     -2.60589780      0.00000000\"\"\"\n",
    "mol.basis='def2-TZVPP'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to specify a single ECP for all atoms present in the system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.ecp='def2-TZVPP'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, if there is an atom in the system that does not have an ECP defined (hydrogen in this case), a warning will be shown once the object is built:\n",
    "\n",
    "`ECP def2-TZVPP not found for  H`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, it is advisable to use the dictionary data structure to specify the ECP for the heavy elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.ecp={'Te': 'def2-TZVPP', 'I': 'def2-TZVPP'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to use ECPs that are not defined in PySCF. Useful resources include the four websites listed earlier, as well as the [Stuttgart/Cologne ECP web page](http://www.tc.uni-koeln.de/PP/index.en.html). For example, if one wants to use the older def-TZVPP ECP for Te only (the newer version, def2-TZVPP, is pre-defined in PySCF), it is possible to get the neccessary information from the [TURBOMOLE basis set library](http://www.cosmologic-services.de/basis-sets/basissets.php):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.ecp={'Te': gto.basis.parse('''\n",
    "ECP\n",
    "Te nelec 46\n",
    "Te F\n",
    "2       1.93927000     -17.86464100\n",
    "Te S\n",
    "2       2.92379400      50.08380500\n",
    "2       1.15275400       1.96814000\n",
    "2       1.93927000      17.86464100\n",
    "Te P\n",
    "2       2.60308600     119.82070200\n",
    "2       0.98544800      -2.03904800\n",
    "2       1.93927000      17.86464100\n",
    "Te D\n",
    "2       1.43501900      37.75721400\n",
    "2       1.93927000      17.86464100\n",
    "'''),'I': 'def2-TZVPP'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PySCF can parse both NWChem-/Dalton- and Gaussian94-formatted ECP files. In addition, PySCF can read a file containing the desired ECP [TODO]:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mol.ecp={'Te': gto.basis.read('ecp/def-TZVPP.dat'),'I': 'def2-TZVPP'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "                        (ne, self.spin))\n",
    "   2316 return nalpha, nbeta\n",
    "\n",
    "RuntimeError: Electron number 87 and spin 0 are not consistent\n",
    "Note mol.spin = 2S = Nalpha - Nbeta, not 2S+1\n",
    "## spin attribute\n",
    "#### DEFAULT: `mol.spin=0`\n",
    "\n",
    "The spin attribute simply sets the number of unpaired electrons. By default, a closed-shell system is assumed: `mol.spin=0` For open-shell systems, the spin must be set. For example, the boron atom has a single unpaired electron: `mol.spin=1` The carbon atom has two unpaired electrons: `mol.spin=2` And finally, nitrogen has three unpaired electrons: `mol.spin=3` The same concept applies to open-shell molecules. For instance, the hydroxyl radical (HO•) should have `mol.spin=1`, since it has one unpaired electron.\n",
    "\n",
    "## unit attribute\n",
    "#### DEFAULT: `mol.charge='Angstrom'`\n",
    "\n",
    "The unit attribute can be either set to 'Angstrom' or 'Bohr'.\n",
    "\n",
    "## verbose attribute\n",
    "#### DEFAULT: `mol.verbose=0`\n",
    "\n",
    "The verbose attribute controls the print level for the molecule object. Setting `mol.verbose=0` will print little to no information, while setting `mol.verbose=4` prints useful information about the basis and number of basis functions [TODO]. Users who want to see detailed information should set `mol.verbose=10`.\n",
    "\n",
    "The final step is to build the molecule object (if any attributes are specified or modified after the molecule is built, this command should be executed again):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7f81c806c450>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol.build()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following Code cell combines the content covered above into a useable sample input for setting up a molecule. The more advanced attributes are listed (commented) at the bottom of the Code cell and the advanced user can find more information about them [here](http://sunqm.github.io/pyscf/gto.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "del mol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7f81c8021f50>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol=gto.Mole()\n",
    "\n",
    "mol.atom=\"\"\"\n",
    "O  -1.551007  -0.114520   0.000000\n",
    "H  -1.934259   0.762503   0.000000\n",
    "H  -0.599677   0.040712   0.000000\n",
    "O   1.350625   0.111469   0.000000\n",
    "H   1.680398  -0.373741  -0.758561\n",
    "H   1.680398  -0.373741   0.758561\"\"\"\n",
    "mol.basis='aug-cc-pVDZ'\n",
    "mol.cart=False\n",
    "mol.charge=0\n",
    "mol.ecp={}\n",
    "mol.spin=0\n",
    "mol.unit='Angstrom'\n",
    "mol.verbose=0\n",
    "\n",
    "mol.build()\n",
    "\n",
    "#ADVANCED ATTRIBUTES – SEE MANUAL\n",
    "#mol.groupname=\n",
    "#mol.incore_anyway=\n",
    "#mol.irrep_id=\n",
    "#mol.irrep_name=\n",
    "#mol.max_memory=\n",
    "#mol.nucmod=\n",
    "#mol.output=\n",
    "#mol.symm_orb=\n",
    "#mol.symmetry=\n",
    "#mol.symmetry_subgroup=\n",
    "#mol.topgroup="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C -2.4029304540002028 0.8689804209456048 -0.2449102395260002\n",
    "C -2.281163443796311 -0.1857634832845413 -1.1231569840665816\n",
    "C -1.0633014911752363 -0.8316364479270776 -1.1800813918528203\n",
    "C -0.014035524047731065 -0.4175648127113501 -0.3738617367804197\n",
    "C 1.2171179376059484 -1.0333512082930716 -0.3986919042759907\n",
    "C 2.267534690772756 -0.6294233235517346 0.39906088314243876\n",
    "C 2.116114451804474 0.43069218042086665 1.2695738382525337\n",
    "C 0.8841753597507057 1.0510997563391844 1.2981420832215316\n",
    "C -0.1637706632591607 0.643156143629797 0.4978932715232942\n",
    "N -1.3702126078038088 1.2404690985952442 0.5168228585183626\n",
    "H -3.3587893785117267 1.4044315562729146 -0.17510154691435847AllChem.MMFFOptimizeMolecule(mol)\n",
    "H -3.1202218930192935 -0.48743088253379613 -1.7418765147321889\n",
    "H -0.9157830536727101 -1.6780309196047807 -1.8639238501786064\n",
    "H 1.3218059425234456 -1.8742982433684525 -1.095731472195006\n",
    "H 3.213437359895837 -1.1631305784166692 0.3265725386987853\n",
    "H 2.9205623001017886 0.772365506775381 1.909000241093106\n",
    "H 0.7494604668030355 1.8894352366908078 1.980269926037916 H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\"---\n",
    "\n",
    "# Mean-field calculations\n",
    "\n",
    "After the molecule object has been built, it is necessary to create a mean-field object [Hartree–Fock (HF) or Kohn–Sham Density Functional Theory (KS-DFT)].\n",
    "\n",
    "In order to create a mean-field object, it is necessary to import the scf (self-consistent field) submodule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscf import scf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For most users, the relevant types of SCF calculations will be:\n",
    "\n",
    "1. RHF (Restricted Hartree–Fock)\n",
    "2. UHF (Unrestricted Hartree–Fock)\n",
    "3. ROHF (Restricted Open-Shell Hartree–Fock)\n",
    "4. RKS (Restricted Kohn–Sham)\n",
    "5. UKS (Unrestricted Kohn–Sham)\n",
    "6. ROKS (Restricted Open-Shell Kohn–Sham)\n",
    "\n",
    "Initializing the mean-field object for these six types of calculations is demonstrated below. These classes require at least one argument, namely, a molecule object, mol.\n",
    "\n",
    "1. `mf=scf.RHF(mol)`\n",
    "2. `mf=scf.UHF(mol)`\n",
    "3. `mf=scf.RHF(mol)` (where mol is defined to be open-shell, mol.spin!=0)\n",
    "4. `mf=scf.RKS(mol)`\n",
    "5. `mf=scf.UKS(mol)`\n",
    "6. `mf=scf.RKS(mol)` (where mol is defined to be open-shell, mol.spin!=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\"---\n",
    "\n",
    "# Creating a mean-field (mf) object with scf.RHF()\n",
    "\n",
    "Since RHF, ROHF, and UHF have nearly identical attributes, the tutorial will proceed with demonstrating an RHF calculation. The first step is to create the mean-field object (RHF in this case):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "mf=scf.RHF(mol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\" H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\" H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\" H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\" H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\" H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\"The input object is the molecule object (mol) that was discussed in detail [above](#molecule).\n",
    "\n",
    "## conv_tol attribute\n",
    "#### DEFAULT: `mf.conv_tol=1e-09`\n",
    "In order for an SCF calculation to converge, PySCF requires two criteria to be met. The first is controlled by the conv_tol attribute, namely, the difference in the SCF energy (in Hartrees) between two sucessive cycles.\n",
    "\n",
    "## conv_tol_grad attribute\n",
    "#### DEFAULT: `mf.conv_tol_grad=numpy.sqrt(mf.conv_tol)`\n",
    "The second criterion for convergence is the conv_tol_grad attribute, namely, the root-mean-square [TODO] of the orbital gradient. This is a vector of length nocc\\*nvirt.\n",
    "\n",
    "## direct_scf_tol attribute\n",
    "#### DEFAULT: `mf.direct_scf_tol=1e-13`\n",
    "The direct_scf_tol attribute is the threshold for discarding integrals.\n",
    "\n",
    "## init_guess attribute\n",
    "#### DEFAULT: `mf.init_guess='minao'`\n",
    "A good initial guess is vital for the efficient completion of the SCF procedure. There are four initial guess options in PySCF:\n",
    "1. `mf.init_guess='minao'`\n",
    "2. `mf.init_guess='atom'`\n",
    "3. `mf.init_guess='1e'`\n",
    "4. `mf.init_guess='chkfile'`\n",
    "\n",
    "The first option generates an initial guess for the density matrix based on the ANO basis, and then projects this onto the basis set specified in `mol.basis`.\n",
    "\n",
    "The second option generates an initial guess based on a superposition of atomic densities computed in the basis set specified in `mol.basis`.\n",
    "\n",
    "The third option sets the initial density matrix to zero. This is generally a very poor guess and should only be used for debugging purposes or for checking values between different packages.\n",
    "\n",
    "The last option is somewhat advanced and reads in an existing density matrix from disk. Further information on this option can be found [here](http://sunqm.github.io/pyscf/scf.html#pyscf.scf.hf.SCF.init_guess_by_chkfile).\n",
    "\n",
    "## max_cycle attribute\n",
    "#### DEFAULT: `mf.max_cycle=50`\n",
    "The max_cycle attribute simply sets the maximum number of SCF cycles that should be performed before the calculation terminates. For systems that are notoriously difficult to converge, the value should be increased to 100 or even 1000.\n",
    "\n",
    "## max_memory attribute\n",
    "#### DEFAULT: `mf.max_memory=4000`\n",
    "The max_memory attribute determines the maximum amount of memory (in Megabytes) that PySCF is allowed to utilize during the SCF procedure. This should be set by the user based on the memory limitations of the computer/server utilized.\n",
    "\n",
    "## verbose attribute\n",
    "#### DEFAULT: `mf.verbose=0`\n",
    "The verbose attribute controls the print level for the mean-field object. Setting `mol.verbose=0` will print only the SCF energy, while setting `mol.verbose=4` prints useful information about SCF settings as well as the SCF energy per iteration, HOMO/LUMO energies, and convergence metrics. Users who want to see detailed information should set `mol.verbose=10`, which provides additional information such as molecular orbital energies.\n",
    "\n",
    "The final step is to run the SCF calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-152.08859934750234"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C -2.4029304540002028 0.8689804209456048 -0.2449102395260002\n",
    "C -2.281163443796311 -0.1857634832845413 -1.1231569840665816\n",
    "C -1.0633014911752363 -0.8316364479270776 -1.1800813918528203\n",
    "C -0.014035524047731065 -0.4175648127113501 -0.3738617367804197\n",
    "C 1.2171179376059484 -1.0333512082930716 -0.3986919042759907\n",
    "C 2.267534690772756 -0.6294233235517346 0.39906088314243876\n",
    "C 2.116114451804474 0.43069218042086665 1.2695738382525337\n",
    "C 0.8841753597507057 1.0510997563391844 1.2981420832215316\n",
    "C -0.1637706632591607 0.643156143629797 0.4978932715232942\n",
    "N -1.3702126078038088 1.2404690985952442 0.5168228585183626\n",
    "H -3.3587893785117267 1.4044315562729146 -0.17510154691435847AllChem.MMFFOptimizeMolecule(mol)\n",
    "H -3.1202218930192935 -0.48743088253379613 -1.7418765147321889\n",
    "H -0.9157830536727101 -1.6780309196047807 -1.8639238501786064\n",
    "H 1.3218059425234456 -1.8742982433684525 -1.095731472195006\n",
    "H 3.213437359895837 -1.1631305784166692 0.3265725386987853\n",
    "H 2.9205623001017886 0.772365506775381 1.909000241093106\n",
    "H 0.7494604668030355 1.8894352366908078 1.980269926037916 H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\"After the SCF kernel has finished running, the mean-field object is updated with several useful output attributes:\n",
    "1. `mf.mo_coeff` - Molecular orbital (MO) coefficients (matrix where rows are atomic orbitals (AO) and columns are MOs)\n",
    "2. `mf.mo_energy` - MO energies (vector with length equal to number of MOs)\n",
    "3. `mf.mo_occ` - MO occupancy (vector with length equal to number of MOs)\n",
    "4. `mf.e_tot` - Total SCF energy in units of Hartrees\n",
    "5. `mf.converged` - Status of SCF convergence (True indicates converged and False indicates unconverged)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following Code cell combines the content covered above into a useable sample input for setting up a mean-field object. The more advanced attributes are listed (commented) at the bottom of the Code cell and the advanced user can find more information about them [here](http://sunqm.github.io/pyscf/scf.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-152.08859934750382"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf=scf.RHF(mol)\n",
    "\n",
    "mf.conv_tol=1e-12\n",
    "mf.conv_tol_grad=1e-8\n",
    "mf.direct_scf_tol=1e-13\n",
    "mf.init_guess='atom'\n",
    "mf.max_cycle=100\n",
    "mf.max_memory=8000\n",
    "mf.verbose=0\n",
    "\n",
    "mf.kernel()\n",
    "\n",
    "#ADVANCED ATTRIBUTES – SEE MANUAL\n",
    "#mf.chkfile=\n",
    "#mf.conv_check=\n",
    "#mf.damp=\n",
    "#mf.diis=\n",
    "#mf.diis_file=\n",
    "#mf.diis_space=\n",
    "#mf.diis_space_rollback=\n",
    "#mf.diis_start_cycle=\n",
    "#mf.direct_scf=\n",
    "#mf.level_shift="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to run an unrestricted Hartree–Fock calculation, the steps taken above can be followed, except the mean-field object should be initiated as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "mf=scf.UHF(mol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Creating a mean-field (mf) object with scf.RKS()\n",
    "\n",
    "The most widely-used SCF method is undoubtedly density functional theory (DFT), and running a DFT calculation in PySCF is nearly as straightforward as running a Hartree–Fock calculation.\n",
    "\n",
    "In order to gain access to additional functions that do not pertain to standard SCF/HF, it is necessary to import the dft (density functional theory) module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscf import dft"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to create the method object (RKS in this case):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "mf=scf.RKS(mol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\" H   0  0  0  0  0  0  0  0  0  0  0  0\"\"\"RHF and RKS (as well as UHF and UKS) share many attributes, so we will only cover the ones that pertain specifically to DFT. These additional attributes primarily concern the exchange-correlation functional and the grid used for its numerical integration.\n",
    "\n",
    "## xc attribute\n",
    "#### (NO DEFAULT)\n",
    "\n",
    "The xc attribute sets the density functional approximation. A vast number of exchange-correlation functionals are available for use in PySCF, essentially all of those that are available in the installed version of libxc. For example, the functionals available in libxc 4.0.4 are listed [here](https://gitlab.com/libxc/libxc/wikis/Functionals-list-4.0.4).\n",
    "\n",
    "While the parsing rules for the xc attribute are explained in the [PySCF manual](http://sunqm.github.io/pyscf/dft.html#customizing-xc-functional), we will go through a main relevant examples here:\n",
    "\n",
    "1. Standard (pre-coded) exchange-correlation functional\n",
    "\n",
    "  * A classic example of this scenario is the beloved B3LYP functional. To run a B3LYP calculation, simply specify `mf.xc='B3LYP'`.\n",
    "  \n",
    "  * As long as commas and operators are not used, PySCF will assume that the string indicates a full exchange-correlation functional. Other examples may include `mf.xc='PBE'`, `mf.xc='TPSS'`, `mf.xc='M06-L'`, etc.\n",
    "\n",
    "2. Combination of a single exchange functional with a single correlation functional\n",
    "\n",
    "  * This is a popular choice, particularly for using various combinations of the older exchange and correlation functionals from the 1980s and 1990s.\n",
    "  \n",
    "  * In PySCF, inserting a comma into the xc attribute string indicates a separation into exchange (left of the comma) and correlation (right of the comma). For instance, if once wants to run B88 exchange with PBE correlation, it is as simple as `mf.xc='B88,PBE'`.\n",
    "  \n",
    "  * Since the PBE exchange-correlation functional is a combination of PBE exchange and PBE correlation, it follows that `mf.xc='PBE'` and `mf.xc='PBE,PBE'` are exactly equivalant ways of using the PBE xc functional.\n",
    "  \n",
    "  * Other popular combinations of a single exchange functional with a single correlation functional include `mf.xc='B88,LYP'`, `mf.xc='revPBE,PBE'`, and `mf.xc='OPTX,LYP'`. These combinations are historically known as BLYP, revPBE, and OLYP, and can equivalently be called as `mf.xc='BLYP'`, `mf.xc='revPBE'`, and `mf.xc='OLYP'`.\n",
    "  \n",
    "  * This option is most useful when one is interested in an unconventional combination, for example, the combination of PBE exchange with TPSS correlation, `mf.xc='PBE,TPSS'`.\n",
    "\n",
    "3. Single exchange functional or single correlation functional\n",
    "\n",
    "  * Using a single exchange functional, for example, PBE exchange only, is as simple as `mf.xc='PBE,'`. Similarly, PBE correlation only can be specified with `mf.xc=',PBE'`. In the first case (PBE exchange only), the fact that there is a comma indicates to the parser that that a pre-coded xc functional is not being utilized. Therefore, anything to the left of the comma (PBE) is taken to be the exchange functional, while anything to the right of the comma (nothing) is taken to be the correlation functional. The reverse situation applies to the second case with PBE correlation only.\n",
    "\n",
    "4. Custom exchange-correlation functional\n",
    "\n",
    "  * Naturally, it is possible to fully customize the xc attribute. Once again, anything to the left of the comma is assumed to be referring to exchange, and anything to the right of the comma is assumed to be reffering to correlation.\n",
    "  \n",
    "  * In order to start with something simple, we will demonstrate how to specify the popular [PBE0](https://aip.scitation.org/doi/abs/10.1063/1.478522) xc functional by building it from its components (the easy way to do this is simply `mf.xc='PBE0'`). PBE0 is comprised of 25% exact (Hartree–Fock) exchange, 75% PBE exchange, and 100% PBE correlation. This can be accomplished by the following line: `mf.xc='0.25*HF+0.75*PBE,PBE'`.\n",
    "  \n",
    "  * A more complicated example is that of B3LYP. B3LYP is comprised of 20% exact (Hartree–Fock) exchange, 8% LDA exchange, 72% B88 exchange, 19% VWN correlation, and 81% LYP correlation. Unfortunately, there are 6 existing parameterization of the VWN LDA correlation functional, and this has led to much confusion through the years.\n",
    "  \n",
    "  * PySCF, TURBOMOLE, GAMESS, ORCA, and PSI4 all use the VWN5 parameterization, and the following two assignments are equivalent:\n",
    "    * `mf.xc='B3LYP'`\n",
    "    * `mf.xc='0.2*HF+0.08*LDA+0.72*B88,0.19*VWN5+0.81*LYP'`\n",
    "  * Q-Chem and NWChem use the VWN1RPA parameterization:\n",
    "    * `mf.xc='0.2*HF+0.08*LDA+0.72*B88,0.19*VWN_RPA+0.81*LYP'`\n",
    "  * Gaussian uses the VWN3 parameterization:\n",
    "    * `mf.xc='0.2*HF+0.08*LDA+0.72*B88,0.19*VWN3+0.81*LYP'`\n",
    "\n",
    "Running B3LYP is as simple as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "mf.xc='B3LYP'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## nlc attribute\n",
    "#### DEFAULT: `mf.nlc=''`\n",
    "\n",
    "The nlc attribute is used to specify the inclusion of the VV10 nonlocal correlation functional. This can simply be accomplished by `mf.nlc='VV10'`. Certain density functionals such as [ωB97X-V](http://pubs.rsc.org/en/Content/ArticleLanding/2014/CP/c3cp54374a#!divAbstract), [ωB97M-V](https://aip.scitation.org/doi/abs/10.1063/1.4952647), and [B97M-V](https://aip.scitation.org/doi/abs/10.1063/1.4907719) require this attribute to run correctly. More information on the VV10 nonlocal correlation functional can be found [here](https://aip.scitation.org/doi/abs/10.1063/1.3521275).\n",
    "\n",
    "For example, the following correctly specifies the ωB97M-V density functional:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "mf.xc='wB97M-V'\n",
    "mf.nlc='VV10'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The internal VV10 parameters (b and C) are automatically set by libxc.\n",
    "\n",
    "## grids attribute/class\n",
    "\n",
    "Although PySCF allows for the comprehensive customization of the DFT integration grid, non-expert users are advised to rely on the default settings as they have been fine-tuned to maximize accuracy and minimize computational cost.\n",
    "\n",
    "### grids.level attribute\n",
    "#### DEFAULT: `mf.grids.level=3`\n",
    "\n",
    "This is the simplest grid attribute to manipulate in order to increase or decrease the number of grid points. While the default value is 3, the range of possible values starts at `mf.grids.level=0` and continues to `mf.grids.level=3`.\n",
    "\n",
    "### grids.atom_grid attribute\n",
    "#### (NO DEFAULT)\n",
    "\n",
    "This attribute can be used to set the total number of radial shells and the number of angular grid points per shell. The format is (radial,angular). It is possible to choose the same setting for all atoms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "mf.grids.atom_grid=(99,590)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or to use Python's dictionary data structure to specify a different setting for different atoms/elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "mf.grids.atom_grid={'H': (50,194),'O': (75,302)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of radial shells can be set to any number larger than 0, while the number of angular grid points per shell must adhere to a valid Lebedev quadrature:\n",
    "\n",
    "{6,14,26,38,50,74,86,110,146,170,194,230,266,302,350,434,590,770,974,1202,1454,1730,2030,2354,2702,3074,3470,3890,4334,4802,5294,5810}\n",
    "\n",
    "Commonly used combinations include:\n",
    "* very coarse: `mf.grids.atom_grid=(23,170)`\n",
    "* coarse: `mf.grids.atom_grid=(50,194)`\n",
    "* moderate: `mf.grids.atom_grid=(75,302)`\n",
    "* fine: `mf.grids.atom_grid=(99,590)`\n",
    "* extremely fine: `mf.grids.atom_grid=(500,974)`\n",
    "\n",
    "Setting this attribute will override the `mf.grids.level` attribute."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### grids.atomic_radii attribute\n",
    "#### DEFAULT: `mf.grids.atomic_radii=dft.radi.BRAGG_RADII` [TODO]\n",
    "\n",
    "This attribute chooses the atomic radius values that will be used to adjust the spacing of the radial shells. If `grids.radii_adjust=None`, setting this parameter will have no affect. There are two built-in options in PySCF:\n",
    "\n",
    "1. `mf.grids.atomic_radii=dft.radi.BRAGG_RADII`\n",
    "2. `mf.grids.atomic_radii=dft.radi.COVALENT_RADII`\n",
    "\n",
    "The first one corresponds to the values found in this 1964 [paper](https://aip.scitation.org/doi/10.1063/1.1725697) by J. C. Slater, while the second one is primarily constructed from the values found in this [paper](http://pubs.rsc.org/en/content/articlelanding/2008/dt/b801115j#!divAbstract).\n",
    "\n",
    "### grids.becke_scheme attribute\n",
    "#### DEFAULT: `mf.grids.becke_scheme=dft.gen_grid.original_becke` [TODO]\n",
    "\n",
    "This attribute selects the partition function used to determine the grid point weights. PySCF allows for two options:\n",
    "\n",
    "1. `mf.grids.becke_scheme=dft.gen_grid.original_becke`\n",
    "2. `mf.grids.becke_scheme=dft.gen_grid.stratmann`\n",
    "\n",
    "The first scheme is from the famous 1988 [paper](https://aip.scitation.org/doi/abs/10.1063/1.454033) by Axel Becke, while the second scheme devised by Stratmann, Scuseria, and Frisch is found in this [paper](https://www.sciencedirect.com/science/article/pii/0009261496006008).\n",
    "\n",
    "### grids.prune attribute\n",
    "#### DEFAULT: `mf.grids.prune=dft.gen_grid.nwchem_prune` [TODO]\n",
    "\n",
    "Pruning integration grids offers the benefit of reducing the number of grid points (as well as the computation time) without significantly affecting the total energy. PySCF has several build in pruning options:\n",
    "\n",
    "1. `mf.grids.prune=dft.gen_grid.nwchem_prune`\n",
    "2. `mf.grids.prune=dft.gen_grid.sg1_prune`\n",
    "3. `mf.grids.prune=dft.gen_grid.treutler_prune`\n",
    "4. `mf.grids.prune=None`\n",
    "\n",
    "### grids.radi_method attribute\n",
    "#### DEFAULT: `mf.grids.radi_method=dft.radi.treutler_ahlrichs` [TODO]\n",
    "\n",
    "Setting up the radial shells.\n",
    "\n",
    "1. `mf.grids.radi_method=dft.radi.treutler_ahlrichs`\n",
    "2. `mf.grids.radi_method=dft.radi.delley`\n",
    "3. `mf.grids.radi_method=dft.radi.mura_knowles`\n",
    "4. `mf.grids.radi_method=dft.radi.gauss_chebyshev`\n",
    "\n",
    "### grids.radii_adjust attribute\n",
    "#### DEFAULT: `mf.grids.radii_adjust=dft.radi.treutler_atomic_radii_adjust` [TODO]\n",
    "\n",
    "There are two options for adjusting the atomic radii.\n",
    "\n",
    "1. `mf.grids.radii_adjust=dft.radi.treutler_atomic_radii_adjust`\n",
    "2. `mf.grids.radii_adjust=dft.radi.becke_atomic_radii_adjust`\n",
    "3. `mf.grids.radii_adjust=None`\n",
    "\n",
    "### grids.verbose attribute\n",
    "#### DEFAULT: `mf.grids.verbose=0`\n",
    "\n",
    "The verbose attribute controls the print level for the integration grid. Setting `mf.grids.verbose=0` will print no information about the grid, while setting `mf.grids.verbose=4` prints useful information such as the total number of grid points, as well as information regarding the schemes used to set up the grid.\n",
    "\n",
    "## nlcgrids attribute/class\n",
    "\n",
    "The nlcgrids attribute functions identically to the grids attribute. The only difference is that nlcgrids controls the grid settings for the nonlocal correlation functional, whereas grids controls the grid for the local exchange-correlation functional. The VV10 nonlocal correlation functional requires a much coarser grid than needed for the local xc functional, which is a good thing b/c the scaling of nlc is quadratic in grid points (loops through r and r'). Example of specifying a coarse grid here.\n",
    "\n",
    "## small_rho_cutoff attribute\n",
    "#### DEFAULT: `mf.small_rho_cutoff=1e-07`\n",
    "\n",
    "This attribute is used to discard grid points that contribute negligibly to the total electron count as computed by the chosen integration grid. Given a vector, **r**, of length ng (number of grid points), that contains the value of the electron density at each grid point, and a vector, **w**, of length ng that contains the weights associated with each grid point, grid points are discarded if **r**∘**w** ≤ `mf.small_rho_cutoff`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "radial grids: Gauss-Chebyshev [JCP 108, 3226 (1998); DOI:10.1063/1.475719) radial grids\n",
      "becke partition: Becke, JCP 88, 2547 (1988); DOI:10.1063/1.454033\n",
      "pruning grids: None\n",
      "grids dens level: 3\n",
      "symmetrized grids: False\n",
      "User specified grid scheme (50, 194)\n",
      "tot grids = 58200\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-152.82239954988563"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf=scf.RKS(mol)\n",
    "\n",
    "mf.conv_tol=1e-12\n",
    "mf.conv_tol_grad=1e-8\n",
    "mf.direct_scf_tol=1e-13\n",
    "mf.init_guess='atom'\n",
    "mf.max_cycle=100\n",
    "mf.max_memory=8000\n",
    "mf.verbose=0\n",
    "\n",
    "mf.xc='B3LYP'\n",
    "mf.grids.atom_grid=(50,194)\n",
    "mf.grids.becke_scheme=dft.gen_grid.original_becke\n",
    "mf.grids.prune=None\n",
    "mf.grids.radi_method=dft.radi.gauss_chebyshev\n",
    "mf.grids.radii_adjust=None\n",
    "mf.grids.verbose=4\n",
    "mf.small_rho_cutoff=1e-10\n",
    "\n",
    "mf.kernel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:qiskit-env]",
   "language": "python",
   "name": "conda-env-qiskit-env-py"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}