Argon_SQM0.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2019 MolSSI Summer School QM project: semiempirical model of Argon\n",
    "\n",
    "## 1. Introduction\n",
    "\n",
    "In this project, we will simulate a cluster of Argon atoms using quantum mechanics (QM). First-principles (a.k.a. ab initio) QM simulations are complicated and expensive, and a quick implementation would rely on a substantial amount of pre-existing software infrastructure (e.g. PySCF or Psi4). Instead, we will implement a much simpler semiempirical QM simulation that has been designed and parameterized to reproduce first-principles QM data using a minimal model. We can then limit our external dependencies to the standard numerical functionality of Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import NumPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As is typically the case in quantum chemistry, we will input a set of atomic coordinates $\\vec{r}_i$ and compute the ground state energy of Argon atoms with those coordinates. All physical quantities in this project will be in Hartree atomic units, where the bohr is the unit of length and the hartree is the unit of energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write down a set of atomic coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More complicated and featureful software would be able to compute other properties besides the ground state energy like atomic forces. It would also have more convenient options for specifying inputs and returning outputs.\n",
    "\n",
    "Note that this QM project is primarily meant as a programming exercise. It contains a theoretical specification of a model, and we will implement the components of the model in software as they are specified. If you have a strong background in quantum chemistry, then you should find the theory reasonably familiar. However, it is only really necessary for you to understand each mathematical expression one-at-a-time to effectively implement this software.\n",
    "\n",
    "## 2. Model Hamiltonian\n",
    "\n",
    "As is standard in quantum chemistry, we will assume that the total energy of our system is defined to be the ground state energy of a quantum many-body Hamiltonian $\\hat{H}$. In second quantization notation, we can write it as\n",
    "\n",
    "$$ \\hat{H} = E_{\\mathrm{ion}} + \\sum_{p,q} \\sum_{\\sigma \\in \\{ \\uparrow , \\downarrow \\} } h_{p,q} \\hat{a}_{p,\\sigma}^{\\dagger} \\hat{a}_{q,\\sigma} + \\tfrac{1}{2}\\sum_{p,q,r,s} \\sum_{\\sigma,\\sigma' \\in \\{ \\uparrow , \\downarrow \\} } V_{p,q,r,s} \\hat{a}_{p,\\sigma}^{\\dagger} \\hat{a}_{r,\\sigma'}^{\\dagger} \\hat{a}_{s,\\sigma'} \\hat{a}_{q,\\sigma} , $$\n",
    "\n",
    "where $\\hat{a}_{p,\\sigma}^{\\dagger}$ and $\\hat{a}_{p,\\sigma}$ are the electron raising and lowering operators for an atomic orbital index $p$ and spin $\\sigma$. We will not be using $\\hat{H}$ itself in our calculations, but we will make use of the coefficient tensors $h_{p,q}$ and $V_{p,q,r,s}$. In first-principles calculations, each element of $h_{p,q}$ and $V_{p,q,r,s}$ would require the evaluation of a complicated integral. In our semiempirical model, we will set most of them to zero and assign a simple analytical form to the rest of them. The notation being used here is mostly consistent with modern quantum chemistry notation, but some objects, particularly $V_{p,q,r,s}$, have multiple conventions in practice.\n",
    "\n",
    "### A. Model design & parameters\n",
    "\n",
    "This semiempirical model combines some standard concepts and methods used in physics and chemistry. First, it will use a minimal number of electronic degrees of freedom. Because Argon is a noble gas, it interacts primarily through London dispersion forces that are mediated by quantum dipole fluctuations. The lowest energy dipole transition is from the occupied $3p$ states to the unoccupied $4s$ state, and we will include these 4 atomic orbitals per atom. Similarly, we will use a multipole expansion to simplify electronic excitations and retain only the monopole and dipole terms, which also restricts electronic polarization to 4 degrees of freedom per atom. We will use $\\{s, p_x, p_y, p_z\\}$ to label both atomic orbitals and multipole moments on each atom. The nuclear charge of Argon is 18, but our model combines the nucleus and the 12 neglected electrons ($1s^2$, $2s^2$, $2p^6$, and $3s^2$) into an ionic point charge with $Z = 6$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# specify basic properties of the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The index of an atomic orbital specifies which atom it is located on and what type it is. We will often extract these individual pieces of information using $\\mathrm{atom}(p)$ to denote the atom's index and $\\mathrm{orb}(p)$ to denote the orbital type. This is the first of many instances in this project where we could either represent something as a pre-tabulate list or a function. We will always make the simpler choice, in this case functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform basic book-keeping operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will discuss the model parameters in more detail as they are used, but it is a good idea to first collect them all in a common data structure, a Python dictionary, for convenient access throughout the notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# REMINDER: atomic units w/ energies in hartree, distances in bohr\n",
    "model_parameters = { 'r_hop' : 3.1810226927827516, # hopping length scale\n",
    "                     't_ss' : 0.03365982238611262, # s-s hopping energy scale\n",
    "                     't_sp' : -0.029154833035109226, # s-p hopping energy scale\n",
    "                     't_pp1' : -0.0804163845390335, # 1st p-p hopping energy scale\n",
    "                     't_pp2' : -0.01393611496959445, # 2nd p-p hopping energy scale\n",
    "                     'r_pseudo' : 2.60342991362958, # pseudopotential length scale\n",
    "                     'v_pseudo' : 0.022972992186364977, # pseudopotential energy scale\n",
    "                     'dipole' : 2.781629275106456, # dipole strength of s-p transition\n",
    "                     'energy_s' : 3.1659446174413004, # onsite energy of s orbital\n",
    "                     'energy_p' : -2.3926873325346554, # onsite energy of p orbital\n",
    "                     'coulomb_s' : 0.3603533286088998, # Coulomb self-energy of monopole\n",
    "                     'coulomb_p' : -0.003267991835806299 } # Coulomb self-energy of dipole"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are no parameters related to orbital overlap because all atomic orbitals are assumed to be orthogonal. The parameter values have been pre-optimized for this project, but the fitting process and reference data are both listed at the end of the project if you'd like to learn more about them.\n",
    "\n",
    "### B. Slater-Koster tight-binding model\n",
    "\n",
    "We will describe the kinetic energy of electrons using a simplified [Slater-Koster tight-binding method](https://en.wikipedia.org/wiki/Tight_binding). Because of the symmetry of atomic orbitals and the translational invariance of the kinetic energy operator, there are 4 distinct, distance-dependent \"hopping\" energies that characterize the interatomic kinetic energy between s and p orbitals:\n",
    "\n",
    "![s-p hopping diagram](hopping_cases.png)\n",
    "\n",
    "All other atomic orientations can be related to these cases by a change of coordinates. While it is compatible with very general functional forms, we will use a Gaussian form to simplify the model and its implementation. The distance-dependence of this simple version is controlled by a single hopping length scale $r_{\\mathrm{hop}}$ and the strength of each type of hopping energy,\n",
    "\n",
    "$$ t_{o,o'}(\\vec{r}) = \\exp(1-r^2/r_{\\mathrm{hop}}^2) \\times \\begin{cases}\n",
    " t_{ss} , & o = o' = s \\\\\n",
    " [\\vec{o}' \\cdot (\\vec{r}/r_{\\mathrm{hop}})] t_{sp}, & o = s \\ \\& \\ o' \\in \\{p_x, p_y, p_z\\} \\\\\n",
    " -[\\vec{o} \\cdot (\\vec{r}/r_{\\mathrm{hop}})] t_{sp} , & o' = s \\ \\& \\ o \\in \\{p_x, p_y, p_z\\} \\\\\n",
    " (r^2/r_{\\mathrm{SK}}^2)\\,(\\vec{o} \\cdot \\vec{o}')  t_{pp2} - [\\vec{o} \\cdot (\\vec{r}/r_{\\mathrm{SK}})] [\\vec{o}' \\cdot (\\vec{r}/r_{\\mathrm{SK}})] (t_{pp1} + t_{pp2}), & o,o' \\in \\{p_x, p_y, p_z\\}\n",
    " \\end{cases} $$\n",
    " \n",
    "where $o$ and $o'$ are the orbital types of the 1st and 2nd atoms and $\\vec{r}$ is a vector pointing from the 2nd atom to the 1st atom. We are assigning direction vectors to the p orbitals, $\\vec{p}_x \\equiv (1,0,0)$, $\\vec{p}_y \\equiv (0,1,0)$, and $\\vec{p}_z \\equiv (0,0,1)$, to simplify the notation. This project has multiple case-based formulas, and we will implement them using a code structure similar to each formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# implement the hopping energy function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### C. Coulomb interaction\n",
    "\n",
    "For the purpose of electrostatics, we will describe all inter-atomic Coulomb interactions with point charges. We need access to both the interaction kernel and its derivatives to define the multipole expansion,\n",
    "\n",
    "$$ V_{o,o'}(\\vec{r}) = \\begin{cases}\n",
    " 1/r , & o = o' = s \\\\\n",
    " (\\vec{o}' \\cdot \\vec{r}) / r^3, & o = s \\ \\& \\ o' \\in \\{p_x, p_y, p_z\\} \\\\\n",
    " -(\\vec{o} \\cdot \\vec{r}) / r^3 , & o' = s \\ \\& \\ o \\in \\{p_x, p_y, p_z\\} \\\\\n",
    " (\\vec{o} \\cdot \\vec{o}') / r^3 - 3 (\\vec{o} \\cdot \\vec{r}) (\\vec{o}' \\cdot \\vec{r}) / r^5, & o,o' \\in \\{p_x, p_y, p_z\\}\n",
    " \\end{cases} $$\n",
    "\n",
    "The Coulomb energy has the same case-based structure and angular dependences as the hopping energy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# implement the Coulomb energy function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The semiempirical model approximations strongly distort the physics of inter-atomic Pauli repulsion, and we compensate for these errors with a short-range ionic pseudopotential, which is a common tool in physics for building effective models of ionic cores:\n",
    "\n",
    "$$ V_o^{\\mathrm{pseudo}}(\\vec{r}) = V_{\\mathrm{pseudo}} \\exp(1 - r^2/r_{\\mathrm{pseudo}}^2) \\times \\begin{cases}\n",
    "1, & o = s \\\\\n",
    "-2 \\vec{o}\\cdot\\vec{r}/r_{\\mathrm{pseudo}}, & o \\in \\{p_x, p_y, p_z\\} \\end{cases} . $$\n",
    "\n",
    "Note that the vector here points from the atom that the pseudopotential is centered on to the orbital that is interacting with it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# implement the pseudopotential energy function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These interaction kernels enable us to define and calculate the ion-ion energy in $\\hat{H}$,\n",
    "\n",
    "$$ E_{\\mathrm{ion}} = Z^2 \\sum_{i < j} V_{s,s}(\\vec{r}_i - \\vec{r}_j) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# implement the ionic energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the vector of electron-ion interactions (omitting on-site Coulomb terms),\n",
    "\n",
    "$$ V^{\\mathrm{ion}}_p = \\sum_{i \\neq \\mathrm{atom}(p)}\n",
    "   V_{\\mathrm{orb}(p)}^{\\mathrm{pseudo}}(\\vec{r}_{\\mathrm{atom}(p)} - \\vec{r}_i)\n",
    "   - Z V_{\\mathrm{orb}(p),s}(\\vec{r}_{\\mathrm{atom}(p)} - \\vec{r}_i)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# implement the potential vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the matrix of electron-electron interaction matrix elements,\n",
    "\n",
    "$$ V^{\\mathrm{ee}}_{p,q} = \\begin{cases}\n",
    " V_{\\mathrm{orb}(p),\\mathrm{orb}(q)}(\\vec{r}_{\\mathrm{atom}(p)} - \\vec{r}_{\\mathrm{atom}(q)}) , & \\mathrm{atom}(p) \\neq \\mathrm{atom}(q) \\\\\n",
    " V^{\\mathrm{self}}_{\\mathrm{orb}(p)} \\delta_{p,q} , & \\mathrm{atom}(p) = \\mathrm{atom}(q)\n",
    "\\end{cases} . $$\n",
    "\n",
    "On-site Coulomb interactions between electrons are described by electronic self-energy parameters $V_o^{\\mathrm{self}}$ for multipole moment $o$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# implement the interaction matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### D. Multipole decomposition\n",
    "\n",
    "To define $V_{p,q,r,s}$ based on the Coulomb matrix elements $V_{p,q}^{\\mathrm{ee}}$, we need to define a mapping from products of atomic orbitals to a linear combination of terms in the multipole expansion of electronic charge. Because the atomic orbitals are normalized, their monopole coefficient with themselves is 1 (i.e. they have unit charge). Because of orbital orthogonality, there is no monopole term for either intra-atomic or inter-atomic transitions between atomic orbitals. We will ignore dipole transitions between atoms, which corresponds to the [neglect of diatomic differential overlap](https://en.wikipedia.org/wiki/NDDO) (NDDO) approximation that is commonly used in semiempirical quantum chemistry. All that remains is the intra-atomic s-p transition, and we will use a model parameter, $D$, to define its dipole strength. These transformation rules between atomic orbitals and multipole moments can be summarized pictorially,\n",
    "\n",
    "![s-p multipole diagram](multipole_cases.png)\n",
    "\n",
    "or mathematically as a 3-index tensor $\\chi_{p,q,r}$ where $p$ and $q$ are the atomic orbital indices and $r$ is the multipole moment index,\n",
    "\n",
    "$$ \\chi_{p, q, r} = \\begin{cases} 1, & \\mathrm{orb}(p) = \\mathrm{orb}(q) \\ \\& \\ \\mathrm{orb}(r) = s \\ \\& \\ \\mathrm{atom}(p) = \\mathrm{atom}(q) = \\mathrm{atom}(r) \\\\ \n",
    " D , & \\mathrm{orb}(q) = \\mathrm{orb}(r) \\in \\{p_x, p_y, p_z\\} \\ \\& \\ \\mathrm{orb}(p) = s \\ \\& \\ \\mathrm{atom}(p) = \\mathrm{atom}(q) = \\mathrm{atom}(r) \\\\\n",
    " D , & \\mathrm{orb}(p) = \\mathrm{orb}(r) \\in \\{p_x, p_y, p_z\\} \\ \\& \\ \\mathrm{orb}(q) = s \\ \\& \\ \\mathrm{atom}(p) = \\mathrm{atom}(q) = \\mathrm{atom}(r) \\\\\n",
    " 0, & \\mathrm{otherwise} \\end{cases} . $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# implement the chi tensor, give direct access to its on-atom value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The multipole expansion plays the same role as an auxiliary basis set in modern resolution-of-identity (RI) methods that are used to accelerate large quantum chemistry simulations. The primary purpose of these methods is to decompose $V_{p,q,r,s}$ into a low-rank factored form, which is\n",
    "\n",
    "$$ V_{p,q,r,s} = \\sum_{t,u} \\chi_{p,q,t} V_{t,u}^{\\mathrm{ee}} \\chi_{r,s,u} $$\n",
    "\n",
    "in our semiempirical model. Unlike the previous vectors and matrices that we have constructed, the $\\chi$ and $V$ tensors are sparse, meaning that most of their entries are zero. Without a more clever implementation, it can be computationally slow and wasteful to store and compute with these zero values. Modern quantum chemistry methods contain many ways of identifying and utilizing sparsity. For example, we could compute tensor elements of $\\chi$ on-the-fly rather than storing the full tensor. This is analogous to \"integral direct\" methods in quantum chemistry. For simplicity, we will still precompute and store $\\chi$, which fails to utilize its sparsity. However, we will avoid the explicit construction of $V_{p,q,r,s}$ by utilizing its factored form.\n",
    "\n",
    "### E. 1-body Hamiltonian\n",
    "\n",
    "The 1-body Hamiltonian coefficients $h_{p,q}$ combine many of the components that we have already discussed and implemented along with the on-site orbital energies, $E_s$ and $E_p$, which are the final two parameters from the semiempirical model,\n",
    "\n",
    "$$ h_{p,q} = \\begin{cases}\n",
    " t_{\\mathrm{orb}(p),\\mathrm{orb}(q)}(\\vec{r}_{\\mathrm{atom}(p)} - \\vec{r}_{\\mathrm{atom}(q)})\n",
    "  , & \\mathrm{atom}(p) \\neq \\mathrm{atom}(q) \\\\\n",
    " E_{\\mathrm{orb}(p)} \\delta_{\\mathrm{orb}(p),\\mathrm{orb}(q)} + \\sum_{r} \\chi_{p,q,r} V_{r}^{\\mathrm{ion}} , & \\mathrm{atom}(p) = \\mathrm{atom}(q) \n",
    "  \\end{cases}. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# assemble the Hamiltonian matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now fully specified the many-body Hamiltonian $\\hat{H}$, and we can now move on to approximating and calculating some of its physical properties.\n",
    "\n",
    "## 3. Hartree-Fock theory\n",
    "\n",
    "Even with a simple model for $\\hat{H}$, we cannot calculate its ground state energy exactly for more than a few Argon atoms. Instead, we will use the [Hartree-Fock approximation](https://en.wikipedia.org/wiki/Hartree–Fock_method), which restricts the ground-state wavefunction to a single Slater determinant. We will find the Slater determinant with the lowest total energy, but a more general wavefunction will usually have an even lower energy. In the next section, we will use many-body perturbation theory to improve our estimate of the total energy.\n",
    "\n",
    "The central objects of Hartree-Fock theory are the 1-electron density matrix $\\rho_{p,q}$ and the Fock matrix $f_{p,q}$. These two matrices depend on each other, which defines a nonlinear set of equations that we must solve iteratively. Iteratively solving these equations is usually referred to as the self-consistent field (SCF) cycle. For this to converge, we must start from a reasonable initial guess for $\\rho_{p,q}$. We will initialize it to the density matrix for isolated Argon atoms,\n",
    "\n",
    "$$ \\rho_{p,q}^{\\mathrm{atom}} = \\begin{cases} \n",
    "   1, & p = q \\ \\& \\ \\mathrm{orb}(p) \\in \\{ p_x, p_y, p_z \\} \\\\\n",
    "   0, & \\mathrm{otherwise}\n",
    "\\end{cases} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct an initial atomic guess for the SCF cycle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of spin symmetry, we use the same $\\rho_{p,q}$ for each electron spin type, which reduces the sum over spin types in $\\hat{H}$ to degeneracy pre-factors. Thus, half of the electrons are spin-up and half are spin-down.\n",
    "\n",
    "The Fock matrix is defined by the density matrix,\n",
    "\n",
    "$$ \\begin{align} f_{p,q} & = h_{p,q} + \\sum_{r,s} ( 2 V_{p,q,r,s} - V_{r,q,p,s} )\\rho_{r,s} \\\\\n",
    "                         & = h_{p,q} + \\sum_{r,s,t,u} ( 2 \\chi_{p,q,t} \\chi_{r,s,u} - \\chi_{r,q,t} \\chi_{p,s,u} )\n",
    "                             V_{t,u}^{\\mathrm{ee}} \\rho_{r,s} ,\n",
    "   \\end{align} $$\n",
    "\n",
    "which is just a sum of two tensor contractions between the $V$ tensor and the density matrix. Because the $V$ tensor has low rank, we can save a lot of computational effort by using its factored form to construct the Fock matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct the Fock matrix using einsum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is our first example of using \"einsum\" in NumPy to trivialize the implementation of a tensor equation using the [Einstein summation convention](https://en.wikipedia.org/wiki/Einstein_notation) and some extra conventions for the order of indices (alphabetical order).\n",
    "\n",
    "The Fock matrix defines the matrix of molecular orbitals $\\phi_{i,j}$ as its eigenvectors,\n",
    "\n",
    "$$ \\sum_{q} f_{p,q} \\phi_{q,r} = E_r \\phi_{p,r}, $$\n",
    "\n",
    "and the density matrix is defined by the occupied orbitals (we occupy the lowest-energy orbitals with all of the electrons available),\n",
    "\n",
    "$$ \\rho_{p,q} = \\sum_{i \\in \\mathrm{occ}} \\phi_{p,i} \\phi_{q,i}. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct the density matrix using NumPy linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate a $\\rho_{p,q}$ and $f_{p,q}$ that are consistent with each other, we mix in a fraction of the new density matrix with the old density matrix and re-calculate them until convergence. A very simple pseudocode for one loop of this cycle is\n",
    "\n",
    "$$ \\begin{align} f_{\\mathrm{new}} &= f(\\rho_{\\mathrm{old}}) \\\\\n",
    "                 \\rho_{\\mathrm{new}} &= w_{\\mathrm{mix}} \\rho(f_{\\mathrm{new}})\n",
    "                 + (1 - w_{\\mathrm{mix}}) \\rho_{\\mathrm{old}} \\end{align} $$\n",
    "\n",
    "This is a crude strategy for an SCF cycle, but it works for weakly interacting atoms like Argon:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# implement a self-contained SCF cycle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have a converged density matrix, the Hartree-Fock total energy is defined as\n",
    "\n",
    "$$ \\begin{align}\n",
    "  E_{\\mathrm{HF}} &= E_{\\mathrm{ion}} + E_{\\mathrm{SCF}} \\\\\n",
    "  E_{\\mathrm{SCF}} &= 2 \\sum_{p,q} h_{p,q} \\rho_{p,q}\n",
    "  + \\sum_{p,q,r,s} V_{p,q,r,s} ( 2 \\rho_{p,q} \\rho_{r,s} - \\rho_{p,s} \\rho_{r,q} ) = \\sum_{p,q} (h_{p,q} + f_{p,q}) \\rho_{p,q}\n",
    "  \\end{align} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the electronic contribution to the total energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bottleneck of this Hartree-Fock implementation is the tensor contraction in the formation of the Fock matrix. Because we are not exploiting the sparsity of $\\chi$, it scales as $O(n^4)$ operations for $n$ atoms because of 3 free tensor indices and 1 summed tensor index in the contractions with the dense $\\chi$ tensor.\n",
    "\n",
    "## 4. 2nd-order Moller-Plesset (MP2) perturbation theory\n",
    "\n",
    "The Hartree-Fock approximation does not describe the London dispersion interaction between Argon atoms. We must use many-body perturbation theory to improve our approximation to the ground state of $\\hat{H}$ and include this physical effect. Thankfully, we only need to use the first correction beyond Hartree-Fock theory, which is [second-order Moller-Plesset (MP2) perturbation theory](https://en.wikipedia.org/wiki/Møller–Plesset_perturbation_theory). It is simple, but computationally expensive. The computational bottleneck is the transformation of the tensor $V_{p,q,r,s}$ from the atomic orbital basis to the molecular orbital basis,\n",
    "\n",
    "$$ \\tilde{V}_{a,i,b,j} = \\sum_{p,q,r,s} \\phi_{p,a} \\phi_{q,i} \\phi_{r,b} \\phi_{s,j} V_{p,q,r,s}. $$\n",
    "\n",
    "where we are restricting some indices to occupied orbitals ($i$ and $j$) and some indices to the unoccupied \"virtual\" orbitals ($a$ and $b$). This use of specific orbital labels to denote a restriction of indices is standard notation in quantum chemistry. To make this easier to program, we can partition the molecular orbitals and their energies into occupied and virtual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# partition the orbitals into occupied and virtual"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the factored form of $V_{p,q,r,s}$ to reduce our intermediate memory usage, but we do need to explicitly store $\\tilde{V}_{a,i,b,j}$,\n",
    "\n",
    "$$ \\begin{align}\n",
    "   \\tilde{\\chi}_{a,i,p} &= \\sum_{q,r} \\phi_{q,a} \\phi_{r,i} \\chi_{q,r,p} \\\\\n",
    "   \\tilde{V}_{a,i,b,j} &= \\sum_{p,q} \\tilde{\\chi}_{a,i,p} V_{p,q}^{\\mathrm{ee}} \\tilde{\\chi}_{b,j,q}\n",
    "   \\end{align} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# transform the interaction tensor from atomic orbitals to molecular orbitals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this simple implementation, we are not taking full advantage of the sparisty of the initial $\\chi$ tensor. Unlike the $V$ tensor, the $\\tilde{V}$ tensor is constructed explicitly and stored in memory. The MP2 correction to the total energy is then just a large summation over previously computed and stored quantities,\n",
    "\n",
    "$$ E_{\\mathrm{MP2}} = - \\sum_{a,b \\in \\mathrm{virt}} \\sum_{i,j \\in \\mathrm{occ}} \\frac{2 \\tilde{V}_{a,i,b,j}^2 - \\tilde{V}_{a,i,b,j} \\tilde{V}_{a,j,b,i}}{E_a + E_b - E_i - E_j} . $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the MP2 correction to the total energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bottleneck of MP2 calculations is the formation of $\\tilde{V}$ from $\\tilde{\\chi}$ and $V^{\\mathrm{ee}}$, which scales as $O(n^5)$ operations for $n$ atoms because of 4 free tensor indices and 1 summed tensor index once the contraction has been optimally arranged into intermediate operations. The 4 nested for loops that accumulate the MP2 energy correction have a lower $O(n^4)$ scaling of cost but with a very large prefactor that will usually cause it to dominate the overall cost.\n",
    "\n",
    "## (Extra Credit) 5. Alternative Fock matrix construction\n",
    "\n",
    "The previous implementation of Fock matrix construction,\n",
    "\n",
    "$$ f_{p,q} = h_{p,q} + \\sum_{r,s,t,u} ( 2 \\chi_{p,q,t} \\chi_{r,s,u} - \\chi_{r,q,t} \\chi_{p,s,u} )\n",
    "                             V_{t,u}^{\\mathrm{ee}} \\rho_{r,s}, $$\n",
    "\n",
    "did not utilize the sparsity of the $\\chi$ tensor. As a result, its computational cost scaling is not as low as it could be, and it will have poor performance for large system sizes. We can improve the scaling by making use of the fact that $\\chi_{p,q,r}$ is only non-zero when all of its indices are on the same atom:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_fock_matrix_fast(hamiltonian_matrix, interaction_matrix, density_matrix, model_parameters):\n",
    "    '''Returns the Fock matrix defined by the input Hamiltonian, interaction, & density matrices.'''\n",
    "    ndof = np.size(hamiltonian_matrix,0)\n",
    "    fock_matrix = hamiltonian_matrix.copy()\n",
    "    # Hartree potential term\n",
    "    for p in range(ndof):\n",
    "        for orb_q in orbital_types:\n",
    "            q = ao_index(atom(p), orb_q) # p & q on same atom\n",
    "            for orb_t in orbital_types:\n",
    "                t = ao_index(atom(p), orb_t) # p & t on same atom\n",
    "                chi_pqt = chi_on_atom(orb(p), orb_q, orb_t, model_parameters)\n",
    "                for r in range(ndof):\n",
    "                    for orb_s in orbital_types:\n",
    "                        s = ao_index(atom(r), orb_s) # r & s on same atom\n",
    "                        for orb_u in orbital_types:\n",
    "                            u = ao_index(atom(r), orb_u) # r & u on same atom\n",
    "                            chi_rsu = chi_on_atom(orb(r), orb_s, orb_u, model_parameters)\n",
    "                            fock_matrix[p,q] += 2.0 * chi_pqt * chi_rsu * interaction_matrix[t,u] * density_matrix[r,s]\n",
    "    # Fock exchange term\n",
    "    for p in range(ndof):\n",
    "        for orb_s in orbital_types:\n",
    "            s = ao_index(atom(p), orb_s) # p & s on same atom\n",
    "            for orb_u in orbital_types:\n",
    "                u = ao_index(atom(p), orb_u) # p & u on same atom\n",
    "                chi_psu = chi_on_atom(orb(p), orb_s, orb_u, model_parameters)\n",
    "                for q in range(ndof):\n",
    "                    for orb_r in orbital_types:\n",
    "                        r = ao_index(atom(q), orb_r) # q & r on same atom\n",
    "                        for orb_t in orbital_types:\n",
    "                            t = ao_index(atom(q), orb_t) # q & t on same atom\n",
    "                            chi_rqt = chi_on_atom(orb_r, orb(q), orb_t, model_parameters)\n",
    "                            fock_matrix[p,q] -= chi_rqt * chi_psu * interaction_matrix[t,u] * density_matrix[r,s]\n",
    "    return fock_matrix\n",
    "\n",
    "fock_matrix1 = calculate_fock_matrix(hamiltonian_matrix, interaction_matrix, density_matrix, chi_tensor)\n",
    "fock_matrix2 = calculate_fock_matrix_fast(hamiltonian_matrix, interaction_matrix, density_matrix, model_parameters)\n",
    "print('difference between fock matrix implementations =', np.linalg.norm(fock_matrix1 - fock_matrix2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This algorithm is much more complicated to implement and it is doing something that should be avoided at all costs in Python: deeply nested for loops. While this algorithm has good scaling, we will soon observe a very large cost prefactor in its performance. We can compare the performance of these algorithms by setting up a sequence of Argon clusters of increasing size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from timeit import default_timer as timer\n",
    "\n",
    "def build_fcc_cluster(radius, lattice_constant):\n",
    "    vec1 = np.array([lattice_constant, lattice_constant, 0.0])\n",
    "    vec2 = np.array([lattice_constant, 0.0, lattice_constant])\n",
    "    vec3 = np.array([0.0, lattice_constant, lattice_constant])\n",
    "    max_index = int(radius/np.linalg.norm(vec1))\n",
    "    atomic_coordinates = np.array([ i*vec1 + j*vec2 + k*vec3 for i in range(-max_index, max_index+1)\n",
    "                                                             for j in range(-max_index, max_index+1)\n",
    "                                                             for k in range(-max_index, max_index+1)\n",
    "                                                             if np.linalg.norm(i*vec1 + j*vec2 + k*vec3) <= radius ])\n",
    "    return atomic_coordinates\n",
    "\n",
    "for radius in np.arange(10.0,40.0,10.0):\n",
    "    coord = build_fcc_cluster(radius, 9.9)\n",
    "    \n",
    "    time1 = timer()\n",
    "    density_matrix = calculate_atomic_density_matrix(coord)\n",
    "    hamiltonian_matrix = calculate_hamiltonian_matrix(coord, model_parameters)\n",
    "    interaction_matrix = calculate_interaction_matrix(coord, model_parameters)\n",
    "    chi_tensor = calculate_chi_tensor(coord, model_parameters)\n",
    "    time2 = timer()\n",
    "    print(len(coord))\n",
    "    print('base matrix construction time =',time2-time1)\n",
    "    fock_matrix = calculate_fock_matrix(hamiltonian_matrix, interaction_matrix, density_matrix, chi_tensor)\n",
    "    time3 = timer()\n",
    "    print('Fock matrix construction time (\"slow\") =',time3-time2)\n",
    "    fock_matrix = calculate_fock_matrix_fast(hamiltonian_matrix, interaction_matrix, density_matrix, model_parameters)\n",
    "    time4 = timer()\n",
    "    print('Fock matrix construction time (\"fast\") =',time4-time3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could improve the implementation of Fock matrix construction in Python by organizing much of the work as dense matrix-matrix multiplications on small blocks of the matrix, which involves more sophisticated bookkeeping. Another way to improve its implementation is simply to port it as-is to C or C++, which have very low performance overheads in nested for loops. Thus an implementation that superficially looks very similar can have dramatically improved performance.\n",
    "\n",
    "## (Extra Credit) 6. Semiempirical model parameterization\n",
    "\n",
    "For a semiempirical model to be useful, it must be parameterized to fit a set of reference data, usually from experiments or high-accuracy quantum chemistry simulation data. The parameters in this project have already been fit to data, which was carried out using the fitting procedure implemented in this section. Data fitting is an optimization problem, and we can make use of pre-existing Python optimization tools by encapsulating the fitting process into a function that inputs a list of parameters and outputs a fitting error that we want to minimize. We begin with our reference data, which is the set of occupied orbital energies, Hartree-Fock binding energies, and MP2 binding energies from all-electron calculations of Argon dimers at various separation distances:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# distance, HF binding energy, MP2 binding energy, 6 occupied HF orbital energies\n",
    "reference_dimer = [ [ 5.67, 0.0064861, -0.0037137, -0.620241, -0.595697, -0.595697, -0.584812, -0.584812, -0.561277 ],\n",
    "                    [ 6.24, 0.0022351, -0.0020972, -0.608293, -0.593456, -0.593456, -0.587742, -0.587742, -0.573604 ],\n",
    "                    [ 6.62, 0.0010879, -0.0014491, -0.603106, -0.592621, -0.592621, -0.588897, -0.588897, -0.578874 ],\n",
    "                    [ 6.99, 0.0005255, -0.0010120, -0.599442, -0.592082, -0.592082, -0.589651, -0.589651, -0.582576 ],\n",
    "                    [ 7.18, 0.0003642, -0.0008494, -0.598041, -0.591886, -0.591886, -0.589922, -0.589922, -0.583988 ],\n",
    "                    [ 7.56, 0.0001737, -0.0006037, -0.595889, -0.591599, -0.591599, -0.590315, -0.590315, -0.586154 ],\n",
    "                    [ 7.94, 0.0000818, -0.0004342, -0.594388, -0.591408, -0.591408, -0.590569, -0.590569, -0.587662 ],\n",
    "                    [ 8.32, 0.0000377, -0.0003159, -0.593344, -0.591281, -0.591281, -0.590731, -0.590731, -0.588708 ],\n",
    "                    [ 9.26, 0.0000047, -0.0001502, -0.591934, -0.591118, -0.591118, -0.590927, -0.590927, -0.590120 ],\n",
    "                    [ 10.21, 0.000000, -0.0000757, -0.591378, -0.591059, -0.591059, -0.590993, -0.590993, -0.590675 ] ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then write a function that encapsulates all of the setup and solving from earlier in the notebook to generate model predictions of the reference data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_reference_data(atomic_coordinates, model_parameters):\n",
    "    '''Returns the occupied orbital, HF, & MP2 energies for the input list of atomic coordinates'''\n",
    "    ndof = len(atomic_coordinates)*orbitals_per_atom\n",
    "\n",
    "    energy_ion = calculate_energy_ion(atomic_coordinates)\n",
    "    density_matrix = calculate_atomic_density_matrix(atomic_coordinates)\n",
    "\n",
    "    hamiltonian_matrix = calculate_hamiltonian_matrix(atomic_coordinates, model_parameters)\n",
    "    interaction_matrix = calculate_interaction_matrix(atomic_coordinates, model_parameters)\n",
    "    chi_tensor = calculate_chi_tensor(atomic_coordinates, model_parameters)\n",
    "\n",
    "    density_matrix, fock_matrix = scf_cycle(hamiltonian_matrix, interaction_matrix, density_matrix, chi_tensor)\n",
    "\n",
    "    occupied_energy, virtual_energy, occupied_matrix, virtual_matrix = partition_orbitals(fock_matrix)\n",
    "    energy_scf = calculate_energy_scf(hamiltonian_matrix, fock_matrix, density_matrix)\n",
    "    energy_hf = energy_ion + energy_scf\n",
    "    energy_mp2 = calculate_energy_mp2(fock_matrix, interaction_matrix, chi_tensor)\n",
    "\n",
    "    return occupied_energy, energy_hf, energy_mp2\n",
    "\n",
    "print(calculate_reference_data(np.array([[0.0,0.0,0.0]]),model_parameters))\n",
    "print(calculate_reference_data(np.array([[0.0,0.0,0.0],[6.99,0.0,0.0]]),model_parameters))\n",
    "print(calculate_reference_data(np.array([[0.0,0.0,0.0],[20.0,0.0,0.0]]),model_parameters))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use this function inside yet another function that loops over the reference data, performs simulations on the reference geometries, and returns the root-mean-squared (RMS) deviation between the reference data and model predictions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_error(reference_dimer, model_parameters):\n",
    "    '''Returns the RMS error between the model & reference data for the input set of model parameters'''\n",
    "    atom_occupied, atom_hf, atom_mp2 = calculate_reference_data(np.array([[0,0,0]]), model_parameters)\n",
    "    rms_error = 0.0\n",
    "\n",
    "    for data_list in reference_dimer:\n",
    "        mol_coord = np.array([ [0,0,0], [data_list[0],0,0] ])\n",
    "        mol_occupied, mol_hf, mol_mp2 = calculate_reference_data(mol_coord, model_parameters)\n",
    "        mol_hf -= 2.0*atom_hf\n",
    "        mol_mp2 -= 2.0*atom_mp2\n",
    "        rms_error += (mol_hf - data_list[1])**2 + (mol_mp2 - data_list[2])**2\n",
    "        rms_error += np.linalg.norm(mol_occupied - data_list[3:])**2\n",
    "\n",
    "#    print(np.sqrt(rms_error)/len(reference_dimer))\n",
    "    return np.sqrt(rms_error)/len(reference_dimer)\n",
    "\n",
    "print(\"rms error =\", model_error(reference_dimer,model_parameters))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can use the SciPy optimizer to calculate the optimal model parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.optimize\n",
    "\n",
    "def objective_wrapper(param_list, reference_dimer, model_parameters, lock):\n",
    "    '''Wraps model_error so that the model parameters are in a list'''\n",
    "    new_parameters = {}\n",
    "    i = 0\n",
    "    for key in sorted(model_parameters.keys()):\n",
    "        if lock != None and key in lock:\n",
    "            new_parameters.update( { key : model_parameters[key] } )\n",
    "        else:\n",
    "            new_parameters.update( { key : param_list[i] } )\n",
    "            i += 1\n",
    "    return model_error(reference_dimer, new_parameters)\n",
    "\n",
    "def fit_model(reference_dimer, model_parameters, lock=None):\n",
    "    '''Returns optimized model parameters that best fit the reference data'''\n",
    "\n",
    "    param_list = []\n",
    "    for key in sorted(model_parameters.keys()):\n",
    "        if lock != None and key in lock:\n",
    "            pass\n",
    "        else:\n",
    "            param_list.append( model_parameters[key] )\n",
    "    min_opt = { 'maxiter' : 2000 }\n",
    "    result = scipy.optimize.minimize(objective_wrapper, param_list,\n",
    "                                     (reference_dimer, model_parameters, lock),\n",
    "                                     method='Nelder-Mead', options=min_opt)\n",
    "\n",
    "    opt_parameters = {}\n",
    "    i = 0\n",
    "    for key in sorted(model_parameters.keys()):\n",
    "        if lock != None and key in lock:\n",
    "            opt_parameters.update( { key : model_parameters[key] } )\n",
    "        else:\n",
    "            opt_parameters.update( { key : result.x[i] } )\n",
    "            i += 1\n",
    "    return opt_parameters\n",
    "\n",
    "new_parameters = fit_model(reference_dimer, model_parameters)\n",
    "print(\"new rms error =\", model_error(reference_dimer,new_parameters))\n",
    "print(new_parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is this a good semiempirical model? Not really. It was designed to be simple in form and capture the correct qualitative behavior, but simplicity and accuracy are usually conflicting design principles. The model is not flexible enough to fit the data to very high accuracy, but its parameters also cannot be uniquely specified by the reference data because we aren't probing the system enough. We can identify this problem by examining the eigenvalues of the Hessian matrix for the objective function that we are minimizing to fit the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_hessian(reference_dimer, model_parameters, lock=None):\n",
    "    param_list = []\n",
    "    for key in sorted(model_parameters.keys()):\n",
    "        if lock != None and key in lock:\n",
    "            pass\n",
    "        else:\n",
    "            param_list.append( model_parameters[key] )\n",
    "    num_param = len(param_list)\n",
    "    hessian_matrix = np.zeros((num_param,num_param))\n",
    "    dx = 0.01\n",
    "    for i in range(num_param):\n",
    "        for j in range(num_param):\n",
    "            param_list[i] += dx\n",
    "            param_list[j] += dx\n",
    "            hes_pp = objective_wrapper(param_list, reference_dimer, model_parameters, lock)\n",
    "            param_list[i] -= dx\n",
    "            param_list[j] -= dx\n",
    "\n",
    "            param_list[i] -= dx\n",
    "            param_list[j] -= dx\n",
    "            hes_mm = objective_wrapper(param_list, reference_dimer, model_parameters, lock)\n",
    "            param_list[i] += dx\n",
    "            param_list[j] += dx\n",
    "\n",
    "            param_list[i] += dx\n",
    "            param_list[j] -= dx\n",
    "            hes_pm = objective_wrapper(param_list, reference_dimer, model_parameters, lock)\n",
    "            param_list[i] -= dx\n",
    "            param_list[j] += dx\n",
    "\n",
    "            param_list[i] -= dx\n",
    "            param_list[j] += dx\n",
    "            hes_mp = objective_wrapper(param_list, reference_dimer, model_parameters, lock)\n",
    "            param_list[i] += dx\n",
    "            param_list[j] -= dx\n",
    "\n",
    "            hessian_matrix[i,j] = (hes_pp - hes_mp - hes_pm + hes_mm)/(4*dx**2)\n",
    "    return hessian_matrix\n",
    "\n",
    "hessian_matrix = model_hessian(reference_dimer, new_parameters)\n",
    "vals, vecs = np.linalg.eigh(hessian_matrix)\n",
    "print('model Hessian eigenvalues:\\n',vals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Small eigenvalues correspond to directions in parameter space where we can change parameters without changing how well the model fits the data. This problem is not specific to our simple QM project - semiempirical model building has persistent problems with model indeterminacy as they add more parameters and details to improve accuracy. Because these models have historically be used to reproduce a limited set of observables (e.g. heat of formation), they've been limited in what reference data they can use for model fitting.\n",
    "\n",
    "We can use this same model and fitting process for other noble gases such as Neon by applying it to a different set of reference data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neon_ref_data = [\n",
    "[ 5.29301 , 0.000356767 , -0.000275745 , -0.857733 , -0.851311 , -0.851311 , -0.849117 , -0.849117 , -0.843045 ],\n",
    "[ 5.48204 , 0.000222991 , -0.000229153 , -0.856225 , -0.851119 , -0.851119 , -0.849443 , -0.849443 , -0.844581 ],\n",
    "[ 5.67108 , 0.000138635 , -0.000190295 , -0.855025 , -0.850968 , -0.850968 , -0.849687 , -0.849687 , -0.845800 ],\n",
    "[ 5.86011 , 8.54899e-05 , -0.000157984 , -0.854072 , -0.850849 , -0.850849 , -0.849869 , -0.849869 , -0.846765 ],\n",
    "[ 6.04915 , 5.21333e-05 , -0.000131269 , -0.853315 , -0.850757 , -0.850757 , -0.850006 , -0.850006 , -0.847529 ],\n",
    "[ 6.23819 , 3.13588e-05 , -0.000109318 , -0.852714 , -0.850684 , -0.850684 , -0.850108 , -0.850108 , -0.848135 ],\n",
    "[ 6.42722 , 1.85686e-05 , -9.13551e-05 , -0.852237 , -0.850626 , -0.850626 , -0.850185 , -0.850185 , -0.848614 ],\n",
    "[ 6.61626 , 1.07945e-05 , -7.66555e-05 , -0.851859 , -0.850582 , -0.850582 , -0.850244 , -0.850244 , -0.848993 ],\n",
    "[ 6.80529 , 6.11193e-06 , -6.45702e-05 , -0.851559 , -0.850547 , -0.850547 , -0.850287 , -0.850287 , -0.849293 ],\n",
    "[ 6.99433 , 3.28677e-06 , -5.45556e-05 , -0.851322 , -0.850520 , -0.850520 , -0.850321 , -0.850321 , -0.849531 ],\n",
    "[ 7.18336 , 1.55256e-06 , -4.61875e-05 , -0.851134 , -0.850499 , -0.850499 , -0.850346 , -0.850346 , -0.849719 ],\n",
    "[ 7.3724 , 4.57149e-07 , -3.91539e-05 , -0.850986 , -0.850483 , -0.850483 , -0.850366 , -0.850366 , -0.849868 ],\n",
    "[ 7.56144 , -2.50044e-07 , -3.3231e-05 , -0.850869 , -0.850471 , -0.850471 , -0.850380 , -0.850380 , -0.849986 ],\n",
    "[ 7.75047 , -7.01547e-07 , -2.82527e-05 , -0.850776 , -0.850461 , -0.850461 , -0.850392 , -0.850392 , -0.850079 ],\n",
    "[ 7.93951 , -9.693e-07 , -2.40853e-05 , -0.850703 , -0.850454 , -0.850454 , -0.850400 , -0.850400 , -0.850153 ],\n",
    "[ 8.12854 , -1.09914e-06 , -2.06105e-05 , -0.850645 , -0.850448 , -0.850448 , -0.850407 , -0.850407 , -0.850211 ],\n",
    "[ 8.31758 , -1.12749e-06 , -1.7719e-05 , -0.850600 , -0.850444 , -0.850444 , -0.850412 , -0.850412 , -0.850257 ],\n",
    "[ 8.50662 , -1.08659e-06 , -1.53105e-05 , -0.850564 , -0.850440 , -0.850440 , -0.850416 , -0.850416 , -0.850294 ],\n",
    "[ 8.69565 , -1.00439e-06 , -1.32963e-05 , -0.850536 , -0.850438 , -0.850438 , -0.850419 , -0.850419 , -0.850322 ],\n",
    "[ 8.88469 , -9.0295e-07 , -1.16011e-05 , -0.850513 , -0.850436 , -0.850436 , -0.850421 , -0.850421 , -0.850345 ],\n",
    "[ 9.07372 , -7.97673e-07 , -1.01637e-05 , -0.850495 , -0.850434 , -0.850434 , -0.850423 , -0.850423 , -0.850363 ],\n",
    "[ 9.26276 , -6.977e-07 , -8.93627e-06 , -0.850482 , -0.850433 , -0.850433 , -0.850425 , -0.850425 , -0.850377 ],\n",
    "[ 9.4518 , -6.07201e-07 , -7.88179e-06 , -0.850470 , -0.850432 , -0.850432 , -0.850426 , -0.850426 , -0.850389 ],\n",
    "[ 9.64083 , -5.27071e-07 , -6.97177e-06 , -0.850462 , -0.850432 , -0.850432 , -0.850427 , -0.850427 , -0.850398 ],\n",
    "[ 9.82987 , -4.56455e-07 , -6.18389e-06 , -0.850455 , -0.850431 , -0.850431 , -0.850428 , -0.850428 , -0.850405 ],\n",
    "[ 10.0189 , -3.93899e-07 , -5.50014e-06 , -0.850450 , -0.850431 , -0.850431 , -0.850428 , -0.850428 , -0.850410 ],\n",
    "[ 10.2079 , -3.38041e-07 , -4.90551e-06 , -0.850445 , -0.850431 , -0.850431 , -0.850429 , -0.850429 , -0.850415 ],\n",
    "[ 10.397 , -2.87888e-07 , -4.38722e-06 , -0.850442 , -0.850431 , -0.850431 , -0.850429 , -0.850429 , -0.850418 ],\n",
    "[ 10.586 , -2.42876e-07 , -3.93422e-06 , -0.850439 , -0.850431 , -0.850431 , -0.850429 , -0.850429 , -0.850421 ],\n",
    "[ 10.775 , -2.02724e-07 , -3.53699e-06 , -0.850437 , -0.850430 , -0.850430 , -0.850429 , -0.850429 , -0.850423 ]\n",
    "]\n",
    "\n",
    "neon_param = fit_model(neon_ref_data, model_parameters)\n",
    "print(\"neon rms error =\", model_error(neon_ref_data,neon_param))\n",
    "print(neon_param)\n",
    "\n",
    "#completed output of Neon parameter fitting saved for posterity:\n",
    "neon_param = {'coulomb_p': -0.010255409806855187,\n",
    "              'coulomb_s': 0.4536486561938202,\n",
    "              'dipole': 1.6692376991516769,\n",
    "              'energy_p': -3.1186533988406335,\n",
    "              'energy_s': 11.334912902362603,\n",
    "              'r_hop': 2.739689713337267,\n",
    "              'r_pseudo': 1.1800779720963734,\n",
    "              't_pp1': -0.029546671673199854,\n",
    "              't_pp2': -0.0041958662271044875,\n",
    "              't_sp': 0.000450562836426027,\n",
    "              't_ss': 0.0289251941290921,\n",
    "              'v_pseudo': -0.015945813280635074}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}