Argon_SQM.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": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "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": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "coordinates =\n",
      " [[0. 0. 0.]\n",
      " [3. 4. 5.]]\n",
      "# of atoms = 2\n"
     ]
    }
   ],
   "source": [
    "atomic_coordinates = np.array([ [0.0,0.0,0.0], [3.0,4.0,5.0] ])\n",
    "number_of_atoms = len(atomic_coordinates)\n",
    "\n",
    "print('coordinates =\\n',atomic_coordinates)\n",
    "print('# of atoms =',number_of_atoms)"
   ]
  },
  {
   "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": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "all orbitals = ['s', 'px', 'py', 'pz']\n",
      "p orbitals = ['px', 'py', 'pz']\n"
     ]
    }
   ],
   "source": [
    "ionic_charge = 6\n",
    "orbital_types = ['s', 'px' ,'py', 'pz']\n",
    "orbitals_per_atom = len(orbital_types)\n",
    "\n",
    "p_orbitals = orbital_types[1:]\n",
    "print('all orbitals =', orbital_types)\n",
    "print('p orbitals =', p_orbitals)"
   ]
  },
  {
   "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": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "index 0 atom 0 orbital s\n",
      "index 1 atom 0 orbital px\n",
      "index 2 atom 0 orbital py\n",
      "index 3 atom 0 orbital pz\n",
      "index 4 atom 1 orbital s\n",
      "index 5 atom 1 orbital px\n",
      "index 6 atom 1 orbital py\n",
      "index 7 atom 1 orbital pz\n",
      "index test:\n",
      "0 0\n",
      "1 1\n",
      "2 2\n",
      "3 3\n",
      "4 4\n",
      "5 5\n",
      "6 6\n",
      "7 7\n"
     ]
    }
   ],
   "source": [
    "def atom(ao_index):\n",
    "    '''Returns the atom index part of an atomic orbital index.'''\n",
    "    return ao_index // orbitals_per_atom\n",
    "\n",
    "def orb(ao_index):\n",
    "    '''Returns the orbital type of an atomic orbital index.'''\n",
    "    orb_index = ao_index % orbitals_per_atom\n",
    "    return orbital_types[orb_index]\n",
    "\n",
    "def ao_index(atom_p,orb_p):\n",
    "    '''Returns the atomic orbital index for a given atom index and orbital type.'''\n",
    "    p = atom_p*orbitals_per_atom\n",
    "    p += orbital_types.index(orb_p)\n",
    "    return p\n",
    "\n",
    "for index in range(number_of_atoms*orbitals_per_atom):\n",
    "    print('index',index,'atom',atom(index),'orbital',orb(index))\n",
    "\n",
    "print('index test:')\n",
    "for index in range(number_of_atoms*orbitals_per_atom):\n",
    "    print(index, ao_index(atom(index),orb(index)))"
   ]
  },
  {
   "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": 5,
   "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": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hopping tests:\n",
      "s s = 0.03365982238611262 0.03365982238611262\n",
      "s px = 0.029154833035109226 -0.029154833035109226\n",
      "px px = 0.0804163845390335 -0.0804163845390335\n",
      "s py = -0.0\n",
      "py px = -0.0\n",
      "py py = -0.01393611496959445 -0.01393611496959445\n"
     ]
    }
   ],
   "source": [
    "vec = { 'px':[1,0,0], 'py':[0,1,0], 'pz':[0,0,1] }\n",
    "\n",
    "def hopping_energy(o1, o2, r12, model_parameters):\n",
    "    '''Returns the hopping matrix element for a pair of orbitals of type o1 & o2 separated by a vector r12.'''\n",
    "    r12_rescaled = r12 / model_parameters['r_hop']\n",
    "    r12_length = np.linalg.norm(r12_rescaled)\n",
    "    ans = np.exp( 1.0 - r12_length**2 )\n",
    "    if o1 == 's' and o2 == 's':\n",
    "        ans *= model_parameters['t_ss']\n",
    "    if o1 == 's' and o2 in p_orbitals:\n",
    "        ans *= np.dot(vec[o2], r12_rescaled) * model_parameters['t_sp']\n",
    "    if o2 == 's' and o1 in p_orbitals:\n",
    "        ans *= -np.dot(vec[o1], r12_rescaled)* model_parameters['t_sp']\n",
    "    if o1 in p_orbitals and o2 in p_orbitals:\n",
    "        ans *= ( (r12_length**2) * np.dot(vec[o1], vec[o2]) * model_parameters['t_pp2']\n",
    "                 - np.dot(vec[o1], r12_rescaled) * np.dot(vec[o2], r12_rescaled)\n",
    "                 * ( model_parameters['t_pp1'] + model_parameters['t_pp2'] ) )\n",
    "    return ans\n",
    "\n",
    "hop_vector = np.array([-model_parameters['r_hop'],0,0])\n",
    "print('hopping tests:')\n",
    "print('s s =', hopping_energy('s','s',hop_vector,model_parameters), model_parameters['t_ss'])\n",
    "print('s px =', hopping_energy('s','px',hop_vector,model_parameters), model_parameters['t_sp'])\n",
    "print('px px =', hopping_energy('px','px',hop_vector,model_parameters), model_parameters['t_pp1'])\n",
    "print('s py =', hopping_energy('s','py',hop_vector,model_parameters))\n",
    "print('py px =', hopping_energy('py','px',hop_vector,model_parameters))\n",
    "print('py py =', hopping_energy('py','py',hop_vector,model_parameters), model_parameters['t_pp2'])"
   ]
  },
  {
   "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": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coulomb tests:\n",
      "s s = 1.0\n",
      "s px = -1.0\n",
      "px px = -2.0\n",
      "s py = 0.0\n",
      "py px = 0.0\n",
      "py py = 1.0\n"
     ]
    }
   ],
   "source": [
    "def coulomb_energy(o1, o2, r12):\n",
    "    '''Returns the Coulomb matrix element for a pair of multipoles of type o1 & o2 separated by a vector r12.'''\n",
    "    r12_length = np.linalg.norm(r12)\n",
    "    if o1 == 's' and o2 == 's':\n",
    "        ans = 1.0 / r12_length\n",
    "    if o1 == 's' and o2 in p_orbitals:\n",
    "        ans = np.dot(vec[o2], r12) / r12_length**3\n",
    "    if o2 == 's' and o1 in p_orbitals:\n",
    "        ans = -np.dot(vec[o1], r12) / r12_length**3\n",
    "    if o1 in p_orbitals and o2 in p_orbitals:\n",
    "        ans = ( np.dot(vec[o1], vec[o2]) / r12_length**3\n",
    "               - 3.0 * np.dot(vec[o1], r12) * np.dot(vec[o2], r12) / r12_length**5 )\n",
    "    return ans \n",
    "\n",
    "coulomb_vector = np.array([-1,0,0])\n",
    "print('Coulomb tests:')\n",
    "print('s s =', coulomb_energy('s','s',coulomb_vector))\n",
    "print('s px =', coulomb_energy('s','px',coulomb_vector))\n",
    "print('px px =', coulomb_energy('px','px',coulomb_vector))\n",
    "print('s py =', coulomb_energy('s','py',coulomb_vector))\n",
    "print('py px =', coulomb_energy('py','px',coulomb_vector))\n",
    "print('py py =', coulomb_energy('py','py',coulomb_vector))"
   ]
  },
  {
   "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": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pseudopotential tests:\n",
      "s = 0.022972992186364977 0.022972992186364977\n",
      "px = 0.04594598437272995\n",
      "py = -0.0\n"
     ]
    }
   ],
   "source": [
    "def pseudopotential_energy(o, r, model_parameters):\n",
    "    '''Returns the energy of a pseudopotential between a multipole of type o and an atom separated by a vector r.'''\n",
    "    ans = model_parameters['v_pseudo']\n",
    "    r_rescaled = r / model_parameters['r_pseudo']\n",
    "    r_length = np.linalg.norm(r_rescaled)\n",
    "    ans *= np.exp( 1.0 - r_length**2 )\n",
    "    if o in p_orbitals:\n",
    "        ans *= -2.0 * np.dot(vec[o], r_rescaled)\n",
    "    return ans\n",
    "\n",
    "pseudo_vector = np.array([-model_parameters['r_pseudo'],0,0])\n",
    "print('pseudopotential tests:')\n",
    "print('s =', pseudopotential_energy('s',pseudo_vector, model_parameters), model_parameters['v_pseudo'])\n",
    "print('px =', pseudopotential_energy('px',pseudo_vector, model_parameters))\n",
    "print('py =', pseudopotential_energy('py',pseudo_vector, model_parameters))"
   ]
  },
  {
   "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": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E_ion = 5.091168824543142\n"
     ]
    }
   ],
   "source": [
    "def calculate_energy_ion(atomic_coordinates):\n",
    "    '''Returns the ionic contribution to the total energy for an input list of atomic coordinates.'''\n",
    "    energy_ion = 0.0\n",
    "    for i, r_i in enumerate(atomic_coordinates):\n",
    "        for j, r_j in enumerate(atomic_coordinates):\n",
    "            if i < j:\n",
    "                energy_ion += (ionic_charge**2)*coulomb_energy('s', 's', r_i - r_j)\n",
    "    return energy_ion\n",
    "\n",
    "energy_ion = calculate_energy_ion(atomic_coordinates)\n",
    "print('E_ion =', energy_ion)"
   ]
  },
  {
   "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": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "V_ion = [-0.8 -0.1 -0.1 -0.1 -0.8  0.1  0.1  0.1]\n"
     ]
    }
   ],
   "source": [
    "def calculate_potential_vector(atomic_coordinates, model_parameters):\n",
    "    '''Returns the electron-ion potential energy vector for an input list of atomic coordinates.'''\n",
    "    ndof = len(atomic_coordinates)*orbitals_per_atom\n",
    "    potential_vector = np.zeros(ndof)\n",
    "    for p in range(ndof):\n",
    "        potential_vector[p] = 0.0\n",
    "        for atom_i,r_i in enumerate(atomic_coordinates):\n",
    "            r_pi = atomic_coordinates[atom(p)] - r_i\n",
    "            if atom_i != atom(p):\n",
    "                potential_vector[p] += ( pseudopotential_energy(orb(p), r_pi, model_parameters) \n",
    "                                         - ionic_charge * coulomb_energy(orb(p), 's', r_pi) )\n",
    "    return potential_vector\n",
    "\n",
    "potential_vector = calculate_potential_vector(atomic_coordinates, model_parameters)\n",
    "np.set_printoptions(precision=1)\n",
    "print('V_ion =', 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": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "V_ee =\n",
      " [[ 3.6e-01  0.0e+00  0.0e+00  0.0e+00  1.4e-01 -8.5e-03 -1.1e-02 -1.4e-02]\n",
      " [ 0.0e+00 -3.3e-03  0.0e+00  0.0e+00  8.5e-03  1.3e-03 -2.0e-03 -2.5e-03]\n",
      " [ 0.0e+00  0.0e+00 -3.3e-03  0.0e+00  1.1e-02 -2.0e-03  1.1e-04 -3.4e-03]\n",
      " [ 0.0e+00  0.0e+00  0.0e+00 -3.3e-03  1.4e-02 -2.5e-03 -3.4e-03 -1.4e-03]\n",
      " [ 1.4e-01  8.5e-03  1.1e-02  1.4e-02  3.6e-01  0.0e+00  0.0e+00  0.0e+00]\n",
      " [-8.5e-03  1.3e-03 -2.0e-03 -2.5e-03  0.0e+00 -3.3e-03  0.0e+00  0.0e+00]\n",
      " [-1.1e-02 -2.0e-03  1.1e-04 -3.4e-03  0.0e+00  0.0e+00 -3.3e-03  0.0e+00]\n",
      " [-1.4e-02 -2.5e-03 -3.4e-03 -1.4e-03  0.0e+00  0.0e+00  0.0e+00 -3.3e-03]]\n"
     ]
    }
   ],
   "source": [
    "def calculate_interaction_matrix(atomic_coordinates, model_parameters):\n",
    "    '''Returns the electron-electron interaction energy matrix for an input list of atomic coordinates.'''\n",
    "    ndof = len(atomic_coordinates)*orbitals_per_atom\n",
    "    interaction_matrix = np.zeros( (ndof,ndof) )\n",
    "    for p in range(ndof):\n",
    "        for q in range(ndof):\n",
    "            if atom(p) != atom(q):\n",
    "                r_pq = atomic_coordinates[atom(p)] - atomic_coordinates[atom(q)]\n",
    "                interaction_matrix[p,q] = coulomb_energy(orb(p), orb(q), r_pq)\n",
    "            if p == q and orb(p) == 's':\n",
    "                interaction_matrix[p,q] = model_parameters['coulomb_s']\n",
    "            if p == q and orb(p) in p_orbitals:\n",
    "                interaction_matrix[p,q] = model_parameters['coulomb_p']                \n",
    "    return interaction_matrix\n",
    "\n",
    "interaction_matrix = calculate_interaction_matrix(atomic_coordinates, model_parameters)\n",
    "print('V_ee =\\n', 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": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi =\n",
      " [[[1.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  2.8 0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  2.8 0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  2.8 0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]]\n",
      "\n",
      " [[0.  2.8 0.  0.  0.  0.  0.  0. ]\n",
      "  [1.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]]\n",
      "\n",
      " [[0.  0.  2.8 0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [1.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]]\n",
      "\n",
      " [[0.  0.  0.  2.8 0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [1.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]]\n",
      "\n",
      " [[0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  1.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  2.8 0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  2.8 0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  2.8]]\n",
      "\n",
      " [[0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  2.8 0.  0. ]\n",
      "  [0.  0.  0.  0.  1.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]]\n",
      "\n",
      " [[0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  2.8 0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  1.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]]\n",
      "\n",
      " [[0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  2.8]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      "  [0.  0.  0.  0.  1.  0.  0.  0. ]]]\n"
     ]
    }
   ],
   "source": [
    "def chi_on_atom(o1, o2, o3, model_parameters):\n",
    "    '''Returns the value of the chi tensor for 3 orbital indices on the same atom.'''\n",
    "    if o1 == o2 and o3 == 's':\n",
    "        return 1.0\n",
    "    if o1 == o3 and o3 in p_orbitals and o2 == 's':\n",
    "        return model_parameters['dipole']\n",
    "    if o2 == o3 and o3 in p_orbitals and o1 == 's':\n",
    "        return model_parameters['dipole']\n",
    "    return 0.0\n",
    "\n",
    "def calculate_chi_tensor(atomic_coordinates, model_parameters):\n",
    "    '''Returns the chi tensor for an input list of atomic coordinates'''\n",
    "    ndof = len(atomic_coordinates)*orbitals_per_atom\n",
    "    chi_tensor = np.zeros( (ndof,ndof,ndof) )\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_r in orbital_types:\n",
    "                r = ao_index(atom(p), orb_r) # p & r on same atom\n",
    "                chi_tensor[p,q,r] = chi_on_atom(orb(p), orb(q), orb(r), model_parameters)\n",
    "    return chi_tensor\n",
    "\n",
    "chi_tensor = calculate_chi_tensor(atomic_coordinates, model_parameters)\n",
    "print('chi =\\n',chi_tensor)"
   ]
  },
  {
   "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": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hamiltonian matrix =\n",
      " [[ 2.3e+00 -1.4e-01 -1.9e-01 -2.4e-01  6.5e-04  5.3e-04  7.1e-04  8.9e-04]\n",
      " [-1.4e-01 -3.2e+00  0.0e+00  0.0e+00 -5.3e-04  2.9e-04  2.2e-03  2.7e-03]\n",
      " [-1.9e-01  0.0e+00 -3.2e+00  0.0e+00 -7.1e-04  2.2e-03  1.6e-03  3.6e-03]\n",
      " [-2.4e-01  0.0e+00  0.0e+00 -3.2e+00 -8.9e-04  2.7e-03  3.6e-03  3.2e-03]\n",
      " [ 6.5e-04 -5.3e-04 -7.1e-04 -8.9e-04  2.3e+00  1.4e-01  1.9e-01  2.4e-01]\n",
      " [ 5.3e-04  2.9e-04  2.2e-03  2.7e-03  1.4e-01 -3.2e+00  0.0e+00  0.0e+00]\n",
      " [ 7.1e-04  2.2e-03  1.6e-03  3.6e-03  1.9e-01  0.0e+00 -3.2e+00  0.0e+00]\n",
      " [ 8.9e-04  2.7e-03  3.6e-03  3.2e-03  2.4e-01  0.0e+00  0.0e+00 -3.2e+00]]\n"
     ]
    }
   ],
   "source": [
    "def calculate_hamiltonian_matrix(atomic_coordinates, model_parameters):\n",
    "    '''Returns the 1-body Hamiltonian matrix for an input list of atomic coordinates.'''\n",
    "    ndof = len(atomic_coordinates)*orbitals_per_atom\n",
    "    hamiltonian_matrix = np.zeros( (ndof,ndof) )\n",
    "    potential_vector = calculate_potential_vector(atomic_coordinates, model_parameters)\n",
    "    for p in range(ndof):\n",
    "        for q in range(ndof):\n",
    "            if atom(p) != atom(q):\n",
    "                r_pq = atomic_coordinates[atom(p)] - atomic_coordinates[atom(q)]\n",
    "                hamiltonian_matrix[p,q] = hopping_energy(orb(p), orb(q), r_pq, model_parameters)\n",
    "            if atom(p) == atom(q):\n",
    "                if p == q and orb(p) == 's':\n",
    "                    hamiltonian_matrix[p,q] += model_parameters['energy_s']\n",
    "                if p == q and orb(p) in p_orbitals:\n",
    "                    hamiltonian_matrix[p,q] += model_parameters['energy_p']\n",
    "                for orb_r in orbital_types:\n",
    "                    r = ao_index(atom(p), orb_r)\n",
    "                    hamiltonian_matrix[p,q] += ( chi_on_atom(orb(p), orb(q), orb_r, model_parameters)\n",
    "                                                 * potential_vector[r] )\n",
    "    return hamiltonian_matrix\n",
    "\n",
    "hamiltonian_matrix = calculate_hamiltonian_matrix(atomic_coordinates, model_parameters)\n",
    "print('hamiltonian matrix =\\n',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": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "atomic density matrix =\n",
      " [[0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "orbital_occupation = { 's':0, 'px':1, 'py':1, 'pz':1 }\n",
    "\n",
    "def calculate_atomic_density_matrix(atomic_coordinates):\n",
    "    '''Returns a trial 1-electron density matrix for an input list of atomic coordinates.'''\n",
    "    ndof = len(atomic_coordinates)*orbitals_per_atom\n",
    "    density_matrix = np.zeros( (ndof,ndof) )\n",
    "    for p in range(ndof):\n",
    "        density_matrix[p,p] = orbital_occupation[orb(p)]\n",
    "    return density_matrix\n",
    "\n",
    "density_matrix = calculate_atomic_density_matrix(atomic_coordinates)\n",
    "print('atomic density matrix =\\n', density_matrix)"
   ]
  },
  {
   "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": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "atomic fock matrix =\n",
      " [[ 5.4e+00  2.5e-04  3.3e-04  4.2e-04  6.5e-04  5.3e-04  7.1e-04  8.9e-04]\n",
      " [ 2.5e-04 -5.9e-01  0.0e+00  0.0e+00 -5.3e-04  2.9e-04  2.2e-03  2.7e-03]\n",
      " [ 3.3e-04  0.0e+00 -5.9e-01  0.0e+00 -7.1e-04  2.2e-03  1.6e-03  3.6e-03]\n",
      " [ 4.2e-04  0.0e+00  0.0e+00 -5.9e-01 -8.9e-04  2.7e-03  3.6e-03  3.2e-03]\n",
      " [ 6.5e-04 -5.3e-04 -7.1e-04 -8.9e-04  5.4e+00 -2.5e-04 -3.3e-04 -4.2e-04]\n",
      " [ 5.3e-04  2.9e-04  2.2e-03  2.7e-03 -2.5e-04 -5.9e-01  0.0e+00  0.0e+00]\n",
      " [ 7.1e-04  2.2e-03  1.6e-03  3.6e-03 -3.3e-04  0.0e+00 -5.9e-01  0.0e+00]\n",
      " [ 8.9e-04  2.7e-03  3.6e-03  3.2e-03 -4.2e-04  0.0e+00  0.0e+00 -5.9e-01]]\n"
     ]
    }
   ],
   "source": [
    "def calculate_fock_matrix(hamiltonian_matrix, interaction_matrix, density_matrix, chi_tensor):\n",
    "    '''Returns the Fock matrix defined by the input Hamiltonian, interaction, & density matrices.'''\n",
    "    fock_matrix = hamiltonian_matrix.copy()\n",
    "    fock_matrix += 2.0*np.einsum('pqt,rsu,tu,rs',\n",
    "                                 chi_tensor, chi_tensor, interaction_matrix, density_matrix, optimize=True)\n",
    "    fock_matrix -= np.einsum('rqt,psu,tu,rs',\n",
    "                             chi_tensor, chi_tensor, interaction_matrix, density_matrix, optimize=True)\n",
    "    return fock_matrix\n",
    "\n",
    "fock_matrix = calculate_fock_matrix(hamiltonian_matrix, interaction_matrix, density_matrix, chi_tensor)\n",
    "print('atomic fock matrix =\\n', fock_matrix)"
   ]
  },
  {
   "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": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "corrected density matrix =\n",
      " [[ 5.4e-08 -4.2e-05 -5.6e-05 -7.0e-05 -4.1e-08 -8.9e-05 -1.2e-04 -1.5e-04]\n",
      " [-4.2e-05  1.0e+00 -1.3e-08 -1.6e-08  8.9e-05 -7.5e-09 -1.0e-08 -1.2e-08]\n",
      " [-5.6e-05 -1.3e-08  1.0e+00 -2.2e-08  1.2e-04 -1.0e-08 -1.3e-08 -1.7e-08]\n",
      " [-7.0e-05 -1.6e-08 -2.2e-08  1.0e+00  1.5e-04 -1.2e-08 -1.7e-08 -2.1e-08]\n",
      " [-4.1e-08  8.9e-05  1.2e-04  1.5e-04  5.4e-08  4.2e-05  5.6e-05  7.0e-05]\n",
      " [-8.9e-05 -7.5e-09 -1.0e-08 -1.2e-08  4.2e-05  1.0e+00 -1.3e-08 -1.6e-08]\n",
      " [-1.2e-04 -1.0e-08 -1.3e-08 -1.7e-08  5.6e-05 -1.3e-08  1.0e+00 -2.2e-08]\n",
      " [-1.5e-04 -1.2e-08 -1.7e-08 -2.1e-08  7.0e-05 -1.6e-08 -2.2e-08  1.0e+00]]\n"
     ]
    }
   ],
   "source": [
    "def calculate_density_matrix(fock_matrix):\n",
    "    '''Returns the 1-electron density matrix defined by the input Fock matrix.'''\n",
    "    num_occ = (ionic_charge//2)*np.size(fock_matrix,0)//orbitals_per_atom\n",
    "    orbital_energy, orbital_matrix = np.linalg.eigh(fock_matrix)\n",
    "    occupied_matrix = orbital_matrix[:,:num_occ]\n",
    "    density_matrix = occupied_matrix @ occupied_matrix.T\n",
    "    return density_matrix\n",
    "\n",
    "density_matrix = calculate_density_matrix(fock_matrix)\n",
    "print('corrected density matrix =\\n', density_matrix)"
   ]
  },
  {
   "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": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SCF density matrix:\n",
      " [[ 5.8e-08 -4.4e-05 -5.8e-05 -7.3e-05 -4.5e-08 -9.2e-05 -1.2e-04 -1.5e-04]\n",
      " [-4.4e-05  1.0e+00 -1.4e-08 -1.7e-08  9.2e-05 -8.1e-09 -1.1e-08 -1.3e-08]\n",
      " [-5.8e-05 -1.4e-08  1.0e+00 -2.3e-08  1.2e-04 -1.1e-08 -1.4e-08 -1.8e-08]\n",
      " [-7.3e-05 -1.7e-08 -2.3e-08  1.0e+00  1.5e-04 -1.3e-08 -1.8e-08 -2.2e-08]\n",
      " [-4.5e-08  9.2e-05  1.2e-04  1.5e-04  5.8e-08  4.4e-05  5.8e-05  7.3e-05]\n",
      " [-9.2e-05 -8.1e-09 -1.1e-08 -1.3e-08  4.4e-05  1.0e+00 -1.4e-08 -1.7e-08]\n",
      " [-1.2e-04 -1.1e-08 -1.4e-08 -1.8e-08  5.8e-05 -1.4e-08  1.0e+00 -2.3e-08]\n",
      " [-1.5e-04 -1.3e-08 -1.8e-08 -2.2e-08  7.3e-05 -1.7e-08 -2.3e-08  1.0e+00]]\n"
     ]
    }
   ],
   "source": [
    "def scf_cycle(hamiltonian_matrix, interaction_matrix, density_matrix, chi_tensor,\n",
    "              max_scf_iterations = 100, mixing_fraction = 0.25, convergence_tolerance = 1e-4):\n",
    "    '''Returns converged density & Fock matrices defined by the input Hamiltonian, interaction, & density matrices.'''\n",
    "    old_density_matrix = density_matrix.copy()\n",
    "    for iteration in range(max_scf_iterations):\n",
    "        new_fock_matrix = calculate_fock_matrix(hamiltonian_matrix, interaction_matrix, old_density_matrix, chi_tensor)\n",
    "        new_density_matrix = calculate_density_matrix(new_fock_matrix)\n",
    "\n",
    "        error_norm = np.linalg.norm( old_density_matrix - new_density_matrix )\n",
    "        if error_norm < convergence_tolerance:\n",
    "            return new_density_matrix, new_fock_matrix\n",
    "\n",
    "        old_density_matrix = (mixing_fraction * new_density_matrix\n",
    "                              + (1.0 - mixing_fraction) * old_density_matrix)\n",
    "    print(\"WARNING: SCF cycle didn't converge\")\n",
    "    return new_density_matrix, new_fock_matrix\n",
    "\n",
    "density_matrix, fock_matrix = scf_cycle(hamiltonian_matrix, interaction_matrix, density_matrix, chi_tensor)\n",
    "print('SCF density matrix:\\n',density_matrix)"
   ]
  },
  {
   "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": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hartree-Fock energy = -17.901180613943193\n"
     ]
    }
   ],
   "source": [
    "def calculate_energy_scf(hamiltonian_matrix, fock_matrix, density_matrix):\n",
    "    '''Returns the Hartree-Fock total energy defined by the input Hamiltonian, Fock, & density matrices.'''\n",
    "    energy_scf = np.einsum('pq,pq',hamiltonian_matrix + fock_matrix,density_matrix)\n",
    "    return energy_scf\n",
    "\n",
    "energy_scf = calculate_energy_scf(hamiltonian_matrix, fock_matrix, density_matrix)\n",
    "print('Hartree-Fock energy =',energy_ion + energy_scf)"
   ]
  },
  {
   "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": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "occupied orbital energies:\n",
      " [-0.6 -0.6 -0.6 -0.6 -0.6 -0.6]\n",
      "virtual orbital energies:\n",
      " [5.4 5.4]\n"
     ]
    }
   ],
   "source": [
    "def partition_orbitals(fock_matrix):\n",
    "    '''Returns a list with the occupied/virtual energies & orbitals defined by the input Fock matrix.'''\n",
    "    num_occ = (ionic_charge//2)*np.size(fock_matrix,0)//orbitals_per_atom\n",
    "    orbital_energy, orbital_matrix = np.linalg.eigh(fock_matrix)\n",
    "    occupied_energy = orbital_energy[:num_occ]\n",
    "    virtual_energy = orbital_energy[num_occ:]\n",
    "    occupied_matrix = orbital_matrix[:,:num_occ]\n",
    "    virtual_matrix = orbital_matrix[:,num_occ:]\n",
    "\n",
    "    return occupied_energy, virtual_energy, occupied_matrix, virtual_matrix\n",
    "\n",
    "occupied_energy, virtual_energy, occupied_matrix, virtual_matrix = partition_orbitals(fock_matrix)\n",
    "\n",
    "print('occupied orbital energies:\\n', occupied_energy)\n",
    "print('virtual orbital energies:\\n', virtual_energy)"
   ]
  },
  {
   "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": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "interaction tensor =\n",
      " [[[[-3.5e-02 -5.2e-16 -8.7e-17 -8.0e-16 -3.9e-16 -2.2e-14]\n",
      "   [ 3.8e-14 -1.2e-15 -1.5e-16 -1.9e-15 -7.2e-16 -3.5e-02]]\n",
      "\n",
      "  [[-5.2e-16 -2.4e-02  4.4e-18  2.1e-14  1.7e-14  1.1e-15]\n",
      "   [-3.0e-16  2.5e-14 -3.1e-16  1.8e-02  1.5e-02 -2.2e-15]]\n",
      "\n",
      "  [[-8.8e-17  2.6e-18 -2.4e-02  1.4e-14 -1.7e-14  1.4e-16]\n",
      "   [-1.2e-16  9.6e-16  2.0e-14  1.5e-02 -1.8e-02 -6.2e-16]]\n",
      "\n",
      "  [[-8.0e-16  2.1e-14  1.4e-14 -1.7e-03  8.5e-19  3.0e-16]\n",
      "   [ 1.9e-16  1.3e-03  1.1e-03 -2.3e-14 -3.1e-15 -8.9e-16]]\n",
      "\n",
      "  [[-3.9e-16  1.7e-14 -1.7e-14  7.3e-19 -1.7e-03  4.2e-17]\n",
      "   [-6.0e-18  1.1e-03 -1.3e-03 -1.9e-15 -2.1e-14 -4.4e-16]]\n",
      "\n",
      "  [[-2.2e-14  1.1e-15  1.3e-16  3.0e-16  4.2e-17  9.2e-03]\n",
      "   [ 9.2e-03  4.0e-16  4.8e-17 -4.5e-16 -3.4e-16 -3.8e-14]]]\n",
      "\n",
      "\n",
      " [[[ 3.8e-14 -3.0e-16 -1.2e-16  1.9e-16 -5.7e-18  9.2e-03]\n",
      "   [ 9.2e-03  5.2e-16  8.7e-17  8.0e-16  3.9e-16  2.2e-14]]\n",
      "\n",
      "  [[-1.2e-15  2.5e-14  9.6e-16  1.3e-03  1.1e-03  4.0e-16]\n",
      "   [ 5.2e-16 -1.7e-03 -9.5e-19 -2.1e-14 -1.7e-14 -1.1e-15]]\n",
      "\n",
      "  [[-1.5e-16 -3.0e-16  2.0e-14  1.1e-03 -1.3e-03  4.8e-17]\n",
      "   [ 8.6e-17  2.6e-20 -1.7e-03 -1.4e-14  1.7e-14 -1.3e-16]]\n",
      "\n",
      "  [[-1.9e-15  1.8e-02  1.5e-02 -2.3e-14 -1.9e-15 -4.5e-16]\n",
      "   [ 8.0e-16 -2.1e-14 -1.4e-14 -2.4e-02  0.0e+00 -3.0e-16]]\n",
      "\n",
      "  [[-7.2e-16  1.5e-02 -1.8e-02 -3.1e-15 -2.1e-14 -3.4e-16]\n",
      "   [ 3.9e-16 -1.7e-14  1.7e-14  3.5e-18 -2.4e-02 -5.1e-17]]\n",
      "\n",
      "  [[-3.5e-02 -2.2e-15 -6.2e-16 -8.9e-16 -4.4e-16 -3.8e-14]\n",
      "   [ 2.2e-14 -1.1e-15 -1.3e-16 -3.0e-16 -5.1e-17 -3.5e-02]]]]\n"
     ]
    }
   ],
   "source": [
    "def transform_interaction_tensor(occupied_matrix, virtual_matrix, interaction_matrix, chi_tensor):\n",
    "    '''Returns a transformed V tensor defined by the input occupied, virtual, & interaction matrices.'''\n",
    "    chi2_tensor = np.einsum('qa,ri,qrp', virtual_matrix, occupied_matrix, chi_tensor, optimize=True)\n",
    "    interaction_tensor = np.einsum('aip,pq,bjq->aibj', chi2_tensor, interaction_matrix, chi2_tensor, optimize=True)\n",
    "    return interaction_tensor\n",
    "\n",
    "interaction_tensor = transform_interaction_tensor(occupied_matrix, virtual_matrix, interaction_matrix, chi_tensor)\n",
    "print('interaction tensor =\\n', interaction_tensor)"
   ]
  },
  {
   "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": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.001278681955212097\n"
     ]
    }
   ],
   "source": [
    "def calculate_energy_mp2(fock_matrix, interaction_matrix, chi_tensor):\n",
    "    '''Returns the MP2 contribution to the total energy defined by the input Fock & interaction matrices.'''\n",
    "    E_occ, E_virt, occupied_matrix, virtual_matrix = partition_orbitals(fock_matrix)\n",
    "    V_tilde = transform_interaction_tensor(occupied_matrix, virtual_matrix, interaction_matrix, chi_tensor)\n",
    "\n",
    "    energy_mp2 = 0.0\n",
    "    num_occ = len(E_occ)\n",
    "    num_virt = len(E_virt)\n",
    "    for a in range(num_virt):\n",
    "        for b in range(num_virt):\n",
    "            for i in range(num_occ):\n",
    "                for j in range(num_occ):\n",
    "                    energy_mp2 -= ( (2.0*V_tilde[a,i,b,j]**2 - V_tilde[a,i,b,j]*V_tilde[a,j,b,i])\n",
    "                                    /(E_virt[a] + E_virt[b] - E_occ[i] - E_occ[j]) )\n",
    "    return energy_mp2\n",
    "\n",
    "energy_mp2 = calculate_energy_mp2(fock_matrix, interaction_matrix, chi_tensor)\n",
    "print(energy_mp2)"
   ]
  },
  {
   "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": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "difference between fock matrix implementations = 2.413952749762047e-15\n"
     ]
    }
   ],
   "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": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "base matrix construction time = 0.00029010800062678754\n",
      "Fock matrix construction time (\"slow\") = 0.0009618200001568766\n",
      "Fock matrix construction time (\"fast\") = 0.016594663999057957\n",
      "19\n",
      "base matrix construction time = 0.18844490700030292\n",
      "Fock matrix construction time (\"slow\") = 0.005319128000337514\n",
      "Fock matrix construction time (\"fast\") = 4.029375972000707\n",
      "55\n",
      "base matrix construction time = 1.3979908880010044\n",
      "Fock matrix construction time (\"slow\") = 0.20557714399910765\n",
      "Fock matrix construction time (\"fast\") = 33.36372286400001\n"
     ]
    }
   ],
   "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": 24,
   "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": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([-0.590921, -0.590921, -0.590921]), -8.950824066074436, -0.0001599824535679023)\n",
      "(array([-0.599361, -0.59236 , -0.59236 , -0.589434, -0.589434, -0.582434]), -17.901094929513107, -0.0013473598084440934)\n",
      "(array([-0.590921, -0.590921, -0.590921, -0.590921, -0.590921, -0.590921]), -17.901648132148868, -0.00032183740115952005)\n"
     ]
    }
   ],
   "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": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rms error = 0.0002863242502397328\n"
     ]
    }
   ],
   "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": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0002858291458085528\n",
      "0.00028586123483649345\n",
      "0.0697827550400981\n",
      "0.00030252506897538535\n",
      "0.09266776547794396\n",
      "0.000285839299206041\n",
      "0.00162673950344277\n",
      "0.0005039444731676446\n",
      "0.00039321739870742673\n",
      "0.000293413434148544\n",
      "0.0002858325294338677\n",
      "0.0002858299256686513\n",
      "0.0002879805982460839\n",
      "0.10430081594517075\n",
      "0.04342655342749679\n",
      "0.07701863364676502\n",
      "0.03308347908650576\n",
      "0.04894288751405319\n",
      "0.020336187501330087\n",
      "0.036470958773284115\n",
      "0.015697468239554632\n",
      "0.022954426558591038\n",
      "0.00951776159423058\n",
      "0.01728118065478587\n",
      "0.007458276234716499\n",
      "0.010762068311981432\n",
      "0.004457001321395454\n",
      "0.008197206239865541\n",
      "0.003555966632285127\n",
      "0.005048931521369579\n",
      "0.0021002097486725194\n",
      "0.0038991862823617245\n",
      "0.001716339684866278\n",
      "0.00238119828960611\n",
      "0.001021116420578156\n",
      "0.0018736199116184168\n",
      "0.0008677566209157689\n",
      "0.0013808615529079776\n",
      "0.0007320266247238685\n",
      "0.0011549595758308909\n",
      "0.0005389077494111427\n",
      "0.0009496348752075263\n",
      "0.0004834126989815779\n",
      "0.0008272193841429489\n",
      "0.0004347885940322048\n",
      "0.0005939924780815616\n",
      "0.0003579608424494638\n",
      "0.0003967144370623498\n",
      "0.0005153175311144255\n",
      "0.00034080748001077824\n",
      "0.0004915464524267442\n",
      "0.00032274680430153835\n",
      "0.0004928258553203911\n",
      "0.00032125170791825683\n",
      "0.00039992352699921954\n",
      "0.00031498316822861064\n",
      "0.00037734546431012915\n",
      "0.00030294487969242934\n",
      "0.0003509038307268281\n",
      "0.0002993215148894099\n",
      "0.00034037524702006756\n",
      "0.0002937758264125941\n",
      "0.0003383829814423393\n",
      "0.00029570786060610356\n",
      "0.00031651600816884763\n",
      "0.0002933101456789131\n",
      "0.0003077007175488972\n",
      "0.0002897661312866954\n",
      "0.00029777249129600634\n",
      "0.00030156993942649396\n",
      "0.00028908950028411326\n",
      "0.00030227927696454014\n",
      "0.0002890470979077263\n",
      "0.00029913254016923163\n",
      "0.0002882250978582377\n",
      "0.00029889206485682665\n",
      "0.0002873508578196126\n",
      "0.00029256123699877146\n",
      "0.00029320382092501375\n",
      "0.0002876765870159818\n",
      "0.0002934911712134024\n",
      "0.0002874131031683572\n",
      "0.0002920427946760553\n",
      "0.00028654977953983365\n",
      "0.00029002093763197737\n",
      "0.00028655099948069997\n",
      "0.00028920232924446\n",
      "0.00028667698987623104\n",
      "0.0002899418743143548\n",
      "0.00028628013660528994\n",
      "0.00028839884485708803\n",
      "0.00028630137061278193\n",
      "0.0002871854804754191\n",
      "0.00028751929534374303\n",
      "0.0002861858633141591\n",
      "0.00028801166419803795\n",
      "0.00028608918579537\n",
      "0.000287386463374112\n",
      "0.0002861200954665288\n",
      "0.00028643848348880165\n",
      "0.0002868095129607359\n",
      "0.00028594233420875526\n",
      "0.0002866690121873296\n",
      "0.0002859399810239029\n",
      "0.0002866704341024858\n",
      "0.0002858950927734411\n",
      "0.00028641765782592185\n",
      "0.00028588088587811383\n",
      "0.00028641624009784043\n",
      "0.00028589036021280215\n",
      "0.0002861165829463739\n",
      "0.00028617012394359305\n",
      "0.000285830107868257\n",
      "0.00028617760626736995\n",
      "0.00028583805307631993\n",
      "0.0002861391218633577\n",
      "0.0002858456872418108\n",
      "0.00028596625258446075\n",
      "0.000285786483294022\n",
      "0.00028595940387839353\n",
      "0.0002857958954353865\n",
      "0.000286043057608991\n",
      "0.0002857734028365267\n",
      "0.00028592991984393993\n",
      "0.00028579905140658785\n",
      "0.00028586216836608497\n",
      "0.0002857768022286488\n",
      "0.00028566143850532194\n",
      "0.00028558679805493924\n",
      "0.00028582670713986813\n",
      "0.0002864136096183406\n",
      "0.00028568340431056464\n",
      "0.0002858200831791772\n",
      "0.0002857795925312008\n",
      "0.00028573641662857325\n",
      "0.00028563899960159557\n",
      "0.00028561318022118475\n",
      "0.0002858542602085255\n",
      "0.0002857299443816888\n",
      "0.0002858264757201742\n",
      "0.0002857281796660446\n",
      "0.00028561583896659945\n",
      "0.00028566621400895776\n",
      "0.0002855989428899332\n",
      "0.00028562994768037203\n",
      "0.00028554037357436744\n",
      "0.0002854734429327066\n",
      "0.00028562585733218725\n",
      "0.0002856255513736557\n",
      "0.0002854817335221198\n",
      "0.0002854396666687342\n",
      "0.0002853442047175593\n",
      "0.0002854338576433066\n",
      "0.00028548620900513427\n",
      "0.0002853880711011487\n",
      "0.00028551367518553323\n",
      "0.00028544829364040593\n",
      "0.0002853937164491849\n",
      "0.0002853078135892701\n",
      "0.0002852124829710207\n",
      "0.0002852367813746575\n",
      "0.00028524213450491556\n",
      "0.00028527710175050853\n",
      "0.0002852645678353287\n",
      "0.0002851852557168241\n",
      "0.00028525975659840925\n",
      "0.0002851209787486873\n",
      "0.0002849782608466969\n",
      "0.00028509569700246163\n",
      "0.0002851396164783686\n",
      "0.0002850525484185789\n",
      "0.00028506741051266787\n",
      "0.0002848784423106755\n",
      "0.00028465987780527366\n",
      "0.0002848526245085766\n",
      "0.0002849683336876866\n",
      "0.0002848566656700951\n",
      "0.0002848073873687285\n",
      "0.000284756579891335\n",
      "0.0002847247871092842\n",
      "0.00028477541410225766\n",
      "0.0002846804547692784\n",
      "0.000284631342901324\n",
      "0.0002843837591829698\n",
      "0.00028451240775928847\n",
      "0.00028448567562342896\n",
      "0.00028439378225140105\n",
      "0.0002843748842412853\n",
      "0.00028427809766751855\n",
      "0.00028424532940120155\n",
      "0.00028397643387683867\n",
      "0.00028419182790256377\n",
      "0.000284125437355545\n",
      "0.0002839797173377755\n",
      "0.00028396577251686857\n",
      "0.00028384912473272044\n",
      "0.0002838854469268473\n",
      "0.0002837797455403039\n",
      "0.00028356263837462384\n",
      "0.00028369454317445184\n",
      "0.0002838305132093163\n",
      "0.0002837789936762164\n",
      "0.0002836176864206911\n",
      "0.0002833677749336993\n",
      "0.00028319902245887193\n",
      "0.00028336242347666194\n",
      "0.0002833172505539001\n",
      "0.0002830505649807657\n",
      "0.00028263253047077495\n",
      "0.0002832962522422102\n",
      "0.00028352542108589344\n",
      "0.00028298653017408507\n",
      "0.0002828941867459031\n",
      "0.0002827985336566875\n",
      "0.00028275447870096035\n",
      "0.00028256266804456767\n",
      "0.0002823434094327269\n",
      "0.00028260095878040895\n",
      "0.0002828335378755381\n",
      "0.00028240068933039033\n",
      "0.0002827020182065319\n",
      "0.00028249356960814785\n",
      "0.00028212549216251\n",
      "0.0002819876309414745\n",
      "0.00028195374326609964\n",
      "0.0002817108117005136\n",
      "0.0002824538929037682\n",
      "0.0002824892362621826\n",
      "0.0002818306606092591\n",
      "0.0002821991395678424\n",
      "0.0002821985748621463\n",
      "0.0002817619578669644\n",
      "0.0002830486893437605\n",
      "0.00028217016486252\n",
      "0.00028187039782700196\n",
      "0.00028176241442617054\n",
      "0.00028164038625702825\n",
      "0.0002818522133356518\n",
      "0.000281410623758219\n",
      "0.00028135321077492147\n",
      "0.0002821695520514379\n",
      "0.00028136037203993\n",
      "0.0002813382966454836\n",
      "0.00028164275813830905\n",
      "0.0002813598947901092\n",
      "0.00028143769467677123\n",
      "0.00028180554665252646\n",
      "0.0002812159790028522\n",
      "0.0002814594265300247\n",
      "0.00028187890671578844\n",
      "0.00028146388263426887\n",
      "0.0002811603647301007\n",
      "0.0002813728015394235\n",
      "0.00028185682102837043\n",
      "0.0002814315837048917\n",
      "0.00028098850970411895\n",
      "0.0002809771931238924\n",
      "0.00028088123492287187\n",
      "0.00028089875770050674\n",
      "0.00028099889759408984\n",
      "0.0002808453360777781\n",
      "0.0002810472378486727\n",
      "0.000281105761748741\n",
      "0.00028181528054281634\n",
      "0.0002810193675151339\n",
      "0.00028061156493286595\n",
      "0.00028045464716536326\n",
      "0.00028164488875011913\n",
      "0.0002809787056900604\n",
      "0.0002806538110443158\n",
      "0.0002807611816307719\n",
      "0.0002803489479039866\n",
      "0.0002802846193597895\n",
      "0.0002804186017001776\n",
      "0.00028028535933308766\n",
      "0.00028048384717284364\n",
      "0.00028012226109351595\n",
      "0.00027995996188682223\n",
      "0.0002803219268314085\n",
      "0.0002799641641169701\n",
      "0.00028164524982698605\n",
      "0.00028041662267991937\n",
      "0.0002803135591534267\n",
      "0.0002801863246353753\n",
      "0.0002797704966412631\n",
      "0.00027971482871419804\n",
      "0.00027982170825198194\n",
      "0.0002802977233894074\n",
      "0.00027953195264145516\n",
      "0.00027932848008738924\n",
      "0.00027950439433735793\n",
      "0.0002805288090718344\n",
      "0.00027997345617819676\n",
      "0.0002804516117958131\n",
      "0.0002799032955056457\n",
      "0.0002799184413618106\n",
      "0.00027953928359767614\n",
      "0.000279706210231764\n",
      "0.0002793181646553862\n",
      "0.0002793962525902985\n",
      "0.00027960850243150865\n",
      "0.00027936217144380464\n",
      "0.00027909019102907685\n",
      "0.00027892544653636573\n",
      "0.00027916139165079637\n",
      "0.0002794591036264192\n",
      "0.0002792204637682669\n",
      "0.00027884821174751383\n",
      "0.0002788450104394647\n",
      "0.00027866181037043355\n",
      "0.0002785890170994098\n",
      "0.00027887464396145515\n",
      "0.0002792265706737879\n",
      "0.0002798246468976355\n",
      "0.0002790064304158748\n",
      "0.00027854958711293873\n",
      "0.0002786974537220741\n",
      "0.00027906252265572193\n",
      "0.000278379062699901\n",
      "0.0002783020051087558\n",
      "0.0002781054499372958\n",
      "0.00027783996537215065\n",
      "0.0002782177651329795\n",
      "0.00027881833587532816\n",
      "0.00027767213421952695\n",
      "0.0002772101586850758\n",
      "0.00027822024170237944\n",
      "0.0002776206117964392\n",
      "0.0002780063321752732\n",
      "0.000277205554586142\n",
      "0.00027693173845328106\n",
      "0.0002774221744037663\n",
      "0.00027745237698764225\n",
      "0.00027828152357299736\n",
      "0.0002773796281550779\n",
      "0.0002768336971382956\n",
      "0.0002770017768699259\n",
      "0.00027721589490541223\n",
      "0.0002781421107866793\n",
      "0.00027753272823036\n",
      "0.00027647928158857146\n",
      "0.0002765159392765535\n",
      "0.00027850728946415316\n",
      "0.0002772378826294957\n",
      "0.00027649797175420185\n",
      "0.0002763650235874526\n",
      "0.0002770941327868413\n",
      "0.0002765930575086195\n",
      "0.000277339766017624\n",
      "0.0002763704361232321\n",
      "0.00027680800795460027\n",
      "0.00027635002338623156\n",
      "0.00027722493934686186\n",
      "0.0002766471010037171\n",
      "0.0002764063381854946\n",
      "0.00027687063696097894\n",
      "0.00027694317863612324\n",
      "0.000276261348632293\n",
      "0.00027610041334508054\n",
      "0.00027661671182422764\n",
      "0.000277188658994471\n",
      "0.0002763035255280698\n",
      "0.0002757966834722303\n",
      "0.000275918303090162\n",
      "0.00027698391296624765\n",
      "0.00027618987864363835\n",
      "0.0002766541355768464\n",
      "0.00027627269626536653\n",
      "0.0002762813985879997\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0002764367593272621\n",
      "0.0002761311622530208\n",
      "0.0002765202789488732\n",
      "0.0002760926958607468\n",
      "0.00027627228685516634\n",
      "0.0002764896941049556\n",
      "0.0002760675188114396\n",
      "0.00027586198736584243\n",
      "0.0002760457729125543\n",
      "0.0002757693528333413\n",
      "0.0002759529958911693\n",
      "0.0002762208853482661\n",
      "0.00027617516309135827\n",
      "0.0002763164074508023\n",
      "0.00027599952758704846\n",
      "0.00027582397611663936\n",
      "0.00027630430553504724\n",
      "0.00027591709605229457\n",
      "0.00027567874610882404\n",
      "0.00027585225125679986\n",
      "0.00027593357411325744\n",
      "0.0002761096129411254\n",
      "0.0002758776033714542\n",
      "0.00027608025618146043\n",
      "0.00027571211300278677\n",
      "0.0002760367967284226\n",
      "0.0002757697533221263\n",
      "0.0002759925984001933\n",
      "0.0002758133943325657\n",
      "0.0002757547559079017\n",
      "0.00027593341987138377\n",
      "0.00027580124011931224\n",
      "0.00027591437004920757\n",
      "0.0002757445672952358\n",
      "0.0002758245721681214\n",
      "0.0002757099421119947\n",
      "0.00027571772896452284\n",
      "0.0002758005819560303\n",
      "0.0002760114663576806\n",
      "0.0002756603469198102\n",
      "0.0002758912307742935\n",
      "0.0002756618081668262\n",
      "0.0002760554350168171\n",
      "0.000275661679410807\n",
      "0.00027578491412425817\n",
      "0.0002754869139433164\n",
      "0.00027543360640063034\n",
      "0.00027575814641436086\n",
      "0.00027571861983484214\n",
      "0.00027564885520109215\n",
      "0.0002757498603718386\n",
      "0.00027588296251039113\n",
      "0.00027559973180570025\n",
      "0.00027573201036258925\n",
      "0.0002756009679479321\n",
      "0.000275801118636286\n",
      "0.00027556625646983467\n",
      "0.0002757168078876105\n",
      "0.00027576339260797136\n",
      "0.000275622186463327\n",
      "0.00027571414913738575\n",
      "0.0002756048323454528\n",
      "0.0002758110513548541\n",
      "0.0002755394609757711\n",
      "0.0002755381132236313\n",
      "0.0002755161836182721\n",
      "0.0002754816193166913\n",
      "0.0002754608155014672\n",
      "0.0002754749803576593\n",
      "0.0002756005576881733\n",
      "0.00027545948781495155\n",
      "0.0002754578852598128\n",
      "0.00027563941742362154\n",
      "0.0002754753858230007\n",
      "0.00027561124513837893\n",
      "0.00027549984023824114\n",
      "0.00027541577539575573\n",
      "0.00027540783821244794\n",
      "0.0002753849019078124\n",
      "0.0002754682765808759\n",
      "0.00027536929390559216\n",
      "0.00027549464759392384\n",
      "0.0002754937999424714\n",
      "0.00027542972514078633\n",
      "0.00027537851068817464\n",
      "0.00027541534026759285\n",
      "0.00027535979162023397\n",
      "0.00027546524105973084\n",
      "0.00027538851001416355\n",
      "0.0002753663212551491\n",
      "0.0002754106898193036\n",
      "0.0002754111053850104\n",
      "0.00027541964920178093\n",
      "0.0002754005958021008\n",
      "0.00027558376071390656\n",
      "0.0002753407946550086\n",
      "0.0002754352986471054\n",
      "0.00027535353002031945\n",
      "0.00027545528382053486\n",
      "0.0002753445570831107\n",
      "0.0002754195165911517\n",
      "0.00027535104486198935\n",
      "0.0002753758134171907\n",
      "0.00027533539786163495\n",
      "0.0002754827977160674\n",
      "0.00027544667344933974\n",
      "0.0002753344420145131\n",
      "0.0002753753007817874\n",
      "0.00027537779919010455\n",
      "0.0002753902555697155\n",
      "0.0002753316812802285\n",
      "0.00027538946603763986\n",
      "0.00027533690857845287\n",
      "0.00027541832161174986\n",
      "0.0002753069537646642\n",
      "0.00027540221463155743\n",
      "0.0002753355725678959\n",
      "0.0002752878446166935\n",
      "0.0002752889475652667\n",
      "0.0002754553300095675\n",
      "0.0002752969924304467\n",
      "0.0002753125344605801\n",
      "0.00027532130031522725\n",
      "0.0002753043070380681\n",
      "0.0002752964871940391\n",
      "0.0002752801506087207\n",
      "0.00027534342674741887\n",
      "0.0002753305527109656\n",
      "0.0002752669281447952\n",
      "0.0002753288757065554\n",
      "0.00027532728764718354\n",
      "0.00027527590923695863\n",
      "0.00027526046152219713\n",
      "0.00027527887406261434\n",
      "0.0002752949281278697\n",
      "0.00027531892415120693\n",
      "0.0002752620271804035\n",
      "0.00027533312811414234\n",
      "0.0002752814569963356\n",
      "0.00027533767876307066\n",
      "0.00027528506019337736\n",
      "0.0002752537410535093\n",
      "0.0002752606948220646\n",
      "0.00027526514095315986\n",
      "0.00027524485961455224\n",
      "0.00027524652004060166\n",
      "0.0002752710296995446\n",
      "0.0002752917965873116\n",
      "0.00027524267951990956\n",
      "0.00027526469816196997\n",
      "0.0002752809192895379\n",
      "0.0002752465130235003\n",
      "0.00027528112248941314\n",
      "0.00027523412194202724\n",
      "0.0002752147197689245\n",
      "0.00027519287117640635\n",
      "0.0002752490764233626\n",
      "0.0002752445500505373\n",
      "0.00027521232154102226\n",
      "0.0002752196639173446\n",
      "0.000275231613150451\n",
      "0.0002752295971998308\n",
      "0.00027521180143566363\n",
      "0.00027519108549952195\n",
      "0.00027517977222135177\n",
      "0.0002752158219806103\n",
      "0.0002752237727910083\n",
      "0.00027517713913342775\n",
      "0.0002752220135727752\n",
      "0.00027520092219732314\n",
      "0.0002751537289105182\n",
      "0.0002751198062892218\n",
      "0.0002751443415821167\n",
      "0.00027519909562752885\n",
      "0.00027514913127222966\n",
      "0.00027517386258268066\n",
      "0.0002751264448009446\n",
      "0.000275133579289788\n",
      "0.0002751258910895374\n",
      "0.00027511437516427\n",
      "0.00027512125753639284\n",
      "0.0002751291485880087\n",
      "0.00027512277666850044\n",
      "0.00027507203116509455\n",
      "0.00027503462394024285\n",
      "0.0002750898381034163\n",
      "0.00027504400401720227\n",
      "0.000275070348997946\n",
      "0.00027506963450028247\n",
      "0.00027503392945914083\n",
      "0.00027500618381689904\n",
      "0.00027502998357027306\n",
      "0.0002750372104763293\n",
      "0.00027499585669450266\n",
      "0.000274959274384131\n",
      "0.00027497446109701635\n",
      "0.0002749933114276755\n",
      "0.0002749355153710114\n",
      "0.00027488713346264114\n",
      "0.00027491013649121385\n",
      "0.00027493860323690744\n",
      "0.0002749481877434761\n",
      "0.0002748836361093608\n",
      "0.0002748552500268512\n",
      "0.0002748668281526177\n",
      "0.0002748733543919336\n",
      "0.0002748269416980609\n",
      "0.0002747938848968412\n",
      "0.0002748202820287847\n",
      "0.0002747691593094089\n",
      "0.00027471482991737377\n",
      "0.0002748038255789365\n",
      "0.00027478680363347276\n",
      "0.0002747793943024864\n",
      "0.00027470659614930846\n",
      "0.0002746883216379516\n",
      "0.0002746551216314102\n",
      "0.0002746007248891402\n",
      "0.00027473004372285085\n",
      "0.00027466106790792613\n",
      "0.0002746283366416904\n",
      "0.000274604495335296\n",
      "0.0002746530390023204\n",
      "0.00027459217367501576\n",
      "0.0002746376411256079\n",
      "0.0002746586093556107\n",
      "0.00027466031742142354\n",
      "0.0002744758926121756\n",
      "0.00027441296425699036\n",
      "0.0002746232107790045\n",
      "0.000274456914402467\n",
      "0.0002746578201596963\n",
      "0.00027471112713041976\n",
      "0.0002745817837166251\n",
      "0.0002745033455940414\n",
      "0.00027456572136751823\n",
      "0.0002744782641896085\n",
      "0.0002746460673302521\n",
      "0.0002744559250076192\n",
      "0.0002746338191904799\n",
      "0.00027453793468568046\n",
      "0.000274471014081152\n",
      "0.0002744537432204629\n",
      "0.00027438534298692037\n",
      "0.00027438599668386677\n",
      "0.0002745344391818323\n",
      "0.0002743263229446269\n",
      "0.0002742358460357647\n",
      "0.00027436106721694705\n",
      "0.00027443048040700686\n",
      "0.000274287838013543\n",
      "0.00027436104762201385\n",
      "0.0002742661193018242\n",
      "0.000274291534162379\n",
      "0.00027425600081176305\n",
      "0.0002742378082432932\n",
      "0.0002741970465764853\n",
      "0.0002742683465549834\n",
      "0.00027453438996660864\n",
      "0.0002743151594676085\n",
      "0.0002742697041966382\n",
      "0.00027411345252235865\n",
      "0.00027406288488218916\n",
      "0.000274052405983435\n",
      "0.0002740792909985656\n",
      "0.00027413472721381904\n",
      "0.00027431034444306727\n",
      "0.00027423801855018797\n",
      "0.00027430202001557885\n",
      "0.0002742139169893496\n",
      "0.000274126096708596\n",
      "0.00027404387422016703\n",
      "0.00027416480586755533\n",
      "0.0002741620577067187\n",
      "0.00027400006241498376\n",
      "0.0002739474936271275\n",
      "0.00027415312464267337\n",
      "0.00027417413338297314\n",
      "0.0002740386527511561\n",
      "0.00027394453408587744\n",
      "0.0002739474881690475\n",
      "0.0002740826999229829\n",
      "0.000273933383952472\n",
      "0.0002741382803991945\n",
      "0.00027408391805421625\n",
      "0.0002740741830265861\n",
      "0.00027401693276092273\n",
      "0.0002740048163190275\n",
      "0.00027394067740280266\n",
      "0.00027401383613787696\n",
      "0.0002740245246584413\n",
      "0.0002740346157863364\n",
      "0.00027397948434151355\n",
      "0.0002738506387255276\n",
      "0.00027385174552125956\n",
      "0.000273953413250107\n",
      "0.0002738528575385641\n",
      "0.00027393957220441315\n",
      "0.000273996822010122\n",
      "0.0002738855444147146\n",
      "0.00027389279541764313\n",
      "0.0002738364901936765\n",
      "0.00027395836148249046\n",
      "0.00027387989229084666\n",
      "0.00027389033856036827\n",
      "0.0002739451850114679\n",
      "0.00027398903158317486\n",
      "0.00027397429974023303\n",
      "0.0002738276435551606\n",
      "0.00027383075980818796\n",
      "0.00027382103124709144\n",
      "0.00027389786294672254\n",
      "0.0002738765808358053\n",
      "0.0002738988042585295\n",
      "0.00027385912995631583\n",
      "0.0002737854431487661\n",
      "0.00027389688594970155\n",
      "0.0002738852093319765\n",
      "0.0002739280635146623\n",
      "0.00027385110289204745\n",
      "0.00027381335841139174\n",
      "0.00027387520339873885\n",
      "0.00027389544311352934\n",
      "0.00027393257806939244\n",
      "0.00027378458814003755\n",
      "0.00027390373515346814\n",
      "0.00027385118289135417\n",
      "0.00027375014681333353\n",
      "0.0002738396339022618\n",
      "0.00027392268942598897\n",
      "0.0002737713454538678\n",
      "0.00027390345351074305\n",
      "0.0002738370551510506\n",
      "0.0002737712032556759\n",
      "0.00027389202168362384\n",
      "0.00027382314681321054\n",
      "0.00027380190706465733\n",
      "0.0002738924277429411\n",
      "0.0002738139462289195\n",
      "0.0002737841191643569\n",
      "0.0002737920139006638\n",
      "0.000273754541640295\n",
      "0.00027387912408580607\n",
      "0.0002737998029567187\n",
      "0.00027374526344280266\n",
      "0.0002738683324023179\n",
      "0.00027387601576138053\n",
      "0.0002737913301088451\n",
      "0.0002737530913425412\n",
      "0.00027392713752797814\n",
      "0.00027377743863136214\n",
      "0.0002738742089360398\n",
      "0.0002737816956773786\n",
      "0.0002737715312135247\n",
      "0.00027387489734278814\n",
      "0.000273776509903767\n",
      "0.00027374666261731866\n",
      "0.0002737435061997686\n",
      "0.00027373856254430653\n",
      "0.0002737413357940725\n",
      "0.0002737197598757048\n",
      "0.0002736906028613023\n",
      "0.0002737242943492083\n",
      "0.00027370312099630215\n",
      "0.00027372433716013367\n",
      "0.0002736994768226974\n",
      "0.0002737058488639872\n",
      "0.0002736871971988369\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0002736597659063175\n",
      "0.00027369377659882684\n",
      "0.0002736674568152972\n",
      "0.00027367361031271713\n",
      "0.0002736665618408264\n",
      "0.00027366326263702503\n",
      "0.0002736443409324293\n",
      "0.0002737680081466554\n",
      "0.000273648027317669\n",
      "0.0002736278008865485\n",
      "0.00027360295920810895\n",
      "0.00027365280093680363\n",
      "0.0002736118751968655\n",
      "0.0002735999322672958\n",
      "0.00027356158950108615\n",
      "0.00027359763470207973\n",
      "0.0002735782547297761\n",
      "0.00027357048284159103\n",
      "0.0002735589722415451\n",
      "0.00027351307485713934\n",
      "0.00027355066581222996\n",
      "0.0002735500053686029\n",
      "0.0002735141813412994\n",
      "0.00027351972554852486\n",
      "0.0002735045721975772\n",
      "0.00027349728337483254\n",
      "0.0002734919787182776\n",
      "0.00027349341128552233\n",
      "0.0002734575862505252\n",
      "0.00027339902339219564\n",
      "0.00027345189193936646\n",
      "0.0002734429977444443\n",
      "0.0002734129064901657\n",
      "0.00027342897576991033\n",
      "0.00027339038037312474\n",
      "0.0002733449780065026\n",
      "0.000273357991381598\n",
      "0.0002733930693036335\n",
      "0.0002733632070028867\n",
      "0.00027339710662113747\n",
      "0.0002733736144377364\n",
      "0.0002733196543755091\n",
      "0.0002733073478284737\n",
      "0.0002733041582747602\n",
      "0.0002733421667168457\n",
      "0.0002733428232976688\n",
      "0.0002733014296994789\n",
      "0.00027332880893903895\n",
      "0.00027325977709055595\n",
      "0.0002732251377347802\n",
      "0.0002733127780068812\n",
      "0.0002732884416189929\n",
      "0.0002732226417424761\n",
      "0.000273189434991334\n",
      "0.00027319068037124153\n",
      "0.00027320317306520086\n",
      "0.0002732118929303368\n",
      "0.00027320749802783764\n",
      "0.00027318202507679173\n",
      "0.00027324308675835267\n",
      "0.00027311884670074894\n",
      "0.0002730843622966225\n",
      "0.000273129692319282\n",
      "0.00027328428819153953\n",
      "0.00027318640245857365\n",
      "0.000273073976967917\n",
      "0.00027306488440780063\n",
      "0.0002730963318475081\n",
      "0.0002732438296384318\n",
      "0.0002731482790935002\n",
      "0.0002731031162686815\n",
      "0.00027307753892766393\n",
      "0.00027317716354884015\n",
      "0.0002730559114280745\n",
      "0.00027310116249769053\n",
      "0.00027319441850382235\n",
      "0.0002731110329889883\n",
      "0.00027310279850782367\n",
      "0.00027314300427296815\n",
      "0.00027313861587143217\n",
      "0.00027307046334122665\n",
      "0.0002730298002408631\n",
      "0.0002730779998916626\n",
      "0.00027312342376483334\n",
      "0.00027322825042825445\n",
      "0.0002730652649181718\n",
      "0.00027310501518536495\n",
      "0.00027311808601465377\n",
      "0.00027306655803171424\n",
      "0.00027313036824937805\n",
      "0.0002730631471720438\n",
      "0.0002731291650634912\n",
      "0.0002730600754310626\n",
      "0.0002730478764846784\n",
      "0.0002730401020024066\n",
      "0.00027306921734272864\n",
      "0.00027304202655221933\n",
      "0.0002730636641974994\n",
      "0.0002730663154841165\n",
      "0.00027307632943612413\n",
      "0.0002730416393753252\n",
      "0.000273037419523581\n",
      "0.00027306837946782636\n",
      "0.0002730400374754067\n",
      "0.00027310641499666015\n",
      "0.00027303695301872583\n",
      "0.00027301260647186903\n",
      "0.0002730114370282191\n",
      "0.0002730657359350143\n",
      "0.00027303444471673\n",
      "0.00027305140779490436\n",
      "0.00027299885812163055\n",
      "0.0002729884727684915\n",
      "0.00027299795117853944\n",
      "0.0002730372054493138\n",
      "0.00027301764465929303\n",
      "0.0002730295837770682\n",
      "0.00027299923318470294\n",
      "0.000273022159229773\n",
      "0.00027297898599457147\n",
      "0.000272969743789863\n",
      "0.0002730388076024299\n",
      "0.00027300714872270903\n",
      "0.0002729991690011087\n",
      "0.00027296217006424217\n",
      "0.00027294754419767873\n",
      "0.00027297182690548945\n",
      "0.00027295985778878035\n",
      "0.00027296509897421127\n",
      "0.00027296937439246627\n",
      "0.00027295515390842845\n",
      "0.0002729519324077783\n",
      "0.00027294595071057415\n",
      "0.00027297835240016027\n",
      "0.0002729563319612966\n",
      "0.0002729714172749221\n",
      "0.00027293855057488193\n",
      "0.00027311919320965997\n",
      "0.00027304060863974033\n",
      "0.00027296287978787604\n",
      "0.00027294534310478403\n",
      "0.00027295998897737637\n",
      "0.00027295443203271326\n",
      "0.0002729806399478733\n",
      "0.00027294451053971815\n",
      "0.00027297812274428083\n",
      "0.0002729428799672827\n",
      "0.00027291400933152576\n",
      "0.00027289914647203535\n",
      "0.00027295198070055864\n",
      "0.0002729630576284309\n",
      "0.00027293820162306007\n",
      "0.0002729340659300143\n",
      "0.0002729280481160102\n",
      "0.00027294912782838035\n",
      "0.00027293910057263907\n",
      "0.00027307109978231116\n",
      "0.00027293406558356783\n",
      "0.00027295610267398843\n",
      "0.0002729303868373887\n",
      "0.00027292306286242404\n",
      "0.0002729235916782253\n",
      "0.0002729322133674125\n",
      "0.0002729015405907895\n",
      "0.00027291934671934154\n",
      "0.00027293393988663624\n",
      "0.0002729193014192276\n",
      "0.0002728976758902839\n",
      "0.0002730378537450381\n",
      "0.0002728966957752284\n",
      "0.0002730393733793345\n",
      "0.00027303417728580844\n",
      "0.0002729181914529697\n",
      "0.0002729259669932955\n",
      "0.000272928170367056\n",
      "0.0002729011403868288\n",
      "0.0002729264983906798\n",
      "0.0002728997839973091\n",
      "0.00027292448876040413\n",
      "0.0002729166095863069\n",
      "0.00027292112192445835\n",
      "0.00027289758046886473\n",
      "0.0002728976239677\n",
      "0.00027292119348120797\n",
      "0.0002729026066169003\n",
      "0.0002729088038385325\n",
      "0.00027302539201291186\n",
      "0.0002729045129993547\n",
      "0.00027301633805534744\n",
      "0.0002729045310914731\n",
      "0.0002729198515020484\n",
      "0.00027289874109187845\n",
      "0.0002730544706461774\n",
      "0.00027289667394751624\n",
      "0.000272882555899215\n",
      "0.00027301559786126954\n",
      "0.0002730225067528153\n",
      "0.0002728972596417358\n",
      "0.0002728971158343552\n",
      "0.0002728922076018365\n",
      "0.000272889014752175\n",
      "0.00027291680434462346\n",
      "0.0002728902023196756\n",
      "0.000272892183822427\n",
      "0.0002730277605220351\n",
      "0.0002728926179073431\n",
      "0.0002728923995105013\n",
      "0.0002728888335842523\n",
      "0.0002729062541352111\n",
      "0.0002728893900488225\n",
      "0.00027302080527881033\n",
      "0.00027289188616487286\n",
      "0.0002728954092023668\n",
      "0.0002728925738613582\n",
      "0.0002730295688965384\n",
      "0.0002728902356546382\n",
      "0.0002728957530809796\n",
      "0.0002728889863952348\n",
      "0.0002728830692497684\n",
      "0.00027289284313595087\n",
      "0.00027288763388357115\n",
      "0.00027289490035515527\n",
      "0.000272887308379747\n",
      "0.0002728944632896183\n",
      "0.0002728873163560075\n",
      "0.00027289409444191255\n",
      "0.00027288739052795555\n",
      "0.00027288045282515604\n",
      "0.0002730175609791798\n",
      "0.00027288222444031074\n",
      "0.0002728857206401741\n",
      "0.0002728807265237257\n",
      "0.00027288829409429627\n",
      "0.0002728803359149354\n",
      "0.00027302234717282954\n",
      "0.0002728871209563371\n",
      "0.0002730277396166732\n",
      "0.0002728842943634903\n",
      "0.000273020281784662\n",
      "0.0002728847420890296\n",
      "0.00027288164483100727\n",
      "0.0002728773944556318\n",
      "0.00027287457460349943\n",
      "0.0002730174900389045\n",
      "0.00027288357331199653\n",
      "0.0002728847124096775\n",
      "0.0002728837222133718\n",
      "0.0002728757472574934\n",
      "0.00027288626051431185\n",
      "0.00027288075489390474\n",
      "0.0002730184707912455\n",
      "0.00027288147198192\n",
      "0.00027288347154390047\n",
      "0.00027301682446099595\n",
      "0.00027288127201785756\n",
      "0.00027288334444622883\n",
      "0.0002728802953129314\n",
      "0.0002728763118927558\n",
      "0.0002728789675257704\n",
      "0.0002728783151335787\n",
      "0.00027287936680625693\n",
      "0.0002728751553821589\n",
      "0.0002728739964575563\n",
      "0.00027287145132156714\n",
      "0.00027287377688275325\n",
      "0.00027301482295731734\n",
      "0.00027287795313666165\n",
      "0.00027287198546120206\n",
      "0.0002728735843757281\n",
      "0.0002728702685835457\n",
      "0.00027286763293917865\n",
      "0.0002730127475149778\n",
      "0.0002728758682188493\n",
      "0.0002728711053668545\n",
      "0.0002730139397658812\n",
      "0.0002728748149613394\n",
      "0.00027286853602329124\n",
      "0.00027286992072310607\n",
      "0.0002730086137484061\n",
      "0.00027287347430194464\n",
      "0.0002728676985072488\n",
      "0.00027286715003084483\n",
      "0.0002728649493343353\n",
      "0.0002728666060521908\n",
      "0.00027286472842549293\n",
      "0.00027286195274645964\n",
      "0.00027286693353736246\n",
      "0.00027286221286780196\n",
      "0.00027285997745799055\n",
      "0.0002729952171389374\n",
      "0.00027286215864271633\n",
      "0.0002728586328959811\n",
      "0.0002728538200968273\n",
      "0.0002728605460588929\n",
      "0.00027285722632362547\n",
      "0.00027285498241415285\n",
      "0.00027285438670057\n",
      "0.00027285772709520646\n",
      "0.0002728509915222827\n",
      "0.00027298630183375137\n",
      "0.00027284834477743903\n",
      "0.0002728406729392163\n",
      "0.0002729901945412061\n",
      "0.00027285942900480577\n",
      "0.0002728508690220681\n",
      "0.00027284674530979753\n",
      "0.00027285277405349136\n",
      "0.00027284854216934296\n",
      "0.0002728437099380144\n",
      "0.0002728415763712957\n",
      "0.0002728432846851113\n",
      "0.0002728398798077937\n",
      "0.0002728352007879028\n",
      "0.00027283854140958163\n",
      "0.0002728357351183659\n",
      "0.0002728332818479947\n",
      "0.0002728264784823712\n",
      "0.0002728340658559122\n",
      "0.0002728316602569504\n",
      "0.000272828237895536\n",
      "0.0002728245230329864\n",
      "0.00027281638805892294\n",
      "0.00027282240118723126\n",
      "0.0002728191461938421\n",
      "0.0002728267581152798\n",
      "0.00027281741825288216\n",
      "0.00027281456592672965\n",
      "0.0002728076851763028\n",
      "0.0002728100517684955\n",
      "0.0002728079676015276\n",
      "0.0002728052441938327\n",
      "0.0002727967058182581\n",
      "0.00027281047686597435\n",
      "0.0002727978436807588\n",
      "0.00027280323561003426\n",
      "0.00027279377744190243\n",
      "0.00027278344806614416\n",
      "0.00027278889277783397\n",
      "0.00027278656448485536\n",
      "0.00027278783421792356\n",
      "0.0002727842776671065\n",
      "0.00027278602643173355\n",
      "0.0002727767948389014\n",
      "0.00027276648384016094\n",
      "0.0002727709084141\n",
      "0.0002727673662554037\n",
      "0.00027276486636795784\n",
      "0.00027275887107179624\n",
      "0.00027275446687363735\n",
      "0.00027273852418496406\n",
      "0.00027275157771003\n",
      "0.0002727474309890396\n",
      "0.00027275221037404817\n",
      "0.00027273890053841653\n",
      "0.00027273782956709633\n",
      "0.00027273037953130265\n",
      "0.00027273712459141976\n",
      "0.0002727187369534911\n",
      "0.0002727008581381907\n",
      "0.0002727221001289193\n",
      "0.0002727151249037755\n",
      "0.0002727126475193771\n",
      "0.00027273342678128836\n",
      "0.0002726922040654898\n",
      "0.0002726654063114556\n",
      "0.0002726916222895348\n",
      "0.0002726888785328944\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00027272566520195903\n",
      "0.00027269420984874393\n",
      "0.000272700244520402\n",
      "0.0002726931038707595\n",
      "0.00027267739184982584\n",
      "0.0002726752729851602\n",
      "0.0002726736416250534\n",
      "0.00027267033199846163\n",
      "0.00027265290875916055\n",
      "0.0002726393642133087\n",
      "0.0002726439588679087\n",
      "0.0002726424502454442\n",
      "0.00027265721600781834\n",
      "0.0002726321376721374\n",
      "0.00027261665187273303\n",
      "0.00027263576014272065\n",
      "0.000272614480858793\n",
      "0.00027259802971151204\n",
      "0.0002725954948689459\n",
      "0.00027256644706132954\n",
      "0.00027260093171937647\n",
      "0.0002725921866288112\n",
      "0.00027261818988465474\n",
      "0.0002725815505739322\n",
      "0.00027259550388583186\n",
      "0.000272551611193995\n",
      "0.0002725228704379847\n",
      "0.00027254294607710795\n",
      "0.00027256039790754435\n",
      "0.00027253982931447714\n",
      "0.0002725700391827475\n",
      "0.0002725273558801145\n",
      "0.0002725322748194271\n",
      "0.0002725214780486731\n",
      "0.0002726859710241006\n",
      "0.0002725196069106832\n",
      "0.0002725423917094962\n",
      "0.0002724971811207754\n",
      "0.00027249168126558683\n",
      "0.00027267490459713877\n",
      "0.00027254890551495334\n",
      "0.0002724793883474763\n",
      "0.0002724644874803736\n",
      "0.0002725007503238413\n",
      "0.0002724556653144858\n",
      "0.00027242787121308214\n",
      "0.00027242941871327025\n",
      "0.00027243290592145935\n",
      "0.0002724105917776989\n",
      "0.0002723659060268304\n",
      "0.00027239733511483607\n",
      "0.00027238957618177324\n",
      "0.00027240545823171096\n",
      "0.0002724423352086346\n",
      "0.0002723702041728624\n",
      "0.0002723193416661039\n",
      "0.0002722683853340632\n",
      "0.00027236392079702936\n",
      "0.0002723910026821827\n",
      "0.0002723148283243815\n",
      "0.000272373124994954\n",
      "0.0002722986445041968\n",
      "0.000272320210817692\n",
      "0.0002722410646733524\n",
      "0.00027217217502571576\n",
      "0.00027223091513366035\n",
      "0.00027224305788366344\n",
      "0.0002723174707126\n",
      "0.0002722673936699441\n",
      "0.00027230264339898755\n",
      "0.0002722038587964199\n",
      "0.0002721863831375541\n",
      "0.0002721945125471703\n",
      "0.00027218834372569475\n",
      "0.000272231841072325\n",
      "0.00027216744380577097\n",
      "0.00027220774887440417\n",
      "0.0002722293342386843\n",
      "0.0002721217581466857\n",
      "0.0002721105282909498\n",
      "0.00027211986218508377\n",
      "0.0002720932678495015\n",
      "0.00027208445119472776\n",
      "0.00027206333245745784\n",
      "0.0002720267352720254\n",
      "0.0002721047466807899\n",
      "0.0002720655530116579\n",
      "0.0002720911089018783\n",
      "0.0002720452269789125\n",
      "0.0002720124624107462\n",
      "0.0002722050715802899\n",
      "0.000272063714316776\n",
      "0.00027196224366656717\n",
      "0.00027194664543042684\n",
      "0.000272022479863045\n",
      "0.00027190483364353875\n",
      "0.0002718576724234008\n",
      "0.00027207410645080355\n",
      "0.000271965192708902\n",
      "0.00027216935934595095\n",
      "0.0002720177710882497\n",
      "0.0002720581832689473\n",
      "0.0002719712692027138\n",
      "0.00027201682707535713\n",
      "0.00027187164236933164\n",
      "0.0002719799089520177\n",
      "0.0002719563052915699\n",
      "0.0002720579814060398\n",
      "0.00027197034054352514\n",
      "0.0002718555302532855\n",
      "0.0002718499444257369\n",
      "0.00027188200441847577\n",
      "0.00027179845788337565\n",
      "0.000271962509349104\n",
      "0.0002719976518331847\n",
      "0.00027191359022827845\n",
      "0.00027193748577930023\n",
      "0.0002718565760949701\n",
      "0.0002718819349153003\n",
      "0.0002717436065689566\n",
      "0.00027187982496496715\n",
      "0.0002719307923997653\n",
      "0.0002719003171315728\n",
      "0.0002716991722427803\n",
      "0.0002716508188530543\n",
      "0.0002720757755568865\n",
      "0.0002718187791034141\n",
      "0.0002719290780761147\n",
      "0.0002718420820350619\n",
      "0.0002718163847282622\n",
      "0.00027195049207475846\n",
      "0.0002718023039640823\n",
      "0.0002720051836853465\n",
      "0.0002718107913269283\n",
      "0.0002719919414102931\n",
      "0.00027180964767165017\n",
      "0.0002717134801917924\n",
      "0.00027170997578532133\n",
      "0.00027176314531747545\n",
      "0.0002717558940274243\n",
      "0.0002716642869795462\n",
      "0.00027183009886658476\n",
      "0.0002717618273304245\n",
      "0.000271766847130663\n",
      "0.0002717958169614513\n",
      "0.0002717786232101729\n",
      "0.00027162790330378624\n",
      "0.000271644571061214\n",
      "0.00027163420501338093\n",
      "0.0002717351431389348\n",
      "0.00027179774123135716\n",
      "0.00027163624322819427\n",
      "0.0002717099419350702\n",
      "0.0002716938840131431\n",
      "0.0002717422959466617\n",
      "0.00027182284317637726\n",
      "0.00027168579303649544\n",
      "0.0002716667327813307\n",
      "0.00027173285548440847\n",
      "0.0002717719930267003\n",
      "0.00027168653807244146\n",
      "0.0002717372185323739\n",
      "0.00027168624320901753\n",
      "0.00027163248394484953\n",
      "0.0002716010083214019\n",
      "0.000271762040783531\n",
      "0.0002716313984703234\n",
      "0.0002717082825859542\n",
      "0.00027164175115858627\n",
      "0.0002715592413412539\n",
      "0.00027166359508028296\n",
      "0.00027155099409558933\n",
      "0.00027166474314572884\n",
      "0.0002717068394562694\n",
      "0.00027163981563933754\n",
      "0.00027180634003800836\n",
      "0.00027162186092016696\n",
      "0.000271666485864813\n",
      "0.0002716176253096887\n",
      "0.00027163320926527527\n",
      "0.00027159599406038437\n",
      "0.0002715539147186702\n",
      "0.00027161486824563835\n",
      "0.00027160323605676087\n",
      "0.0002716155310290302\n",
      "0.0002717541678303337\n",
      "0.00027160137892382185\n",
      "0.000271756294374033\n",
      "0.00027160021644366317\n",
      "0.0002715460705748954\n",
      "0.0002715222045634295\n",
      "0.000271536479512574\n",
      "0.0002715787191408706\n",
      "0.00027157745148234744\n",
      "0.00027171821635467105\n",
      "0.00027158161267892866\n",
      "0.0002715991855632114\n",
      "0.00027153746078624464\n",
      "0.00027152171178075054\n",
      "0.00027150546606648565\n",
      "0.0002716626737765689\n",
      "0.00027157004175299316\n",
      "0.0002715465130093108\n",
      "0.00027165884192976075\n",
      "0.0002715646153474549\n",
      "0.00027150336011948277\n",
      "0.0002716489170573307\n",
      "0.00027154567084864003\n",
      "0.0002717109798712495\n",
      "0.000271544941555155\n",
      "0.0002716528383135206\n",
      "0.00027154783167981586\n",
      "0.0002716554917919721\n",
      "0.0002715445453490703\n",
      "0.0002715198521348684\n",
      "0.0002715056971587768\n",
      "0.0002715034819596149\n",
      "0.0002714883023644236\n",
      "0.0002714696683469796\n",
      "0.0002715360306961172\n",
      "0.0002715559166904998\n",
      "0.0002715168070875676\n",
      "0.00027149686458837456\n",
      "0.00027147532097878004\n",
      "0.00027149491055428535\n",
      "0.00027149210005018756\n",
      "0.0002714925576240994\n",
      "0.0002714810275822475\n",
      "0.0002714604649348834\n",
      "0.00027144788901154723\n",
      "0.000271471634515059\n",
      "0.00027145324676424115\n",
      "0.00027163162956693324\n",
      "0.00027148470025155103\n",
      "0.00027161239055250956\n",
      "0.0002714853393532933\n",
      "0.0002714771244972902\n",
      "0.0002714903943034278\n",
      "0.00027144235574072984\n",
      "0.0002715877828436004\n",
      "0.00027149309887176606\n",
      "0.0002714711657428602\n",
      "0.000271440989834003\n",
      "0.0002715949820231452\n",
      "0.0002716459579154164\n",
      "0.00027146852645931014\n",
      "0.0002714387407270273\n",
      "0.00027159447052393037\n",
      "0.00027160657794003315\n",
      "0.0002714677839154522\n",
      "0.00027144436618971307\n",
      "0.00027162000086452333\n",
      "0.00027146093563646333\n",
      "0.0002714499828770527\n",
      "0.00027144692699122944\n",
      "0.00027143277081616225\n",
      "0.0002715962175745814\n",
      "0.0002714407744452092\n",
      "0.00027158614457265613\n",
      "0.0002714541480339049\n",
      "0.0002715921933242683\n",
      "0.00027145273682923244\n",
      "0.0002714271450871276\n",
      "0.00027158980401561654\n",
      "0.0002714228093061007\n",
      "0.0002715765241999698\n",
      "0.0002715892450786791\n",
      "0.000271443507231145\n",
      "0.00027158692690005886\n",
      "0.0002714431530301441\n",
      "0.00027160065323831846\n",
      "0.00027144020984008467\n",
      "0.0002714358844271289\n",
      "0.00027144159797417645\n",
      "0.0002714415241208153\n",
      "0.0002715885380820932\n",
      "0.0002714379401167025\n",
      "0.0002714220932907679\n",
      "0.0002715804920680378\n",
      "0.0002714296153641266\n",
      "0.000271442594689245\n",
      "0.0002714330367371627\n",
      "0.0002714346508082892\n",
      "0.00027144132448808184\n",
      "0.0002714324743757505\n",
      "0.0002715964501359512\n",
      "0.0002714333237792043\n",
      "0.0002715866898198485\n",
      "0.00027143403408756\n",
      "0.0002714210682957618\n",
      "0.0002714167559633669\n",
      "0.0002714157341533815\n",
      "0.00027157338954157353\n",
      "0.0002716004893225143\n",
      "0.0002714282297753119\n",
      "0.0002714309214754835\n",
      "0.00027141570636992073\n",
      "0.0002715757129457339\n",
      "0.00027158123453894964\n",
      "0.0002714277978331539\n",
      "0.0002714206610349769\n",
      "0.0002714168167593699\n",
      "0.0002714143572459226\n",
      "0.0002714149526956848\n",
      "0.0002714204235895608\n",
      "0.00027141526996584417\n",
      "0.00027141499731372676\n",
      "0.00027157439179889927\n",
      "0.0002714214034585571\n",
      "0.00027140671183432436\n",
      "0.0002714006339357374\n",
      "0.0002714064897674772\n",
      "0.00027141217799096605\n",
      "0.0002714030742829981\n",
      "0.00027140550757727324\n",
      "0.0002714055210879667\n",
      "0.0002714039125500824\n",
      "0.000271566576561211\n",
      "0.0002714113100856982\n",
      "0.0002714014812392031\n",
      "0.00027140295756101476\n",
      "0.00027140560930229705\n",
      "0.00027139822493004026\n",
      "0.00027139935967289614\n",
      "0.000271395518005673\n",
      "0.0002713948418076729\n",
      "0.00027139422943568686\n",
      "0.0002715619178884227\n",
      "0.0002713878690436283\n",
      "0.00027137984529686325\n",
      "0.00027138914404658854\n",
      "0.0002713991733239026\n",
      "0.00027140257731312836\n",
      "0.0002713834488716184\n",
      "0.00027138487969866206\n",
      "0.0002713874542938168\n",
      "0.0002713837709108752\n",
      "0.00027138664832121256\n",
      "0.00027138520599220205\n",
      "0.0002713781951768445\n",
      "0.0002713791627466379\n",
      "0.0002713802003056328\n",
      "0.00027136811710446303\n",
      "0.0002715244070416172\n",
      "0.00027155102529824384\n",
      "0.00027138538036688465\n",
      "0.00027137293324752486\n",
      "0.00027136992104945636\n",
      "0.0002713694367804616\n",
      "0.000271388121443421\n",
      "0.0002713784592200211\n",
      "0.0002713677623627662\n",
      "0.00027153165257987397\n",
      "0.0002713762738987708\n",
      "0.00027136620255041214\n",
      "0.000271532901975099\n",
      "0.00027137473503148957\n",
      "0.0002713613994016666\n",
      "0.00027135479543946534\n",
      "0.0002713622293834285\n",
      "0.0002713587180119316\n",
      "0.0002713607486527243\n",
      "0.0002713531384992814\n",
      "0.00027151636590941244\n",
      "0.0002713577348266372\n",
      "0.00027136190282695764\n",
      "0.00027151732363423067\n",
      "0.00027136596023113395\n",
      "0.0002713525674956874\n",
      "0.0002715186714274838\n",
      "0.00027135585278260737\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0002713478166335437\n",
      "0.00027134110805349625\n",
      "0.0002713472120634738\n",
      "0.0002713504918691341\n",
      "0.00027150761613404324\n",
      "0.000271358526496306\n",
      "0.0002713503830172564\n",
      "0.0002713602617230501\n",
      "0.0002713410582687678\n",
      "0.0002715028949721528\n",
      "0.0002713594121193372\n",
      "0.000271351808655043\n",
      "0.00027151626577655877\n",
      "0.00027135282466486974\n",
      "0.0002715092354412629\n",
      "0.0002713532537484347\n",
      "0.00027134556466081776\n",
      "0.00027134351390865116\n",
      "0.00027134010717332757\n",
      "0.0002713362063404094\n",
      "0.0002713365402812933\n",
      "0.0002713434019784763\n",
      "0.000271503145419378\n",
      "0.0002713474039701004\n",
      "0.0002713384535618135\n",
      "0.00027134020051802375\n",
      "0.00027133569997114283\n",
      "0.0002715075485602115\n",
      "0.00027133465935871864\n",
      "0.00027133715585883893\n",
      "0.0002713295799700447\n",
      "0.0002714931312364851\n",
      "0.000271326446797184\n",
      "0.0002713177714443299\n",
      "0.0002713260049883849\n",
      "0.0002713278005300429\n",
      "0.00027132587786682534\n",
      "0.0002713223915536341\n",
      "0.00027132113990047417\n",
      "0.0002713196445323079\n",
      "0.00027131948824890193\n",
      "0.0002713112170323598\n",
      "0.00027129924882960137\n",
      "0.0002714855598963602\n",
      "0.0002713274862154314\n",
      "0.0002713110716382713\n",
      "0.0002713078722581049\n",
      "0.00027130456570652003\n",
      "0.00027131033570456406\n",
      "0.00027129996130047273\n",
      "0.00027130069661610663\n",
      "0.0002712972482016542\n",
      "0.00027129107558466575\n",
      "0.00027128961963923963\n",
      "0.0002712750145349227\n",
      "0.00027128741430182564\n",
      "0.00027128452132204624\n",
      "0.0002712786628192288\n",
      "0.0002712770392228862\n",
      "0.00027127451049253475\n",
      "0.0002712652017464033\n",
      "0.0002712684162919198\n",
      "0.00027126244728042145\n",
      "0.00027124872557482123\n",
      "0.0002712506395666842\n",
      "0.00027124970595626224\n",
      "0.000271240359510656\n",
      "0.00027121598764636166\n",
      "0.0002712227341161707\n",
      "0.00027123831482821014\n",
      "0.00027122666261942266\n",
      "0.00027121325979238316\n",
      "0.000271185121094448\n",
      "0.0002712054530818375\n",
      "0.00027136668188332746\n",
      "0.00027125387835897144\n",
      "0.0002711933865353894\n",
      "0.00027119327960855433\n",
      "0.0002712014210282838\n",
      "0.0002711740294359851\n",
      "0.0002711353640540166\n",
      "0.00027116577232536016\n",
      "0.00027115478936382453\n",
      "0.0002711477652843883\n",
      "0.0002711603202009005\n",
      "0.00027114416818188624\n",
      "0.0002711275070391126\n",
      "0.00027109602674626283\n",
      "0.0002711052891009317\n",
      "0.00027110788558880973\n",
      "0.00027110413276766306\n",
      "0.00027108694580564715\n",
      "0.00027104966279207484\n",
      "0.0002710660215980886\n",
      "0.0002710588048742254\n",
      "0.0002710454472687408\n",
      "0.00027116913694479565\n",
      "0.0002712235117547675\n",
      "0.00027112424612836387\n",
      "0.0002710442486780817\n",
      "0.00027118444740165896\n",
      "0.00027105670628824267\n",
      "0.0002710207737128859\n",
      "0.00027115909425003837\n",
      "0.00027104335911867035\n",
      "0.0002710049507715905\n",
      "0.0002711369142575586\n",
      "0.00027101076071903363\n",
      "0.00027100637466362566\n",
      "0.00027101441888958814\n",
      "0.0002711378199753117\n",
      "0.00027106103860501167\n",
      "0.00027100984615662004\n",
      "0.0002709874632053857\n",
      "0.00027112956861503326\n",
      "0.0002709972466044158\n",
      "0.0002711634644739687\n",
      "0.00027103198061459644\n",
      "0.0002709926776499377\n",
      "0.00027099667949473226\n",
      "0.0002709915571310833\n",
      "0.00027113412322808\n",
      "0.00027102001764905124\n",
      "0.00027114774073119807\n",
      "0.00027101345771417804\n",
      "0.0002710149794543558\n",
      "0.00027097540221600514\n",
      "0.0002711299733677276\n",
      "0.0002710118285165865\n",
      "0.0002709964366354181\n",
      "0.0002711471920973757\n",
      "0.00027100212701429313\n",
      "0.0002710087739647352\n",
      "0.00027099012816377623\n",
      "0.00027114890589939757\n",
      "0.0002709986314966007\n",
      "0.00027100206866359315\n",
      "0.00027097419395597736\n",
      "0.00027096413277242625\n",
      "0.00027098098156918003\n",
      "0.00027096768965276675\n",
      "0.0002710010894877232\n",
      "0.0002709860102407608\n",
      "0.0002711366850276025\n",
      "0.0002709888963640866\n",
      "0.0002709781792097331\n",
      "0.00027096951675698825\n",
      "0.0002709950345788296\n",
      "0.0002709811694728167\n",
      "0.00027098134513146307\n",
      "0.00027096148388025955\n",
      "0.0002711257081155577\n",
      "0.0002709664698306034\n",
      "0.0002709548473452279\n",
      "0.0002711171213213303\n",
      "0.00027097020906385676\n",
      "0.00027095688142544016\n",
      "0.00027114111995222754\n",
      "0.00027097086179111134\n",
      "0.0002711338846955443\n",
      "0.0002709713643234558\n",
      "0.0002711331463593193\n",
      "0.0002709710536950316\n",
      "0.000271129180239092\n",
      "0.00027096997313820123\n",
      "0.0002709586584540769\n",
      "0.00027095646202699786\n",
      "0.00027095527779324525\n",
      "0.0002711280985321673\n",
      "0.00027096460581232256\n",
      "0.0002711492505800838\n",
      "0.00027096198156016557\n",
      "0.00027112536209630524\n",
      "0.00027096430161234594\n",
      "0.0002709556393362421\n",
      "0.0002709589525879051\n",
      "0.00027095542018967344\n",
      "0.0002711245641738095\n",
      "0.00027096033608361194\n",
      "0.00027094958627278833\n",
      "0.0002711201860182107\n",
      "0.00027112519838646474\n",
      "0.00027095902338676435\n",
      "0.0002711277450314116\n",
      "0.0002709578338450739\n",
      "0.000270961099127871\n",
      "0.00027095649765892106\n",
      "0.000271125335685441\n",
      "0.0002709573470270719\n",
      "0.00027095187012773984\n",
      "0.0002709573916335104\n",
      "0.00027095112964102997\n",
      "0.00027112418290348123\n",
      "0.00027095546255889497\n",
      "0.0002709572284353567\n",
      "0.0002711241780817409\n",
      "0.00027095517922761677\n",
      "0.00027095697607490036\n",
      "0.000270954117884824\n",
      "0.0002709558538545626\n",
      "0.0002709531311608458\n",
      "0.0002709537009930449\n",
      "0.0002711296964871142\n",
      "0.0002709535030915702\n",
      "0.0002711262694312719\n",
      "0.0002709537845427078\n",
      "0.0002709494770675394\n",
      "0.00027112150392430586\n",
      "0.00027112655689553794\n",
      "0.00027095324067711354\n",
      "0.00027095589248628423\n",
      "0.00027095275264088324\n",
      "0.0002711236314362589\n",
      "0.00027095336284423084\n",
      "0.00027095048474883997\n",
      "0.000270949927210334\n",
      "0.00027094857052770355\n",
      "0.0002711214895761647\n",
      "0.0002709524203666606\n",
      "0.00027112408332413473\n",
      "0.0002709518381975309\n",
      "0.0002711225756343718\n",
      "0.000270951869294393\n",
      "0.00027112385551044825\n",
      "0.00027095162723773724\n",
      "0.00027095208680888455\n",
      "0.0002709485103638419\n",
      "0.00027094766469447953\n",
      "0.00027095231614457214\n",
      "0.00027095046887381494\n",
      "0.000270952140225969\n",
      "0.0002709504112151578\n",
      "0.00027095103042948864\n",
      "0.0002709480204397297\n",
      "0.00027094748042406534\n",
      "0.0002711213350913047\n",
      "0.0002709469969218596\n",
      "0.00027112077967931424\n",
      "0.0002709498327958165\n",
      "0.0002709493477928539\n",
      "0.00027112183207494994\n",
      "0.0002709494164470882\n",
      "0.0002709471422694358\n",
      "0.0002711215181663303\n",
      "0.00027094914532569184\n",
      "0.0002709466142373622\n",
      "0.0002711208718544658\n",
      "0.00027112346285702036\n",
      "0.0002709484942165246\n",
      "0.00027094692918516034\n",
      "0.0002709474188221302\n",
      "0.00027112079526697415\n",
      "0.00027094838238271606\n",
      "0.0002709482092883851\n",
      "0.00027094593337923054\n",
      "0.00027112008750769143\n",
      "0.00027094799597884704\n",
      "0.00027094578904135914\n",
      "0.00027111971385869883\n",
      "0.000270945491112337\n",
      "0.0002711196445072921\n",
      "0.00027112249697825217\n",
      "0.00027094715102236124\n",
      "0.00027094534820136117\n",
      "0.000271120000699541\n",
      "0.0002711215673663437\n",
      "0.0002709469398521038\n",
      "0.00027112088466235305\n",
      "0.0002709468996613519\n",
      "0.0002709456219860872\n",
      "0.0002711209087399348\n",
      "0.000270946612274866\n",
      "0.0002711210046812711\n",
      "0.00027094656788701616\n",
      "0.00027094694189458223\n",
      "0.00027094601297312305\n",
      "0.0002709466826921656\n",
      "0.00027112036974145167\n",
      "0.000270946425692703\n",
      "0.000270945562248089\n",
      "0.0002711202435969946\n",
      "0.0002709463381581653\n",
      "0.00027094639188837065\n",
      "0.00027094610486008057\n",
      "0.00027094540766022483\n",
      "0.00027112023428297683\n",
      "0.0002709460719023182\n",
      "0.00027094489377813946\n",
      "0.0002711192346113634\n",
      "0.0002709456297275544\n",
      "0.00027111990795743656\n",
      "0.00027094586136832695\n",
      "0.0002711208512984319\n",
      "0.00027094566333298433\n",
      "0.0002711199621993157\n",
      "0.00027094574009451283\n",
      "0.0002711210883700098\n",
      "0.0002709455712281121\n",
      "0.0002709450310211099\n",
      "0.0002709448675670396\n",
      "0.0002711196622583846\n",
      "0.0002709453314483488\n",
      "0.00027112004336165456\n",
      "0.0002709454716215842\n",
      "0.00027112014619574086\n",
      "0.0002709454183619324\n",
      "0.00027094578071899106\n",
      "0.00027094530914316287\n",
      "0.00027112024418924334\n",
      "0.000270945363646386\n",
      "0.0002711202799455186\n",
      "0.00027094534950897496\n",
      "0.00027094522327451533\n",
      "0.0002709449282405445\n",
      "0.00027094479367681114\n",
      "0.00027111979081388254\n",
      "0.000271119948279155\n",
      "0.0002709452346811358\n",
      "0.0002709453212962276\n",
      "0.00027111999626957845\n",
      "0.0002709451972811391\n",
      "0.0002711200638546612\n",
      "0.0002709451840022821\n",
      "0.0002709449123738299\n",
      "0.00027094498612987303\n",
      "0.00027094536711363894\n",
      "0.0002709450607475276\n",
      "0.0002709449404195455\n",
      "0.00027094463619084107\n",
      "0.0002711195729573739\n",
      "0.00027094487156708637\n",
      "0.00027094454484926224\n",
      "0.0002709442682365983\n",
      "0.00027094440306240595\n",
      "0.00027094461391390275\n",
      "0.0002711195481879461\n",
      "0.0002709448414656883\n",
      "0.00027094475930736906\n",
      "0.00027094469782482864\n",
      "0.00027111959966507053\n",
      "0.00027094477023167446\n",
      "0.0002711197059970703\n",
      "0.00027094475878621036\n",
      "0.0002709445217094546\n",
      "0.00027111969096865594\n",
      "0.0002709447077985712\n",
      "0.00027094438636217154\n",
      "0.0002709443212099645\n",
      "0.00027094430538367225\n",
      "0.00027094423170994215\n",
      "0.0002711192198643765\n",
      "0.0002709444366584876\n",
      "0.0002711193014798589\n",
      "0.0002709445776177724\n",
      "0.0002709441535349264\n",
      "0.00027111903466358077\n",
      "0.0002709442910974628\n",
      "0.00027094416558992736\n",
      "0.0002709441193011148\n",
      "0.00027094410496344555\n",
      "0.0002709438961821782\n",
      "0.0002709435867237951\n",
      "0.00027094391451857084\n",
      "0.0002711195604224464\n",
      "0.00027094422398748406\n",
      "0.00027094384133778465\n",
      "0.00027111896660442995\n",
      "0.00027094421061760903\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00027094382718461193\n",
      "0.0002711189542893902\n",
      "0.00027094415255945624\n",
      "0.00027094401267103134\n",
      "0.0002709437792304456\n",
      "0.0002709436547788299\n",
      "0.00027094368123745615\n",
      "0.0002709435540284374\n",
      "0.00027094332603835994\n",
      "0.00027094344168238186\n",
      "0.00027094345088352295\n",
      "0.00027094320886723405\n",
      "0.00027094282022671124\n",
      "0.00027094342723318194\n",
      "0.0002709431899387073\n",
      "0.00027094327255597535\n",
      "0.0002709430294906382\n",
      "0.0002709429494775275\n",
      "0.0002709427541079141\n",
      "0.000271117709174621\n",
      "0.0002709427493580236\n",
      "0.0002709424154407709\n",
      "0.00027111803616626523\n",
      "0.0002709433403047918\n",
      "0.00027094261813331573\n",
      "0.00027094255859057314\n",
      "0.0002709425822539194\n",
      "0.00027094282161893265\n",
      "0.0002709423630643805\n",
      "0.00027094207074272237\n",
      "0.00027094215507261716\n",
      "0.00027111729260212035\n",
      "0.0002709429141809682\n",
      "0.0002709425435610889\n",
      "0.00027094205811299454\n",
      "0.0002711173085975865\n",
      "0.00027094205645369484\n",
      "0.0002709418068293904\n",
      "0.00027094180666840256\n",
      "0.00027111686973951707\n",
      "0.0002711172724000679\n",
      "0.0002709424944146696\n",
      "0.0002709417878075586\n",
      "0.0002709414969628333\n",
      "0.0002709418222699712\n",
      "0.0002711172024929267\n",
      "0.0002709423109424491\n",
      "0.0002709417078772268\n",
      "0.0002709420696506375\n",
      "0.00027094146121189617\n",
      "0.0002709411888691188\n",
      "0.00027094125007592977\n",
      "0.00027094138245264954\n",
      "0.0002709411492687678\n",
      "0.00027094092547848617\n",
      "0.00027094141138653705\n",
      "0.00027094098738859716\n",
      "0.0002709410160120688\n",
      "0.0002709409713607366\n",
      "0.00027094067334831976\n",
      "0.00027094031825746206\n",
      "0.0002709403999963892\n",
      "0.0002709405631915836\n",
      "0.00027094059726522524\n",
      "0.00027094021635520074\n",
      "0.00027093989959997686\n",
      "0.0002709402393605632\n",
      "0.0002709401879679027\n",
      "0.0002709400412651354\n",
      "0.0002709405046255602\n",
      "0.0002709405958273421\n",
      "0.0002709400858984519\n",
      "0.00027093972081687113\n",
      "0.0002711155216850546\n",
      "0.0002709398314435118\n",
      "0.00027093974385910765\n",
      "0.0002709401087906215\n",
      "0.00027093952334852824\n",
      "0.0002711150670491932\n",
      "0.0002709403372100073\n",
      "0.0002709401073789215\n",
      "0.00027094036373886206\n",
      "0.00027093994789473936\n",
      "0.0002709396936346321\n",
      "0.0002709397627989058\n",
      "0.0002709397842009776\n",
      "0.0002709403853092899\n",
      "0.0002709397784736185\n",
      "0.00027093983092698693\n",
      "0.0002709394307860791\n",
      "0.0002709393896115321\n",
      "0.0002709392176432865\n",
      "0.0002711147667166618\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-27-b7b6a4e8e33a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     37\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mopt_parameters\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     38\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 39\u001b[0;31m \u001b[0mnew_parameters\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfit_model\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreference_dimer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     40\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"new rms error =\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_error\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreference_dimer\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnew_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     41\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnew_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-27-b7b6a4e8e33a>\u001b[0m in \u001b[0;36mfit_model\u001b[0;34m(reference_dimer, model_parameters, lock)\u001b[0m\n\u001b[1;32m     25\u001b[0m     result = scipy.optimize.minimize(objective_wrapper, param_list,\n\u001b[1;32m     26\u001b[0m                                      \u001b[0;34m(\u001b[0m\u001b[0mreference_dimer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlock\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 27\u001b[0;31m                                      method='Nelder-Mead', options=min_opt)\n\u001b[0m\u001b[1;32m     28\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     29\u001b[0m     \u001b[0mopt_parameters\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/_minimize.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[1;32m    589\u001b[0m                       callback=callback, **options)\n\u001b[1;32m    590\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'nelder-mead'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 591\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_minimize_neldermead\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    592\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'powell'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    593\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_minimize_powell\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36m_minimize_neldermead\u001b[0;34m(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, adaptive, **unknown_options)\u001b[0m\n\u001b[1;32m    566\u001b[0m         \u001b[0mxbar\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnumpy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msim\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    567\u001b[0m         \u001b[0mxr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mrho\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mxbar\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mrho\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msim\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 568\u001b[0;31m         \u001b[0mfxr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    569\u001b[0m         \u001b[0mdoshrink\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    570\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[0;34m(*wrapper_args)\u001b[0m\n\u001b[1;32m    291\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    292\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 293\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    294\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    295\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-27-b7b6a4e8e33a>\u001b[0m in \u001b[0;36mobjective_wrapper\u001b[0;34m(param_list, reference_dimer, model_parameters, lock)\u001b[0m\n\u001b[1;32m     11\u001b[0m             \u001b[0mnew_parameters\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;34m{\u001b[0m \u001b[0mkey\u001b[0m \u001b[0;34m:\u001b[0m \u001b[0mparam_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m}\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m             \u001b[0mi\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mmodel_error\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreference_dimer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     14\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfit_model\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreference_dimer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlock\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-26-27d580213203>\u001b[0m in \u001b[0;36mmodel_error\u001b[0;34m(reference_dimer, model_parameters)\u001b[0m\n\u001b[1;32m      6\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mdata_list\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mreference_dimer\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m         \u001b[0mmol_coord\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mdata_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m         \u001b[0mmol_occupied\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmol_hf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmol_mp2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_reference_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmol_coord\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      9\u001b[0m         \u001b[0mmol_hf\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0;36m2.0\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0matom_hf\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m         \u001b[0mmol_mp2\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0;36m2.0\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0matom_mp2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-25-e9ea21355ab8>\u001b[0m in \u001b[0;36mcalculate_reference_data\u001b[0;34m(atomic_coordinates, model_parameters)\u001b[0m\n\u001b[1;32m     10\u001b[0m     \u001b[0mchi_tensor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_chi_tensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0matomic_coordinates\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m     \u001b[0mdensity_matrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfock_matrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscf_cycle\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhamiltonian_matrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minteraction_matrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdensity_matrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mchi_tensor\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m     \u001b[0moccupied_energy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvirtual_energy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moccupied_matrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvirtual_matrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpartition_orbitals\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfock_matrix\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-17-8da032be3a7d>\u001b[0m in \u001b[0;36mscf_cycle\u001b[0;34m(hamiltonian_matrix, interaction_matrix, density_matrix, chi_tensor)\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0mCONVERGENCE_TOLERANCE\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1e-4\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0miteration\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mMAX_SCF_ITERATIONS\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m         \u001b[0mfock_matrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_fock_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhamiltonian_matrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minteraction_matrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdensity_matrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mchi_tensor\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      8\u001b[0m         \u001b[0mnew_density_matrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_density_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfock_matrix\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-15-91080e61f96a>\u001b[0m in \u001b[0;36mcalculate_fock_matrix\u001b[0;34m(hamiltonian_matrix, interaction_matrix, density_matrix, chi_tensor)\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0mfock_matrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mhamiltonian_matrix\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     fock_matrix += 2.0*np.einsum('pqt,rsu,tu,rs',\n\u001b[0;32m----> 5\u001b[0;31m                                  chi_tensor, chi_tensor, interaction_matrix, density_matrix, optimize=True)\n\u001b[0m\u001b[1;32m      6\u001b[0m     fock_matrix -= np.einsum('rqt,psu,tu,rs',\n\u001b[1;32m      7\u001b[0m                              chi_tensor, chi_tensor, interaction_matrix, density_matrix, optimize=True)\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/numpy/core/einsumfunc.py\u001b[0m in \u001b[0;36meinsum\u001b[0;34m(*operands, **kwargs)\u001b[0m\n\u001b[1;32m   1249\u001b[0m     \u001b[0;31m# Build the contraction list and operand\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1250\u001b[0m     operands, contraction_list = einsum_path(*operands, optimize=optimize_arg,\n\u001b[0;32m-> 1251\u001b[0;31m                                              einsum_call=True)\n\u001b[0m\u001b[1;32m   1252\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1253\u001b[0m     \u001b[0mhandle_out\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/numpy/core/einsumfunc.py\u001b[0m in \u001b[0;36meinsum_path\u001b[0;34m(*operands, **kwargs)\u001b[0m\n\u001b[1;32m    929\u001b[0m         \u001b[0;31m# If we're broadcasting, nix blas\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    930\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0midx_removed\u001b[0m \u001b[0;34m&\u001b[0m \u001b[0mbcast\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 931\u001b[0;31m             \u001b[0mdo_blas\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_can_dot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtmp_inputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout_inds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0midx_removed\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    932\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    933\u001b[0m             \u001b[0mdo_blas\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/numpy/core/einsumfunc.py\u001b[0m in \u001b[0;36m_can_dot\u001b[0;34m(inputs, result, idx_removed)\u001b[0m\n\u001b[1;32m    467\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mc\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minput_left\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0minput_right\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    468\u001b[0m         \u001b[0;31m# can't deal with repeated indices on same input or more than 2 total\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 469\u001b[0;31m         \u001b[0mnl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minput_left\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcount\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minput_right\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcount\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    470\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mnl\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mnr\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mnl\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mnr\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    471\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "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": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.11610654314199755\n",
      "0.11610654589654693\n",
      "0.1563845691212358\n",
      "0.11610316005164277\n",
      "0.06264069799422883\n",
      "0.11610752325292063\n",
      "0.11613172626874084\n",
      "0.11612114365668312\n",
      "0.11611108696412328\n",
      "0.1161068246790123\n",
      "0.11610643020834462\n",
      "0.11610654204537149\n",
      "0.11610746458627959\n",
      "0.06694133826682498\n",
      "0.09897820302291077\n",
      "0.09613303530792087\n",
      "0.09280773781154286\n",
      "0.08892879552107283\n",
      "0.08439960351638082\n",
      "0.07911638071740897\n",
      "0.0729529041015107\n",
      "0.06576261039185995\n",
      "0.05737541993738627\n",
      "0.028051890020415353\n",
      "0.04270419106966426\n",
      "0.030494066045974315\n",
      "0.036227576610483196\n",
      "0.029067806939240725\n",
      "0.021813715194559057\n",
      "0.014295303365823995\n",
      "0.008835291525373901\n",
      "0.03203695896361719\n",
      "0.002913053621202024\n",
      "0.04194831731209078\n",
      "0.008186498536307403\n",
      "0.015190565873969518\n",
      "0.022624427527434633\n",
      "0.036141519975961184\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-71-26fc18985595>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     32\u001b[0m ]\n\u001b[1;32m     33\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 34\u001b[0;31m \u001b[0mneon_param\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfit_model\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mneon_ref_data\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     35\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"neon rms error =\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_error\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mneon_ref_data\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mneon_param\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     36\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mneon_param\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-27-b7b6a4e8e33a>\u001b[0m in \u001b[0;36mfit_model\u001b[0;34m(reference_dimer, model_parameters, lock)\u001b[0m\n\u001b[1;32m     25\u001b[0m     result = scipy.optimize.minimize(objective_wrapper, param_list,\n\u001b[1;32m     26\u001b[0m                                      \u001b[0;34m(\u001b[0m\u001b[0mreference_dimer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlock\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 27\u001b[0;31m                                      method='Nelder-Mead', options=min_opt)\n\u001b[0m\u001b[1;32m     28\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     29\u001b[0m     \u001b[0mopt_parameters\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/_minimize.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[1;32m    589\u001b[0m                       callback=callback, **options)\n\u001b[1;32m    590\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'nelder-mead'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 591\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_minimize_neldermead\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    592\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'powell'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    593\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_minimize_powell\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36m_minimize_neldermead\u001b[0;34m(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, adaptive, **unknown_options)\u001b[0m\n\u001b[1;32m    566\u001b[0m         \u001b[0mxbar\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnumpy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msim\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    567\u001b[0m         \u001b[0mxr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mrho\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mxbar\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mrho\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msim\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 568\u001b[0;31m         \u001b[0mfxr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    569\u001b[0m         \u001b[0mdoshrink\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    570\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[0;34m(*wrapper_args)\u001b[0m\n\u001b[1;32m    291\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    292\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 293\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    294\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    295\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-27-b7b6a4e8e33a>\u001b[0m in \u001b[0;36mobjective_wrapper\u001b[0;34m(param_list, reference_dimer, model_parameters, lock)\u001b[0m\n\u001b[1;32m     11\u001b[0m             \u001b[0mnew_parameters\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;34m{\u001b[0m \u001b[0mkey\u001b[0m \u001b[0;34m:\u001b[0m \u001b[0mparam_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m}\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m             \u001b[0mi\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mmodel_error\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreference_dimer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     14\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfit_model\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreference_dimer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlock\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-26-27d580213203>\u001b[0m in \u001b[0;36mmodel_error\u001b[0;34m(reference_dimer, model_parameters)\u001b[0m\n\u001b[1;32m      6\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mdata_list\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mreference_dimer\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m         \u001b[0mmol_coord\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mdata_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m         \u001b[0mmol_occupied\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmol_hf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmol_mp2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_reference_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmol_coord\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      9\u001b[0m         \u001b[0mmol_hf\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0;36m2.0\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0matom_hf\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m         \u001b[0mmol_mp2\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0;36m2.0\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0matom_mp2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-25-e9ea21355ab8>\u001b[0m in \u001b[0;36mcalculate_reference_data\u001b[0;34m(atomic_coordinates, model_parameters)\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m     \u001b[0mhamiltonian_matrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_hamiltonian_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0matomic_coordinates\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m     \u001b[0minteraction_matrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_interaction_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0matomic_coordinates\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m     \u001b[0mchi_tensor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_chi_tensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0matomic_coordinates\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-46-345c7257ba50>\u001b[0m in \u001b[0;36mcalculate_interaction_matrix\u001b[0;34m(atomic_coordinates, model_parameters)\u001b[0m\n\u001b[1;32m      7\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0matom\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0matom\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m                 \u001b[0mr_pq\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0matomic_coordinates\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0matom\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0matomic_coordinates\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0matom\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m                 \u001b[0minteraction_matrix\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcoulomb_energy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0morb\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0morb\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr_pq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mp\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mq\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0morb\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m's'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m                 \u001b[0minteraction_matrix\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel_parameters\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'coulomb_s'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-38-802ce7b64f6e>\u001b[0m in \u001b[0;36mcoulomb_energy\u001b[0;34m(o1, o2, r12)\u001b[0m\n\u001b[1;32m      9\u001b[0m         \u001b[0mans\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvec\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mo1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr12\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mr12_length\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mo1\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mp_orbitals\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mo2\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mp_orbitals\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m         ans = ( np.dot(vec[o1], vec[o2]) / r12_length**3\n\u001b[0m\u001b[1;32m     12\u001b[0m                - 3.0 * np.dot(vec[o1], r12) * np.dot(vec[o2], r12) / r12_length**5 )\n\u001b[1;32m     13\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mans\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "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
}