{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Strain Derivatives (strain_derivative)\n",
    "## Background and Basic Usage\n",
    "\n",
    "Let us first unpack some data that we will need below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# unpack data if not done yet\n",
    "import numpy as np\n",
    "import os\n",
    "if not os.path.exists(\"data/nacl\"):\n",
    "  import tarfile\n",
    "  tar = tarfile.open(\"data/nacl.tgz\")\n",
    "  tar.extractall(path='data/')\n",
    "  tar.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we explore the `strain_derivative` class, which can be used to compute both linear and nonlinear elastic constants in crystals. This class does not do that much work: it mainly instantiates and uses several of our core classes. We begin with the generic considerations of how the class works. We will begin working with NaCl, given that it will be treated as an example below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us begin by importing and instantiating the class. The constructor is listed for clarity.\n",
    "\n",
    "```python\n",
    "def __init__(self,xtal,pgn,vname,order,target,fdtype='central',sim_trans=None,sim_trans_basis='standard',lattice_dim=3):\n",
    "```\n",
    "\n",
    "The variable *xtal* defines the crystal structure (a POSCAR for now). You will recognize familiar variable names from the `finite_difference` class. We will explore the variables `\"xx,yy\"` at order `(1,1)`: a second derivative. Please note that you must specify a *target*, which is either *energy* or *stress*; and this is what you will be taking the derivative of. \n",
    "\n",
    "At first, we will use the class in the simplest way possible: we will evaluate a specified energy derivative, without any considerations of group theory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from strain_derivative import strain_derivative\n",
    "sd=strain_derivative('data/nacl/c450_k20/static/POSCAR','Oh',['xx','yy'],[1,1],'energy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The constructor will call the representation class which will symmetrize the strain tensor according to the provide point group: a variable called `strain_rep` will store this object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['dim', 'G', 'irrep_matr', 'chi', 'irrvec_dict', 'rank', 'rdim', 'pg', 'irrep_counter']\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "OrderedDict([('A1g', 1), ('Eg', 1), ('T2g', 1)])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print vars(sd.strain_rep).keys()\n",
    "sd.strain_rep.irrep_counter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The strain matrices are also stored in the ordereddict `sd.strain_matr`. We print them out here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{align}\\epsilon_{xx}=\\begin{bmatrix}\n",
       "  1 & 0 & 0 \\\\\n",
       "  0 & 0 & 0 \\\\\n",
       "  0 & 0 & 0 \n",
       "\\end{bmatrix}&&\\epsilon_{yy}=\\begin{bmatrix}\n",
       "  0 & 0 & 0 \\\\\n",
       "  0 & 1 & 0 \\\\\n",
       "  0 & 0 & 0 \n",
       "\\end{bmatrix}&&\\epsilon_{zz}=\\begin{bmatrix}\n",
       "  0 & 0 & 0 \\\\\n",
       "  0 & 0 & 0 \\\\\n",
       "  0 & 0 & 1 \n",
       "\\end{bmatrix}&&\\epsilon_{xy}=\\begin{bmatrix}\n",
       "    0 & \\frac{1}{2} & 0 \\\\\n",
       "  \\frac{1}{2} &   0 & 0 \\\\\n",
       "    0 &   0 & 0 \n",
       "\\end{bmatrix}&&\\epsilon_{xz}=\\begin{bmatrix}\n",
       "    0 & 0 & \\frac{1}{2} \\\\\n",
       "    0 & 0 &   0 \\\\\n",
       "  \\frac{1}{2} & 0 &   0 \n",
       "\\end{bmatrix}&&\\epsilon_{yz}=\\begin{bmatrix}\n",
       "  0 &   0 &   0 \\\\\n",
       "  0 &   0 & \\frac{1}{2} \\\\\n",
       "  0 & \\frac{1}{2} &   0 \n",
       "\\end{bmatrix}\\end{align}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, Math,Markdown\n",
    "from coordinate_tools import table_text as ttext\n",
    "display(Math(r\"\\begin{align}\"+\"&&\".join([r\"\\epsilon_{{{}}}=\".format(s)+ttext(ss,cdelim='&',prec=None).latex().replace(\"0.5\",r\"\\frac{1}{2}\") for s,ss in sd.strain_matr.items()])+\"\\end{align}\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can review my course notes for more details on the above. Other important properties defined in the constructor are the standard basis (sbasis), irreducible basis (ibasis), and the transformed basis (tbasis):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['xx', 'yy', 'zz', 'xy', 'xz', 'yz']\n",
      "('A1g', 'Eg.0', 'Eg.1', 'T2g.0', 'T2g.1', 'T2g.2')\n",
      "['v0', 'v1', 'v2', 'v3', 'v4', 'v5']\n"
     ]
    }
   ],
   "source": [
    "print sd.sbasis\n",
    "print sd.ibasis\n",
    "print sd.tbasis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For multidimensional irreps, the row will always be denoted with a right period followed by an integer: `Eg.0`. If there are multiple instances of any irreps, these will be distinguished by a left period preceded by an integer: `1.A1g`, etc. The transformed basis has been rotated by a matrix provided to the constructor, and this can be used to perform BID; but I will leave this for another part of the docs. For now we will either use the standard basis or the irreducible basis. The irreducible strain vectors can be seen in the standard variable of the representation instance (this class will be documented later)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.57735027  0.57735027  0.57735027  0.          0.          0.        ]\n",
      " [ 0.70710678 -0.70710678  0.          0.          0.          0.        ]\n",
      " [-0.40824829 -0.40824829  0.81649658  0.          0.          0.        ]\n",
      " [ 0.          0.          0.          1.          0.          0.        ]\n",
      " [ 0.          0.          0.          0.          0.          1.        ]\n",
      " [ 0.          0.          0.          0.          1.          0.        ]]\n"
     ]
    }
   ],
   "source": [
    "print np.array(sd.strain_rep.irrvec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a variety of methods which help you impart or measure strains, create directories, etc. However, I will only outline the higher level methods that bring everything together. In order to setup a calculation, the high level method to create all of the strain directories is:\n",
    "\n",
    "```python\n",
    "def generate_fd_runs(self,delta,dname=None,dprefix='run',fpinp=None,zerodisp=False):\n",
    "```\n",
    "\n",
    "You must supply the delta, which has the same conventions as the `finite_difference` class, and everything else is optional. The variable *dname* provides a name for the overall directory that will hold the runs, and a default is chosen is nothing is provided; *dprefix* is the prefix for the subdirectories; *fpinp* is the location of the input files of the first-principles method which will be copied to every directory (e.g. INCAR, KPOINTS, etc if using VASP); *zerodisp* declares whether or not undistorted directories will be created.\n",
    "\n",
    "The next major method is the one that collects all the data after the runs finish. \n",
    "```python\n",
    "def gather_all_fp_data(self,fp_dirs,unstrained=None):\n",
    "```\n",
    "where *fp_dirs* is a list of first-principles directories that contain the data; and *unstrained* is a path to the unstrained data.\n",
    "\n",
    "The next method simply instantiates the finite_difference class and feeds it the data.\n",
    "```python\n",
    "def finite_difference(self,delta,use_available=False): \n",
    "```\n",
    "where these are standarad finite difference variables.\n",
    "\n",
    "The next method instantiates the fit_quadratic class and feeds it the error tails.\n",
    "```python\n",
    "def fit_quadratic(self,dmin=3,dmax=None,sigi=None,stype='consec'):\n",
    "```\n",
    "where these are standard fit_quadratic variables.\n",
    "\n",
    "\n",
    "Finally, this method outputs the data to relevant files.\n",
    "```python\n",
    "def output_data(self,sigi=None,units=None):\n",
    "```\n",
    "\n",
    "These are the basic methods that are called in the main program."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples\n",
    "\n",
    "The above is best illustrated using a series of examples, and we will explore several:\n",
    "\n",
    "1. Linear elastic constants of NaCl.\n",
    "1. Linear and nonlinear elastic constants of ThO2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NaCl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us start by creating the necassary directories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/cam1/git_projects/principia_materia_doc_v0/strain/sample_mkruns\n"
     ]
    }
   ],
   "source": [
    "%cd sample_mkruns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are starting in this directory with a directory containing the first-principles static calculation in a directory called `static`, which contains the static calculation of the undistorted structure, in addition to a yaml input file for strain_derivative.py."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;36minputs_mkruns.txt\u001b[0m@  \u001b[01;36mstatic\u001b[0m@\r\n"
     ]
    }
   ],
   "source": [
    "%ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delta: 0.005,0.03,0.005\r\n",
      "fpinp: static/\r\n",
      "order: 1\r\n",
      "pgname: Oh\r\n",
      "schema: mkruns\r\n",
      "unstrained_path: static/\r\n",
      "vname: A1g\r\n",
      "xtal: static/POSCAR\r\n"
     ]
    }
   ],
   "source": [
    "!more inputs_mkruns.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the variables should be clear. \n",
    "\n",
    "*fpinp* - path to the input files for DFT.\n",
    "*schema* - mkruns tells us to generate the runs.\n",
    "\n",
    "The above directory indicates that we will setup runs for the first derivative wrt to the A1g strain. However, please note that both first and second derivatives can be obtained from the same stencil. Now let us execute the code, and we can see all the runs that it has created. You will see that the highest level directory is called *strain_A1g*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "!strain_derivative.py -i inputs_mkruns.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;36minputs_mkruns.txt\u001b[0m@  \u001b[01;36mstatic\u001b[0m@  \u001b[01;34mstrain_A1g\u001b[0m/\r\n"
     ]
    }
   ],
   "source": [
    "%ls "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;34mdelta_0.5000\u001b[0m/  \u001b[01;34mdelta_1.0000\u001b[0m/  \u001b[01;34mdelta_1.5000\u001b[0m/  \u001b[01;34mdelta_2.0000\u001b[0m/  \u001b[01;34mdelta_2.5000\u001b[0m/\r\n"
     ]
    }
   ],
   "source": [
    "%ls strain_A1g/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;34mrun_m1\u001b[0m/  \u001b[01;34mrun_p1\u001b[0m/\r\n"
     ]
    }
   ],
   "source": [
    "%ls strain_A1g/delta_0.5000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INCAR  info_strain.yaml  \u001b[0m\u001b[01;32mKPOINTS\u001b[0m*  POSCAR  POTCAR\r\n"
     ]
    }
   ],
   "source": [
    "%ls strain_A1g/delta_0.5000/run_m1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INCAR  info_strain.yaml  \u001b[0m\u001b[01;32mKPOINTS\u001b[0m*  POSCAR  POTCAR  \u001b[40;31;01mWAVECAR\u001b[0m@\r\n"
     ]
    }
   ],
   "source": [
    "%ls strain_A1g/delta_0.5000/run_p1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how it has automatically created WAVECAR links in case you will run sequentially (XX add flag for this). Now let us clean up and look at some results. I have listed where you can find the runs in my dft_runs directory as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "%rm -r strain_A1g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/cam1/git_projects/principia_materia_doc_v0/strain/data/nacl/c450_k20\n"
     ]
    }
   ],
   "source": [
    "%cd ../data/nacl/c450_k20/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/cam1/dft_runs/nacl/pbe/c450_k20/\r\n"
     ]
    }
   ],
   "source": [
    "%cat ../readme.txt # this is where the data was copied from"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you can see that we have run a bunch of different strain derivatives. These are all either first or second derivatives. Please note that there is one cross derivative, denoted as *strain_xx_yy*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inputs_energy_2.txt  \u001b[0m\u001b[01;34mstatic\u001b[0m/       \u001b[01;34mstrain_T2g.0\u001b[0m/  \u001b[01;35mUntitled.png\u001b[0m\r\n",
      "inputs_mkruns.txt    \u001b[01;34mstrain_A1g\u001b[0m/   \u001b[01;34mstrain_xx\u001b[0m/\r\n",
      "inputs_stress_1.txt  \u001b[01;34mstrain_Eg.0\u001b[0m/  \u001b[01;34mstrain_xx_yy\u001b[0m/\r\n"
     ]
    }
   ],
   "source": [
    "%ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is our input file for processing these results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "order: 2\r\n",
      "pgname: Oh\r\n",
      "target: energy\r\n",
      "unstrained_path: static/\r\n",
      "xtal: static/POSCAR\r\n",
      "scale_delta: 0.5\r\n"
     ]
    }
   ],
   "source": [
    "%cat inputs_energy_2.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we use the above input file and we add one more variable on the command line which tells us which data set to use. Please note that we have decided to take second order energy derivatives, and thus we must rescale the deltas using *scale_delta* by a factor of 0.5; because the grid was created for a first derivative, not a second derivative. Now let us run the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "!strain_derivative.py -i inputs_energy_2.txt -dpath strain_A1g\n",
    "!strain_derivative.py -i inputs_energy_2.txt -dpath strain_Eg.0\n",
    "!strain_derivative.py -i inputs_energy_2.txt -dpath strain_T2g.0\n",
    "!strain_derivative.py -i inputs_energy_2.txt -dpath strain_xx\n",
    "!strain_derivative.py -i inputs_energy_2.txt -dpath strain_xx_yy -order 1,1 --scale_delta 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you can see the output file. The notiation *-2* indicates that this is a second derivative. The first line is the output of our fit_quadratic function, and we are using default configurations, though you can pass configurations through argparse (see code)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# 20.5300403857 340.287852585 0.00260950832044 [2, 3, 4]\r\n",
      "#gy s=[0,o] l=[1,0] lw=2 lc=bu sfp=1 sfc=re sc=re ss=1.6 xr=0,0.012625 yr=20.26134,20.80196 sts=2 st=c\\z{0.6}\\v{-0.5}1\\N=20.530\r\n",
      "#gf 0,0.0125,-100,20.5300403857+340.287852585*x^2\r\n",
      "0.00250000 20.59600000\r\n",
      "0.00500000 20.46600000\r\n"
     ]
    }
   ],
   "source": [
    "!head -5 strain_A1g/Y_energy_A1g-2.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will want to compare our results to published results and experiment. Let us start with experiment. \n",
    "\n",
    "```\n",
    "Elastic constants of alkali halides at 4.2 degrees k \n",
    "J. T. Lewis A. Lehoczky C. V. Briscoe \n",
    "Physical Review 161 877 1967 \n",
    "```\n",
    "\n",
    "They provide rather obscure units of $10^{11}$ dyne/cm$^2$, so we will use our units module to convert, and we will see a factor of 10 brings us to GPa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.0"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from unit_conversion import pressure as units_pressure\n",
    "1E11*units_pressure(\"dyne/cm^2\",\"GPa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can see their results at T=4K in units of GPa (Table 4).\n",
    "\n",
    "|c$_{11}$|c$_{12}$|c$_{44}$|\n",
    "|---|---|---|\n",
    "|57.33 | 11.23 | 13.31|\n",
    "\n",
    "We can then rotate these results to the irreducible basis (see my class notes).\n",
    "\n",
    "|c$_{A_{1g}}$|c$_{E_{g}}$|c$_{T_{2g}}$|\n",
    "|---|---|---|\n",
    "| 79.79| 46.1 | 13.31|\n",
    "\n",
    "Perhaps it will be easiest to convert the above units into eV/(unit cell) by using the conversion between GPa and eV/A^3 and then multiplying by the volume. Let us do this.\n",
    "\n",
    "```python\n",
    "from periodica import structure\n",
    "xtal=structure('static/POSCAR')\n",
    "conv=xtal.vol*units_pressure(\"GPa\",'eV/A^3')\n",
    "print \"|\".join([\"{:.2f}\".format(conv*s) for s in [79.79,46.1, 13.31]])\n",
    "# 22.50|13.00|3.75\n",
    "print \"|\".join([\"{:.2f}\".format(conv*s) for s in [57.33,11.23,13.31]])\n",
    "# 16.17|3.17|3.75\n",
    "```\n",
    "\n",
    "So in units of eV/unit-cell, we have:\n",
    "\n",
    "|c$_{11}$|c$_{12}$|c$_{44}$|\n",
    "|---|---|---|\n",
    "| 16.17| 3.17 |3.75 |\n",
    "\n",
    "|c$_{A_{1g}}$|c$_{E_{g}}$|c$_{T_{2g}}$|\n",
    "|---|---|---|\n",
    "|22.50|13.00|3.75|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now we know what the answers should roughly be. Let us first look at the irreducible measurements, then we can look at the direct measurement of $c_{12}$.\n",
    "\n",
    "We could directly rangle the raw data, but the comments in the file are actually gracey commands such that we can easily interactively plot this. Our notation for the value of the intercept, inconveniently, is $c_1$ in all cases; not to be confused with a given elastic constant. You can see that the agreement is pretty rough: 20-30%. \n",
    "\n",
    "|c$_{A_{1g}}$|c$_{E_{g}}$|c$_{T_{2g}}$|\n",
    "|---|---|---|\n",
    "|20.53|10.294|3.633|\n",
    "\n",
    "If we rerun things with LDA and the same parameters (not shown, but check my dft_runs directory):\n",
    "\n",
    "|c$_{A_{1g}}$|c$_{E_{g}}$|c$_{T_{2g}}$|\n",
    "|---|---|---|\n",
    "|25.26|15.12|3.39|\n",
    "\n",
    "Somewhat reassuringly, LDA and PBE bound the results in all cases except for the T2g modes (shear). \n",
    "\n",
    "Let us look at all of the raw error tails."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xmgrace -nosafe -noask -hardcopy -hdevice PNG -batch /home/cam1/.ppbatchfile\r\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABXgAAAK8CAMAAABiCycpAAAAQlBMVEX///8AAAD/AAAA/wAAAP///wC8j4/c3NyUANMA////AP//pQByIbxnB0hA4NAAiwDAwMCBgYFCQkLv4+Pfx8fPq6vxjVQBAAAACXBIWXMAAAsTAAALEwEAmpwYAAAADnRFWHRUaXRsZQBVbnRpdGxlZPPPpzIAAAALdEVYdEF1dGhvcgBjYW0xmdSm6wAAABV0RVh0U29mdHdhcmUAR3JhY2UtNS4xLjI1jnpa/wAAIABJREFUeF7t3d2WosqyBtCyxxnrGuzd7/+qR/4hwSqTEChwzotuSSF3dOzyW1mI+PUFAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAARygfihfG5u7Dw+qAYWsY+3kSgE9T3mr3cUI+xh6ZOR0btEdUB7UjxX263e32zSQAH6uNzIchIO9NghbjsZH+iO7poh8YJe8PkwB8qnu9LG3StwvIx1bzoFwMzVnMFrd7dZahXvb2pxt+mATgU5VtUtah2edoH5/JyYPGPR3rB6pJ2oc/TQLwqYZMHJa8/Vp1ebVaDKvadqCfo1oLtw9/mATgUxV9OtbvmNWJWoxWqOPHnftwTGN04UIfvD9NAvCpxleAdcFbXYwwHuwfN9oTusv6iP1hEoCPNY7QLnjH75BNN2rNtWTLl+cOVzD8MAnAddTXJ3yzJP1Gd463+7tWXfTQb9Tq3K3+Rxait+zPKfwwCcBVFPcqc6u3uBZS8Se35oTA6A2yr3p9Oz0/O3x4Yn7idnjb7YdJAK6iaPOtCsfs5O3OylYHT0bn52fLNn2TNC2GVfALkwBcQDG6pGD2u327Sp2Y7NBd//VaZrZnekcjzaeGh/fWXpgE4OzufRAunFTt03Yk2aE55NXMrHO2X1d3N2sYrbn7Pb+ZBODURjdFKJZvFJaaPJ+9WJ18QLgyWgW/PAnAmZWRdOtPU2RkZhWz05F62Vs9eH0SgBMLpdu9Xy2/npkLz1XJW030+iQAJxZJt3J0tjbNzKdXgk33bNyacw2vTwJwYgsx+KpifPVZmpmz08W9hf/F7m291ycBOK92sbnCJHen83w76ULwdhn7+iQA53X79hf624L2qWnuTtenw14LFvK0bC+teH0SgPOq3sIa3R9h9Ewlydxa88wkd4v6/Oxrd3QcXb/WK9tZX54E4MSqN7T6gJvdJye9hrdSPzHJ3fozv30kDwvYRcMtcQb9QvfVSQDOrL4xWR1xxfDtPMPXQzwxjujmsNFpgiHKu4/FDfeDXFrwdneLeDIJwMVUS96H+qsr20isHk53SnUf9W2P7caaJ4dsrc5jVE+Ozmf0l/4OGTvO4oVJAC6nTd4hG4ti+ETasmnudlE5/2b27vPBVfDWy+pyuBVZPUfz5GhpuzAJwPW0t6oZ/2b/ffAmudufdagWsY/V8nCz8yrTq43hZjhDmvZ5n9wbfTYJwCV175l1vg/e56o338aJWXTTFsP7cr16aCFfnwwDXNra4AVgJcELsDPBC7AzwQuwM8ELsDPBC7AzwQuwM8ELsDPBC7AzwQuwq+ruZHef2gUAAAAAAADgUM2718Pm/OsAqtH6+1tmwwDka7+PZfhG2KL6Qpbpd7sW9ffE+gYsgHe4V6vb6tuwmqhtv/61nEZsP+qKeoCool3a3trvI+yidRKxo9HkC7MAyNWdza2+1bX+q03W/kH7ZPNA8AKE9edym0ztk3VyVqFNZcEL8E5d8Lab/dmFdqM+EVEm77kBEFAH7WidO1ncVpeYFbNrHQAIaBazz4K3vuZsfNq3uQQNIMsoVeguGHsavPXZhtHHKrZp34lm3WTSM826yaRnmnWTSa8/66qDrqtoYvZ58Dafs+iTd5v/ir1rnqlNZt1k0jPNusmkZ5p1k0mvOOs2WXER9+bs7dPgLR6ZW0Vvl7zbtO9Es24y6Zlm3WTSM826yaTXn3XVQZdVtu+aPQve5m21UfJu074TzbrJpGeadZNJzzTrJpNef9ZVB11V9/G18UVkk98J2sDtP1i8UftONOsmk55p1k0mPdOsm0x6/VlXHXRRfe5WCfvkAxTd092D4ck3OtGsm0x6plk3mfRMs24y6fVnXXXQNVXnbzv1Xcgqyx8ZHhJ4ePKNTjTrJpOeadZNJj3TrJtMev1ZVx10SUPulqPPpk3uALlX8G7jRLUqdRMnqvVEpa6rddVBVzTkbn0z9DZwuw+ptfeD7M5A3Lc9x7uNE9Wq1E2cqNYTlbqu1lUHXVD1wYjxFXZN4nZvspXtu2ntLdCH997O1L4T1arUTZyo1hOVuq7WVQddUPv9E5VmUfsI1/u9+5BaMR4t++Fzte9EtSp1Eyeq9USlrqt11UGfoSjHV/CORkfDZ2rfiWpV6iZOVOuJSl1X66qD6JypfSeqVambOFGtJyp1Xa2rDqJzpvadqFalbuJEtZ6o1HW1rjqIjvYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHyK2VjgMsEBnvoH1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon18/f2bjnSeP8OHkxwh2sfXv3/pSOPvv//+l45BTXKEaB9///tvaWH7iN3/BC9PSI4Q7eMRsPMlbx27gpdnJEevvN9u93LYvD02i9Hz1digGdA+qoBNx77+1ithwcsTkqPziN3Kvd98hG7RbzaaXWrtwORpPtD/nq5sn42D5Ojcq9VtUS16m80qd6tl73jNW9669W4XyNr38f6rTyqkoxXBy1OSo1F0SXqro7bol7TjBvUp/EjgdqR/js/097//VWvepYR9MgySo9OdzS1vdab2ydo/eCj6h/eubdr36f799+Qsr+DlG5Kj0Z/LbYK3+fOrDuLumZH+1O/Sk3yQv1W2VkvehSvKZsH7v+qsxD9hjOSY6YK33ezPOYwNy+CFJ/kk/6sSt7p+YX5F2Sx4//33729zde/DQlDzQSRHog7a0Tq3X/uO9GcatO/TNYFbpek8SZPg/Tfsm0YyH0dyTDXXK/wQvMNFZtr32eoFb3OuYZ6k08H/ddm8uC8fRnJMNZeRfR+8ozfcRhf2Dhf38jG6N9WqNez0mWZwlLD9Hv9bPC/B1cmKbxRNpn4fvMOZBv/d+mz1W2uVxSXvk+CtzggPw3wmyTFxb04ifB+8o4+zad9Hq64laywteZ8E7+gRH0tyjJVtpH4bvONLe7Xvk/0dzhks3YpsFrztluBFcox1H18bX0Q2PxkzOtOgfR+tvk3DID1zOw3ef6PgTSOajyM5Bn3ujta5Cx+gGN84Z/YkH+S/f9UVuY0qeZMryqbBW53arXf4a8GL5BgUoxvi3Jc+Mrwwon0f7H/J1WLpkncavF3y/v33b37JL59GcnSG3C27y3m/qqVvckve7r5lDe37YMNba5VqyTvengVvnbw+MkxNcrSG3K1vht4GbndD3lHaTu9XNnrMZ+mvJWvM315LB6x06UmORpFc2twkbvcmWzncIH167kH7Pte/6Tndej07HkiD11KXgeRotN8/UWmC9ZG593v33T/FcFHZ5EyD9n2u2XtkVfBOsnW63X66GCqS46mi+2K1emN4OKZ9H2v23e31hQ3jgWnw1qd3awIYyRGjfZ9qfkOy+rZjo/O49c13k6dbzvZ+PMkRon2fqfvoRL+iHT5L0Qw1F/bWi9xul3HyOu3w6SRHiPbxqn/1mYY2ftNPufFhJEeI9vGav5OrfqeXAPN5JEeI9vGa9IqH8QafR3KEaB8v+ZdcfOZUw4eTHCHax0uS99Nc1vDpJEeI9vGS8QeK//osBZIjRPt4SXODnJb1LpIjRPt4UX+lr9hFcgRpHy/76/PCdCRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHyK2VjgMsEBnvoH1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APskRon1APsnRKe+3270cNm+PzWL0fK8o7/1+2gfkkxytR+xW7v3mI3SLfnNQ3G9DOl+0fX/SAeCtrpkc+e7V6raoFr3NZpW71bI3XfM+FsLjoeu1708nfQJ4m+slxypFu7S9Nbla3Nq+dH937snA1drXx67ohQ1dLTlW6s7mPha01YmEsjud0D/on54ugS/WvknuSl7YysWSY63+XG4TvM2fX3XSds981QvhSQ5frX1J7kpe2Mi1kiOuC952sz/nUEtPNFysfWns/pG8sI1LJccb1NE6WueO17hF9dZbUZZXvaohTd0/ghe2cankiCvrqxqeBO9j+N5cddaf6L1S+9LQraU7AW9wpeR4g+YysifBe+/eeRuSt7n4t9fvekZp5tbSnYBVrpQV71Y0MfskeLvHxdC2K7UvzdxauhPwBldKjrh7c3XDD8FbLX3bJe+V2pdmbi3dCXiDKyVHWNleVfY0eNsH7dW+1VD/5OmlkdtKdwPiLpQcYd3H18YXkY3PxYyv7r1e8C4nb7oT8AZXSo6gPnenEdsPfvW3xxG8QMSVkiOmGH0aeBSxw5mGYWP45PCV2pdmbi3dCVhr9HK6UnKEDLlbVsna3zRndG+G/gzEsA6+UvvSzK2lOwFrJC+pKyVHxJC79U3O28Dtbsjb3iWyWwgPt+S9VPtGedtJdwFWmL2qLpUc61WX5o6vbG4St1viVh9Za/dqbho53FOne3AFkx+ORroLkG/+srpUcqzXfv9EpVnMVvdluHff/TPclewxVN5H30txrfYlPx5yF94hfV39uVpyvFMxvhnOcKZ3co+cq7Vv9uMBRCUvq8rVkmNvV2vf9IcDiJu8qhqXS46dXa99w48G8AajvB1cLzl2dc32CV14mzRza9dMjt1oH/DE/9XSzK1JjhDtA8aatB1JM7cmOUK0D5inbaN6Ko3chuQI0T7YwDneZkhztpHulWZuTXKEaB+8WR9O6RO/Q5qzjXSvkVHcDiRHiPbBW03iKX3yMGnOttLdFk3+RR3JEaJ98E5JPqVP7yqN2Va620+Sf1JDcoRoH7xRmk8HJG+as410rxzpv+mPT65FaR+8T5pPfzZN3jRcZ9IDVkr/SX8Eb5T2wfuk+fRng+BN03UmPSAu/Te5O1mU9sHbpPlUS3daJQ3XRrrXZub/JMkRon3wNkk+NdKdXpfmbCvdbQ+zf5HkCNE+eJtJPHXSnX6Qxmwr3W13yb9HcoRoH7zNKG4H6U6DNFxn0gOONvq3SI4Q7YN3SSO3Nd0pDdeZ6e6/leQI0T54mzRya/Uzabo2ksPPRHKEaB+8TZq5tQul7YjkCNE+iGszNc3c2pXidiA5QrQP1uoXsZ00c2vpUdcgOUK0D7KkYdton0xD989Vc1dyxGgfPJOm60x6gODlNdoHE2m4zqQHTKSxe9XclRwx2sfBjoumNFHn0iNe8CG5KzlitI8DHZFPabjOpAdkGqXurv+unUmOEO3jMDtGVBqujXSvN9np33QsyRGifRxllLrvTKk0XWfSAzbxtn/ObyU5QrSPgyS5uyJ500T9TnosQZIjRPs4Rhq7f15M3jRR59Ij2ILkCNE+jpGm7p958KaJOpMewH4kR4j2cYg0dGvVE2m4zqQzcQjJEXJrpeOwqTRzazL2BETGO2gfe0gjdfl2MvUz6aH8SpIjRPt4ryRgn0kzt5bOxe8lOUK0j9XSMP1Ocmgaua1kL34vyRGifbwgjdHvpMcuSyO3lu7E7yU5QrSPqTRHv5MemyHN3Fq6E7+X5AjRvs+Vxuh30mPD0sytpTvxe0mOEO37CGmOfic9dhtp5tbSnfi9JEeI9l1NmqPfSA/dVRq6f+TuqUiOEO07mTQ8X5ROc7w0df8I3lORHCHa93ul4fmidJrfKo1duXsqkiNE+46TJuZr0lnOS+6emeQI0b7tpJGZLZ3wauTuiUmOEO3LlIZjRDr35xG7pyU5QrSvkibim6T/MywRuqckOUIu3L40Bt8k/Z+BT3Th5NjDL25fmnhbSv+3gW/94uQ4gf+7pRF0Iek/FngXwRvxq4M3LRb4LQRviPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHiPYB+SRHr7zfbvdy2Lw9NovR8+1w495saR+QT3J0HrFbaRP1sfkI3aLf7BTNXtWTFe0D8kmO1r1a3RbVorfZbKK17BK2c29XvO2m9gH5JEej6M4dNIvZx8K22+53qRRJv7QPyCc5Gt3Z3PJ2q1azZf3n+EFjdA64pn1APsnR6M/lNsHb/PlVB3H3zFe14O1PMjS0D8gnORJd8Lab/TmHWnWlw/jKB+0DVpAciTpoR+vcfu3bbDRXPvTvuGkfkE9yTJX1VQ3Pgverudq3v5pM+4AVJMdUcxnZN8HbXsrbPtY+IJ/kmCiamP02eOvPWrSD7bmHznQ3gJas+Ma9ubrhh+Ct1rzNI+0D8kmOse4WDD8Eb7XkbR5oH5BPcox0H18bX0S2+DtBH8xLTwJ8T3IM+twdrXOnH6DoCF4gQHL0itENce5PPjLc6ke1D8gnOTpD7pbd5bxf1dJ3fkve0fkH7QPySY7WkLv1R4LbwO1uyNtc3lu069xhX+0D8kmORneD8+4KuyZxuzfZyuY2vbfm1O/onIT2AfkkR6P9/olKs6x9ZO793t2VoWhGq53upZvkADGS46lifAvI7na9ZTk556t9QD7JEaJ9QD7JEaJ9QD7JEaJ9QD7JEaJ9QD7JEaJ9QD7JEaJ9QD7JEaJ9QD7JEaJ9QD7JEaJ9J/InHYCjSI4Q7TuJP530CTiC5AjRvlPoY1f08jtIjhDtO4NJ7kpefgHJEaJ9J5DkruTleJIjRPt+vzR2/0heDic5QrTv90tT94/g5XCSI0T7fr00dGvpTrAvyRGifb9emrm1dCfYl+QI0b5fL83cWroT7EtyhGjfr5dmbi3dCfYlOUK077dLI7eV7ga7khwh2vfrpZFbS3eCfUmOEO379dLMraU7wb4kR4j2/Xpp5tbSnWBfkiNE+369NHNr6U6wL8kRon2/Xxq6f+Quh5McIdr3+6Wp+0fwcjjJEaJ9J5DGrtzlcJIjRPvOQO7y20iOEO07BbnLLyM5QrTvJMQuv4rkCNG+ExG6/BqSI0T7gHySI0T7gHySI0T7gHySI+TWSscBFoiMd9A+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/kCNE+IJ/k6JT32+1eDpu3x2Yxen5wH5qmfUA+ydF6xG7l3m8+QrfoN8cee/aPtQ/IJzka92p1W1SL3mazyt1q2Ttf8xY3wQuESI5a0S5tH6FaRW3RZesoYzs3K14gRnLUurO55e1WneYt6z/HD3r3shS8QIjkqPXncpvgbf78qoO4e6bxWBoLXiBGckx1wdtu9ucc+u1pGGsfkE9yTNWhOorWfu3bqE5JCF4gRnJMlPVVDU+Dt0ye1T5gBckx0VxG9ix4qxMNSfBO9eMAY7LiuaKJ2WfB21z7YMULxEiOsXtzdcOT4G0/UCx4gRjJMVKfwq3+Xgze7kMWgheIkRyDLlnHF5GNQrZcOEWjfUA+ydHrc3e0zh2vbgUv8B6So1OMbohzf/6RYacagCjJ0Rpyt6yytb9pzvz2ZIIXiJEcjSF362sX2sDtbsjb3iWyIXiBGMlRq26yOz572yRu9yZbOdwgvdnqH2sfkE9y1Nrvn6g0J3UfmXu/d9/9U0w+RiF4gRjJ8UxRjj8sPDyc0D4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QI0T4gn+QIubXScYAFIuMdtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzlCtA/IJzl65f12u5fD5u2xWYyerxWPnW7DTtoH5JMcnSpRq6ztNx+hW/Sb/ehkJ+0DVpAcrXu1uq3Ws02o1rlbLXsna95yupP2AWtIjkbRRuljOVtFbXFrG9P93Rh2age0D8gnORrd2dyyOYNbdudx+wf1Rve34AUCJEdjOGtbJ23/BtqQsSOlUw1AgORIdMHbbvbnHMbaE8Bf2gesITkSddCO1rnji8daxXCVmfYB+STHVHMW4dvgHZ/21T4gn+SYas4ifBO81ecqmisf6mcnzwG8QnJMFE3MPg/eoizr6O2fnRjvCdCTFd+4N5crPA/e2r0f1D4gn+QYK9vLxH4I3qL/D5b2Afkkx0j38bXxRWSLvxPcBS+wnuQY9Lk7WucufoBiGF16EuB7kqNXjG6Ic1/8yHCv/+ia9gH5JEdnyN1ynKzT25O1+o+uaR+QT3K0htytb4beBm53Q97uLpHNHu1FZ9Vu7d8Ar5McjepChfEVdk3idm+ylc0deMv2npHD91RoH5BPcjTar5aoNKn6yNz7vfvun6IZrdO5rG+Z3tI+IJ/keKooR++rdbfrLcvJOV/tA/JJjhDtA/JJjhDtA/JJjhDtA/JJjhDtA/JJjhDtA/JJjhDtA/JJjhDtA/JJjhDtA/JJjhDtA/JJjhDtA/JJjpAzte9EtSp1Eyeq9USlrqt11UF0ztS+E9Wq1E2cqNYTlbqu1lUH0TlT+05Uq1I3caJaT1TqulpXHUTnTO07Ua1K3cSJaj1RqetqXXUQnTO170S1KnUTJ6r1RKWuq3XVQXTO1L4T1arUTZyo1hOVuq7WVQfROVP7TlSrUjdxolpPVOq6WlcdRGeb9p1o1k0mPdOsm0x6plk3mfT6s646iM427TvRrJtMeqZZN5n0TLNuMun1Z111EJ1t2neiWTeZ9EyzbjLpmWbdZNLrz7rqIDrbtO9Es24y6Zlm3WTSM826yaTXn3XVQXS2ad+JZt1k0jPNusmkZ5p1k0mvP+uqg+hs074TzbrJpGeadZNJzzTrJpNef9ZVB9HZpn0nmnWTSc806yaTnmnWTSa9/qyrDqJzA8iXRgkAAAAAAAAAAACQpbjfbrcyHT3SckXJ6HSzKKuton92LytKbYbmB21uoYz56MJOZbl7ZxeqmI8u/ATcZ8fs56Wa65Hx1q7mxSyPjko86oV1fcWt6vL99nt6u1xRMjrdLG6Nvf8RK0pth+Y//1tbKGM+Ot+pvN337upCFQuj082mzHs9doiXaq7qPO5TB7NilkfHJR71wrq+pu1fX8f9xKaWK0pGp5vN+rH6gdn3B2RFqY37/sG7VMZsdLZTcUSUzapYGk1/Atqt3fvaeqnmRzsPDN60mOXRSYlHvbA+wL39SS2P+omdWa4oGZ1utmuyxw/IvimxotTm4QHBu1DGfDTd6YiV+byKxdHkJ6DfOijWXqr5qziwwlkxy6OTEo96YX2AvseH/TyklitKRieb/SmpvVcT+aXWitsBK7N5GQuj6U7HvODSKhZHp5tdPw+LtZdqrhxW4VIxy6N9iYe9sK5v+DHYPwqWLVeUjE43+3P/O/9Qryi19lhH7N7thTLmo+lO93372UqrWBxNNrtk2PknoPdSzcnYzpaKWR7tx456YX2A4deO7leOoy1XlIwu77T3z8dyFT+WWr31fkDwpmUsjCabx5xoSKtYHk027+3a/KiTkS/V3I7t+jM6WCpmeXRe4nyEmIX/zB1suaJkdHmn0U/RLpar+KnU+ve33SNtVsbSaLJ5r2KsLPe+mGhNqY/tKnnTd+x3k5bzfPTACIkDAAADlElEQVSw19lSMcuj8xJ3fmF9gKW+H2u5oh9edq2d1zvLVfxUal3jMLqTWRlLo7PNW1lfTLTvid41pVar89v9vv+Vb620nOejh73OlopZHp2XuPML6wMs9f1YyxX98LJrdFfG7GW5ih9Kba7xH0Z3kpaxODrdfGRZFWTV9fW79nVFqZXZxwD2NC/n2ehhr7OlYpZHZyXu/cL6AEt9P9ZyRT++7Cr3nf8Fy1V8X2r7RvHuGbGi1n7HnS9+W1Fqpb7W/6iAmJfzbPSw19lSMcujsxL3fmF9gKW+H2u5oh9fdl/VK2/n34cWq0hHp5tFu9++UfaVlvFkdLrZ71gtfdvxPawo9eF+b043tKM7m5XzdPSw19lSMcujaYm7v7A+wFLfj7Vc0U8vu3qfvX88FqtIR6eb3UnIYXQnK2od3lHZ9zLOFaW2b6sdl7xpOc9HD3udLRWzPJqWuPsL6wOMryHZOQqeWK4oGV3aaf83VpaqmI1ON+v3qjp7NnxFrcPrb9/gXVHq4486cKvk3f3HoPJSze2jPXs5slTM8mhS4v4vrA8wfnEN/28cabmiZHRhpwN+PBaqmI9ONw8L3hW1FgcF74pS+z3Lfbvae6nmZGxnS8Usj05LPOCF9QHGL67JE4dZrigZne90P+BXzHkVC6PLO+2eD8tlfF/raGvP7q4rte3n7o1tvFRz7bDgXSpmeXRS4hEvrE/Q/WrW/q72CyxXlIymO/U/HsWe/31Oq1gcXdxp/3xYLCMdnW72v4buXO2KUo8O3rS656OHBe9SMcuj4xKPeWF9gD652v4fL6mo/SsZTXfqfzx2/alOq1isbbHD++fDilq71dCwKtrHilK7rCi6HNnbSzU3Q/s2c7CmxINeWJ+gafhwjv14k4rK7tfcpM7JZnXt/AHnTdeU2tq70K9VtbY34x29JvcRKXX3vrZeqrl+vHc3e/klHvbC+gTV5fG/q6/jioq+tKTO0WYx/Hjs/d/l3FI7RzR8Ra3VbuUBSbGi1Gor7fKuXqm5PDbFcks88IX1EXa/DcqPxhUNj5I6f0fZJyp1Ta3FQaWvKHX/2/kkXqr5WCcoEQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAOCT/T9J3RvNu/RfNgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "!gracey.py -com strain_A1g/Y_energy_A1g-2.txt PNG strlt=view strc=br str=A1g,0.95,0.7\n",
    "from IPython.display import Image\n",
    "Image(filename='Untitled.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xmgrace -nosafe -noask -hardcopy -hdevice PNG -batch /home/cam1/.ppbatchfile\r\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABXgAAAK8CAMAAABiCycpAAAAQlBMVEX///8AAAD/AAAA/wAAAP///wC8j4/c3NyUANMA////AP//pQByIbxnB0hA4NAAiwDAwMCBgYFCQkLv4+Pfx8fPq6vxjVQBAAAACXBIWXMAAAsTAAALEwEAmpwYAAAADnRFWHRUaXRsZQBVbnRpdGxlZPPPpzIAAAALdEVYdEF1dGhvcgBjYW0xmdSm6wAAABV0RVh0U29mdHdhcmUAR3JhY2UtNS4xLjI1jnpa/wAAIABJREFUeF7t3dt2o0yWqFHL4x99jcjK93/VFscIQjitIAQImPNil4QQFb22/RWJDv76AgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgL3c0/sPVbJtNPNg1Wybbur9sBng2u63W3r/0dFbPZvex4ON+MGqbve/1WHTID00AH1J4w11V9Dqdpspb9119xY9eO9v189PqJJDAzCcwEZb6uHeENRY3Z4Md/XtHwyBfq5ss2OyCeDimmRO6/jo6HBZ9vniwb1/rGqe0z1Yhb2erjY8nUwD0GjyGO6NJ7xzp7whrOMpb3MOHG2Mn1DdXGoAmDMJb3QCO7ndbxhvNk9qgxs/OY5w81A9c/UBgGl4m3c0jHfSbMbvDevDOzmnTc+dn48AwFeSy+FE9vnOw+R0tnssrXa41lA1N4UXOLn27QZ17kcW0vCGy7TJpYOJfsen8IYLvs0t4QVOraqb5jb/9n96F9g/xe2cvhr29DaFSL/jj+Gt22cKL3BmVR/Jyb/3XxG3c3qVdnpvYrgWPDnJje7cuzX8fACAw6uGxDUnrekFgtuM8cFF4R3fdXaLzopDeNsLvN2jw4MAZ1OHf+RvEd7xv6O5rDycYIfwdhcahBc4s+ibFarnbwRrvjwsFR5cEN77eJrbnGAPt8cIN+8ka/14AIDDuy8v3ILwjtc1utvdd5U9YtzvP1xoEF7gzAoKtyC8dfzyXfvFDW196/7kd/LJ4rAjwKkUFC7O62tvJ+vfsjAavgf91l3iHS7wdlvG2wDnUlC4yXltGt6ny8Vf8ZWEqeFA3QlwLN0T4ARu8418RRrecJz5g/7U3XH3SXNb6Z4AJ3D74apAJw1hY3xwEt7JSW78wOjH7o7Hid468dg0eQcFwHk0iYu+ZCF6pDEpbm98cBLef34tZLcxfl0t3GwvDz8lOT40wLnEb6f9evqenOgcdBQenNQxunOfKemku5M/hznXXeEFzqz9YrK2fFV4V0H1fMb67OktZEOT409GDNvipkfvXmjfU/bcXeEFzqw55b01V1RDAZub051mTcMbvoUhfBiuvVQ7PBgJEZ69zvAlvMC59eUNBayq+ANmP2qfF1+3ff7z7k1u28hOuzsefPzkxJNoJ4Dz6T9BFhfw1+xV/Z9qf5wph/Y2Vxse28MV3OF8Nuluf8LbHCJ69sTt1xUAHFv8mlljafaaF9/ilFbpq3XBvZ7uCnBtS8MLwELCC7Ax4QXYmPACbEx4ATYmvAAbE16AjQkvwMaEF2BjwguwqeZLFGqf5wUAAAAAAAAAgIuJ/5xs86cJJ394u9/cmf/bWABkmfyd2Lr581jV0x8fHP4W4nORAchU1XF42+42LU4KW/dnvNOtACxQtae8w53hVvJp0cqHRwHeKYT3Pv7R7eFGp3aqC/BOIby3EN74HLe6ucgA8E5xePsb4zWHVnMV+OasF+BtxvBG57njuW93pxXeZdZvAMgRqsKv4X1oviQwvJtsnfEd6KirHPRIR13loEc66ioHPf9RFz3prF4Jb/9W3v729H/E3jXMdx1napWjrnLQIx11lYMe6airHPSMR12nFafwWnibN/kOG9cZ34GOuspBj3TUVQ56pKOuctDzH3XRk87q1fA257zdrXXGd6CjrnLQIx11lYMe6airHPT8R130pLN6NbzNKW93Y53xHeioqxz0SEdd5aBHOuoqBz3/URc96azG3kZvIpu9GDPuOPdguQMddZWDHumoqxz0SEdd5aDnP+qiJ51V/D7e2Q9QDIR3sMpBj3TUVQ56pKOuctDzH3XRk84qVLb+4SPDvXHrkcZ3oLVa6ioOtNYDLXXZWhc96axCeO/D90HOfwFkODOebP5sB1qrpa7iQGs90FKXrXXRk85q8ppaG9zhC3m7b4ms+vPcauUPUKzjQGu11FUcaK0HWuqytS560lnV4fy2K+7wItujyM3dW3fpN3T3UOM70FotdRUHWuuBlrpsrYuedEr3/hMl41ntra6Hb2Wous3Np4Xre/wlOUca34HWaqmrONBaD7TUZWtd9KRrqOKvgOzPce/3++Sa75HGd6C1WuoqDrTWAy112VoXPYnBkcZ3oLVa6ioOtNYDLXXZWhc9iYHxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo4jxAfmUo8itl24HmCEZ72B8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QD7lKGJ8QL6Ll6Oqw+377Xarq3A/Uo9jqu6tYb+Ljw9Y5NLlaFI73qlvj5hWt6jEozrs9rjZGPa69PiAhS5cjqqOw9t2t2nx8zlvFXZ7PKc17HTh8QGLXbgcVXvKO9wZbkXnwINbOOOt75NHrjw+YLFrlyOE937rkzreGD1OcofdqnRc6X2A3127HCG8txDeZCRVHbbV9XiRoXPt8QHLXLsccXj7G+M1h/F+tFv7wlrc3muPD1jm2uUYixqd547nvp3mjWPDo82LcW16x0evPT5gmWuX4/fw3ps3jkWPVm18x/ecXXt8wDLXLsev4e1eTZte923eyjvscu3xActcuxy/hrf7hNo0vPH7evtrD4N4L4CRVgS/hbd/224S3uact9/n2uO7rv/934z/pXvBT65djl/CO3yTQxreSngv78/ftrWD5p7w8rJrl2MsavQmsiiy7Qtpc/82EF7+NOGN7v8VXl537XKEU9nZD1AILz9LwvtHeHndtcsRKlv//JHh50sN4az42uO7tiS8X8LL665djlDU+/De3NvM15Ol4Q1xvvb4ri0N758/8T34l2uXY/KaWhvc4Qt5+2+J7Ay7DR8Wrn2AgqfwwuuuXY46nN92xR1eZLtHn04bw9tsbHYP3b34+K5NeFnuwuUYXjrrLxs8mlvXw9/+CW8Y6/fsNz42j3/2p3Hh8V3eNLxxg9s3l/11yZefKUdQ3eMPC4ebQfMX1yYbjO+6fgzv3//7+6d9m2/zBl8XfpmjHEWM77om4f0b3/zb/UfLeS9zlKOI8V1Xm9Ve1ODH7e4sV3X5mXIUMb7r6s5oR9Hm7sajwO2pLzxTjiLGd11xbf/MhLf5TPGwEaaUo4jxXVcc3sfpbbT56RYklKOI8V3XJLzhXQ3/N17aFV5+pBxFjO+6fgjv3yi8XlzjB8pRxPiuaxreUXNpt31bQ3TdFxLKUcT4ruuH8A7l/fP3rw9P8BPlKGJ81/Uc3v5jau1XpPvIMP+iHEWM77qewvunv+tMl18pRxHju66n8PYfG3aqy++Uo4jxXVca3v6Prv2v/8Qw/INyFDG+y0r+2GXzdWRtcdvLuy0B5kfKUcT4Lqv99rE+se235PQZ7r+VrOVqLz9QjiLGd01daBP9V+LE5XXZgXnKUcT4SP1tT4P7/Pp+MmYpRxHjY+pP9JXok+9Hh4hyFDE+ppIvaBBeZilHEeNj4m/63t7JPegpRxHjYyJ5Pc3bGpinHEWMj4n2jwv3t//4LAU/UY4ixsdE9wU5Pee7/EQ5ihgfifEtvrLLz5SjiPHx5I/PC/Mb5ShifEA+5ShifEA+5Shy66XbAWZIxjsYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOSJVHW7fb7dbXYX7kToMzfiAfMoxalI73qlvj+hWt6jEozrazfiAfMrRq+o4vG13mxY/n/NWcZ+ND3jRd7ipHL2qPeUd7gy3osYObs54gUzfg+6ucgQhvPfbPbkxqu9hN+MDXjBmd0ivcgShqLcQ3mRAVT3ZZnzAbybdbcurHEEc3v7GeM1hvD+NsfEBv0i625RXOYKxqFFax3PfTvMGM+EFXpdm90E5Ir+H9968vUx4gdel1f0W3olfw9tcaEjCOzVuB2j1rdWKn/wa3u6TbM54gZdNTnUHyhH8Ft66uym8wMvS5raUI/glvMM3OQgv8KL//kub21KOYCxq9CayKLL3mUs0xgf84L9GmtyOcgThVHb2AxTCC7yojW7DpYbfhMrWP39k2KUG4N/G6DbS5raUIwhFvQ/fB3mb+Xoy4QV+MoluI21uSzmCyWtqbXCHL+TtvyWyI7zAnCG6obpfwvurOpzfdsUdXmR7pHb6tynG28YHtOai20qj++2Ta8Hw0ll/TffR3Loe/vZPNfkYhfACiR+i20ir+y28/1Dd4w8Lh5sTxgf8o7qNNLu+nayU8cHV/ZLdr6fyfilHIeODS/u9uo20u8pRxvjgsrro/p7dxiS7ylHI+OCShui+VN2evzL8LsYH17MgugnlKGJ8cC3l0W0oRxHjg+sYoltY3S/lKGR8cBFvi25DOYoYH1zAW6PbUI4ixgcnN0T3fdX9Uo5CxgdntkZ0G8pRxPjgrNaKbkM5ihgfnNKa1f1SjkLGB+ezbnQbylHE+OBcVj7V7SlHEeOD8xiiu3J1v5SjkPHBSWwW3YZyFDE+OIFNo9tQjiLGBwc3RHe76n4pRyHjgyPbI7oN5ShifHBUe0W3oRxFjA8Oac/qfilHIeODw9k5ug3lKGJ8cCBDcvet7pdyFDI+OIiPiW5DOYoYHxzAByW3oxxFjA8+3MdFt6EcRW69dDvwCT4uupLxDsYHn+rjohtRjiLGByv4TjfkGqL7kdX9Uo5Cxgdv9j1IH3jZh0e3oRxFjA/easzuwvQeILoN5ShifPBOk+5ml/cg1f1SjkLGB2+UdDenvMeJbkM5ihgfvE+a3e+Xyjsk9zDV/VKOQsYH75NW9/vX8B4xug3lKGJ88DZpdFvpTsExk9tRjiLGB2+TNreV7tQ5cnQbylHE+OBt0ua20p2+jvY62izlKGJ88DZpc1vJPieIbkM5ihgfvEua3F60x0mq+6UchYwP3iZNbmt48DzRbShHEeODt0mb2zruW8b+RTmKGB+8Tdrc1vmi21COIsYHb5M2t3W25HaUo4jxwfuk0f1+elPDWShHEeOD90mr+y28zDE+eKM0u2ftrnKUMT54m//+++8i3VWOMsYHbzG8d+Ea3VWOMsYHpcY3jHVvX7hAdpWjkPFBiSS6g1NHt6EcRYwPlppv7jUoRxHjgyWuHN2GchQxPsh19eg2Ll6Oqg6377fbra7C/U5VPzbfw917a9jv4uODLGNyLx3dxqXL0aR2vFPfHjGtblGJu62tcWty/9LjgwyiG7twOao6Dm/b3abFk3Pe9tS2OentS/t4TmvY6cLjg5dJburC5ajaU97hznArOgd+6Hs7BroOFx26B6Z3gZTozrh2OUJ478N13PFGe2f4z36/Kh1Xeh+IiO4Prl2OEN7xBbSwKXLvLzXU9XiRoTOzL+CS7i+uXY44vP2N8ZpDrL8A3FxyeBQ6au/Mvnyq038a6lNEzRXdH1y7HGN4o/Pc+M1jvap/91jzYlyb3vGRa4/vQC7x+f/9Se6rrl2O18IbX/at2viO7y679vgOY8yu9K5Hc3NcuxyvhLc7zY0uL8SfqLj2+I5i0l3lXYHk5rp2OV4Ib/NRtaa88aZwt7/2MIh24mMk3VXe9xLd12hF8EJ4W/V0Y7h77fEdQ5rdb+V9G9Fd6trleDW80Tluf1d4DyOt7rfwvsPYXNFd4trlGHsbvYls9h8B9XSr8B5HGt1WuhM5QnNFd6lrlyOc6P77AxTpVuE9jrS5rXQnXhMlV3OLXLscoad1CO/zlYZm6+Rby0Kv4618orS5rXQnfqe573TtcoTwjmWdvHNsNHx0rRPifO3xHULa3Fa6E/8mue927XJMXlNr0zp8Ie/wLZHdg1VX2uHDwrUPUBxGmtxeuhs/Et01XLscdTi/7Yo7vMh27z6ddu8er7pvg2w2NndDdy8+vkNIk9tKd2Ke6K7lwuVoPxdxG18nezS3roe//dO/Yax5G9ntXoeNj3vjn/1pXHh8R5E2t5XuxJOxuaK7BuUIqnv0ulpf1+ivTTQb7/EuX8Z3AGlzW+lOxEJzRXctylHE+D5e2txWuhOdKLmauyrlKGJ8ny+N7rfuztHcTSlHEeP7fGl1v4V3SnJ3oBxFjO8A0uzq7khz96IcRYzvCHT3meTuSzmKGN8h6G5Mcz+AchQxvoOQ3ZbkfgrlKGJ8B3Lt6GruR1GOIsbHx4uSq7mfQjmKGB87++d5vOZ+KuUoYnzs6B9XriX3sylHEeNjN2N2k/Rq7gEoRxHjYy+T7vblldyjUI4ixsdOku5+f2vukShHEeNjH2l2v7vyprvxoZSjiPGxj7S633OvsPGxlKOI8bGLNLqtdCc+l3IUMT4211xSSJvbSnfkcylHEeNjS+PrZ2lzW+nefC7lKGJ8bGRsbvsKWtrcVvoUPpdyFDE+Vhcld3jXQprc3uRpfDLlKGJ8rGmmuZ00ua3JHnw05ShifKzjx+R20ua20p34XMpRxPh4u1+a20qb20p34nMpRxHj441eSW4nbW4r3YnPpRxFjI83iIP7a3M7aXS/dfdQlKOI8VEiKe5Lye2k1f0W3kNRjiK3Xrod/m1pcEdpdnX3GCTjHYyPbIXFHejukSlHEeMjw5uS29PdA1OOIsbHa97b3J7sHpZyFDE+frVGcgPRPSTlKGJ8/CgK7jrN5biUo4jx8WQaXM1lhnIUMT4igsuLlKOI8fEluGRTjiLGd2lJcCWXVylHEeO7KsGlhHIUMb6rEVzeQTmKGN9VTIOruZRRjiLGd36Cy/spRxHjOy/BZT3KUcT4TigpruTyfspRxPhORXDZiHIUMb5TmAZXcVmdchQxvuNKattK94F1KEcR4zueNLaNdB9Yl3IUMb5DUVs+hHIUMb5DEFw+jHIUMb7PNa1tJ90H9qEcRYzvw6Sl7aW7wb6Uo4jxfYo0tf+pLR9MOYoY3/7UluNRjiLGt59pcBWXI1GOIsa3saS2rXQf+HjKUcT4tpCWtpfuBoehHLGqDrfvt9utrsL9TlU/Nt/Hu8a3nrSznXQvOCLlCJrUjnfq2yO61S0qcbe1NW41vjdLO9tJ94KDU45BVcfhbbvbtHhyzntvToGbk96hvMb3Dmlne+lucBrKMajaU97hznArOgd+6HsbAm18y6Wd7aR7wRkpRySE9z5cxx1vtHeG/xTepdLOdtK94OSUIxIFNYR3ZkJ3lxqypJ3tpbvBZShHJA5vf2O85hDrLwB/Gd8/pJntpbvBFSlHZAxvdJ4bv3msV4V3mRlfI43rk/QJcHHKEXktvPFl3wuPL43rk/QJwODC5Xj2Snib95zdwqWGqcmep5TWtZPuBSQu14rXvRDe6n5v0zs8Gj92Lmlcn6RPAF514nLkeyG8rXrceLbxpXF9lj4DyHe2chR5NbzVeMp74PGlRX2SPgF4lwOX4/3G3kZvIpu9GFMfK7xpUv8hfSqwgmOUYyPx+3j/9QGKsHXuwT2lHf2X9LnAVj6tHLsKla1nPzI8Gj+6tuP40o7+S/pcYFc7luPzhPCGsk6/nqw3fnRtxfGl8XxNehTgA61YjuOZvKbWpnX4Qt7hWyK7B6vxNPjn8aVFfL/0vxE4ip/LcUF1OL/tiju8yHbvvoH33j1e1UN3/7ulNXyj/r8DOB3hHbSfi7iNr6o9mlvXw9/+qbrNzdvIbvc6fFXDv8I77AOQEN4fVffodbW+tff7fXLN1/iAfMpRxPiAfMpRxPiAfMpR5Jzj+043AG91znJs5nzj+x6kDwBvc75ybOps4xuzK72worOVY2MnG9+ku8oLazlZObZ2rvEl3VVeWMm5yrG5U40vze638sI6TlWO7Z1qfGl1v4UX1nGqcmzvTONLo9tKdwLe4Ezl2MGZxpc2t5XuBLzBmcqxgzONL21uK90JeIMzlWMHZxpf2txWuhPwBmcqxw5ONL40ub10N6DcicqxhzONL01uK90JeIMzlWMHZxpf2txWuhPwBmcqxw7ONL60ua10J+ANzlSOHZxpfGlzW+lOwBucqRw7ONX40uh+6y6s41Tl2N6pxpdW91t4YR2nKsf2zjW+NLu6C+s4Vzk2d7Lx6S5s4mTl2NrZxqe7sIWzlWNj5xuf7ML6zleOTZ1zfKIL6zpnOTZjfEA+5Shy66XbAWZIxjsYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5BPOYoYH5Dv4uWo6nD7frvd6irc77fWj6338W51bw37XXx8wCKXLkeT2vFOfXvEtLpFJe62tsatyf1Ljw9Y6MLlqOo4vG13mxZPznnr5tS2ak56u/uP57SGnS48PmCxC5ejak95hzvDregcOFyJeAS6S2100aF7YHoX4AXXLkcI7/3WJ3W80Rgu5T72a7dW6bjS+wC/u3Y5Qnj7ssabHsZLu/3DdT1eZOi3T+4BvOLa5YjD298YrzlM9OFtX1iL2zu3L8C/XbscY3ij89zx3DfWPdy8GNemN2wfbwG86trleDW89/H9Y1Ub33AJYrgB8LJrl+PV8PbvNRvuhF2uPT5gmWuX48XwVtNNVXj7b3/tYRDvBTDSiuDF8NbPH2fr97n2+IBlrl2O18J7T7rbnPIKL7DYtcsx9jZ6E9nTPwLiL9LpCS9Q4NrliN/HO/cBisZMd4UXKHHtcoTK1iG80ysN1fRbczqh15PNAK+4djlCeMd36g5fh9ML3Y16HOJ87fEBy1y7HJPX1NrCDl/I279zN3S3+V6y4cPC4V0O1x4fsMy1y1GH89uuuMOLbPfu02nNO3bD++6ajc3u0bvLrj0+YJkLl+PeF7W/bPCIbF0Pf/unf8NY//cm+r3aDIc/+9O48PiAxZQjqO7RddyZl9TaPaYvvRkfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85ihgfkE85itx66XaAGZLxDsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOIsYH5FOOSFWH2/fb7VZX4X6/tX5svYf7xgfkU45Rk9rxTn17RLe6RSXutrbCVuMD8ilHr6rj8LbdbVo8Oeetm1PgqjnpHbYYH5BPOXpVe8o73BluRefA4UrEI9BDj40PyKccQQjv/dZfxx1vNIZLvo/9hq3GB+RTjiCEdyxr2PQQLjAIL1BAOYI4vP2N8ZrDhPACJZQjGMMbneeGxkaih+PNAC9RjuDV8N69qwEooRzBq+Ht32vW6N7XO4r3AhhpxU9eDG8VbTI+IJ9yBC+Gt44+zmZ8QD7lCF4L7z3+GLHxAfmUIxh7G72J7OlaTPxFOsYHLKEcwW8foGhMu2t8wALKEYTK1rMfGf5qzoWn3xRpfEA+5QhCeMd36oavw2mF7vY9Nj4gn3IEk9fU2sIOX8jbv3M3dHf4MnTjA/IpR1CH89uuuMOLbPfuG3ir57c/Gx+QTzl6976o/ansI7J1Pfztn6rb3P/9iWgv4wMWUI6fVPfodbXpS2qB8QH5lKOI8QH5lKOI8QH5lKOI8QH5lKPIkcZ3oLVa6ioOtNYDLXXZWhc9icGRxnegtVrqKg601gMtddlaFz2JwZHGd6C1WuoqDrTWAy112VoXPYnBkcZ3oLVa6ioOtNYDLXXZWhc9icGRxnegtVrqKg601gMtddlaFz2JwZHGd6C1WuoqDrTWAy112VoXPYnBkcZ3oLVa6ioOtNYDLXXZWhc9icE64zvQUVc56JGOuspBj3TUVQ56/qMuehKDdcZ3oKOuctAjHXWVgx7pqKsc9PxHXfQkBuuM70BHXeWgRzrqKgc90lFXOej5j7roSQzWGd+BjrrKQY901FUOeqSjrnLQ8x910ZMYrDO+Ax11lYMe6airHPRIR13loOc/6qInMVhnfAc66ioHPdJRVznokY66ykHPf9RFT2KwzvgOdNRVDnqko65y0CMddZWDnv+oi57EIPwxIICXpSkBAAAAAAAAAAAAslT17Xa7p1v3NL+iZOv0bnVv7lXjo1tZsNRu0/OTVjezjOetMzvd75tPdmYVz1tnfgK5DOvAAAAFnUlEQVTqp+ds56U1t1vie5t6Xsz81miJe/1inV91a6Zc3z5ntvMrSrZO71a3ztb/RyxYar/p+ed/bTPLeN76vNP9Vm891ZlVzGyd3u2WWbfbdvHSmpt17vepg6fFzG+Nl7jXL9b5dWP/+trvJzY1v6Jk6/Rud/7Y/MBs+wOyYKmdevvwzi3jaevTTtUeKXtaxdzW9Cegv7f5XHsvrfkxzh3Dmy5mfutkiXv9Yl1A3f+k3vf6iX0yv6Jk6/Ruf072+AHZthILltrd3CG8M8t43prutMeZ+fMqZrcmPwHjvZ2y9tKav6odV/i0mPmtkyXu9Yt1AeOMd/t5SM2vKNk6uTtektr6bCJ/qa3qtsOZ2fMyZramO+3zC5euYnbr9O4wz92y9tKaG7utcG4x81vHJe72i3V+4cdg+xTMm19RsnV6d7z2v/EP9YKlth7nEZtPe2YZz1vTnept59lLVzG7Nbk7lGHjn4DRS2tOtm1sbjHzW8dte/1iXUD4Z8fwT469za8o2Tq/09Y/H/Or+HWpzUvvO4Q3XcbM1uTuPhca0lXMb03u1v25+V4XI19ac79t05/RYG4x81ufl/i8hTIz/zO3s/kVJVvnd4p+ijYxv4rfltr++23zpD0tY25rcrduMna/b/1moiVLfdxvypu+Yr+ZdDk/b93t92xuMfNbn5e48S/WBczNfV/zK/rl16638fnO/Cp+W2q7xrB1I0/LmNv6dPd2b99MtO2F3iVLbc7Ob3W9/Tvfeulyft662+/Z3GLmtz4vceNfrAuYm/u+5lf0y69dZ3hnzFbmV/HLUrv3+IetG0mXMbt1evfRsiZkzfvrN53rgqU2nj4GsKXn5fy0dbffs7nFzG99WuLWv1gXMDf3fc2v6Ndfu0a98f8F86v491L7F4o3b8SCtY47bvzmtwVLbbTv9d8rEM/L+Wnrbr9nc4uZ3/q0xK1/sS5gbu77ml/Rr792X81v3sb/HppdRbp1erfq99s2ZV/pMn7YOr077tic+vbbt7BgqQ913V1u6Ldu7Gk5P27d7fdsbjHzW9Mlbv6LdQFzc9/X/Ip++7Vr99n6x2N2FenW6d3hImTYupEFaw2vqGz7Ns4FS+1fVtuvvOlyft662+/Z3GLmt6ZL3PwX6wLi95BsnIIfzK8o2Tq30/YvrMyt4mnr9G77WtVgy4EvWGv4/ds2vAuW+vh/2uA25d38x6Dx0pr7W1vOMjK3mPmtyRK3/8W6gPiXK/z/xp7mV5Rsndlphx+PmVU8b53e3S28C9Za7RTeBUsd97xvO9XRS2tOtm1sbjHzW6dL3OEX6wLiX67JA7uZX1Gy9Xmneod/Yj6vYmbr/E6b92F+Gf9ea3Rvy+kuW2o/z80H23lpza3dwju3mPmtkyXu8Yt1BcM/zfp/q32A+RUlW9Odxh+Pasv/fU5XMbt1dqft+zC7jHTr9O74z9CNV7tgqXuHN13dz1t3C+/cYua3xkvc5xfrAsZy9fPfX7Ki/j+SrelO44/Hpj/V6Spm1zY74e37sGCtw9lQOCvaxoKlDq2oho5s7aU1d5u2HWawZIk7/WJdQTfwcI19f5MV3Yd/5ibrnNxt3ju/w3XTJUvtbb3Qr0Vr7b+MN/qd3EbJUjefa++lNbe3t57mKH+Ju/1iXUHz9vjPmmu8ompcWrLO6G4Vfjy2/t/l3KUO9hj4grU2u913KMWCpTb30ilv6pU13/etWO4Sd/zFuoTNvwblV/GKwq1knZ+x7AMtdclaq52WvmCp23+dT+KlNe/rAEsEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAuLL/BzFeHywSWvUyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "!gracey.py -com strain_Eg.0/Y_energy_Eg.0-2.txt PNG strlt=view strc=br str=Eg,0.95,0.7\n",
    "Image(filename='Untitled.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xmgrace -nosafe -noask -hardcopy -hdevice PNG -batch /home/cam1/.ppbatchfile\r\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABXgAAAK8CAMAAABiCycpAAAAQlBMVEX///8AAAD/AAAA/wAAAP///wC8j4/c3NyUANMA////AP//pQByIbxnB0hA4NAAiwDAwMCBgYFCQkLv4+Pfx8fPq6vxjVQBAAAACXBIWXMAAAsTAAALEwEAmpwYAAAADnRFWHRUaXRsZQBVbnRpdGxlZPPPpzIAAAALdEVYdEF1dGhvcgBjYW0xmdSm6wAAABV0RVh0U29mdHdhcmUAR3JhY2UtNS4xLjI1jnpa/wAAIABJREFUeF7t3dt2m8q2QFHJbbf1DCT5/1/d5ioo8KU0BQjo/eHEYMQqz2OPTbCk3G4AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACwr6IsyyLdmfg8pBo2lh7wuatMdgGwqLy30pKOlNX9/ijtwgOKqt31aDMAX+mS+U15i+nnFh5QH9FSXoCfVO2lbHMVm36uVU1zuvCA4l7Vdxmay163GwC+V/SXrfVF61I0k8vdxQdUfZjr8nYfArCsGqJa38Ydf6Y16/HCA4rhgrg+uv8YgEWPmwjlYniT+wyLDxg9wUF4AX6vXPrtWn0nN93XW3pAmmkAvlYtNHN2o2Fs4QHFQosBTq5+xm37HINcS8389pdlCw8o5y0GOLeiqpubPhHhV6rHb80e6ojfFl+ntviA4uvLY4BzKroLzvrO7KyK36uWmlmfp+xfp5Z8fuEBn//5zP8qwMEV/X2BpVuzXT4nhk+2cZ1Vs9k9vFxifBdh4QHtq4bdagAu5XERuvCc3K61E92n+mvaWXlH++qWD1FdeED/Zg3KC1zJ6CkFxfyNwuo3D0s9Ptt187GjMd6V3L5YekDb49l/GeC0ylk4czTvc5NEc7In/fTSA5oaT3cBnFiweZObCa1JWWdPLVt4QHNQesMC4LTSMOaav0ptEt75p+d75tfFAGe2kMEs83e4+SG88wcs/loP4LTCF5uzjk5aOw/v/AHNUbFFABzIfX7LdeS+YHZI8vifwzv7D5bu8QIXMn3CV5rEpLmN2SHpxep9dMaFa9n5nuBTKwCOZfIsg9n75KTP4a0lR8ybOb5jO7+WXXjAUosBzqt5Y7ImjsXwj/GM/nmIn1RDWIeXwI1/fTa/r/B4wKB/twiAa2he0vB5yVm/NWSXxPrD6UGJPtV1Rvtm1vcsuvI+Lnn718UtPOBxE8L78QJX05X3cWe2KIY3zlnWPKIuaTG6fK2vnPsbBn3D+6QuPaB5vVr9AG/HC1xP9w4K4/x9H942mn04O3Vchwp/Xs5WZflI6sIDht57V0jgktLfmf0Q3uZdztNfsxWT7eaIZDP9vVxNdQEaP4UXgBcTXoCNCS/AxoQXYGPCC7Ax4QXYmPACbEx4ATYmvAAbE16ATdXvTlZ5MS8AAAAAAAAAAAAAAAAAAAAAAABfuAPkS1NCjnXGd6CzrnLSI511lZMe6ayrnPT8Z33qQfTWGd+BzrrKSY901lVOeqSzrnLS85/1qQfRW2d8BzrrKic90llXOemRzrrKSc9/1qceRG+d8R3orKuc9EhnXeWkRzrrKic9/1mfehC9dcZ3oLOuctIjnXWVkx7prKuc9PxnfepB9NYZ34HOuspJj3TWVU56pLOuctLzn/WpB9FbZ3wHOusqJz3SWVc56ZHOuspJz3/Wpx5E70jjO9BaLXUVB1rrgZb63FqfehC9I43vQGu11FUcaK0HWupza33qQfSONL4DrdVSV3GgtR5oqc+t9akH0TvS+A60VktdxYHWeqClPrfWpx5E70jjO9BaLXUVB1rrgZb63FqfehC9I43vQGu11FUcaK0HWupza33qQfSONL4DrdVSV3GgtR5oqc+t9akH0TM+IJ9yhBgfkE85QowPyKccIcYH5FOOEOMD8ilHiPEB+ZQjxPiAfMoRYnxAPuUIMT4gn3KEGB+QTzlCjA/IpxwhxgfkU44Q4wPyKUeI8QH5lCPE+IB8yhFifEA+5QgxPiCfcoQYH5BPOUKMD8inHCHGB+RTjpB7J90PsEAyXsH4gHzKEWJ8QD7lCDE+IJ9yhBgfkE85QowPyKccIcYH5FOOEOMD8ilHiPEB+ZQjxPiAfMoRYnxAPuUIMT4gn3KEGB+QTzlCjA/IpxwhxgfkU44Q4wPyKUeI8QH5lCPE+IB8yhFifEA+5QgxPiCfcoQYH5BPOUKMD8inHCHGB+RTjhDjA/IpR4jxAfmUI8T4gHzKEWJ8QD7lCDE+IJ9yhBgfkE85QowPyKccIcYH5FOOEOMD8ilHiPEB+ZQjxPiAfMoRYnxAPuUIMT4gn3KEGB+QTzkGRXW/38t0b60oq6r7RNkcVPSfMT4gn3L06qJ+qtL9dZD7HLcfFvd7X17jA/IpR6esivaiNy1v+cjsZ5ybDz93dTuMD8inHJ2ut/ehqf3+0Y6i/3i4JWF8QD7laPU3Ex4Xs8P2cL37+KTwAgHKMVVObzUUk1+3lf2W8AIByjHV3cV9bI4H9Jnh5tOPOhsfkE85Jor6V2yjzfp3bUVZjm5EfJa3eFwVGx8L/n76k+6EB+UYK5Pn8X6GtmqfZtb1uHkW7+gg47uo/xa1n/vbbvyTXr6kHA/19ez4d2lNZ+vIthe6jfpuw+iquH3u72DYz8nVZa2va5vIDh80n/rXZ/g/5eVBK75Q31JIJtL/Eq3ObberuQD2yrWrG6o69Pb2p/3gXxPitr7KyxeUY6q7xu0MG1UX2+Lzjzq9Xrl2cf/96z8YwvuZ3Ft9n+Fvs/GnTm9/ECSUY2p0bXurw9t90D2RrP212qi8xndNf4aL2VF4/9QfPGLrkpevKUdi8gSy4Yq3C28X3McLi43vmv72sR2H9/ZfF99Wfde3vfqFlHIkJi9dG94epw3v6JVr/QftH1zMn+FSdhzev/UTyfqN5lPCyzLlSExeujY8c6x95fAQ3keB2z+4rHF4b017B8LLl5QjMXnp2vCuOG1phZdUEt6x8T3eP+3zztzzpaUcne7Stuiucbv+9vcauj/7e76Ve7y0vg9v/+Gf+uK3bu+/+im/44O4JuVotfcSbkX3L03UL1mr/+ze87x/kfCw2Y/N+K7u6/D+fdxp+NNe+/6pD/7yeK5EOVr108juZdW/KO3xrmSfF7flcIHbvHlDvel5vLS+Dum/0dMduieZ1de8w+e5MuXoleXjn1L79Pj48R45zf7JpvFd3dfhfVzwPq59//PUXhrKEWJ8V/dleP8+XkrxCG/9euJ+L1emHCHGd3Vfhbd754bGI7yj+75cmnKEGN/VfRXef6N7Co9bu8JLSzlCjO/qvgjv3/G93PrZDP3upYO5HuUIMb6rWw7v4010Gv/6pzX8c8FLQzlCjO/qFsObdHcob/+ekVyecoQY39UthXfW3f6fpfCSYTrKEWJ8V7cQ3kl32w//eEt0JpQjxPiubh7eSXfbf/Fy9C69UFOOEOO7unl4/43u4/5rr3T9G0AklCPE+C6ued+b9BkMI02Em3+LojU+kgtTjhDju7bmn7T8r72f0Jp2t70Y7t+VrKG93JQjyPiurHlz82lOk+52uyflddsB5QgyPn7hT3unoSuv55ShHDHGx88mr5v442YDyhFkfPwoeYMGv2FDOYKMj580/97aiJevoRxBxsdP/k5/n+ZFbNyUI8j4+En9lIbhCWd/u5dUcHHKEWJ8/GjyFDM3eKkpR4jx8bP2ZRayy4NyhBgfv+L1wkwoR4jxAfmUI8T4gHzKEWJ8QD7lCDE+IJ9yhBgfkE85QowPyKccIcYH5FOOEOMD8ilHiPEB+ZQjxPiAfMoRYnxAPuUIMT4gn3KEGB+QTzlC7p10P8ACyXgF4wPyKUeI8QH5lCPE+IB8yhFifEA+5QgxPiCfcoQYH5BPOUKMD8inHCHGB+RTjhDjA/IpR4jxAfmUI8T4gHzKEWJ8QD7lCDE+IJ9yhBgfkE85QowPyKccIcYH5FOOEOMD8ilHiPEB+ZQjxPiAfMoRYnxAPuUIMT7gV/73v/89NpQjxPiAH/2v9dihHCHGB3yvq+4ou8oRZHzANxaiW1OOEOMDvrB0qdtRjhDjAxb00V2q7k05gowPSH0f3ZpyhBgfMPZzdGvKEWJ8QK+P7g/VvSlHkPEBjZ+j+/H4UDlCjA/4TXR77aZyhBgfrGB0bfj+fq7uI7t9epUjxPjgxSaBens/R/eWdLf5wpQjxPjgpdJCvbNfXOo2Jl9T+3UpR4jxwSvNCvWu+uj+WN3Z11RTjkFR3e/3Mt1bK8qqGn1ivGl88EJpoN6zvBnRraVf04fwPtTZ/VSl++sgj3M83TQ+eJ20Tx/vVt4hub+N7m3xi/r4UI5WWRXtRW9a3vJ+L77eND54nTRPH+8U3meiW0u/ooZytLrefl7zJvunO5JN4YXXSevUSA/axXPJbaVfUEM5Gv3dgzIJ6/fXu8ILL5TWqZEetLlIdGvpF9RQjolyequhmP66Ldm8CS+8UFqnRnrQtoLRvX3xRQnvVDW9ov3hRoPwwguldWqkB20nHt1G+gU1lGOsqH/FNtqsf9dWlGV3mZts1owPXiWNUyc9bBsvqu7tiy9LOUbK5E5C+Vna9mlmTY+TzZrxwcukdWqkB23hZdGtpV9QQzkG5bSpt+bWQl3i+hP17mSzZnzwMmmdGulBa3vdpW4n/YIaytGr7yHU5R3t6n+XVrS7k832iIluL/CEtE6N9KA19dF9XXVr7dehFd/oLmo7w0bVXOMmm80R7R9A3LS4nfSgtawT3Vr6FX14yXBqdDF7q0vbfVA2zU02myO6PUBc2qePTbo7JHeF6DbSr8m7k81MnjE29LUP72SzOaL9A3iBNFAfq4d39eg25l+TckxNXro2vB9OW9pks2Z88EJJodbt7hbJ7cy+JuWYmrx0bXh6WftS4WSzZnzwSpNCrdnd7aLbSr4k5ZiavHSt6C9/2+vgZLNmfPBSo+qu192No9sbfT3K0equZYvuorbrb39zofsz2bwZH7zcytnd+lJ3mXI0upsHRfdPS5Td+/IW3e7u/kOyeTM+WMXK0d25ujfl6NRPI7uXVf9WDY+3IavuVVk9QptsGh8cxPtEt6YcnbIsx68Wfnw8eVOcdNP44N0NyX2T6NaUI8T44J29Y3RryhFifPCu3jO5LeUIMT54R+8c3ZpyhBgfvJ03j25NOUKMD97KAaJbU44Q44P3cZDq3pQjyPjgPRwnujXlCDE+2Fuf3MNU96YcQcYHezpidGvKEWJ8sJdjJrelHCHGB3s4cnRryhFifLC1o0e3phwhxgcbGpp75OjWlCPE+GAbj+YePbo15QgxPljbKLlnaG5DOUKMD9Z0wuY2lCPE+GAt50xuSzlCjA/WcObo1pQjxPjg1c4e3ZpyhBgfvM6Q3FNHt6YcIcYHr3Gd6NaUI8T4IO5KyW0pR4jxQdDloltTjhDjg4BLRremHCHGB0/pk3vF6t6UI8j4INeouZeMbk05QowPfk9ye8oRYnzwO5o7phwhxgc/k9yUcoQYH/xAdBcoR4jxwTdE9wvKEWJ8sKhPruouUo4Q44PUqLmi+wXlCDE+eJDc31KOEOODhuZmUY6QeyfdD5chuTkk4xWMjyvT3GcpR4jxcVGSG6IcIcbH9WhunHKEGB9XMkqu5oYoR4jxcRGa+1LKEWJ8nJ7krkA5QoyPU9PclShHiPFxWpK7IuUIMT7OSXTXpRwhxsf5iO76lCPE+DiTPrmquzblCDE+TmLUXNFdn3KEGB+HJ7k7UI4Q4+PINHcvyhFifByT5O5LOUKMj2MZB1dz96McIcZ3IB/pjktJiiu5+1KOEOM7iI9e+okLENw3pBwhxncIQ3avll7FfVfKEWJ8RzDp7kXKK7nvTTlCjO8Aku6evryaewDKEWJ87y/N7seJyyu5R6EcIcb3/tLqfpwzvJp7KMoRYnxvL41uIz3o0EbJ1dyjUI4Q43t7aXMb6UFHpblHpRwhxvf20uY20oOOR3KPTTlCjO/tpc1tpAcdxzi4TXMP/LVcmXKEGN+7S5PbSQ97f0lxm8vcA385V6ccIcb39ka1fUgPemvz4LaO+xWhHDHG9/YmeeqlB72r5eK2Dvol0VCOEON7e0mfWulB7+e75DbSL+kAXxMPyhFifG8v7VMjPeh9jIP7VXNr6Vf08c5fFDPKEWJ87y/t08dbJioJ7jfJbaVf0sdbflV8RTlCjO/9pX36eLdEZQW3k35FjfQg3pdyhBjfAaR9eqNAZQe3l35JjfQg3pdyhBjfEbxjn0bJzWxuI/mSWulBvC/l6BXV/X4v0721oqyq8Seqx9CM7xDeKE/j4D7X3MbkK+qlB/G+lKNTZ/dTle6vgzzN8eeRw8fGdxC7xykJ7vPJbYxqO5YexttSjlZZFe1Fb1re8n4vJjuKu/Ae015demFwB2lyG+lBvC/laHW9HUe13Z/uuN1d8fIrawS3lza3kR7E+1KORn8zoUw6O7vevVXl+BjjY2Ya3Jc3t5E2t5EexPtSjolyequhmP26ragmcTY+RtYPbi9tbiM9iPelHBPV9AJ3dqOhuE+vio2P25bBHaTR/dDdQ1GOsaL+Fdtos/5dW1GWj6ve+vPCSycJ7jbJbaXV/RDeQ1GOkTK5sfCZ2Kp9mlnX47K+ESG8pM1NP7uBNLu6eyjKMSjHiW1U7S3e+hPN7vpGg/Be2+7BHejukSlHr76lUJd3tKv/1Vr/1N32RsQkvFPDfk5lWttWesz2dPdQtOIb3TVuZ9iomkve7nXDrngvI41tLT1mT7J7WMoxMXlZ2m34uKwTXHTPNBPeK3jf2iZE95CUY2ryBLLhircJb3sPePo3BeM7m8MEl0NTjqnx5exteHsc4T23aW1b6THwQsoxNXnp2vD0snL8bAe3Gk4jjW0tPQZWoBxTk5euFX1iJ9fBwnsCasuulKPVXdoW3TVu19/+XsPkLXmF98AEl7egHI3uXkLxeMZYc8eh6HZP3jpHeA8nyW0tPQS2pByN+mlk97Lq36rh8a5k1b0qq+lblgnvMaSpbaVHwR6Uo1OW5fjVwqPfpY3eI2fG+N5M2tlOehjsSzlCjO8tpJ1tpUfB21COEOPbT9rZVnoUvCPlCDG+7aWpraXHwHtTjhDj25LachbKEWJ8WxBczkY5QozvxZLGJtKj4aCUI8T4XiHt60z6ADg45QgxvueldW2lR8EZKUeI8WVKO9tKj4KTU44Q4/uNtLOd9DC4DOUIMb4vpZntpIfBFSlHiPHV0rjOpA+Ai1OOkAuPL43rTPoAoHfhcrzC5caX1rWVHgV863LleK0Tjy+N60z6AOC3TlyOLZxjfGlSv5M+Fsh3jnLs5mjjSzP6nfSxwKscrRxv5o3Hl3b0O+ljgVW9cTmO4D3Gl3b0G+lDgR28RzkOa5PxpfH8pfQ0wLvYpBznlTu+NI4vkv5ngLeWWw7G/ndPE7iK9D8LHJvwRuSHNz0DcEHCG2J8QD7lCDE+IJ9yhBgfkE85QowPyKccIcYH5FOOEOMD8ilHiPEB+ZQjxPiAfMoRYnxAPuUIMT4gn3KEnHN8H+kO4KXOWY7NnG98H730E8DLnK8cm7p30v1HNWRXemENZ0vGPk42vkl3lRfWcrJybO1c40u6q7ywknOVY3OnGl+a3Q/lhXWcqhzbO9X40up+CC+s41Tl2N6ZxpdGt5EeBLzAmcqxgzONL21uIz0IeIEzlWMHZxpf2txGehDwAmcqxw7ONL60uY30IOAFzlSOHZxofGlyO+lhQNyJyrGHM40vTW4jPQh4gTOVYwdnGl/a3EZ6EPACZyrHDs40vrS5jfQg4AXOVI4dnGl8aXMb6UHAC5ypHDs41fjS6H7oLqzjVOXY3qnGl1b3Q3hhHacqx/bONb40u7oL6zhXOTZ3svHpLmziZOXY2tnGp7uwhbOVY2PnG5/swvrOV45NnXN8ogvrOmc5NmN8QD7lCDE+IJ9yhBgfkE85QowPyKccIcYH5FOOEOMD8ilHiPEB+ZQjxPiAfMoRYnxAPuUIMT4gn3KEGB+QTzlCjA/IpxwhxgfkU44Q4wPyKUeI8QH5lCPE+IB8ytErqvv9XqZ7a0VZVd0nys+D+o9rxgfkU45Ond1PVbq/DvJQ2tlBxgfkU45WWRXtRW9a3vJ+L/qPq9lBxgfkU45Wl9LPy9lk/2hH8Tiob7HxAfmUo9HfTCiT8I6vd2/19W63sz/e+IB8yjFRTm81FJNftz1uMAgvEKAcE9XoArfZXJyP8AIRyjFW9HcTus3612hFWaZPMnvk2PiAfMoxUibP4y0/w9s+g2zS49H9COMD8inHoJwltmrvKdSfmO4ettrn9Q5GBwE8aMWy+pZCMpD+Xm4x2V2MrouND8inHBPdNW5n2KgmzyobPfHB+IB8yjExvbYdPh49c/dWjp9wZnxAPuWYmjyBbMjtKLz9y9daxgfkU46pyUvXhrfHeYR32l3jA56gHFOTl64NTy8bXjlcTJ9YZnzAE5RjavLStaK//O2vgx/d7YpsfEA+5Wj1dxK6a9yuv/29hu7PR3f7N0M3PiCfcjS6ewlFF9T6JWv1n0W3e9hKn/5sfEA+5Wg0TS2bNzrvNodn8FZl1d337f79iebQ7nHGB+RTjk5ZluPfmz0+nr9HzojxAfmUI8T4gHzKEWJ8QD7lCDE+IJ9yhBgfkE85QowPyKccIcYH5FOOEOMD8ilHiPEB+ZQjxPiAfMoRYnxAPuUIMT4gn3KEGB+QTzlCjA/IpxwhxgfkU44Q4wPyKUeI8QH5lCPE+IB8yhFifEA+5QgxPiCfcoQYH5BPOUKMD8inHCHGB+RTjhDjA/IpR4jxAfmUI8T4gHzKEWJ8QD7lCDE+IJ9yhBgfkE85QowPyKccIcYH5FOOEOMD8ilHiPEB+ZQjxPiAfMoRYnxAPuUIuXfS/QALJOMVjA/IpxwhxgfkU44Q4wPyKUeI8QH5lCPE+IB8yhFifEA+5QgxPiCfcoQYH5BPOUKMD8inHCHGB+RTjhDjA/IpR4jxAfmUI8T4gHzKEWJ8QD7lCDE+IJ9yhBgfkE85QowPyKccIcYH5FOOEOMD8ilHiPEB+ZQjxPiAfMoRYnxAPuUIMT4gn3KEGB+QTzlCjA/IpxwhxgfkU44Q4wPyKUeI8QH5lCPE+IB8yhFifEA+5QgxPiCfcoQYH5BPOUKMD8inHCHGB+RTjhDjA/IpR4jxAfmU4wtFdb/fy3RvrSirqv+E8QH5lGNZnd1PVbq/DvIox8YH5FOORWVVtBe9aXnL+70YbRofkE85FnW9/bzmTfYnO4wPyKccS/qbCWXS2eR61/iAZyjHd8rprYZi9us24wPyKcd3qukFbnqjwfiAZyjHN4r6V2yjzfp3bUVZelYDEKIcXyuTGwvlZ3jbp5kNPTY+IJ9yfKWcJvbW3GmoS1x/ot9tfEA+5fhCfUuhLu9oV/+rteKxu32ZxeBxLMCIVvxed43bGTaq4ZLX+IB8yvGd0bXtrQ5v90E5JNj4gHzK8a3JE8iG3AovEKEc35q8dG14exzhBSKU41uTl64NTy97vHLY+IB8yvGtyUvXiv7y93EdbHxAPuVY1F3aFt01btff/l7D4y15jQ/IpxxLunsJRfcvTdQvWav/LLrdj/sPxgfkU44l9dPI7mXVv1XD413JqntVVqP7vsYH5FOOZWVZjl8t/Ph48h45xgc8QTlCjA/IpxwhxgfkU44Q4wPyKUeI8QH5lCPE+IB8yhFifEA+5QgxPiCfcoQYH5BPOUKMD8inHCHGB+RTjhDjA/IpR4jxAfmUI8T4gHzKEXKk8R1orZa6igOt9UBLfW6tTz2I3pHGd6C1WuoqDrTWAy31ubU+9SB6RxrfgdZqqas40FoPtNTn1vrUg+gdaXwHWqulruJAaz3QUp9b61MPonek8R1orZa6igOt9UBLfW6tTz2I3pHGd6C1WuoqDrTWAy31ubU+9SB6RxrfgdZqqas40FoPtNTn1vrUg+itM74DnXWVkx7prKuc9EhnXeWk5z/rUw+it874DnTWVU56pLOuctIjnXWVk57/rE89iN464zvQWVc56ZHOuspJj3TWVU56/rM+9SB664zvQGdd5aRHOusqJz3SWVc56fnP+tSD6K0zvgOddZWTHumsq5z0SGdd5aTnP+tTD6K3zvgOdNZVTnqks65y0iOddZWTnv+sTz2I3jrjO9BZVznpkc66ykmPdNZVTnr+sz71IHp3gHxpSgAAAAAAAAAAAIAsRXW/38t0756WV5TsnW4WZb1VDJ/dyhNLbXfNH7S6hWXM9y4cVJabT3ZhFfO9C98B1ewx2/nVmps9461NzRezvHe0xL1+sM6vuNdTru7vM9vlFSV7p5vFvbX1F/HEUrtd8+//tS0sY753flB5r7ae6sIqFvZON9tlVs2+XfxqzfU693vVwWwxy3vHS9zrB+v82rHfbvt9x6aWV5TsnW6214/1N8y23yBPLLVVbR/epWXM9s4OKvZI2WwVS3vT74Bua/O5dn615s9x7hjedDHLeydL3OsH6wKq7ju13Os7dmZ5Rcne6WZ3Tfb5DbJtJZ5YavvhDuFdWMZ8b3rQHlfm81Us7k2+A4atnbL2qzXfih1XOFvM8t7JEvf6wbqAYca7fT+klleU7J1sDrektr6ayF9qo7jvcGU2X8bC3vSgfX7g0lUs7p1u9vPcLWu/WnNttxUuLWZ577DE3X6wzu/xbbB9CpYtryjZO90c7v1v/E39xFIbn9cRm097YRnzvelB1bbz7KSrWNybbPZl2Pg7YPCrNSf7Nra0mOW9w769frAu4PHXjv6vHHtbXlGyd/mgrb8/llfx41LrX73vEN50GQt7k819bjSkq1jem2xW3bX5Xjcjf7Xmbt+m36MPS4tZ3jtf4nwPMQv/M7ez5RUle5cPGn0XbWJ5FT8ttfn72+ZJmy1jaW+yWdUZK8utn0z0zFI/t+vypr+x30y6nK/37vZztrSY5b3zJW78g3UBS3Pf1/KKfvix62x8vbO8ip+W2qzxsXcjs2Us7Z1t3svmyUTb3uh9Zqn11fm9qrZ/5lsnXc7Xe3f7OVtazPLe+RI3/sG6gKW572t5RT/82LX6Z8ZsZXkVPyy1fY7/Y+9G0mUs7p1ufrasDln9/PpN5/rEUmuzlwFsab6cr/bu9nO2tJjlvbMlbv2DdQFLc9/X8opSLAOIAAADDklEQVR+/LGrVRt/Bcur+H6p3S+KN2/EE2sdDtz4yW9PLLXWPNd/r0DMl/PV3t1+zpYWs7x3tsStf7AuYGnu+1pe0Y8/drf6J2/jvw8triLdO90suuO2TdktXcYXe6ebw4H1pW+3fwtPLPVTVbW3G7q9G5st58u9u/2cLS1meW+6xM1/sC5gae77Wl7RTz92zTFbf3ssriLdO93sb0I+9m7kibU+fqOy7dM4n1hq92u1/cqbLufrvbv9nC0tZnlvusTNf7AuYPwcko1T8IXlFSV7lw7a/hcrS6uY7Z1uNr+r6m058CfW+vj52za8Tyz18/80wa3Lu/m3Qe1Xa+4+2nKWI0uLWd6bLHH7H6wLGP9wPf6/saflFSV7Fw7a4dtjYRXzvdPN3cL7xFqLncL7xFKHI8ttpzr41ZqTfRtbWszy3ukSd/jBuoDxD9fkE7tZXlGyd35QtcNfMeerWNi7fNDmfVhexvdrHW1tOd3nltrNc/PBtn615sZu4V1azPLeyRL3+MG6gv6vZt3f1d7A8oqSvelBw7dHseX/PqerWNy7eND2fVhcRrp3ujn8NXTj1T6x1L3Dm67u6727hXdpMct7x0vc5wfrAoZydfPfX7Ki7o9kb3rQ8O2x6Xd1uorFtS1OePs+PLHW/mrocVW0jSeW2rei6DuytV+tud217TAfnlniTj9YV9AO/HGPfX+TFZX9X3OTdU426+fO73Df9JmldrZe6O2ptXZvxjv6mdxGZKmbz7XzqzU3H289zUH+Enf7wbqC+unx7zXX8YqKYWnJOkebxePbY+v/Xc5dam+PgT+x1vqwcodSPLHUeiud8qZ+s+Zy34rlLnHHH6xL2PxtUH40XtHjo2Sd77HsAy31mbUWOy39iaVu/3Y+iV+teV8HWCIAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADAlf0fAykOV1iIUMsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "!gracey.py -com strain_T2g.0/Y_energy_T2g.0-2.txt PNG strlt=view strc=br str=T2g,0.95,0.7\n",
    "Image(filename='Untitled.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xmgrace -nosafe -noask -hardcopy -hdevice PNG -batch /home/cam1/.ppbatchfile\r\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABXgAAAK8CAMAAABiCycpAAAAQlBMVEX///8AAAD/AAAA/wAAAP///wC8j4/c3NyUANMA////AP//pQByIbxnB0hA4NAAiwDAwMCBgYFCQkLv4+Pfx8fPq6vxjVQBAAAACXBIWXMAAAsTAAALEwEAmpwYAAAADnRFWHRUaXRsZQBVbnRpdGxlZPPPpzIAAAALdEVYdEF1dGhvcgBjYW0xmdSm6wAAABV0RVh0U29mdHdhcmUAR3JhY2UtNS4xLjI1jnpa/wAAHfRJREFUeF7t3dt22koWQFHw6HGeJZL8/6+2dUUqCuzyBlkq5nxpuywH1YazmsjgnE4AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADwu5q2bZt0cfb5xTazlP+Gzy9c0jUA1trzIFvS5jJ8cVnT9nz5zGvmG7rFO0EGYDaWNV/eZv7iXN7PFPcfN+saD8dm/ggA1i7Dc9T+aW/6tc+WXrqrDP3T3ulyw9jd/luWlyAuaYgByGmmJ6nd89X0Uu7pMqW0K+/80fS0dvkM19NdgG+6zLXsrs8uv3LqnunOH03hnT849U95pwOy2QYg43p1ILly0C9dn8NOvV0etYiw6wwA5dqH1wqmsK6eF89Pc7srxPMqAN8z/9Asp5mqvLqkcB5z7EID8M66l9IOL0Uo9fAJbztVOQ1v/0T3+qM3gHfTXLrm/ugVBpfrj9luNXNts+Ednvk+fv8bQJWa8Xlpd8m1MIGXR9cKPv/cxUvIrk9ux0+6m2un9789+GMAqtNMTcxdcx27uDJ/cajmnVYP7xq+vmvietz4Z/TfPb8N48GVYoDaXJ+1Zl6TO7Z2ZfzS9GQ1W97plzWMQe2aPqd1/DMW37r6MkDt5lce9JdbV1/61P1CsdT1q2NerwsrQ5mHw7sDk8ouv/MHVzkADqu9G87v6Cp68zR51nf5+mF3XHOZPlp946M/BaAyi+edP/H4KsHiie78Y7Tp1ROr1nppGfBGosV7/PazZV27yxTN9RtW4X38pwBUJVq85e+/uZX/eV2/JLzAu1r17yceJrO9+dPnxK5aK7zAGzk/ukbbv/QgdXPIg+9vb16tcJ5SLLzAu1q/kitt6Dq5g5tDHjxjvnnNxOI3QZ4XN3z7zBigWquXJdz8npz0Nbyd5Ig0rStpT9vljS1/8HbzzBigXv0vJuur18z/Zs/iX5H4ymUuZuYXN0y/BmKy7O7qx3IPr1cAVKZ/D8Tnk8/FP73efbg+KDGluqvtVMzumkVf3utlg8W74nqL6wzjp3cOBKjbWN7rJdemmX9xTl7/HV16m+vz3fm9acP71bqP5l/HO1hfTe5Mqddd4N2Mv3JhWcnH4R3aOvV1NL0j7Rry5e/qbbqrDDeXIvrFNgk0wFtIf2b2RXj7X1+e/pituf6orLN6Dnu+OXrU/0HpIsAb+iq8ADyZ8AJsTHgBNia8ABsTXoCNCS/AxoQXYGPCC7Ax4QXYmPACbKr77WSX9Xt+AQAAAAAAAAAAAAAAAAAAAAAAYHYGKJemhBIbjq/Om9ryttxUkJsKWtzUhrdaow3HV+dNbXlbbirITQUJ77NsOL46b2rL23JTQW4qSHifZcPx1XlTW96WmwpyU0HC+ywbjq/Om9ryttxUkJsKEt5n2XB8dd7UlrflpoLcVJDwPsuG46vzpra8LTcV5KaChPdZ6hxfnbuqdFt2dRzC+yx1jq/OXVW6Lbs6DuF9ljrHV+euKt2WXR2H8D5LneOrc1eVbsuujkN4n6XO8dW5q0q3ZVfHIbzPUuf46txVpduyq+MQ3mepc3x17qrSbdnVcQjvsxgfUE45QowPKKccIcYHlFOOEOMDyilHiPEB5ZQjxPiAcsoRYnxAOeUIMT6gnHKEGB9QTjlCjA8opxwhxgeUU44Q4wPKKUeI8QHllCPE+IByyhFifEA55QgxPqCccoQYH1BOOUKMDyinHCHGB5RTjhDjA8opR4jxAeWUI8T4gHLKEWJ8QDnlCDE+oJxyhBgfUE45QowPKKccIcYHlFOOEOMDyilHiPEB5ZQjxPiAcsoRYnxAOeUIMT6gnHKEGB9QTjlCjA8opxwhxgeUU44Q4wPKKUeI8QHllCPkPErXATIk4xmMDyinHCHGB5RTjhDjA8opR4jxAeWUI8T4gHLKEWJ8QDnlCDE+oJxyhBgfUE45QowPKKccIcYHlFOOEOMDyilHiPEB5ZQjxPiAcsoRYnxAOeUIMT6gnHKEGB9QTjlCjA8opxwhxgeUU44Q4wPKKUeI8QHllCPE+IByyhFifEA55QgxPqCccoQYH1BOOUKMDyinHCHGB5RTjhDjI/Xnz/rzv5+SJd6ecoQYH2t//v33d/n53/96/6SXJeUIMT6WPrP73yq83ecD5WVBOUKMj6s+u6vwfi78/ft3WFZerpQjxPi4+kzrn1V4/44f90X+Ny+DcoQYH2ur8F5j6ykvK8oRYnysLcP75795ufsZ2+qHbrw35QgxPtaWgf27SK3wsqQcIcbH2iq8d9ZBOUKM7531L1f4t+7pvcC6xsuScoQY3/v6869rbvcyhmVRH4Q3XeKNKUeI8b2tP+NLFrofmy3Keye8f/PLvCvlCDG+d/Vnega7fuXuvfD+84SXJeUIMb53df2VDN27067rd8KbX+VtKUeI8b2pP9frC3+Wrxq7E96/3rfGinKEGN+b+nvv0kE2vPN1CRgoR4jxvam7L1LIhvefl5KxphwhxvemisL7V3dJKEeI8b2pkvD+0V1SyhFifG8q09fB7Rd0l1vKEWJ8b+q/e79f9ya8ukuGcoQY35tav19t0eA0vKvuSjAj5QgxvjfVvV9tzu3y9+Qk4V111794yUQ5QozvXfW/mKwP6Z9/y4sOSXiXTV4dx3tTjhDje1fdU94ust2vhrz/28mu/8hwf/D1C7w55Qgxvrc1lve/dXfX/6rlurt3XoDGO1KOEON7X+M/5r64fvC3+5FbtzT+9oaku57wMlOOEON7a39Xvx8Hvk05Js3lfD636WrvMk+p/TzosjjI+IByyjHqsvvpkq6f+i9dP1ofZHxAOeUYtJdmeNJ7W97mPIX3cnOQ8QHllGMwpnRu7MJ5esbbXA9qpq+N/wvwfcrRmy7btrfhvbTTYvd8t9NerwXfHA3wJeVYaW8uNXw+zZ3DO60JLxChHCuX+SLCqDlnngYLLxChHEvNdDVh1i1kwjt/tFwGuOd///vf9RPlWGhvXsfbdpcX0vAurkcYH/Cl/w2uC8oxa/vX6K6e8nYXGm7Du7geMbyud7Y8CuA0VVcr8pq2e/nCeiDDlYckvM3iebHxAQ/cPNcdKMfKZfWu4fHNwUl4L4sXPhgfcMcY3ZvqnpQjcX2XWvfJWNh1ePvLvhPjAzKm6Oaqe1KO1PX3MkwXfZNrMlONB8YHpB5Ht6Mca8tnt7nwrrtrfMDK19HtKMfa7VvXVjFukjdYGB8wmaL7RXVPypG6eevaaRnea3fHH8EZH9D7dnQ7yjEYQzq9VGzV3zm81+5Ovwzd+ICy6HaUo9cO75xori8gW1xxmMLbveRhfcXX+ODtFVf3pByjvqlt/4vOx08XL+edwjv++xP9oeOXjA/e2g+i21GOUdu2y6u7t1d6s4wP3tZPnuqOlCPE+OAtTdH9SXVPyhFkfPB+YtHtKEeI8cF7iUe3oxwhxgfvY4pusLon5QgyPngTT4tuRzlCjA/ewFOj21GOEOODyk3RfV51T8oRZHxQs1dEt6McIcYHtXpVdDvKEWJ8UKVXVvekHEHGB/V5bXQ7yhFifFCXFz/VHSlHiPFBPabovri6J+UIMj6oxGbR7ShHiPFBBTaNbkc5QowPDm6K7nbVPSlHkPHBkf1GdDvKEWJ8cFS/Fd2OcoQYHxzSb1b3pBxBxgfH87vR7ShHiPHBsfzyU92RcoQYHxzHFN1fru5JOYKMD45hP9HtKEeI8cHezcndSXQ7yhFifLBne4xuRzlCjA/2ap/JHShHiPHBHu05uh3lCDE+2J2dR7ejHCHGB7tygOh2lCPE+GA/DlLdk3IEGR/sw3Gi21GOEOOD3zYl9zDVPSlHkPHBbzpidDvKEWJ88FuOmdyBcoQYH/yGI0e3oxwhxgebO3h0O8oRYnywqQqi21GOEOOD7VRS3ZNyBBkfbKOe6HaUI8T44NWm5FZT3ZNyBBkfvFKN0e0oR4jxwavUmdyBcoQYH7xCzdHtKEeI8cGz1R7djnKEGB88z5zcqqPbUY4Q44PneJ/odpQjxPgg7p2SO1COEOODoLeLbkc5QowPAt4yuh3lCDE++JEpue9Y3ZNyBBkfFHvz6HaUI8T4oMjbJ3egHCHGB98mujPlCDmP0nVgRXQHkvEMxgdfmZv77tFdUI4Q44MHrs0V3RXlCDE+yFokV3NvKUeI8cENzf2ScoQYH6xI7rcoR4jxwUx0v005QowPeqJbRDlCjI+3NydXdL9POUKMj3e2aK7oFlGOEOPjTUluiHKEGB/vR3PjlCPE+Hgni+RqbohyhBgfb0Jzn0o5QoyP6knuCyhHiPFRNc19EeUIMT6qJbkvpBwhxkedRPe1lCPE+KiP6L6ecoQYHzWZkqu6r6YcIcZHJRbNFd3XU44Q4+PwJPcXKEeI8XFkmvtblCPE+Dgmyf1dyhFifByO5u6AcoQYH0ciuXuhHCHGxy/7SBfu0NxdUY4Q4+MXfUzSL6wskqu5e6EcIcbHr5mzez+9mrtXyhFifPyWVXdvyiu5+6YcIcbHL0m6uyiv5h6AcoQYH78jze5HX17JPQrlCDE+fkda3U+aeyDKEWJ8/Io0uj3NPQ7lCDE+fkXa3F56EPulHCHGx69Im9tLD2K/lCPE+NjWeBk3bW4vPZb9Uo4Q42Mj84/O7nZXeQ9EOUKMj5dbJXf88Vma3N7629gz5QgxPl7ptrijtLm91RHsmnKEGB+vcTe5g7S5vfQg9ks5QoyPp/uiub20ub30IPZLOUKMjyf6TnJHaXQ/dPdQlCPE+HiCZXC/bm4vre6H8B6KcoQYHxFJcb+V3FGaXd09FOUIMT5+5sfBnenukSlHiPFRLFrcie4emHKEGB8FnpXciewelnKEGN+B/F6elsF9UnNnv7crApQjxPgO4peeGybFfW5yOTDlCDG+Q5izu116BZdHlCPE+I5g1d0Xl3cdXMUlTzlCjO8Aku6+prxJcCWXh5QjxPj2L83ux5PLK7iUU447msv5fG7T1d7lOjTj27+0uh/PC6/g8kPKkddl99MlXT/1X5o/Nr7dS6PbSw8qJ7lEKEdWe2mGJ7235W3OwnskaXN76UHftwiu5vJjypE19nbZ2MnZM95DSZvbSw/6wrq2g/QYKKAcOdO13fY2vJd2uXjzZfYmbW4vPeiONLad9Bj4AeV4pL251NBcVjU2vr1LkztKD7uhtryScjxyOTfrhea8fhpsfLuXJreXHnQluGxBOR5ouh+xrXQLwnsoaXN7yTHr2g6SQ+CZlOO+9uZ1vG135UF4DyVtbm/8Whrbzuqb4TWU4562fyHv6ilvd6FBeA8mbW4vTW76TfBaynFH03YvX1i/qmG48rAK79q8zn6k0f0kuGxOK77vsnrX8GX42DPeY0mr+9GHNz0KtqQcj6zepdaMLy0T3sMYntam2U1/tAabU46Hlu9SGy76rv+mYHw7NV9MyJQ3PZb9qf5OUo6Hls9uhfcAVsGdLyjo7pHUe08tdqQcD92+dc2lhh1KcttJjqj3P+bazPdUZfdWsifleOjmrWsn4d2LNLWD9KiFmv4zrtaiutdKHd/NppQja3wtQzO+qGHVX+H9RWlnR+lhHNQqUFOkju92U8qR0w7vnGiuLyBbXHEQ3q2lmR2lh3F4SaA66SFHlO7pQznyupeRndvL9KsamtXLeYV3G2lnB+lR1CQN1EcV4U239Ek57mjbdnl19/ZK78D4nizt7CA9ijqleeqlBx1PuqMP4Y0yvmdIOztKD6N2aZ566UGHk26opxwhxvdzaWcH6VG8kbROvfSgw0k31FOOEOMrlHZ2kB7Fe0rr1EsPOpx0Qz3lCKlzfM9+rKedHaWH8ebSOI3Sw44m3U+vznJspr7xzY+M9Aul0syO0sNgtgjTVXrQ0aT7GdRXjk3VNr7VYyP94l1pXG+k3wA5q4ffJD3ocNIN9Worx8YqG1/y4Ei/vJLG9Ub6DfCF5PE3SA86nHRDvcrKsbW6xpc+ODKP+bSug/Qo+In08ddLDzqcdEO9usqxuarGlz42Phb/Rk5G+u0QlT7+Pirobm5TwhtU1fjSx8ZHLrzpN8HzpI+/jxrCm91VVeXY3oHHlxY1/XcaBul3wSulj78qHoDpnj6EN+oY40sTm5c+NnrpHwUvVeXjL9lUt61jlGO39ja+NKaPpN+bPjh66UHwWlU+/Fab6re1t3IczNbjS+P5Tekfk5U8OAbpQfBqNT74Fv9JDfvauhyVecr40kxGpH/2960eGlfpYbCB+h53yX9RTynH2/rfOQ3fq6Vn8ESL2l6lBwE/tfjPSXgjnhTe9I/9HWlze+lBwBMIb0hN40ub20sPAp6gpnL8gprGlza3lx4EPEFN5fgFVY0vje6H7sJrVFWO7VU1vrS6H8ILr1FVObZX1/jS7OouvEZd5dhcZePTXdhEZeXYWm3j013YQm3l2Fh945NdeL36yrGpOscnuvBadZZjM8YHlFOOEOMDyilHiPEB5ZQjxPiAcsoRYnxAOeUIMT6gnHKEGB9QTjlCjA8opxwhxgeUU44Q4wPKKUeI8QHllCPE+IByyhFifEA55QgxPqCccoQYH1BOOUKMDyinHCHGB5RTjhDjA8opR4jxAeWUI8T4gHLKEWJ8QDnlCDE+oJxyhBgfUE45QowPKKccIcYHlFOOEOMDyilHiPEB5ZQj5DxK1wEyJOMZjA8opxwhxgeUU44Q4wPKKUeI8QHllCPE+IByyhFifEA55QgxPqCccoQYH1BOOUKMDyinHCHGB5RTjhDjA8opR4jxAeWUI8T4gHLKEWJ8QDnlCDE+oJxyhBgfUE45QowPKKccIcYHlFOOEOMDyilHiPEB5ZQjxPiAcsoRYnxAOeUIMT6gnHKEGB9QTjlCjA8opxwhxgeUU44Q4wPKKUeI8QHllCPE+IByyhFifEA55QgxPqCccoQYH1BOOUKMDyinHCHGB5RTjhDjA8opR4jxAeWUI8T4gHLKEWJ8QDnlCDE+oJxyhBgfUE45QowPKKccIcYHlFOOEOMDyilHiPEB5ZQjxPiAcsoRYnxAOeUIMT6gnHKEGB9QTjlCjA8opxwhxgeUU44Q4wPKKUeI8QHllCPE+IByyhFifEA55Zg1l/P53H6x2vafNtOnxgeUU45JV9RPl0ernxX+bHBzPk/lNT6gnHKM2kszPL1dlTdZvQzFbc/T2IwPKKcco7G357mpmdVm+uJ88cH4gHLKMZiu4l6fzGZW5y8KLxCgHGttepG3N662U3CFFwhQjrXxKm5iXG3Ow4/VrnU2PqCccqw03Q/TbsyrbV/e5vqs2PiAcsqx1N6+jve0Wu1fxbs4yPiAcspx1T2fvb5EN7/aXW1YPCseXuU7m9cBlrTijqZt+8g+Xu3fUeGda0CAcqxdMu8aXqw2n83t0uuda8DPKcdadykhXbuuDj9WW5Q3cyzAF5QjccmFd1odg3t9Y3HuWIDHlCOxfuvaZFhdvHNt+mA+AuC7lCPx6K1rc3ivBZ6PAPgu5Ug8euua8ALPoByj8bUMzfjyhekXQK5Xpxc3XFzjBX5OOQbt8HOz5jKUtR1+fJasTr8Cff71kMYH/IByDLoXjJ3by/SmtGZ4bpus9guX9nJ975rxAeWUY9K2139K7dP4cbI6vJPt+pnxAeWUI8T4gHLKEWJ8QDnlCDE+oJxyhBgfUE45QowPKKccIcYHlFOOEOMDyilHiPEB5ZQjxPiAcsoRYnxAOeUIMT6gnHKEGB9QTjlCjA8opxwhxgeUU44Q4wPKKUeI8QHllCPE+IByyhFifEA55QgxPqCccoQYH1BOOUKMDyinHCHGB5RTjhDjA8opR4jxAeWUI8T4gHLKEWJ8QDnlCDE+oJxyhBgfUE45QowPKKccIcYHlFOOEOMDyilHiPEB5ZQjxPiAcsoRYnxAOeUIMT6gnHKEGB9QTjlCjA8opxwhxgeUU44Q4wPKKUeI8QHllCPE+IByyhFifEA55QgxPqCccoQYH1BOOUKMDyinHCHGB5RTjhDjA8opR4jxAeWUI8T4gHLKEWJ8QDnlCDE+oJxyhBgfUE45Qs6jdB0gQzKewfiAcsoRYnxAOeUIMT6gnHKEGB9QTjlCjA8opxwhxgeUU44Q4wPKKUeI8QHllCPE+IByyhFifEA55QgxPqCccoQYH1BOOUKMDyinHCHGB5RTjhDjA8opR4jxAeWUI8T4gHLKEWJ8QDnlCDE+oJxyhBgfUE45QowPKKccIcYHlFOOEOMDyilHiPEB5ZQjxPiAcsoRYnxAOeUIMT6gnHKEGB9QTjlCjA8opxwhxgeUU44Q4wPKKUeI8QHllCPE+IByyhFifEA55QgxPqCccoQYH1BOOUKMDyinHCHGB5RTjpA6x1fnrirdll0dx2JXdW5wM3WOr85dVbotuzoO4X2WOsdX564q3ZZdHYfwPkud46tzV5Vuy66OQ3ifpc7x1bmrSrdlV8chvM9S5/jq3FWl27Kr4xDeZ6lzfHXuqtJt2dVxCO+zbDi+Om9qy9tyU0FuKkh4n2XD8dV5U1velpsKclNBwvssG46vzpva8rbcVJCbChLeZ9lwfHXe1Ja35aaC3FSQ8D7LhuOr86a2vC03FeSmgoT3WTYcX503teVtuakgNxUkvM+y4fjqvKktb8tNBbmpIOF9ljNAuTQlAAAAAAAAAAAAQJHmcj6f23T1IPInn6ymn7a9SzOv7NW3dtevLD/bq9vzzq8ud1PZfdV+fnq5OWiPfrKd49xZe9Ccu8f55XzIaeVPPllND+oePJ/2H6v0xO+studDvJD95rzzq+vd1HVfHWY36YnnV9PtpJ/zwDDL0+mQ48qffLKaHtRchv9j3v3/06Qnnl/93M75COFNzzu/muymrvvq0j0X7J40rg7aox9t5zB31i5cxr85tLd/r9i//Mknq+lBx/ir3un2xPOrTf8kcXHATqXnnV9NdlPVfTVdQvn8f5a9t+lH2znMnbUL86P8CP/xpvInn6wmnzaH2ee3dtc5RHgz551fXeymrvtquvbZZi6e7sxPtnOcO2sPro/y/T8abuRPPllND7pcDvJ3ofTE768eIby5886vLnZT1301/7V99/+p/Wg7h7mzduH6d4np7xEHkj/5ZDU9qL/+f4SHSHri91ePEd7b886vLnZT4X3VOUB4cyeeX523c5g7axeyTzSOIn/yyWryads/Pnb/0D/dnPiD1SPcdbnzzq+uPqruvhrWF5/sUf7E86vTdo5zZ+3CvWEeQv7kH4f3U9M/Rnb/k+XbE7+3eoS7Lnfe+dXVbqq7r/pP9r6f/InnVxfbOcidtQt3hnkM+ZP/Mryn7u9K1+W9yp54dvUId13uvPOr6W7quq9O3Yb2/vfx/InnV9fbOcKdtQt3hnkM+ZP/TnhPzXn3282feG71CHdd7rzzqze7qeq+6raz+zLlTzy/mmznAHfWLuSHeRD5k/9WeLv/a9754z9/4rnVI9x1ufPOr97upqb76nM3+/+7eP7E86vpdvZ/Z+1CfpgHkT/574W32f3jI3/iudUj3HW5886v3u6mpvvq1O6/u/kTz6/ebGf/d9YuLF/Fc7h55U8+Wc0ftHwY7VT+xHOrt6nan9x551czu6novjrE7zPKnXh+NbOd3d9Zu3B9lB9wXvmTT1bzB6Wf7VD+xHOrmVTtTu6886uZ3dRzX2VCtUOZE8+v5raz+ztrF5rFMFdfOIL8ySer+YPSz3Yof+K51Uyqdid33vnVzG5uV3Ymt4/carP7FzT0bk88v5rdzu7vrH2YfsPF/l9cmJE/+WQ1f9ABrqzkTzyzmknV/mTOO796u5tq7qtrqHa+ofTE86vZ7RzgztqFaYj7f3FhRnLy4/8kq+tPp/c0pj+M3aFv7W5YSlO1Qz/YTW331TVUe/9FXj/YzoHurH0Yprh+9/VhrE6+nd41k2xp+Wl3TP/ZER4e39pd//E1W/tVvJva7qvuJa6T+ft2qng7h7qz9qF74d1hL4gvT/76QpZkS4tP+0fLYf51ku/srh0f+/u/A0t3U9t91b2rq6I7a72dY91ZO3HoXxq/PPnrR8mWFp92/zDU8kv79q3dHUbpbtxXv6h0O8e6swAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAD26f9kV9czs6Y+fgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "!gracey.py -com strain_xx_yy/Y_energy_xx_yy.txt PNG strlt=view strc=br str=c12,0.95,0.7\n",
    "Image(filename='Untitled.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do not need to use energy derivatives: we can also use stress. Let us take first derivatives of stress."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pgname: Oh\r\n",
      "target: stress\r\n",
      "unstrained_path: static/\r\n",
      "xtal: static/POSCAR\r\n",
      "order: 1\r\n"
     ]
    }
   ],
   "source": [
    "%cat inputs_stress_1.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the command line, we add the variable *sigi*, which tells us which component of the stress to take the derivative wrt. Of course, you get all of them for free, and you can select whatever you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "!strain_derivative.py -i inputs_stress_1.txt -dpath strain_A1g    -sigi A1g\n",
    "!strain_derivative.py -i inputs_stress_1.txt -dpath strain_Eg.0   -sigi Eg.0\n",
    "!strain_derivative.py -i inputs_stress_1.txt -dpath strain_T2g.0  -sigi T2g.0\n",
    "!strain_derivative.py -i inputs_stress_1.txt -dpath strain_xx     -sigi xx\n",
    "!strain_derivative.py -i inputs_stress_1.txt -dpath strain_xx     -sigi yy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note that this same stencil can also be used to obtain the second derivatives of stress."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!strain_derivative.py -i inputs_stress_1.txt -dpath strain_Eg.0  -sigi A1g -order 2 -scaled 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us inspect one example, for the case of *A1g*. Here we see that the stress does not agree well with the energy. We would need to study the convergence parameters to see why (see my dft_runs directory). The example of ThO2 below shows nice agreement between energy and stress."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xmgrace -nosafe -noask -hardcopy -hdevice PNG -batch /home/cam1/.ppbatchfile\r\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABXgAAAK8CAMAAABiCycpAAAAQlBMVEX///8AAAD/AAAA/wAAAP///wC8j4/c3NyUANMA////AP//pQByIbxnB0hA4NAAiwDAwMCBgYFCQkLv4+Pfx8fPq6vxjVQBAAAACXBIWXMAAAsTAAALEwEAmpwYAAAADnRFWHRUaXRsZQBVbnRpdGxlZPPPpzIAAAALdEVYdEF1dGhvcgBjYW0xmdSm6wAAABV0RVh0U29mdHdhcmUAR3JhY2UtNS4xLjI1jnpa/wAAIABJREFUeF7t3dt2m0oWQFHLI+M8S0rn/3+1DboApZJ12YCqijkf+liYGNhyrxCM5K8vAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA+JRDuuCx/eHHdNH0cWYFAE4Ou93kwdjwibH9z1rHtKqTL7M/7n66+7PWsAiAk760w8Npd/PdPOY+sR9/mZ8vuu/+ezz/F4CL8wlu+vgiPavtdIXN1PQ4+jL76xrjrw3ATyG/psX8eTBq6r3A3p7uTk+c98MqrjYA3OqSefl4P07t5KLtxSTTg58/OHziODpVzp8eA2zaOLyHcSSPmSsNd64zfP2c1w5fZnoKnfkqANs2Du9ELrF3rvsed6Pcpj9my391gO26l8bclYb7FxpG4Z18we7Bbb8BGnHounhzf+0j98KbudLQncveLPzRL7wf3twfAWjA/tg1995F2PvuhTfzdc6nr+lr0o79jQvCC2zN/nzf1sv/tr8T3tyVhl33tbvT6h/H60YOp+1NwnttrfACzdpfspe7HHBK5dT1k3fCm7nS0H+d/lJGX99zeS+3oA1fplvr/KHwAu0aOpm5gWvU26vrJ++EN3Pe3K15KepQ3tOFhnF4R1UWXqBZw4t0u7cFm3zq63RRNjV8Mhve3JWGSUQvEe7uJLssOH/UnXTfBhqgLblOPikf3tyVhsma5wfDa91Gn+zL2y0/7O7cgAZQvUDe8uHNnadO1jxfSx6u5k4+ef0J3J23dgCoXjaez8mGN3sGPV2zD+/lAu/p8fDJ6/ugn/IM0J5cPJ+UDW/uSkMuvKcT27FhhU72iwO0YPf+iWW2jbkrDclW+geT5vZGq19WAmjR7tdLqWkcO9dP5sKbvdKQC+/oNonu8fhuiU7uawO0oSvccI6aNngc3IvrJ3NxPOYzfhjnPf1j6eOvd16/DFCN8a2zXzfvkzM6L70aPpkp5p0LBJM3fEzPsbNfRneBdvVvTNZXbj/cabDPnrcmMuFNrzRcf9Q2elXcIY3qzZfpbiXTXaBd3cnorrvKOtSu+3C6UlYmvMmVhv7ybf/R5Pf7JFVPv4zrDEDrzuUdarffX9845zf9n0sCmVxp6M6mzwv2lwsMN91NwuuVE8AGnF8tNq7dw/Du+/dO77o6/m1r6ZWG8blrt5VDd159czI7Dm/3ZSe/vw2gUeOfmXUehjfv5m129pPH3Y/mfo3q4fhgBYBmvRleAN4lvAArE16AlQkvwMqEF2BlwguwMuEFWJnwAqxMeAFWJrwAq+reMOHotbsAAAAAAAAAAAAAAAAAAAAAANCgP3/SJRCzA3hdmhJeseL42tzUmtuyqaC3NvXnrRPetzb1ns9sasWttmjF8bW5qTW3ZVNB72zqve6+tak3fWZTK261RSuOr81Nrbktmwp6Y1NvdvedTb3rM5tacastWnF8bW5qzW3ZVNDrm3q3u29s6m2f2dSKW23RiuNrc1Nrbsumgl7f1LvdfWNTb/vMplbcaotWHF+bm1pzWzYV9PKm3j7hfX1T7/vMplbcaotWHF+bm1pzWzYV9Oqm3u/uy5sK+MymVtxqi9ocX5tH1ehhlXtUge4WfFQRwjuXNsfX5lE1eljFHlWku+UeVYjwzqXN8bV5VI0eVqlHFepusUcVI7xzaXN8bR5Vo4dV6lGFulvsUcUI71zaHF+bR9XoYRV6VLET3lKPKkh459Lm+No8qkYPq8yjCna30KOKEt65tDm+No+q0cMq8qii3S3zqMKEdy7GBzfC3W3U9/ChcoQYH6R0N+P74vRQOUKMDxK6e+ua3Ut6lSPE+GBKd29NutuXVzlCjA8mdPdW0t2uvMoRYnwwpru30uz+UI4Y44Mx3b2VVvdbeKOMD0ac8N5Ko9tTjhDjg4HuZqTN7SlHiPHBle7mpM3tKUeI8cGF7malze0pR4jxwZnuZqXJPVGOEOODk1m7O3pXg+qlze0pR4jxQW++7l7jlH6iUqPcDpQjxPigM1t3J3lKP1mlyRFdKEeI8cHXUt1to7zJIZ0oR4jxwXLdbba8yhFjfLBgd5sob3pI38IbZXww3xs0pH36biK8t4elHEHGB0ue8DZZ3i/lCDI+Nm+27qZ9OklXqtLNESlHiPGxdfN1t+Hwjo7t9FA5QoyPjZuxu02HtzM6FuUIMT62bfHuNlXegXKEGB+bNmd375Q3XakNyhFifGzZvN0VXkYOx93ueJgsuT4yPjZs5u4KL4Of7HaOl8f7nwXCCzO+cOIsbW4vXakNyvHA8bg/xfZc3mN3Aiy8MPsJb7a86SqNUI7f7c+9/Tnn3Z8WdB8LL8zfXeHlrDvf7RzGtRVeWKC7mfKmK7RCOX53vbQrvDCySHdvypt+uhnK8SThhcFC3c28q0GblONJu2FSwsvmLdXdr5t3NWiTcjznMNxPJrxs3mInvGdNR7ejHM85nm9q6EyuOkxcV4GWLd3dFmnFG/ZDa53xsnW6G6YcTzkOFxqEl43T3TjleMZh3F3hZdN0dwbK8YTLy9fOhJcN0905KMdjSXeFlw3T3Vkox0P70Q0NPeFls3R3HsrxyNDdS26Fl63S3ZkoxwNDd69vhi68bJTuzkU5fre/veV5P3pXdONjQ3R3Nsrxu/Pvn+icTnOnFTY+tkN356McIcbHZujujJQjxPjYCt2dk3KEGB8bobuzUo4Q42MbdHdeyhFifGyC7s5MOUKMjy3Q3bkpR4jxsQG6OzvlCDE+NkB3Z6ccIcZH+5zwzk85QoyP5unuApQjxPhone4uQTlCjI/G6e4ilCPE+Gib7i5DOUKMj6bp7kKUI8T4aJnuLkU5QoyPhunuYpQjxPhol+4uRzlCjI9m6e6ClCPE+Pj6+zddcnH/MzXQ3QUpR4jx8fXvX7rk5O+///6XLquIE94lKUeI8fH3v/9yJ7Y/2f2v5vDq7qKUI8T4+Ans7Slvn92aw6u7y1KOEOOjC2y67OtvfyZcb3h1d2HKEWJ8m/e/u2e295ZXQHeXphwhxrd5//UXFdKlnXrDq7uLU44Q49u6v//9rzvnzRX2zuLy6e7ylCPE+Lbu3393rvLWG17dXYFyhBjfxv3t2tqd8mbuKLsJ7/+6qxL/So+x7q5BOUKMb+P+1xW3u3/h9o6ym/D+++/f39PdvT8yoS6D7q5COUKMb+NOwe1qelvSJLz/hnXTJBdEd9ehHCHGt239Ce/pWsNtSacL/3dpc3bdUujuSpQjxPi27fJDte4cdvqZ08JRYa9r/C97XaIMursW5Qgxvk3rf7TWyZ7y3glvd0V4WFwU3V2NcoQY36Z195Kd5E5574R39FFhdHc9yhFifFv2d7hmkHsrspvwnh+VGl7dXZFyhBjflvVv0zBIr9xOw/tvFN400UXQ3TUpR4jxbdl//7o7ck+68iZ3lE3D213a7Vf4W+YJr+6uSjlCjG/D/pfcLZae8k7Deynv33//bm/5/TzdXZdyhBjfhg0/Wut0p7zjxzfh7ctb6kuGdXdlyhFifNt1vZfs5PbHa+mCIs90T3R3bcoRYnzb9W96Tbc/nx0vSMNb5qluT3dXpxwhxrdZNz8j68I7aev08fnVxSXS3fUpR4jxbdbN727vb2wYL5iGt7+82ystwLr7AcoRYnxbdfuGZP3bjo2u4/Zvvpt8+qyoq726+wnKEWJ823R56cT1jHZ4LcVp0enG3v4k97LKuLwFXXbQ3Y9QjhDj41n/+isN5/ymr3L7GN39DOUIMT6e83dy1+/0FuAP0t0PUY4Q4+M56R0P4wefo7ufohwhxsdT/iU3n5VxqUF3P0Y5QoyPpyQ/Tyvjtgbd/RzlCDE+njJ+QfHfQl5LobsfpBwhu7N0OUyc3iDnzPnuhknGHIyPJ13v9C0ju7r7WcoRYnw87W9JrxfW3c9SjhDjo0q6+2HKEWJ81Eh3P005QoyPCunuxylHiPFRH939POUIMT6qo7sFUI4Q46M2ulsC5QgxPiqju0VQjhDjoy66WwblCDE+6qK7ZVCOEOOjKk54C6EcIcZHTXS3FMoRYnxURHeLoRwhxkc9dLccyhFifFRDdwuiHCHGRy10tyTKEWJ8VEJ3i6IcIcZHHXS3LMoRYnxUQXcLoxwhxkcNdLc0yhFifFRAd4ujHCHGR/l0tzzKEWJ8FE93C6QcIcZH6XS3RMoRYnwUTneLpBwhxkfZdLdMyhFifBRNdwulHCHGR8l0t1TKEWJ8FEx3i6UcIcZHuXS3XMoRYnyU6o/uFkw5QoyPQulu0ZQjxPgok+yWTTlCjI8i6W7hlCPE+CiR7pZOOUKMjwLpbvGUI8T4KI/ulk85QoyP4uhuBZQjxPgoje7WQDlCjI/C6G4VlCPE+CiL7tZBOUKMj6LobiWUI8T4KInu1kI5QoyPguhuNZQjxPgoh+7WQzlCjI9i6G5FlCPE+Piw78sHulsT5QgxPj7o++JLdyujHCHGx8dcs9vR3booR4jx8SmT7v6UN/08JVOOEOPjQ5Lu9pcbqIZy/O5w3O2Oh+Hh7ufhfvi08fEZaXa/lbcqyvGrn+x2jteHP9HdXx8aH5+SVvdbeKuiHL85dme3++6k9/Sw62532ns95zU+PiKNbi9diXIpxy/251Pbn3PeLrX73Xlal/8aHx+SNreXrkS5lOMXl6u5h92uu8x76P93/IHx8Rlpc3vpSpRLOX5xvZZ7Cu/pf7/6EF8+cVkD1pQ2t5euRLmU4xmX8J4fXq85GB8fkSb3LF2NYinHM/rQDue5w7mv8fERaXJ76UqUSzmecOjvahBeipE2t5euRLmU4wmn28iEl2Kkze2lK1Eu5Xhsf8psNrxTl8/DstLm9tKVKIlWvOp4urshG97LElhVGt1v3a2Kcjx0ON9VJrwU409a3W/hrYpyPHJ5+droJrLhxjLj4xP+/Lktb7oOJVOOB67dHZ3negEFH9W/67nu1kw5frcf3hDn6+glw5Tg/NsmdLdiyvGrobuHy+28X92p72Wp8bG64bf8yG61lOM3Q3f7N0M/B3f0hrzGx9qS364mulVSjl/sk9vuTsUd/ZDN+Fib32rZBOX4xfn3T3ROF3V/mns8jn/3j/GxLt1tg3K8Zn+4/FytZ3ysSncboRwhxseadLcVyhFifKxId5uhHCHGx3p0tx3KEWJ8rEZ3G6IcIcbHWnS3JcoRYnys44/uNkU5QoyPVehuY5QjxPhYg+y2RjlCjI8V6G5zlCPE+Fie7rZHOUKMj6W5vNsi5QgxPhamu01SjhDjY1my2yblCDE+FqW7jVKOEONjQS4zNEs5QoyP5ehuu5QjxPhYjOw2TDlCjI+l6G7LlCPE+FiI7jZNOUKMj0W4vNs45QgxPpagu61TjhDjYwGy2zzlCDE+5qe77VOOEONjdrq7AcoRYnzMTXe3QDlCjI+Z6e4mKEeI8TErtzNshHKEGB9z0t2tUI4Q42NGsrsZyhFifMxHd7dDOUKMj9no7oYoR4jxMRfd3RLlCDE+ZqK7m6IcIcbHPHR3W5QjxPiYhe5ujHKEGB8zcPvu5ihHiPERp7vboxwhxkeY7G6QcoQYH1G6u0XKEWJ8BOnuJilHiPERo7vbpBwhxkeI7m6UcoQYHwFuZ9gs5QgxPt6nu9ulHCHGx9tkd8OUI8T4eJfubplyhBgfb9LdTVOOEOPjLS7vbpxyhBgf79DdrVOOEOPjDbK7ecoRYny8zOkuyhGzO0uXwz26u2mSMQfj40Wyy5dyBBkfL3G6S085QoyPV+guJ8oRYny8QHY5U44Q4+NpTne5Uo4Q4+NZustAOUKMjyfJLiPKEWJ8PEd3GVOOEOPjKbrLhHKEGB9PcHmXhHKEGB+P6S4p5QgxPh6SXW4oR4jx8Yjucks5QoyPB3SXDOUIMT5+p7vkKEeI8fEbP1YjTzlCjI9f6C53KEeI8XGf7HKPcoQYH/c43eU+5QgxPu7QXX6hHCHGR57s8hvlCDE+cpzu8jvlCDE+MnSXB5QjxPi4Ibs8pBwhxkdKd3lMOUKMj4Ts8gTlCDE+Jpzu8hTlCDE+xnSX5yhHiPExIrs8STlCjI8rp7s8TTlCjI8L3eV5yhFifJzILq9QjhDjo6e7vEQ5QoyPjuzyGuUIMT6c7vI65QgxPnSX1ylHiPEhu7xOOUKMb+uc7vIO5Qgxvo3TXd6iHCHGt22yy3uUI8T4tszpLu9SjhDj2zDd5W3KEWJ82yW7vE85Qoxvs3SXAOUIMb6NcpmBEOUIMb5t0l1ilCPE+DZJdglSjhDj2yCnu4QpR4jxbY/uEqccIca3ObLLDJQjxPg2xukus1COEOPbFt1lHsrx2P44fHjY7XaH/fWx8W2K7DIT5XikS+3l4/3u2HV4dy2v8W2I011moxy/2x9H4e27+2N3La/xbYfuMh/l+N2+P+U9PzjuDv1/D+cAG992yC5zUo6HhvBePthfT3mNbyN0l1kpx0PX8I4LfDr1Nb5tkF1mphwPZcN7/sj4NqDLru4yK+V4SHi3TXaZn3I8NA7v5dKu8G6F7LIE5XhofKJ7vpthCO/UeT2aobvMRCteNL6d7HTKu7/OzfjaJrssQzkeGsLblfd4OByP7mrYBKe7LEU5HhqFt3/98HHvPt5N0F0WoxwPjcObLkk/Qztkl+Uox0O34b1eaTC+dukuC1KOh27Ce7y+VYPxtcplBhalHA+l4R1u5zW+Vukuy1KOh5LwDj9Z+zK+RskuC1OOhy53714eDb+Pwvia5HSXxSnH77r7xzqnn6btD7vjKMLG1yLdZXnK8YL94XI3w4XxtUZ2WYNyhBhfW7rs6i7LU44Q42uK7LIS5QgxvobILqtRjhDja4fush7lCDG+Vsgua1KOEONrhO6yKuUIMb4myC4rU44Q42uAe8hYnXKEGF/9ZJf1KUeI8dVOdvkE5Qgxvop8pwu+dJcPUY4Q46vE98VkqezyIcoRYnxVuGZ3ml7d5VOUI8T4ajDp7rW8ssvnKEeI8VUg6e6pvO4h45OUI8T4ypdm97srr+zyUcoRYnzlS6v7Q3b5MOUIMb7ipdHt6S6fpRwhxle8tLm9dCVYl3KEGF/x0ub20pVgXcoRYnzFS5vbS1eCdSlHiPGVLk3uWboarEo5QoyveGlye+lKsC7lCDG+4qXN7aUrwbqUI8T4ipc2t5euBOtSjhDjK17a3F66EqxLOUKMr3xpdL91l49TjhDjK19a3W/h5eOUI8T4yvcnza7u8nHKEWJ8peveD0d3KY1yhBhf4U7vQ6a7FEY5QoyvaMPbP8ouRVGOEOMrWPpLJkSXYihHiPEVK80uFEQ5QoyvULJL0ZQjxPiKJLsUTjlCjK9EskvplCPE+Moju5RPOUKMrzi6SwWUI8T4CiO7VEE5QoyvLLpLHZQjxPiKIrtUQjlCjK8gTnephnKEGF85dJd6KEeI8ZVCdqmJcoQYXxm67Oou9VCOEOMrgexSG+UIMb7Pk13qoxwhxvdpskuNlCPE+D5MdqmScoQY30fJLpVSjpDdWbqcFbjKQH0kYw7G9zGyS8WUI8T4PkR2qZpyhBjfR8gulVOOEOP7ANmlesoRYnyrk10aoBwhxrcy2aUJyhFifKuSXRqhHCHGtyLZpRnKEWJ8q5FdGqIcIca3EtmlKcoRYnyrkF0aoxwhxrcG2aU1yhFifMuTXdqjHCHGtzRXGWiRcoQY37JklzYpR4jxLUl2aZVyhBjfcmSXdilHiPEtRXZpmXKEGN8yZJe2KUeI8S2gr67s0jLlCDG+2ckuG6AcIcY3M9VlE5QjxPhmJbtshHKEGN+MZJfNUI4Q45uN7LIhyhFifDORXTZFOUKMbxayy8YoR4jxzUF22RrlCDG+ONlle5QjxPiiXGVgi5QjxPhiZJdtUo4Q4wvoqyu7bJFyhBjf22SXDVOOEON7j+qybcoRYnxvOFVXdtkw5QgxvpepLihHjPG9SHXhSzmCjO8lsgs95QgxvhfILpwpR4jxPculXRgoR4jxPUd2YUw5QozvGaoLU8oRYnyPyS6klCPE+B6RXbilHCHG9yuXdiFLOUKM7xeyC3coR4jx3aO6cJ9yhBhfnuzCb5QjxPgyVBceUI4Q40udqiu78BvlCDG+KdWFZyhHiPGNqC48STke2x+HDw+73e54uD42vguXGOB5yvFIl9rh4+P+6+u4u5bY+HqqCy9Rjt/tj6Pw7s/F3e0u57zG9+USA7xMOX637095zw+O5+AOS4xPdeF1yvHQKLPCO+USA7xl4+V4xji8p4+Et6e68KZNl+M540sNp4u8u93+vGS741NdeN92y/G0IbzdKe+x6++lu5sdn+xCxFbL8YJRePddeY/dLWVnmxyf6kLQJsvxmlF4+6sN13vJvvpT4LHhE806VVd24TXba0XQJLzdOe/w+onN/b2lujCHrZXjDePwHo+nyw2Xx22O7ztdcKK6MJM2yzGrUXj7H6uNy9ve+L4vkuUuMcB82ivH7IbwHk7B7cp7/vlaa+O7ZneS3nN0VRdm0lo5FpB95dr5B2yNjW/S3VN5L9FVXZhPY+VYwm14Rx+cP9OGpLvf36ILi2irHIvYTHjT7H735U3XAsKaKscyxtd4Tx/t27zGm1b3++YnbMAsmirHMo7DWzOc3wL98vaQbY0vjW4vXQmYQUvlWMLh/DqTc2l/InwYfrTW1vjS5vbSlYAZtFSOVRwOh+GtGpoaX9rcXroSMIOWyvEBLY0vbW4vXQmYQUvl+IB2xvcnTe5Zuh4Q1045PqKN8Z1u1k2T20tXBWbQRjk+pv7xXV4j8celBlhN/eX4qMrHd41uJ21uL/kDwBwqL8enVTy+SXQ7aXN7oz8AzKXicpSg0vFdojt9PXAa3W/dhWVUWo5S1Di+bHQ7aXW/hReWUWM5ClLb+O5Gt5dmV3dhGbWVozBVje/36nZ0F1ZRVTnKU834Hke3p7uwhmrKUaYKxndJ7uPq9mQXlldBOUpW+PhejO6F6MKyCi9H6Qoe31vJBdZQcDlqUOj4RBeKVmg5alHg+EQXildgOWpS1viuzRVdKFpZ5ahOMeMbmiu6ULxiylGnAsY3Sq7mQh0KKEfNPjs+zYU6fbYc1fvU+CQXavapcjTiA+PTXKjeB8rRklXHJ7nQiFXL0Z5VxjcOruZCA1YpR7uWHV9SXMmFRixbjuYtNT7BhZYtVY6NmHt80+AqLrRp7nJszCzjS2rbS9cBGjJLObYrML40tSfpWkCDAuXg9fGlnT1J1wLa9mo5mLg/vjSut9I/AWzF/XLw2J9dWtNH0q8AbJDwRvwW3nRdgDPhDTE+4HXKEWJ8wOuUI8T4gNcpR4jxAa9TjhDjA16nHCHGB7xOOUKMD3idcoS0Ob42j6rRw3JU9RgdVZsHuJo2x9fmUTV6WI6qHsI7lzbH1+ZRNXpYjqoewjuXNsfX5lE1eliOqh7CO5c2x9fmUTV6WI6qHsI7lzbH1+ZRNXpYjqoewjuXNsfX5lE1eliOqh7CO5cVx9fmptbclk0F2VSQ8M5lxfG1uak1t2VTQTYVJLxzWXF8bW5qzW3ZVJBNBQnvXFYcX5ubWnNbNhVkU0HCO5fdWbp8AWts42zFTa25LZsKsqmgblMrJqNhK46vzU2tuS2bCrKpoNGmVtxqi1YcX5ubWnNbNhVkU0HCO5fLvxsAXpCmBAAAAAAAAAAAAHjJ/rjb7Q7p0krkdz5Zmj489I7765JSPXV0/ZLxo1Ld7nd+6fhoGnuuDj8Pjzcrleidw6nnySrBftd9nx93VU4rv/PJ0nSl7pvnR/mxSnf8ztLDroob2W/2O790ejRtPVfVHE264/ml6eGkj/nFaZZfX1WOK7/zydJ0pf3x9Bdz8X/TpDueX/pzOLsawpvud35pcjRtPVfH7lywO2mcrFSitw6nmierCMfzvxwOt/+uKF9+55Ol6Up1/FPv63bH80v3/UniaIVCpfudX5ocTVPP1eUSys/fLKW36a3DqebJKsL1u7yG//Om8jufLE0e7qs5zqeOrlNFeDP7nV86Opq2nqvLtc9D5uJpYd45nHqerBIM3+XlfzfcyO98sjRd6Xis5N9C6Y7fX1pDeHP7nV86Opq2nqvrP9uL/7/aW4dTzZNVhOHfEpd/R1Qkv/PJ0nSl/vp/Dd8i6Y7fX1pHeG/3O790dDQNPledCsKb2/H80uvhVPNkFSF7olGL/M4nS5OHh/77o/hv/a+bHf9laQ1PXW6/80snHzX3XJ2Wjx6UKL/j+aWXw6nnySrCvWFWIb/zv4f3x77/Hin+J8u3O35vaQ1PXW6/80snR9Pcc9U/KP148jueXzo6nEqerCLcGWYd8jv/MLxf3b+VhsWlyu54dmkNT11uv/NL06Np67n66g6o9H+P53c8v3R6ODU8WUW4M8w65Hf+mfB+7XfFH25+x3NLa3jqcvudX3pzNE09V93hFF+m/I7nlyaHU8GTVYT8MCuR3/mnwtv91Vz4939+x3NLa3jqcvudX3p7NC09Vz9HU/6/xfM7nl+aHk75T1YR8sOsRH7nnwvvvvjvj/yO55bW8NTl9ju/9PZoWnquvg7ldze/4/mlN4dT/pNVhPFdPNXNK7/zydL8SuNvo0Lldzy39DZV5cntd35p5mgaeq6qeD+j3I7nl2YOp/gnqwjDd3mF88rvfLI0v1L6qED5Hc8tzaSqOLn9zi/NHE07z1UmVAXK7Hh+ae5win+yirAfDXPyiRrkdz5Zml8pfVSg/I7nlmZSVZzcfueXZo7mdklhcseRW7ov/oaG3u2O55dmD6f4J6sMl3e4KP/mwoz8zidL8ytVcGUlv+OZpZlUlSez3/mlt0fTzHM1hKrwA0p3PL80ezgVPFlFuAyx/JsLM5KdP/8nWTp9eHlNY/rD2AI9dXSnRWmqCvTG0bT2XA2hKv2NvN44nIqerDKcpjh99XU1Jjt/uLxqJjmk8cNunf7fjM3KAAAA5UlEQVRRDd8eTx1d//GQrXK9fDStPVfdLa4X1z9XqJcPp6onqwzdjXfVXhAf7/xwI0tySKOH/XdLNb+d5JmjO5y/98t/Al89mtaeq+5VXQ09WdPDqevJKkTVbxo/3vnho+SQRg+7Xww1/lTZnjq6arx6NJ6rD3r1cOp6sgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADK9H+AP9TQQGNapQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "!gracey.py -com strain_A1g/Y_stress_A1g-2.txt PNG strlt=view strc=br str=A1g,0.95,0.7\n",
    "Image(filename='Untitled.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ThO2\n",
    "### Linear Elastic (GGA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "%cd ../.."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operationally, doing ThO2 is no different, so I will just provide the results (see my dft_runs directory for all the data). Units of GPa.\n",
    "\n",
    "```bash\n",
    "# /home/cam1/dft_runs/tho2/gga/c500_k20/irreducible_elastic\n",
    "ls -d strain_{A1g,Eg.0,T2g.0}| robo_run.py 'pwd;head -2 Y_ener*2.txt|tail -1|pcol.py 3' fold=t | pcol.py 2|transpose.sh |format.sh 2\n",
    "ls -d strain_{A1g,Eg.0,T2g.0}| robo_run.py 'pwd;head -2 Y_stre*2.txt|tail -1|pcol.py 3' fold=t | pcol.py 2|transpose.sh |format.sh 2\n",
    "```\n",
    "\n",
    "|Target|c$_{A_{1g}}$|c$_{E_{g}}$|c$_{T_{2g}}$|\n",
    "|---|---|---|---|\n",
    "|Energy|572.51|246.68|79.43|\n",
    "|Stress|570.25| 244.79| 78.29|\n",
    "|Spunar|562.7|246.5|70.9|\n",
    "\n",
    "Here we see excellent agreement between energy and stress; and the published data of Spunar:\n",
    "\n",
    "```\n",
    "Theoretical investigation of structural and thermo-mechanical properties of thoria up to 3300 k temperature\n",
    "B. Szpunar J. A. Szpunar\n",
    "Solid State Sciences 36 35 2014\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nonlinear Elastic (GGA)\n",
    "That is great, but we will also need the third and four derivatives of the energy with respect to the identity strain. Let us check out the results of the third derivative, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xmgrace -nosafe -noask -hardcopy -hdevice PNG -batch /home/cam1/.ppbatchfile\r\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABXgAAAK8CAMAAABiCycpAAAAQlBMVEX///8AAAD/AAAA/wAAAP///wC8j4/c3NyUANMA////AP//pQByIbxnB0hA4NAAiwDAwMCBgYFCQkLv4+Pfx8fPq6vxjVQBAAAACXBIWXMAAAsTAAALEwEAmpwYAAAADnRFWHRUaXRsZQBVbnRpdGxlZPPPpzIAAAALdEVYdEF1dGhvcgBjYW0xmdSm6wAAABV0RVh0U29mdHdhcmUAR3JhY2UtNS4xLjI1jnpa/wAAIABJREFUeF7t3dt2o0YagFGRNSvXQKff/1VHHIpDUbKFS8iU2PtibCGFrvzjfCEI5NsNAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA3q5u7ur1tsSmflsTbZukdtJtjjZ1L9u+CuBSmqq6x7Ct2kUPm6rtNlWLTXXbbetePG+bpXZy6/+i1cvvL+tJL3Bl91j2X+s5h11jx03Dc8unw+tXEjvpHt3TvQpv91h5gau7H4OO380HuFNbF8e303fLGgepnQzZXYX3fvjbPdkf9s5bAS5mLuNU1EU8p2PTOa3Lp4PETu7f9o8W4a3DX7feDHAty3O243Fol8XF00NH19uiaiZ2Mj+aXzuf/203uwC4jEQzl5umCL8kvNN3210AXEZ3vnX8NpwlWB2OhnIuCro91ZDYyejBOYVmswuAy1icbg2HoatYVnONQ1vXx7SdxE5GD8KbvDICoBjd9bZVm+rbM7q/uD/6nI5V4/BOJ343b71NtjsZPQhvYhcAxajbrrldF3/Wsu6v7OI4JzMZ3r6t3btjqe4mdjJKh7eNbrMAKEk9lq47zfqzmvXRbJv5nrOptesH/VW5TXQnWrDZySgZ3gf7AChCHbK4OM066WIYi17S6aO5eGbV8MUzm/shluKdjBJ/Qbf7n/5LAuD3zQePiUtjhxauRS/pDdGczhHMp3Nv68PfobwPkhntZFTFqxqy+3g3ACdXzwGrt58b1n+YWCR+TWc8Wp2iOb1VFkV4CO+DZMY7GVRxeO/GW4njzQBFaF7Rr+7arr6Fq8J2vazb8F334J7coa6p8m530kuFd4x0YjvA+b3iwHG4pnZ4dyxsCycEmm7z0Nnx4LqP67a8qZ10HgR2dSANUJKfhHc6yzpdJtZ/jT7CsTstUS/uSavDk/1fPr1s9GAnD8O7vNcNoCg/ydc6vIs7zbpoLl7Xm8I5P9cdrUYxfbyT7WsH3U7ibQAleNS1r/S/omd6o22xh3RRhz4uj1C3gX68k+0uR5t9AJSh+vJUaZWQeM1UxnZTyenZZXi3R6uPd7J8ZuXLhQOc1/p+tThli95Oopesyrg4YzBopzyuzsludvN4J8tnVh5tBzi59WW2ccqmcwoL0UviZq6uV2jmnUfh/eJq3fVOHoW3idsNUIrufOvw+Qh1O9Wwjg99v7CIa1zD5VNdQRe3EQ8xrcNHMzzeyaPwtolL0gCKMN4vtvxV7N236xd9aVHUdSTn8wyd9dHv8E33J0/fpXeyCm/4d0S47BegSGN55/LV3e1l69d8ab5FYlXD9dnj26LD450U4x89PPlgJ8MBeXgw3F5Rd0fkjneBko2ffLDM3a7wDgezTbP8SMe637R8Tbexu16hO7SenlgeE293cmu6eHdLC2eWh5VWD84+ABQkfs9sZ3j7HfT3qY26vi6enQxXAC82LP6azU5SNjsA+BC7wwtAHuEFeDPhBXgz4QV4M+EFeDPhBXgz4QV4M+EFeDPhBXgz4QV4q+7Tydpv7t0FAAAAAAAAAAAAAAAAAAAAAACArQpgvzgl7HHc+A7bsx3PDtuzHc8O23PROz7sz7iG48Z32J7teHbYnu14dtiei97xYX/GNRw3vsP2bMezw/Zsx7PD9lz0jg/7M67huPEdtmc7nh22ZzueHbbnond82J9xDceN77A92/HssD3b8eywPRe948P+jGs4bnyH7dmOZ4ft2Y5nh+256B0f9mdcw3HjO2zPdjw7bM92PDtsz0Xv+LA/4xoKHF95Sy5vxQUuubwVF7hk4f1C99u9q/G3e9fNZNzQPzm/uMDxlbfk8lZc4JLLW3GBSxbeh+5lvXe1rqo+tM18o0lf27pq7//bDk92ChxfeUsub8UFLrm8FRe4ZOF9aIxqM9zSV7Xj8e6weehutzm8vMDxlbfk8lZc4JLLW3GBSxbeR+pwD3V/iNvMge0390fDty7L4WxDgeMrb8nlrbjAJZe34gKXLLyPjEe6Y3in7tZDaadPtpi/Gb8WpLwll7fiApdc3ooLXLLwPtKEt85Wb6HdN3dnGqYqz88WOL7yllzeigtccnkrLnDJwvtIXQ1vqzXzadzOUNz5DEM452B8wA8ox1p3HUM9vYs2ms40hPDOZyTGDQDPU45If6Hu6jxDONMgvMBrKEesO9vQThfq9sbMpsK7Nj4NsKYV3+iOeedbJG7TmYZkeMcNAM9Tjkh9b26X3kV5xzMNwgu8hnLcpo9kmG9OW5d3dW1vT3iBDMpxGy5lqIasjsG9lze+e2JxEdl8YZnxAfspx215xLt81yw8G840fMgNFMCvU46Vqa1zZKczDfMHOcybjA/YTzlWEuFd3EwRTvzO97UZH7CfcqyFkwjt1NbF3RQhuPMH8hofsJ9yrI0fgT6fVZhPK9xCcec32YwP+AHliNyT2zbtfO/a+mMb2u5+4sUNxcYH7KccG901DvG2WfjtawPjA/ZTjizGB+ynHFmMD9hPObIYH7CfcmQxPmA/5chifMB+ypHF+ID9lCOL8QH7KUcW4wP2U44sxgfspxxZjA/YTzmyGB+wn3JkMT5gP+XIYnzAfsqRxfiA/ZQji/EB+ylHFuMD9lOOLMYH7KccWYwP2E85shgfsJ9yZDlsfP/EG4DPcVg5ruGY8f0ziDcDH+KYclxGNYq35xizK73weY5IxvUcML5Fd5UXPtMB5biSA8YnvPDxDijHlbx+fKvuKi98pNeX41JePz7hhc/3+nJcyuvHJ7zw+V5fjkt5+fii7iovfKKXl+NaXj8+3YXP9/pyXMrrxye88PleX45Lef34hBc+3+vLcSkHjE934eMdUI4rOWB8wgsf74ByXMkR49Nd+HRHlONCjhmf7MJnO6Ycl3HY+FQXPthh5bgG4wP2U44sxgfspxxZjA/YTzmyGB+wn3JkMT5gP+XIYnzAfsqRxfiA/ZQji/EB+ylHFuMD9lOOWNNWVdXU46O6uT9qm/Bk3T8ZHhkf8BPKsXYv672rdVUN5W2q9v5NW7Xjs/037fjkzfiAn1COtTGq9+Pc7ssQ2vuUhqPc6eHYYeMDfkI5VuohuCG1/eHvberw/DCcbTA+YD/lWBkLG8IbzueOm8OTi2/Gr1zYnz/xluDxM1yccqw0IbUhvMN8hvBOVZ6CbHzcbn//xlsGf/7++1+8DXrKsVJXw9tqzXAat62Gr2FjOMMQzjkYH7c///6bOrC9Z/df4eUB5Vjrrh6rp3fRukPednrHbb6QbD4jMW7guu6B3R7y9tkVXh5Rjkh/oe5U2O4IuG27S8puwktSF9h42+1PfyQsvDygHLG+tdOFukOIh29T4V0bn+ZC/nt4ZPtoO5ekFd/oUjvfItF1eH01780RL7N/+5MK8daO8PKQckTqe3O79I7lbdvhELj7XnjZ+PPvf90xb6qwDzaDcvTqpte1NtwUPJa3f1stlFd42fj774OzvMLLF5TjNlzKUIUrd8fg9qkdryobLzKbLiKbLywzvov707W1O+RNXFG2Ce9/3VmJv2KMcvTmI97lu2b9/w6BHe6rcAMFsf+64nbXL2yvKNuE9++/f/8MV/feJULNhSjHytTW4Zv1BQ3TBznc5m/Gr1zUENyuptuSRuH9O782TjKXoxwrX4Y3nIcIZyC6LeNXrqk/4B3ONWxLut74X2hz8rVcjHKshdQOH8EbOjx+Pm8I7vyBvMZ3beFNte4Ydv3MsHFR2OkV/yXPS3AtyrE2JjacVRg/Aj28rTYUd36TzfiurX9rrZM85H0Q3u6M8LyZa1KOSN39pp92unetvR8CD2+tTQ/nR8Z3cd21ZIPUIe+D8C6+47KUY6O7xmHxcLzAd364fGR8V/ZnPmeQ+iiyTXjHR8KLcuQxvivrP6ZhFp+5XYf37yK8caK5HOXIYnxX9u/f7orcQVfe6IqydXi7U7v9C/444EU58hjfhf0XXS0WH/KuwxvK++fv3+0lv1yNcmQxvgub31rrdIe8y8eb8PbldcswPeXIYnzXNV1LNti+vRZvcKTLRDmyGN91/V2f0+2PZ5cb4vA61GWmHFmM77I275F14V21df14vLsYOsqRxfgua/O72/sLG5Yb1uHtT+/2BBjlyGN8V7X9QLL+Y8cW53H7D9+Nnh4523t5ypHF+K4p3DoxHdHO91IMm4YLe/uD3PCSZXmddrg65chifDzrb3+mYcxvfJcbF6McWYyP5/xZXfW7vgSY61GOLMbHc+IrHpYPuB7lyGJ8POVvdPGZUw0XpxxZjI+nRO+nuazh6pQji/HxlOUNxX/cS4FyZDE+njJ8QM7I8S7KkcX4eNJ0pa/sohyZjI+n/XG/MIFyZDE+YD/lyGJ8wH7KkaUaxdsBEiTjFYwP2E85shgfsJ9yZDE+YD/lyGJ8wH7KkcX4gP2UI4vxAfspRxbjA/ZTjizGB+ynHFmMD9hPObIYH7CfcmQxPmA/5chifMB+ypHF+ID9lCOL8QH7KUcW4wP2U44sxgfspxxZjA/YTzmyGB+wn3JkMT5gP+XIYnzAfsqRxfiA/ZQji/EB+ylHFuMD9lOOWNNWVdXU46O6WT661f2T4ZHxAT+hHGv3st67WlfV0Nq6avttY3n7h7fpofEBP6Eca2NU78e53ZchtPcpDVunh8OX7rvwDcDTlGOlHoLbpbY7odAf/t66DvepnR+Gsw3GB+ynHCvjkW4Ib3g0nnoIDxffjF8BnqccK01466z/OmU4+bD/ZtwA8DzlWLkf2vbHtsO5hWVpq+UZhnDOwfiAH1COte7qsTq8ixaFd76QbD4jMW4AeJ5yRPoLdac30cKFY8ILvI5yxLqzDe3U2/HCsYfhXRufBljTim90x7zjkW47XcD7KLzjBoDnKUekvqe2S+9U3rZp2rZvrvACL6Ect+4DGXpda8NNwaG83ZttbT1cxyu8wEsox22oaxVumeiD2x3prl/Qb5zD63Iy4MeU47Y84l2+a7Z4wRBlN1AAL6EcK1Nb58jeukPd/vB3+iCH2/zN+BXgecqxkgzvdDnvfOI3nIcwPmA/5VgLJxHGg9zO9OG8U3DnD+Q1PmA/5VgbKzufVVgmeCzu/Cab8QE/oByRurtytw33rtXNdBdbr7ugd76qzPiAn1COje4ah/nb1VN3w+W+gfEB+ylHlvLG90+8AXi78spxKoWN759BvBl4r8LKcTZFjW/MrvTCbyuqHOdT0vgW3VVe+FUlleOEShqf8MJZlFSOEypofKvuKi/8poLKcUYFjU944TQKKscZFTQ+4YXTKKgcZ1TO+KLuKi/8onLKcUoFjU934TQKKscZFTQ+4YXTKKgcZ1TQ+IQXTqOgcpxRSePTXTiLkspxQiWNT3jhLEoqxwkVNT7dhZMoqhznU9j4ZBdOobBynE1541Nd+H3lleNUjA/YTzmyGB+wn3JkMT5gP+XIYnzAfsqRpRrF2wESJOMVjA/YTzmyGB+wn3JkMT5gP+XIYnzAfsqRxfiA/ZQji/EB+ylHFuMD9lOOLMYH7KccWYwP2E85shgfsJ9yZDE+YD/lyGJ8wH7KkcX4gP2UI4vxAfspRxbjA/ZTjizGB+ynHFmMD9hPObIY38wvjodnKUcW4xv9M4g3AynKkcX4emN2pReeoxxZjK+z6K7ywhOUI4vxdYQX9lGOWNNWVdvMD6v7wzo8qu9PVvOTxtdZdVd54XvKEWmruqttu3xYh4fDN/223v+q/z1jfPWnEl7YSXjXxqiG8k4P+y8hwFOWnwzvDuOOiyK8sJPwrtTVOJCqT+3iYfe/7XiWoZnONjw1vjiuOeJ9n0DUXeWFbz1VjuuYktr2R7XTw+GbkOHFN+PXV4k7+6R4N++lu7DTq8tRuOmds6Zv6/rhsG31snePLy7uV+K/9jDCCzu9uxwntyxtsziy7c85zGcYwjmHk4wvTu5G/Be8lvDCTucox2mswzsf4vZPzBeSTU+ceHxxfAfxq15Cd2GfE5fjN1ThgoXyw7sQ13cj/gt2El7Yp4xyvE0zXM3wYeGdxcXdiP+C5+gu7FJaOQ5Wh0Pe9snwroVXFyau7yB+1ZdkF770Ga04yr28VdM0Vbu+jOFheMPzHyKO7yh+WZLqwrM+rRw/Uje94e604dMY+ut4rxfehTi+o/hlb6DofJwPLsfzmuHgP2S10z+ablzrHi8uIrtNF5ZdY3xxfAfxq47hHAaf6Brl+MbyiHcwHtJOMe4fzwfA0/bLjS+u7yB+1cuEt+ykl89yuXI8pR4vbpiOcftD3PkAeP5m/HpBcXxH8cuyLLqrvHySC5fjC+HahunjIYcQz9eahY8nM75OHN9B/KofEF4+lHIkDJ+Q0xlTO34eZPRhkd3z41dGcX0H8auetOqu8vJBlGMjnGcYvu9SO51jGIo7v8lmfI/F8e3Er/mG8PKplCPSLLJ6G+6oaOff/dPdV7G8/MH4nvHT/Aovn0o51lbXNvTqZhni6AXGt8fOAEfdVV4+h3JkMb79ovx24peMdJdPpRxZjO/n4vr+bxNg4eVTKUcW43uFR/kVXj6VcmQxvlfaBFh3+VDKkcX4Xm+RXuHlQylHFuM7TlTe+GkomHJkMb6jjdmdDoLj56FEypHF+N6gO9idwivAfALlyGJ8bya/fATlyGJ8v0SAKZpyZDG+XxXlV4AphXJkMb5T+GF9XSnBb1GOLMZ3Knvq6yo1fpFyZDG+81nU93F+w9XB0suvUI4sxndSy/omArzorvLyC5Qji/GdWpTfOcDCy+9SjizGV4S4vqvuKi/vpxxZjK8oIb3Cyy9TjizGVx7h5fcpRxbjK1HUXeXl7ZQji/EVKerueOYX3kY5shhfkVLhVV/eSDmyVKN4O6cWhbfbNMdXfjmQZLyC8ZVp092B+vImypHF+Mr0ILwD9eVwypHF+Ar1RXcH6suRlCOL8RXry+yOFvWVX15JObIYX8m+rm6gvryecmQxvqvYld/nis6FKUcW47uUp+r7zDkMrk45shjfBX1Z3/CWnfTyFeXIYnxX9aC+i+4qL48pRxbju7RFfcf8Ci9PUY4sxseyvqvuKi8PKUcW42MgvOyhHFmMjwXh5UnKkcX4WIi6+88/icse4KYcmYyPpbi7vfhFoBx5jI+lKLzTG2/x67g65chifCzF4e2IL1vKkcX4WNl2txPaK76MlCOL8bGSDm9niq/8ohyZjI+1R93tLeKrvtemHFmMj9jD7Abyi3LkMT4SvqhuIL7Xphyxpq2qtlltmR/V9yerxZPGRw7tvSrliLRVfU9t1YbHXWqn0tb99v4lA+Mjl/hekXKsjVGdytt2B8AhvEN370Obsmx8vEBor/hehnKs1NU4kGo8qq27b0N42/G7ZtpifLzIFF/5vQLlWJmS2i6OaufMhmnN34xf4QW+j+8Tb9tRAuVYmSLbTG1NbZs2GR8v9zC+316oRjGUY2UZ2fk4d9oWn3MwPo6RiG+4L0N6P4FyrHwZ3nnTdOxrfBxmHd9Fd5W3fMqxUoVTu8LLGUzxXf9yi/hllEY5VppwNYPwchbb7ipv8ZRjpQ6HvIurd78M79r4NLya8JZOK75yL2/VNE3VJq5qSIV33ADHEt7Pohx3ddPrTzKMn8aQuo5XePktUXf9Gs3SKcet72gnZLWzeDR9O11ENl9YZny8R9zdUfwyCqEct/UR72Bx/8QcXjdQ8Gui8D5xkxtnphwpdbi4oTNVdvogh5tbhnmzOLwD8S2VcqRMl/MOD+brG8K1ZuFp4+NNUt3tiW+JlCNh8c7abRneENz5A3mNjzd5GN6O9pZGOTZW5xkWl/beQnHnN9mMj7f5ors97S2JckSaRVZv3eHu6nLn7r6K5eUPxsf7fJXdngPfYijH2urahpT1C4yPt/qiugPtLYNyZDE+zkd8z085shgfpxTaK74npRxZjI/TEt8TU44sxsepTfGV33NRjizGx+mJ7wkpRxbjoxDftPfbyyV4KeXIYnwU5FF8v71AmFdTjizGR1kS7Q33w0nvGylHFuOjPOv4LrqrvG+jHFmMjyKF9v5PeH+HcmQxPoo1lHfVXeV9F+XIYnwUTXh/iXJkMT4KJ7y/QjmyGB9li7qrvG+iHFmMj8JF3d1e5MsRlCOL8VG4TXjnC804jnJkMT4KF4U3dYcFr6ccWYyP0q2729PewylHFuOjdInw3uK723g15chifBQv1d2O9h5IObIYHx8gld2B+B5EObIYH58hVd1BaK/4vpJyZKlG8Xb4JNr7MpLxCsbHVYjvKylHFuPjQkJ7xTebcmQxPi5Ge19CObIYHxekvdmUI4vxcU0OfPMoRxbj47K0N4NyZDE+Lk17f0g5shgfVye9P6EcWYwPvmzv41viLk05shgf3B6ecnj8IRBXpxxZjA8GY3sX8Q2feSa9W8qRxfhgEto7xHfRXeWNKUcW44OVub3C+wXlyGJ8sNGXd9Vd5Y0oRxbjgxTh/ZpyZDE+SBPeryhHFuODpKi7yrumHFmMD9Ki7m4v8r005chifJC2CW+4zIybcmQyPkiLwvvo5raLUo4sxgcPrLvbkd6JcmQxPnhgG96b9gbKEWvaqmqb8Khuqqpq6ulh2z0Mj4wPHkt09/bw83QuRjkibXWvbFO1w6O6Gozlrfvt/UsGxgePbbPbG9t75fgqx9oY1bG8dX902x309luH7t6HNmbZ+OAbm+oOQnuvGl/lWLkf4Q7fDKlth0Pbe3n71LbjWYZmOttgfPBTV26vcqxMSW271NbhyLYaehyqvPhm/Ar8xFXTqxwr0ztnTdfW6U21/tH4v6uXGR9kumR7lWNlGd7FxQshvGFTOOdgfJDvgulVjpWH4e0ezJumY1/jg1e4WnuVY2V8Fy0O7/BWm/DCYYb0XiW+yrHShEt2V+EdLyMTXjhSaO8F4qscK3U45G2X4W2HKaXCuzY+DfzQx7ZXK77S3arWNE3VLiZTj0fBqfCOG4CX+dT2LijHrftAht5wd9rwaQz9dbyjcIew8MKbfHp6lePWd7SzeDdtEdlw+9riIrL5wjLjg4N8dHuV47Y+4h3M90rM3V1snLJsfHCYDz7loBwp9fR5ZLd2PuUwfZDDbf4mPAccYGxvIr4PPn2nEMqRMl3OO3e3rrvN4Vqz6UMcxq/AQUJ7l/F98HmT5VCOhPmdtbm73aBCcOcP5DU+eIN1e8PHqxecXuXYWJ5nmC+760/qDsWd32QzPniXKb2L7hZbXuWINIusht8/0Rm2dPdVLC9/MD54n6G9wvt5Vtc2pKxfYHzwTnF3Sy2vcmQxPngz4cX44M2EF+OD94q6W2h5lSOL8cGbfUJ3lSOP8cGbReFN3dV2fsqRxfjgzTbhne+sKIdyZDE+eLdVd9O3FJ+ecmQxPni3KLyd4tqrHFmMD95u091eUelVjizGB79gm91eOe1VjizGB79jU91eKelVjizGBydTRHuVI4vxwekM6T11e5Uji/HBGZ09vcqRxfjgpE7dXuXIYnxwWic+5aAcWYwPzuys7VWOLMYHJ3fK9CpHFuOD8ztfe5Uji/FBCc52ykE5shgfFOJU7VWOLMYH5ThPepUjSzWKtwOn9NvtlYxXMD4ozClOOShHFuOD8vx+epUji/FBkX65vcqRxfigUA/Tm/6k39dSjizGB+VKtPfB77Z4NeXIYnxQsiG9U3vDr3I7PL3KkcX4oHCL9C66e3B5lSOL8UH5QnuFtxDGB5+gL++qu8eWVzmyGB98COEth/HBxxDeUhgffIqou4eWVzmyGB98jPd1VznyGB98DOEthfHBx4jCu7ix4uWUI4vxweeIu3tcepUji/HB51iF95b8LIcXUY4sxgcfZN3d2+azHF5GObIYH3yUVXZ7h6RXObIYH3ya7eUMr2+vcmQxPriCV59yUI4sxgcX8dL2KkesaauqbeaH1f1hHR7V9yer+Unjgwt5XXqVI9JWdVfbdnhUV4OxvHW/vX/JwPjgUl7UXuVYG6Maytsf3XYHvf3Wobv3rWOWjQ8u5yWnHJRj5X6EO3wzpLYZDm2b8fRCO55laKazDcYH15PfXuVYmZLa9ke1U1+HzaHKi2/Gr8ClZKZXOVamd86aqa3D5v74dxHe8WXGB1eV017lWFmGd754YTy3O59hCOccjA8u7OenHJRjJR3e8eqyedN07Gt8cGk/TK9yrFThgoVFeOt2vI5XeIGNn7RXOVaacMnuFN7+lokhx8ILJOxPr3Ks1OGQtw2VbZruVrb+QSq8a+PTwMV8216t+Ep3q9q9tVW7nEy38ZYOb3gJcG3fpndJOe7qpjfcnTYc3w7X8QbD8a/wAl94Pr3Kces72pkvY1hGtlP3D6eLyOYLy4wPmD2bXuW4rY94B+v7J8YOu4EC+MZz6VWOlHr6PLJR/3j6IIebW4aBB5452ascKdPlvKPx1rX5WrPwtPEBsSG9m/YufqeQciRM76zV43FvO4wpBHf+QF7jAxI26V3/Fk3l2FicZxiPfJvwKyiG4s5vshkf8MCyveGXxof0KkekWWS1u4ysaufPxhmuK1te8GB8wCNTehfdHcqrHGuraxuG6x2Wj+MXGB/whSG9wvtaxgd86R7eVXf78ipHFuMDviG8r2Z8wHeE98WMD/hG1N2uvMqRxfiA72y6qxx5jA/4jvC+mPEB3xHeFzM+4Ftxd5Ujj/EB3xLe1zI+4HtRd5Ujj/EBz1hmVzkyGR/wJJ/H+yrGB+ynHFmMD9hPObIYH7CfcmQxPmA/5chifMB+ypHF+ID9lCOL8QH7KUcW4wP2U44sxgfspxxZjA/YTzmyVKN4O0CCZLyC8QH7KUcW4wP2U44sxgfspxxZjA/YTzmyGB+wn3JkMT5gP+XIYnzAfsqRxfiA/ZQji/EB+ylHFuMD9lOOLMYH7KccWYwP2E85shgfsJ9yZDE+YD/lyGJ8wH7KkcX4gP2UI4vxAfspRxZtRHJvAAAH0UlEQVTjA/ZTjizGB+ynHFmMD9hPObIYH7CfcmQxPmA/5Yg1bVW1zWpTXYXH9f3J6cHN+ICfUI5IW9X3+FbtYlM9tbbut/cvGRgfsJ9yrI1RXZW3DeEdunsf2vSk8QH7KcfK/eB2+Kaaj2q7cw9DeNvxazOdbTA+YD/lWJmS2k5HtXV1C+ENVV58M34FeJ5yrExnc5uprW0dts7bppcZH7CfcqwswzueXmimrfMZhnDOocTxlbfk8lZc4JLLW3GBS16suLzFH2oT3ro74zB+P19INh37Fji+8pZc3ooLXHJ5Ky5wycL7SBVO7Ybw9u+xCe+vKm/FBS65vBUXuGThfaQJVzOM4R3upBDeX1XeigtccnkrLnDJwvtIHQ55h0vI+hMNX4Z3bXz61IpY5Ep5Ky5wyeWtuLQlF9iKN7qXt2qapmq7ydRTXh+Gd9xQkPKWXN6KC1xyeSsucMmLFZe3+APUTa8/yTB+GkN/HW873kQhvL+qvBUXuOTyVlzgkoV3rRkO/kNWO/2j5X8Z9DGewzt+d9z4DtuzHc8O27Mdzw7bc9E7PuzPKMnyiHcwHNJG4X3rDRSH7dmOZ4ft2Y5nh+256B0f9mcUrV58VMNt6uz0QQ63N9wyfNie7Xh22J7teHbYnove8WF/RtGmy3nDw+EAd77WLDx93PgO27Mdzw7bsx3PDttz0Ts+7M8o2fwJOYMQ3hDc+QN5jxvfYXu249lhe7bj2WF7LnrHh/0Z5YrOM9zm8I7Fnd9kO3B8h+3ZjmeH7dmOZ4ftuegdH/ZnlKpZZDWYL3jo7qtYXv5w3PgO27Mdzw7bsx3PDttz0Ts+7M8o1OrahpT1C1bXPQA8Z1ERAAAAAAAAAAAAYL/xs3vPK73A7dbxN22cwHZtia3N/eHwS5nO4Kklp1/0S9KLSW1tz3K5aWpx263Dpww2p/hpjtf2cGvdtOf5YS5C3X90w/zBDaeTXuBma/dxxIvnf9Nmbamt3c9uFX2O0e8pb8lPrbjXnuUHI7W47dbu18d0zvAPZLy2R1vrxJ2yfGkY4u12kn+cttILjLfW7XnCG68tubXtfjFId+Rwirk/teTmTEt+asXjtpP8YKQWl9jajke880t+zWZtD7ZOv12Xp4V/VU2/juJs0guMt9bhY99PIF5bams4LXKSI5tnlnyblzy94vc8teLe8JsHTyC1uO3W8CsSzyBe24OtZxlwUbYfi34y6QUmtp4mvIm1bbaGX4S3+rCi3/PMksM6zzHnZ1bcux9AnmLBycUltp7pTGm8tvRWx7s/kPhFQOeSXmBq61n++UqtbbN1+q+0U4z9qSXPm09wqiG9tsTW+39anOQHI7G4xNa6OsVJhl68tvTW+hQ/wqWZ/yPipKfH0wtMbT3JP1/JtT3aeprwphaX3nqL32r5Fem1bbd2/+F+kh+M7eJSW5v+7csT/EzctmtLb3Wi4SfOd8AYSS8wtfUsfwOptT3a2m1fPPgt6cWlt9bhLMmvSq9tu7Vb7El+MLaLS20dL2koZ8h192Zr94t2pxfwhPRsTyS9wNTWs/wNpNb2aOvtFP/d/mBxya3r91l+TXJt2639xbAn+cHYLO7R1u4C71OcNk2sbbv1/j/tcJ3hCVZcjvRsTyS9wNTWs/wNpNb2aOs5/rv9weISW5uz/BOWWFti63CFwEl+MOLFfbG1v5Q3PPg9ybXFW9vhcfeTcYKfi2KkZ3si6QWmtp7lbyC1tkdb62nrr0ovbru1+y/KgpowXjtykh+MeHFfbQ05+13ptUVbw8PTXC5dhvRsTyS9wNTWs/wNpNb2aGt7hhMNjxaX3lpQE8b3qE7yg/HUkoNTZCy9tgfh7X4uHPI+LT3bE0kvMLX1LH8DqbU92HqO+/EfLO7B1nKaEG5SOckPxjNLnp3hWoH02jbhnR+FJ/jW8mqsU44tvcDU1pP885VcW3rraT7VJ7W4R1vP0YT02tZbm/6EdDC96Lc8s+TZGX6a02uLtk65Fd495v9/Tzq29AJTW8/wo9pJrS259TTdTS3u4dZzDDq9tvXWk4X3mSXPShny3OHmpAU5p3oxxdUTZ5FeYGrrGX5UO6m1pbbW5zkltl3c463doH//3xjptaW3nuQHI7249NZu0b9fsfTaoq3TQhvnePcI0zrDP01J6QUmtp7kn6/k2hJb5+7+/j9hm8V9tfV2ikvg0mtLbj3LD0ZycQ+2Rhn+Jem1rbdOHT7LmAsRZnqKf5pSogWOXxLLPs3/8U+teO7uGe4PfWrJ4QzJCY7FNmtLr3hwlh+MZ5ZcT0Ne/A38mmdWfAvnGuZzvzxjGN+Jp7ZaYBM+DXa77PNczvLEirtLA05z9vH21JLH/5Ssz/Avitt2bYkVj84S3meWXA0nSs/R3adW3P0o9z8X45M8q7su89TnxZcLrKeVrpfdjBE7x9/G9ysebrIsasn9vyua/gPcT+H7FY9OE94nltz9XLTNSf7ddntmxf3D+5J1d7+mOcs/TA8sFzh/d+Zll7fiJ5bcnGz936/4dL5fcpFD9hk5AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAK/xf/OoOsoOvxTXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "!gracey.py -com tho2/Y_energy_A1g-3.txt PNG ss=1.2 strlt=view strc=br str=A1g,0.95,0.7\n",
    "Image(filename='Untitled.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us tabulate some of these results.\n",
    "\n",
    "$$ \\Large\n",
    "\\begin{align}\n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{A_{1g}}^2}= 159 \\,\\,\\textrm{eV/uc}\n",
    "&& \\hspace{4mm}&&\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}^3}= -900 \\,\\,\\textrm{eV/uc}\n",
    "&& \\hspace{4mm}&&\n",
    "\\frac{\\partial^4 E}{\\partial \\epsilon_{A_{1g}}^4}= 4258 \\,\\,\\textrm{eV/uc}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "If we want to relate these to volumetric strain, we need to perform chain rule (see course notes on strain):\n",
    "\n",
    "$$\\Large\n",
    "\\begin{align}\n",
    "   \\frac{d^3 E}{d \\epsilon_{v}^3}&=\n",
    "   \\frac{d^3 E}{d \\epsilon_{A_{1g}}^3} \\left(\\frac{d \\epsilon_{A_{1g}}}{d \\epsilon_{v}}\\right)^3\n",
    "    + 3 \\, \\frac{d^2 E}{d \\epsilon_{A_{1g}}^2} \\frac{d \\epsilon_{A_{1g}}}{d \\epsilon_{v}} \\frac{d^2 \\epsilon_{A_{1g}}}{d \\epsilon_{v}^2}\n",
    "\\end{align} \n",
    "$$\n",
    "\n",
    "where we have:\n",
    "\n",
    "$$\\Large\n",
    "\\begin{align}\n",
    "\\left.\\frac{d \\epsilon_{A_{1g}}}{ d\\epsilon_{v} }\\right|_{\\epsilon_{v}=0} = \\frac{1}{\\sqrt3}\n",
    "&&\n",
    "\\left.\\frac{d^2 \\epsilon_{A_{1g}}}{ d\\epsilon_{v}^2 }\\right|_{\\epsilon_{v}=0} =\\frac{-2}{3\\sqrt3}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Plugging in, we see that:\n",
    "\n",
    "$$\\Large\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{v}^3}= -279 \\,\\,\\textrm{eV/uc}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BID at second order (GGA)\n",
    "\n",
    "It is useful to show how one can manually perform BID, measuring all the derivatives simultaneously. The runs can be found here.\n",
    "\n",
    "```\n",
    "/home/cam1/dft_runs/tho2/gga/c500_k20/bundled\n",
    "```\n",
    "There is a yaml input file called `inputs_mkruns.txt` with the following:\n",
    "\n",
    "```yaml\n",
    "fpinp: ../static/\n",
    "delta: 0.02,0.08,0.01\n",
    "order: 1\n",
    "pgname: Oh\n",
    "schema: mkruns\n",
    "unstrained_path: ../static/\n",
    "vname: v0\n",
    "xtal: ../static/POSCAR\n",
    "_sim_trans:\n",
    "  v0:\n",
    "    A1g: 1.\n",
    "    Eg.0: 1.\n",
    "    T2g.0: 1.\n",
    "```\n",
    "\n",
    "You can see that we are taking the first derivative with respect to the variable v0, which is defined in the supplied similarity transformation *sim_trans* (the leading underscore means that it will bypass the parser). You can see that this new variable is equal parts A1g, Eg (zeroth row), and T2g (zeroth row). This will clearly probe all irreps, and because there are no repeating irreps, there is no interference. Therefore, we can measure all elastic constants via one stress derivative. Going into the following directory:\n",
    "\n",
    "```\n",
    "/home/cam1/dft_runs/tho2/gga/c500_k20/bundled/strain_v0\n",
    "```\n",
    "We can execute the following to extract all elements.\n",
    "\n",
    "```yaml\n",
    "data_path: ./\n",
    "pgname: Oh\n",
    "target: stress\n",
    "unstrained_path: ../../static/\n",
    "xtal: ../../static/POSCAR\n",
    "order: 1\n",
    "units: GPa\n",
    "sigi: A1g,Eg.0,T2g.0\n",
    "```\n",
    "\n",
    "We can add these results to our table.\n",
    "\n",
    "\n",
    "|Target|c$_{A_{1g}}$|c$_{E_{g}}$|c$_{T_{2g}}$|\n",
    "|---|---|---|---|\n",
    "|Energy|572.51|246.68|79.43|\n",
    "|Stress (LID)|570.25| 244.79| 78.29|\n",
    "|Stress (BID)|570.71|244.73|76.93|\n",
    "|Spunar|562.7|246.5|70.9|\n",
    "\n",
    "Here you can see the advantage of BID: you get everything at sufficient precision from one measurement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear and Nonlinear Elastic (LDA)\n",
    "\n",
    "It is also interesting to explore more complicated derivatives. Let us consider the following two new ones, in addition to one previous one (but now with LDA, just because):\n",
    "\n",
    "```\n",
    "/home/cam1/dft_runs/tho2/lda/c500_k20/irreducible_elastic/strain_A1g-3\n",
    "/home/cam1/dft_runs/tho2/lda/c500_k20/irreducible_elastic/strain_A1g_Eg.0-2\n",
    "/home/cam1/dft_runs/tho2/lda/c500_k20/irreducible_elastic/strain_A1g_xx-2\n",
    "```\n",
    "\n",
    "$$\\Large\n",
    "\\begin{align}\n",
    "\\hspace{-4mm}\n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{A_{1g}}^2}=  170 \\,\\,\\textrm{eV/uc}\n",
    "&& \n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{E_{g}^{(0)}}^2}= 67 \\,\\,\\textrm{eV/uc}\n",
    "&& \n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{T_{2g}}^2}= 25.5 \\,\\,\\textrm{eV/uc}\n",
    "&& \n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{xx}^2}= 101.4 \\,\\,\\textrm{eV/uc}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The latter derivative, with respect to $\\epsilon_{xx}$, can be obtained from the previous ones via the chain rule:\n",
    "\n",
    "$$\\Large\n",
    "\\begin{align}\n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{xx}^2}&=\n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{A_{1g}}^2}\n",
    "\\left(\n",
    "\\frac{\\partial \\epsilon_{A_{1g}}}{\\partial \\epsilon_{xx}} \n",
    "\\right)^2\n",
    "+\n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{E_{g}}^2}\n",
    "\\left[\n",
    "\\left(\n",
    "\\frac{\\partial \\epsilon_{E_{g}^{(0)}}}{\\partial \\epsilon_{xx}} \n",
    "\\right)^2\n",
    "+\\left(\n",
    "\\frac{\\partial \\epsilon_{E_{g}^{(1)}}}{\\partial \\epsilon_{xx}} \n",
    "\\right)^2\n",
    "\\right]\n",
    "\\\\&=\n",
    "\\frac{1}{3}\n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{A_{1g}}^2}\n",
    "+\n",
    "\\frac{2}{3}\n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{E_{g}}^2}\n",
    "=101.3 \\,\\,\\textrm{eV/uc}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We can see that it is nearly identical to the measured value.\n",
    "\n",
    "Now for nonlinear.\n",
    "\n",
    "$$\\Large\n",
    "\\begin{align}\n",
    "\\hspace{-3mm}\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}^3}= -957 \\,\\,\\textrm{eV/uc}\n",
    "&& \\hspace{0mm} &&\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}\\partial \\epsilon_{E_{g}^{(0)}}^2}= -108 \\,\\,\\textrm{eV/uc}\n",
    "&& \\hspace{0mm} &&\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}\\partial \\epsilon_{xx}^2}= -391 \\,\\,\\textrm{eV/uc}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Once again, we can obtain the latter derivative from the preceding two via the chain rule:\n",
    "\n",
    "$$\\Large\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}\\partial \\epsilon_{xx}^2}=\n",
    " \\frac{1}{3} \\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}^3} \n",
    " +\\frac{2}{3}\\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}\\partial \\epsilon_{E_{g}^{(0)}}^2}\n",
    " = -391 \\,\\,\\textrm{eV/uc}\n",
    "$$\n",
    "\n",
    "The result obtained from the chain rule gives essentially identical results as compared to the direct measurement. \n",
    "\n",
    "If we want to convert to volumetric strain, we have:\n",
    "\n",
    "$$\\Large\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{v}\\partial \\epsilon_{xx}^2}=\n",
    "\\frac{1}{\\sqrt3}\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}\\partial \\epsilon_{xx}^2}\n",
    "=-226 \\,\\,\\textrm{eV/uc}\n",
    "$$\n",
    "\n",
    "The above result does not agree with what you would get if you manually computed the second derivative at a series of different volumes, given that such a procedure would change the reference of the strain. However, this is actually the quantity we want, because this is what is measured in experiment. In this case, you need to *rebase* the derivative. In other words, you need to renormalize strain such that it is measured from the expanded coordinate system (see course notes for derivation). We will denote the rebased strain as $\\mathcal{E}$. In this particular case, we have: \n",
    "\n",
    "$$\\Large\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}\\partial \\mathcal{E}_{xx}^2}\n",
    "=\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}\\partial \\epsilon_{xx}^2}\n",
    "+ \\frac{2}{\\sqrt3} \n",
    "\\frac{\\partial^2 E}{\\partial \\epsilon_{xx}^2}\n",
    "= -273.9 \\,\\,\\textrm{eV/uc}\n",
    "$$\n",
    "\n",
    "Using the chain rule, we can then move to volumetric strain:\n",
    "\n",
    "$$\\Large\n",
    "\\frac{\\partial^3 E}{\\partial \\epsilon_{v}\\partial \\mathcal{E}_{xx}^2}=\n",
    "\\frac{1}{\\sqrt3} \\frac{\\partial^3 E}{\\partial \\epsilon_{A_{1g}}\\partial \\mathcal{E}_{xx}^2}\n",
    "= -158.1 \\,\\,\\textrm{eV/uc}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}