LAMMPS-Tutorials-02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LAMMPS Tutorials 02. Using variables and iterating in LAMMPS!\n",
    "\n",
    "### Author: Mark Tschopp, mark.a.tschopp.civ at mail.mil\n",
    "\n",
    "Please contact me if you have a problem with this tutorial, so I can modify in Github.  I have added FAQs, and will update my versions of LAMMPS in the future to keep the scripts current.\n",
    "\n",
    "The latest version of this [Jupyter Notebook](http://ipython.org/notebook.html) tutorial is available at https://github.com/mrkllntschpp/lammps-tutorials.\n",
    "\n",
    "The original tutorials are given here: https://icme.hpc.msstate.edu/mediawiki/index.php/LAMMPS_tutorials.  A number of these tutorials are out of date and have been ported over into the current iPython Jupyter Notebook tutorials on github.\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Description:\n",
    "<a id=\"Sec1\"></a>\n",
    "This is part two of a tutorial for running LAMMPS simulation on a Windows machine. For this example, the molecular dynamics simulation calculates the equilibrium lattice constant and the corresponding cohesive energy for aluminum. This tutorial will introduce the use of variables via command line and using MATLAB/Python for running LAMMPS. \n",
    "\n",
    "<img src=\"https://icme.hpc.msstate.edu/mediawiki/images/e/ef/Fcc_stereo.gif\" alt=\"Face-centered Cubic Lattice Structure\" title=\"FCC Lattice Structure\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complete Tutorial One \n",
    "\n",
    "If you have not done so already, complete the first tutorial available [here](https://github.com/mrkllntschpp/lammps-tutorials). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Step 1: Modify the input file \n",
    "\n",
    "This input script was run using the 24Jan2020 version of LAMMPS. Changes in some commands may require revision of the input script. In the below input script, the initial lattice constant (4) was replaced with the variable **latconst** (\\${latconst}). Either modify the initial script `calc_fcc.in` or copy the text below and paste it into a new text file. Name the modified file: `calc_fcc_ver1.in`. If pasting, use the `Paste Special` command with `Unformatted Text`. Note, this requires the variable **latconst** to be manually passed in from the command line when executing lammps, e.g. `lmp_win_no-mpi -in calc_fcc_ver1.in -var latconst 4`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting calc_fcc_ver1.in\n"
     ]
    }
   ],
   "source": [
    "%%writefile calc_fcc_ver1.in\n",
    "######################################\n",
    "# LAMMPS INPUT SCRIPT\n",
    "# Find minimum energy fcc configuration\n",
    "# Mark Tschopp\n",
    "# This requires the variable latconst to be input via the command line\n",
    "# e.g., lmp_win_no-mpi -var latconst 4 < calc_fcc_ver1.in\n",
    "\n",
    "######################################\n",
    "# INITIALIZATION\n",
    "clear \n",
    "units metal \n",
    "dimension 3 \n",
    "boundary p p p \n",
    "atom_style atomic \n",
    "atom_modify map array\n",
    "\n",
    "######################################\n",
    "# ATOM DEFINITION\n",
    "lattice fcc ${latconst}\n",
    "region box block 0 1 0 1 0 1 units lattice\n",
    "create_box 1 box\n",
    "\n",
    "lattice fcc ${latconst} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1  \n",
    "create_atoms 1 box\n",
    "replicate 1 1 1\n",
    "\n",
    "######################################\n",
    "# DEFINE INTERATOMIC POTENTIAL\n",
    "pair_style eam/alloy \n",
    "pair_coeff * * Al99.eam.alloy Al\n",
    "neighbor 2.0 bin \n",
    "neigh_modify delay 10 check yes \n",
    " \n",
    "######################################\n",
    "# DEFINE COMPUTES \n",
    "compute eng all pe/atom \n",
    "compute eatoms all reduce sum c_eng \n",
    "\n",
    "#####################################################\n",
    "# MINIMIZATION\n",
    "reset_timestep 0 \n",
    "fix 1 all box/relax iso 0.0 vmax 0.001\n",
    "thermo 10 \n",
    "thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms \n",
    "min_style cg \n",
    "minimize 1e-25 1e-25 5000 10000 \n",
    "\n",
    "variable natoms equal \"count(all)\" \n",
    "variable teng equal \"c_eatoms\"\n",
    "variable length equal \"lx\"\n",
    "variable ecoh equal \"v_teng/v_natoms\"\n",
    "\n",
    "print \"Total energy (eV) = ${teng};\"\n",
    "print \"Number of atoms = ${natoms};\"\n",
    "print \"Lattice constant (Angstoms) = ${length};\"\n",
    "print \"Cohesive energy (eV) = ${ecoh};\"\n",
    "\n",
    "print \"All done!\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Step 2: Run in LAMMPS \n",
    "OK. Get to the command prompt and run the script that we just wrote.  It should look like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LAMMPS (24 Jan 2020)\n",
      "OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:94)\n",
      "  using 1 OpenMP thread(s) per MPI task\n",
      "OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:94)\n",
      "  using 1 OpenMP thread(s) per MPI task\n",
      "Lattice spacing in x,y,z = 4 4 4\n",
      "Created orthogonal box = (0 0 0) to (4 4 4)\n",
      "  1 by 1 by 1 MPI processor grid\n",
      "Lattice spacing in x,y,z = 4 4 4\n",
      "Created 4 atoms\n",
      "  create_atoms CPU = 0 secs\n",
      "Replicating atoms ...\n",
      "  orthogonal box = (0 0 0) to (4 4 4)\n",
      "  1 by 1 by 1 MPI processor grid\n",
      "  4 atoms\n",
      "  replicate CPU = 0 secs\n",
      "WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:178)\n",
      "Neighbor list info ...\n",
      "  update every 1 steps, delay 0 steps, check yes\n",
      "  max neighbors/atom: 2000, page size: 100000\n",
      "  master list distance cutoff = 8.28721\n",
      "  ghost atom cutoff = 8.28721\n",
      "  binsize = 4.1436, bins = 1 1 1\n",
      "  1 neighbor lists, perpetual/occasional/extra = 1 0 0\n",
      "  (1) pair eam/alloy, perpetual\n",
      "      attributes: half, newton on\n",
      "      pair build: half/bin/atomonly/newton\n",
      "      stencil: half/bin/3d/newton\n",
      "      bin: standard\n",
      "Setting up cg style minimization ...\n",
      "  Unit style    : metal\n",
      "  Current step  : 0\n",
      "WARNING: Energy due to 1 extra global DOFs will be included in minimizer energies\n",
      "Per MPI rank memory allocation (min/avg/max) = 4.495 | 4.495 | 4.495 Mbytes\n",
      "Step PotEng Lx Ly Lz Press Pxx Pyy Pzz c_eatoms \n",
      "       0   -13.417787            4            4            4     29590.11     29590.11     29590.11     29590.11   -13.417787 \n",
      "      10   -13.439104         4.04         4.04         4.04    5853.9553    5853.9553    5853.9553    5853.9553   -13.439104 \n",
      "      20       -13.44    4.0500047    4.0500047    4.0500047 2.1131131e-10 2.2023682e-10 2.142121e-10 1.9948502e-10       -13.44 \n",
      "      29       -13.44    4.0500047    4.0500047    4.0500047 4.7215918e-10 4.8063841e-10 4.5787837e-10 4.7796075e-10       -13.44 \n",
      "Loop time of 0.029984 on 1 procs for 29 steps with 4 atoms\n",
      "\n",
      "0.0% CPU use with 1 MPI tasks x 1 OpenMP threads\n",
      "\n",
      "Minimization stats:\n",
      "  Stopping criterion = energy tolerance\n",
      "  Energy initial, next-to-last, final = \n",
      "        -13.4177872966     -13.4399999527     -13.4399999527\n",
      "  Force two-norm initial, final = 3.54599 5.81382e-14\n",
      "  Force max component initial, final = 3.54599 5.80057e-14\n",
      "  Final line search alpha, max atom move = 1 5.80057e-14\n",
      "  Iterations, force evaluations = 29 46\n",
      "\n",
      "MPI task timing breakdown:\n",
      "Section |  min time  |  avg time  |  max time  |%varavg| %total\n",
      "---------------------------------------------------------------\n",
      "Pair    | 0.0049961  | 0.0049961  | 0.0049961  |   0.0 | 16.66\n",
      "Neigh   | 0          | 0          | 0          |   0.0 |  0.00\n",
      "Comm    | 0.0050011  | 0.0050011  | 0.0050011  |   0.0 | 16.68\n",
      "Output  | 0          | 0          | 0          |   0.0 |  0.00\n",
      "Modify  | 0          | 0          | 0          |   0.0 |  0.00\n",
      "Other   |            | 0.01999    |            |       | 66.66\n",
      "\n",
      "Nlocal:    4 ave 4 max 4 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "Nghost:    662 ave 662 max 662 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "Neighs:    280 ave 280 max 280 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "\n",
      "Total # of neighbors = 280\n",
      "Ave neighs/atom = 70\n",
      "Neighbor list builds = 0\n",
      "Dangerous builds = 0\n",
      "Total energy (eV) = -13.4399999527352;\n",
      "Number of atoms = 4;\n",
      "Lattice constant (Angstoms) = 4.05000466178543;\n",
      "Cohesive energy (eV) = -3.3599999881838;\n",
      "All done!\n",
      "Total wall time: 0:00:00\n"
     ]
    }
   ],
   "source": [
    "!\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -var latconst 4.0 < calc_fcc_ver1.in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try modifying the starting lattice constant from the command line, *e.g.*, to start with a lattice constant of 3.0 Angstroms, type 'lmp_serial -var latconst 3.0 < calc_fcc_ver1.in'. This will not affect the final values (above), but it will start the simulation with a lower lattice constant for the FCC structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LAMMPS (24 Jan 2020)\n",
      "OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:94)\n",
      "  using 1 OpenMP thread(s) per MPI task\n",
      "OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:94)\n",
      "  using 1 OpenMP thread(s) per MPI task\n",
      "Lattice spacing in x,y,z = 3 3 3\n",
      "Created orthogonal box = (0 0 0) to (3 3 3)\n",
      "  1 by 1 by 1 MPI processor grid\n",
      "Lattice spacing in x,y,z = 3 3 3\n",
      "Created 4 atoms\n",
      "  create_atoms CPU = 0 secs\n",
      "Replicating atoms ...\n",
      "  orthogonal box = (0 0 0) to (3 3 3)\n",
      "  1 by 1 by 1 MPI processor grid\n",
      "  4 atoms\n",
      "  replicate CPU = 0 secs\n",
      "WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:178)\n",
      "Neighbor list info ...\n",
      "  update every 1 steps, delay 0 steps, check yes\n",
      "  max neighbors/atom: 2000, page size: 100000\n",
      "  master list distance cutoff = 8.28721\n",
      "  ghost atom cutoff = 8.28721\n",
      "  binsize = 4.1436, bins = 1 1 1\n",
      "  1 neighbor lists, perpetual/occasional/extra = 1 0 0\n",
      "  (1) pair eam/alloy, perpetual\n",
      "      attributes: half, newton on\n",
      "      pair build: half/bin/atomonly/newton\n",
      "      stencil: half/bin/3d/newton\n",
      "      bin: standard\n",
      "Setting up cg style minimization ...\n",
      "  Unit style    : metal\n",
      "  Current step  : 0\n",
      "WARNING: Energy due to 1 extra global DOFs will be included in minimizer energies\n",
      "Per MPI rank memory allocation (min/avg/max) = 4.502 | 4.502 | 4.502 Mbytes\n",
      "Step PotEng Lx Ly Lz Press Pxx Pyy Pzz c_eatoms \n",
      "       0     9.965675            3            3            3    3458818.8    3458818.8    3458818.8    3458818.8     9.965675 \n",
      "      10    8.2278955         3.03         3.03         3.03    3326311.7    3326311.7    3326311.7    3326311.7    8.2278955 \n",
      "      20    6.5407146         3.06         3.06         3.06      3152679      3152679      3152679      3152679    6.5407146 \n",
      "      30    4.9138013         3.09         3.09         3.09    2971490.2    2971490.2    2971490.2    2971490.2    4.9138013 \n",
      "      40    3.3576717         3.12         3.12         3.12    2772390.3    2772390.3    2772390.3    2772390.3    3.3576717 \n",
      "      50    1.8866582         3.15         3.15         3.15    2554805.1    2554805.1    2554805.1    2554805.1    1.8866582 \n",
      "      60   0.51216154         3.18         3.18         3.18    2331874.9    2331874.9    2331874.9    2331874.9   0.51216154 \n",
      "      70  -0.76371356         3.21         3.21         3.21      2122446      2122446      2122446      2122446  -0.76371356 \n",
      "      80   -1.9484596         3.24         3.24         3.24    1938269.5    1938269.5    1938269.5    1938269.5   -1.9484596 \n",
      "      90   -3.0533151         3.27         3.27         3.27    1778431.7    1778431.7    1778431.7    1778431.7   -3.0533151 \n",
      "     100   -4.0877221          3.3          3.3          3.3    1637176.2    1637176.2    1637176.2    1637176.2   -4.0877221 \n",
      "     110   -5.0566619         3.33         3.33         3.33    1502862.8    1502862.8    1502862.8    1502862.8   -5.0566619 \n",
      "     120   -5.9590126         3.36         3.36         3.36    1368469.7    1368469.7    1368469.7    1368469.7   -5.9590126 \n",
      "     130   -6.7911412         3.39         3.39         3.39    1233014.6    1233014.6    1233014.6    1233014.6   -6.7911412 \n",
      "     140   -7.5508303         3.42         3.42         3.42    1101952.7    1101952.7    1101952.7    1101952.7   -7.5508303 \n",
      "     150   -8.2406817         3.45         3.45         3.45    981949.61    981949.61    981949.61    981949.61   -8.2406817 \n",
      "     160   -8.8658333         3.48         3.48         3.48    874574.93    874574.93    874574.93    874574.93   -8.8658333 \n",
      "     170   -9.4330444         3.51         3.51         3.51     781176.4     781176.4     781176.4     781176.4   -9.4330444 \n",
      "     180   -9.9493468         3.54         3.54         3.54    700382.06    700382.06    700382.06    700382.06   -9.9493468 \n",
      "     190   -10.420951         3.57         3.57         3.57    629924.78    629924.78    629924.78    629924.78   -10.420951 \n",
      "     200    -10.85267          3.6          3.6          3.6    567324.29    567324.29    567324.29    567324.29    -10.85267 \n",
      "     210   -11.247853         3.63         3.63         3.63    510216.16    510216.16    510216.16    510216.16   -11.247853 \n",
      "     220   -11.608415         3.66         3.66         3.66    456609.65    456609.65    456609.65    456609.65   -11.608415 \n",
      "     230   -11.935085         3.69         3.69         3.69     404837.7     404837.7     404837.7     404837.7   -11.935085 \n",
      "     240   -12.227243         3.72         3.72         3.72    352930.13    352930.13    352930.13    352930.13   -12.227243 \n",
      "     250   -12.483171         3.75         3.75         3.75    300228.65    300228.65    300228.65    300228.65   -12.483171 \n",
      "     260   -12.701267         3.78         3.78         3.78    248104.16    248104.16    248104.16    248104.16   -12.701267 \n",
      "     270   -12.882122         3.81         3.81         3.81    200021.24    200021.24    200021.24    200021.24   -12.882122 \n",
      "     280   -13.028802         3.84         3.84         3.84    158291.05    158291.05    158291.05    158291.05   -13.028802 \n",
      "     290   -13.146194         3.87         3.87         3.87    124258.04    124258.04    124258.04    124258.04   -13.146194 \n",
      "     300   -13.239431          3.9          3.9          3.9    96679.419    96679.419    96679.419    96679.419   -13.239431 \n",
      "     310   -13.312518         3.93         3.93         3.93    73797.678    73797.678    73797.678    73797.678   -13.312518 \n",
      "     320   -13.368221         3.96         3.96         3.96     54117.68     54117.68     54117.68     54117.68   -13.368221 \n",
      "     330   -13.408037         3.99         3.99         3.99    35673.589    35673.589    35673.589    35673.589   -13.408037 \n",
      "     340    -13.43198         4.02         4.02         4.02    17614.659    17614.659    17614.659    17614.659    -13.43198 \n",
      "     350       -13.44         4.05         4.05         4.05     2.726913     2.726913     2.726913     2.726913       -13.44 \n",
      "     360       -13.44    4.0500047    4.0500047    4.0500047 -1.5619633e-10 -1.7003086e-10 -1.512873e-10 -1.4727082e-10       -13.44 \n",
      "     366       -13.44    4.0500047    4.0500047    4.0500047 -2.7892201e-10 -3.025746e-10 -2.5169922e-10 -2.8249221e-10       -13.44 \n",
      "Loop time of 0.064966 on 1 procs for 366 steps with 4 atoms\n",
      "\n",
      "24.1% CPU use with 1 MPI tasks x 1 OpenMP threads\n",
      "\n",
      "Minimization stats:\n",
      "  Stopping criterion = energy tolerance\n",
      "  Energy initial, next-to-last, final = \n",
      "         9.96567496491     -13.4399999527     -13.4399999527\n",
      "  Force two-norm initial, final = 174.865 2.58036e-14\n",
      "  Force max component initial, final = 174.865 2.56996e-14\n",
      "  Final line search alpha, max atom move = 1 2.56996e-14\n",
      "  Iterations, force evaluations = 366 382\n",
      "\n",
      "MPI task timing breakdown:\n",
      "Section |  min time  |  avg time  |  max time  |%varavg| %total\n",
      "---------------------------------------------------------------\n",
      "Pair    | 0.014991   | 0.014991   | 0.014991   |   0.0 | 23.07\n",
      "Neigh   | 0          | 0          | 0          |   0.0 |  0.00\n",
      "Comm    | 0          | 0          | 0          |   0.0 |  0.00\n",
      "Output  | 0          | 0          | 0          |   0.0 |  0.00\n",
      "Modify  | 0          | 0          | 0          |   0.0 |  0.00\n",
      "Other   |            | 0.04998    |            |       | 76.93\n",
      "\n",
      "Nlocal:    4 ave 4 max 4 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "Nghost:    662 ave 662 max 662 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "Neighs:    400 ave 400 max 400 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "\n",
      "Total # of neighbors = 400\n",
      "Ave neighs/atom = 100\n",
      "Neighbor list builds = 1\n",
      "Dangerous builds = 0\n",
      "Total energy (eV) = -13.4399999527351;\n",
      "Number of atoms = 4;\n",
      "Lattice constant (Angstoms) = 4.05000466178543;\n",
      "Cohesive energy (eV) = -3.35999998818377;\n",
      "All done!\n",
      "Total wall time: 0:00:00\n"
     ]
    }
   ],
   "source": [
    "!\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -var latconst 3.0 < calc_fcc_ver1.in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Step 3: Iterate and Run in MATLAB \n",
    "\n",
    "Now, the exciting part! In this section, MATLAB will be used to run a LAMMPS simulation in Windows. Follow these steps: \n",
    "1. Open MATLAB \n",
    "1. Change to the directory that contains the LAMMPS executable (`lmp_serial.exe`), the input script (`calc_fcc_ver1.in`), and the potential file (`Al99.eam.alloy`). \n",
    "1. Type `system('lmp_serial -var latconst 4 < calc_fcc_ver1.in');` \n",
    "1. The simulation should show up on the MATLAB screen. When it completes, it will display \"All done\"! \n",
    "\n",
    "Now that you can run a simulation from MATLAB, it is easy to loop over a number of different lattice constants to calculate the energy. For instance, modify the input script again to only display the energy for the lattice constant input through the command line. The script should look like this: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting calc_fcc_ver2.in\n"
     ]
    }
   ],
   "source": [
    "%%writefile calc_fcc_ver2.in\n",
    "######################################\n",
    "# LAMMPS INPUT SCRIPT\n",
    "# Find minimum energy fcc configuration\n",
    "# Mark Tschopp\n",
    "# This requires the variable latconst to be input via the command line\n",
    "# e.g., lmp_win_no-mpi -var latconst 4 < calc_fcc_ver1.in\n",
    "\n",
    "######################################\n",
    "# INITIALIZATION\n",
    "clear \n",
    "units metal \n",
    "dimension 3 \n",
    "boundary p p p \n",
    "atom_style atomic \n",
    "atom_modify map array\n",
    "\n",
    "######################################\n",
    "# ATOM DEFINITION\n",
    "lattice fcc ${latconst} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1  \n",
    "region box block 0 1 0 1 0 1 units lattice\n",
    "create_box 1 box\n",
    "create_atoms 1 box\n",
    "replicate 1 1 1\n",
    "\n",
    "######################################\n",
    "# DEFINE INTERATOMIC POTENTIAL\n",
    "pair_style eam/alloy \n",
    "pair_coeff * * Al99.eam.alloy Al\n",
    "neighbor 2.0 bin \n",
    "neigh_modify delay 10 check yes \n",
    " \n",
    "######################################\n",
    "# RUN 0\n",
    "run 0\n",
    "\n",
    "variable natoms equal \"count(all)\" \n",
    "variable teng equal \"pe\"\n",
    "variable length equal \"lx\"\n",
    "variable ecoh equal \"v_teng/v_natoms\"\n",
    "\n",
    "print \"Total energy (eV) = ${teng};\"\n",
    "print \"Number of atoms = ${natoms};\"\n",
    "print \"Lattice constant (Angstoms) = ${length};\"\n",
    "print \"Cohesive energy (eV) = ${ecoh};\"\n",
    "print \"%% ecoh = ${ecoh};\"\n",
    "\n",
    "print \"All done!\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `run 0` command can be used to calculate the energy of the system without actually running an iteration. Name this `calc_fcc_ver2.in` and combine this with a simple MATLAB script like this: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```OCTAVE\n",
    "% MATLAB Script for running LAMMPS multiple times\n",
    "\n",
    "count = 0;\n",
    "for i = 3.0:0.10:5.0\n",
    "    command_line = ['lmp_win_no-mpi -var latconst ' num2str(i) ' < calc_fcc_ver2.in'];\n",
    "\n",
    "    % this next line executes the command line\n",
    "    system(command_line)\n",
    "\n",
    "    % all that is left is to mine the 'log.lammps' file for the energy\n",
    "    fid = fopen('log.lammps');\n",
    "        tline = fgetl(fid);\n",
    "        while ~feof(fid)\n",
    "           matches = strfind(tline, '%%');\n",
    "           num = length(matches);\n",
    "           if num > 0 && matches == 1\n",
    "                teval = strrep(tline,'%%','');\n",
    "                eval(teval)\n",
    "           end\n",
    "           tline = fgetl(fid);\n",
    "        end\n",
    "    fclose(fid);\n",
    "\n",
    "   % store the values in a matrix\n",
    "   count = count + 1;\n",
    "   X(count) = i; Y(count) = ecoh;\n",
    "end\n",
    "\n",
    "plot(X,Y,'-*r'), axis square\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This script increases the lattice constant (variable `i`) from 3.0 to 5.0 Angstroms in increments of 0.10 Angstroms (`3.0:0.10:0.50`) AND calculates the energy for each lattice constant. As an easy extension, MATLAB can now be used to run hundreds of simulations, extract energies, and plot the lattice constant-energy curve in a matter of a minute or so. Notice the line `print \"%% ecoh = ${ecoh};\"` in the `calc_fcc_ver2.in` script.  The '\\%\\%' is used as a flag for when MATLAB mines the `log.lammps` file. This MATLAB script iteratively reads every line in the `log.lammps` file and checks to see if the flag \\%\\% is contained within each line. When MATLAB finds a line with the flag, it removes the flag, and evaluates (`eval`) the rest of the line in MATLAB (which converts the string to a variable in the MATLAB stack). Notice that vector `X` in MATLAB contains the lattice constant that is passed to LAMMPS and vector `Y` in MATLAB now contains the potential energy `ecoh` that was computed in LAMMPS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Run Iteratively in Python\n",
    "\n",
    "Now let's try it in Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.00 -log log_000.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.10 -log log_001.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.20 -log log_002.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.30 -log log_003.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.40 -log log_004.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.50 -log log_005.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.60 -log log_006.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.70 -log log_007.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.80 -log log_008.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.90 -log log_009.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.00 -log log_010.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.10 -log log_011.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.20 -log log_012.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.30 -log log_013.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.40 -log log_014.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.50 -log log_015.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.60 -log log_016.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.70 -log log_017.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.80 -log log_018.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.90 -log log_019.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 5.00 -log log_020.lammps\n",
      "[3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0]\n",
      "--- 0.4868137836456299 seconds ---\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import subprocess \n",
    "import shlex\n",
    "import time\n",
    "\n",
    "start_time = time.time()\n",
    "j = 0\n",
    "x = []\n",
    "for i in list(np.arange(3.00, 5.05, 0.10)):\n",
    "    command_line = f'\\\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\\\" \\\n",
    "-in calc_fcc_ver2.in -var latconst {i:0.02f} -log log_{j:03d}.lammps'\n",
    "    x.append(round(i,3))\n",
    "    j += 1\n",
    "    print(command_line)\n",
    "    args = shlex.split(command_line)\n",
    "#    p = subprocess.call(args) # Success! but much slower...\n",
    "    p = subprocess.Popen(args) # Success!\n",
    "    \n",
    "print(x)\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow!  That ran quickly.  And, it created a bunch of log files.  I purposedly separated the different parts of this process in Python.  Now, we are going to read the cohesive energy for each lattice constant into a variable and store it in a vector.  While I removed the `;` in the flagged line for execution in Python, I could have just modified the script so that it was Python executable, as opposed to MATLAB executable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- 0.06231498718261719 seconds ---\n"
     ]
    }
   ],
   "source": [
    "# how do I read each line in script, find the flag, and assign ecoh to a vector Y for plotting\n",
    "\n",
    "import glob, os\n",
    "import time\n",
    "\n",
    "start_time = time.time()\n",
    "target = \"%%\"\n",
    "y = []\n",
    "for file in glob.glob(\"log_*.lammps\"):\n",
    "    with open(file) as f:\n",
    "        for line in f:\n",
    "            if target in line and \"print\" not in line:\n",
    "                exec(line.replace(\"%% \",\"\").replace(\";\",\"\"))\n",
    "                y.append(ecoh)\n",
    "\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK.  Now that I have both of those, let's put everything together and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.00 -log log_000.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.10 -log log_001.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.20 -log log_002.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.30 -log log_003.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.40 -log log_004.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.50 -log log_005.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.60 -log log_006.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.70 -log log_007.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.80 -log log_008.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 3.90 -log log_009.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.00 -log log_010.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.10 -log log_011.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.20 -log log_012.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.30 -log log_013.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.40 -log log_014.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.50 -log log_015.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.60 -log log_016.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.70 -log log_017.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.80 -log log_018.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 4.90 -log log_019.lammps\n",
      "\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\" -in calc_fcc_ver2.in -var latconst 5.00 -log log_020.lammps\n",
      "--- 11.010871171951294 seconds ---\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd7wcVfnH8c83PQESkIQQIBCaQCgKBEVa6E1qgpjQAxq6IAqI2Iii/sSKFCkiTZpIpAkIEkroAYEkFMGESElCQDqEtOf3x5nlbpZb9ube3dl79/t+vfZ1Z2anPDuZzDNzzsw5igjMzKz+dMk7ADMzy4cTgJlZnXICMDOrU04AZmZ1ygnAzKxOOQGYmdUpJ4AORNJhkiZWaN1/kPT9Sqy7GiRdKukn7bzOrSU9357rtE+TdI+kr7Vh+fclrdGeMdWLbnkHYLUhIo7KO4ZaExH3A+vkHYc1kHQPcGVEXFyYFhFL5xdRx+Y7ALNOSpIv8KxZTgA1RtJ3JP1H0nuSnpG0bxPzDZEUxf/Ji2+ls+KiByT9RtLbkqZJ2iKb/rKk1yUdWrTsJ0UokraV9Iqkb2XzzZQ0prHtFG1rYtF4SDpG0gvZ7/ixpDUlPSTpXUnXSerRzD74uqRni/bBJtn09bJtvy1pqqS9ShZdTtKt2XKPSFqzaJ3rSrpT0v8kPS9p/6Lvds+2856kVyV9u3g/FP27XF8S5+8knZ0N95P0x2xfvSrpJ5K6NvH7uhT9O7+Z7Y/PlPy7Hirpv5LekHR6K5c9QtJ/gbuz6YdImpHN/31JL0naUdKKkj6UtHzR+jeVNEdS91bGfbuk40rmf0rSiGx4C0mPSXon+7tFE/vmR5KuLBr/5DiXdCawNXCOUrHPOdk8IWmton+Hy7PfMEPS9yR1yb47TNJESb+U9Jak6ZJ2ayyOeuEEUHv+QzrI+wFnAFdKGrSE6/oi8DSwPHAVcA2wGbAWcBDpP1JTt88rZjGsDBwBnCtpuVZse1dgU2Bz4BTgQuBAYDCwATC6sYUkfQX4EXAI0BfYC3gzOyHdDPwDWAE4HvizpOIimtGkfbYc8CJwZrbOpYA7s32wQjbfeZLWz5b7I3BkRCyTxXZ3I6FdDewuqW+2zq7A/tk6AS4DFpD27cbAzkBT5drfAPYBhgMrAW8B55bMsxWp+GkH4AeS1mvFssOB9YBdJA0FziPt+0E0/JsSEbOAe7LfUXAQcE1EzG9l3FdR9G+abXc14NYsSdwKnE06Fn+dTV+eVoiI04H7geMiYumIOK6R2X6f/cY1sjgPAcYUff9F4HmgP/AL4I+S1Jo4OpWI8KeGP8CTwN7Z8GHAxGx4CBBAt6J57wG+VjTvC0XfbZjNP7Bo2pvA57PhS4GfZMPbAh+VrPt1YPPS7ZTGlY0HsGXR+OPAqUXjvwJ+28TvvQM4oZHpWwOzgC5F064GflQU/8VF3+0OPJcNfxW4v2R9FwA/zIb/CxwJ9C2ZZ1vglaLxicAh2fBOwH+y4YHAx0DvonlHAxOa+I3PAjsUjQ8C5pPq5Ar/rqsUff8oMKoVy65R9P0PgKuLxvsA84Adi/bNA9lw12wff2EJ4l4G+ABYLfvuTOCSbPhg4NGSdT0EHNbIcfsjUhl/Yb7Cb+rW2LFXdLytlcX/MTC06LsjgXuKjtMXS/ZFACvm/f88r4/vAGpMdrv+ZFbM8TbpirT/Eq5udtHwRwARUTqtqTuANyNiQdH4h83MW862y93uYNJdUKmVgJcjYlHRtBlkV7OZWUXDxfGuBnyxsE+z/Xog6S4HYCQpYcyQdK+kLzURW/FV7gE0XP2vBnQHZhat/wLS3UZjVgPGF837LLCQlEjK+S0tLfty0fBKxeMR8SEp8RfcCAxVeopmJ+CdiHi0tXFHxHukq/xR2byjgD8XxTCjZF2l/3btoT/Qo2RbTR4j2b6A1h3XnYoTQA2RtBpwEXAcsHxELAtMARq7Rf0g+9unaNqKjcxXCR9UcLsvA2s2Mv01YHChPDezKvBqmeu8NyKWLfosHRFHA0TEYxGxN+mE/TfguibW8xdgW0mrAPvSkABeJl159i9af9+IWL+J9bwM7FYST6+IKPe3tLRscRO/M4FVCiOSepOKYch++9zs9x5IulK/og3bvhoYnSXQ3sCEbPprpORRrKl/u5aOreaaL36DdEdSvK1yj5G65ARQW5YiHeBzAJQqXjdobMaImEM6sA+S1FXS4TR+4qyEJ4ERkvpklW9HtOO6Lwa+nVVGStJaWWJ8hHRyOEVSd0nbAnuS6jVacgvwWUkHZ8t2l7SZUqVyD0kHSuoXqdz7XdJV7adk+/we4E/A9Ih4Nps+k1Q38StJfbPK0jUlDW8inj8AZ2a/C0kDJO1d1t5p/bLXA3tmlbA9SHUkpRcUl5OKR/YCrqRpLW3776ST7zjg2qK7tb+T9v8BWWXuV4GhpH+XUk8C20haVVI/4LSS72eTyvc/JSIWkpLZmZKWyeI8qYXfVNecAGpIRDxDKh9/iHSgbwg80MwiXwdOJt3Srw88WOkYM78hlSPPJlV+/rn52csXEX8hlR9fBbxHuiL/TETMI52gdiNd6Z1HKo9/rox1vkeqlB1FuhqdBfwf0DOb5WDgJUnvAkeRKkKbchWwIw1X/wWHkIofniFVjl5PKiNvzO+Am4B/SHoPeJhUOVmOVi0bEVNJFebXkO4G3iPV53xcNM8DwCLgiYh4aUm3HREfAzdQsn8i4k1gD+BbpGP1FGCPiHijkXjvBK4lPbzwOJ9OEr8D9sue4jm7kRiPJ10oTCPV2VwFXNLMb6pryipDzKwOZE99vQ2sHRHTi6bfDVwVRS9YWefnOwCzTk7Snllx3VLAL4HJwEtF328GbEK68rY64gRg1vntTSr6eg1Ym/RIaQBIugy4CzgxKyqzOuIiIDOzOuU7ADOzOtWhGovq379/DBkyJO8wzMw6lMcff/yNiBhQOr1DJYAhQ4YwadKkvMMwM+tQJJW+iQ24CMjMrG45AZiZ1SknADOzOuUEYGZWp5wAzMzqVH0kgJkzYfhwmDWr5XnNzOpEfSSAH/8YJk6EcePyjsTMrGZ07gTQuzdIcP75sGhR+iul6WZmda5zJ4Bp0+CAA6Bb9r5b795w4IEwfXrzy5mZ1YHOnQAGDYK+fWFh1sHT3LlpfMVq9ZxoZla7OncCAJg9G446CgYMgNVWc0WwmVmmQ7UFtERuuCH97doVLr4YLr8833jMzGpE578DKBg5MhUB3X573pGYmdWE+kkAW20F/fvDX/+adyRmZjWhfhJAt26wzz5wyy3pTsDMrM7VTwKAVAz0/vtw1115R2Jmlrv6SgDbbw/9+rkYyMyMeksAPXrAnnvCTTfB/Pl5R2Nmlqv6SgAAI0bA//4H996bdyRmZrmqvwSwyy7Qp0/D+wFmZnWq/hJAnz6w++4wfnxqIM7MrE7VXwKAVAw0axY8+GDekZiZ5Sa3BCBpsKQJkp6VNFXSCVXb+Je/nCqEXQxkZnUszzuABcC3ImI9YHPgWElDq7Llvn1hp51SAoioyibNzGpNbgkgImZGxBPZ8HvAs8DKVQtg5EiYMQMef7xqmzQzqyU1UQcgaQiwMfBII9+NlTRJ0qQ5c+a030b32iu1EOpiIDOrU7knAElLA38FToyId0u/j4gLI2JYRAwbMGBA+214+eVh223TW8EuBjKzOpRrApDUnXTy/3NEVP9SfORI+Pe/4Zlnqr5pM7O85fkUkIA/As9GxK9zCWKffVIn8W4byMzqUJ53AFsCBwPbS3oy++xe1QgGDYIttnACMLO6lFuXkBExEVBe2//EyJFw0knw4ouw1lp5R2NmVjW5VwLnbt99018/DWRmdcYJYMgQ2HRTJwAzqztOAJCKgR55BF55Je9IzMyqxgkAUuNw4LsAM6srTgAA66wD66/vBGBmdcUJoGDkSLj/fnj99bwjMTOrCieAghEjUgcxN96YdyRmZlXhBFCw0Uaw5pp+KczM6oYTQIGU7gL++U946628ozEzqzgngGIjR8KCBXDLLXlHYmZWcU4AxTbbDFZZxcVAZlYXnACKdemSmoa44w54//28ozEzqygngFIjR8LcuXDbbXlHYmZWUU4ApbbaCgYMcDGQmXV6TgClunZNHcXcemu6EzAz66ScABozcmSqA7jzzrwjMTOrGCeAxmy3HfTr57aBzKxTcwJoTI8esNdeqVmI+fPzjsbMrCKcAJoycmR6I/iee/KOxMysIpwAmrLzzrDUUi4GMrNOywmgKb17w+67w/jxsHBh3tGYmbU7J4DmjBgBs2enPoNnzco7GjOzduUE0Jwvfzk1D/HUUzBuXN7RmJm1KyeApvTuDX37pk5iAM4/PzUZ3bt3vnGZmbUTJ4CmTJsGBxyQHgkF6NkTDjwQpk/PNy4zs3bSrbkvJa0CjAK2BlYCPgKmALcCt0XEorZsXNIlwB7A6xGxQVvW1e4GDUp3AAsWpPGPP07jK66Yb1xmZu2kyTsASX8CLgHmAf8HjAaOAe4CdgUmStqmjdu/NFtXbZo9G446Kr0U1r07vPpq3hGZmbWb5u4AfhURUxqZPgW4QVIPYNW2bDwi7pM0pC3rqKjCOwD33Qc33QT77ZdvPGZm7ai5OoBdsyKgRkXEvIh4sQIxLUbSWEmTJE2aM2dOpTfXuK23Th3GX3JJPts3M6uA5hLAysBDku6TdLSk/tUKqlhEXBgRwyJi2IABA/IIIT39M2ZMahZi2rR8YjAza2dNJoCI+CapiOf7wEbA05Juk3SIpGWqFWDNOOSQlAguuyzvSMzM2kWzj4FGcm9EHA0MBn4LfBOYXY3gasrgwbDTTnDppQ3vBpiZdWBlvQcgaUNgHHAu6amg77bHxiVdDTwErCPpFUlHtMd6K2bMGPjvf+Huu/OOxMyszZp8CkjS2qR3AEYDC4FrgJ0jot0KwSNidHutqyr22QeWXRb+9CfYcce8ozEza5PmHgO9A7ga+GpETK5SPLWtV6/0dvAll8Dbb6dkYGbWQTVXCbxGRJweEZMlrSZpRwBJveuyErhgzJjUWfw11+QdiZlZm7RYByDp68D1wAXZpFWAv1UyqJq26aaw4YapGMjMrAMrpxL4WGBL4F2AiHgBWKGSQdW0wjsBjz4KU6fmHY2Z2RIrJwF8HBHzCiOSugFRuZA6gIMOgm7dfBdgZh1aOQngXknfBXpL2gn4C3BzZcOqcQMGwJ57whVXwPz5eUdjZrZEykkA3wHmAJOBI4G/A9+rZFAdwpgx8PrrcNtteUdiZrZEmu0PACBr8/+i7GMFu+0GAwemR0L32ivvaMzMWq25/gBulrSnpO6NfLeGpHGSDq9seDWsW7fUPtCtt6Y7ATOzDqa5IqCvk3oCe07SY5L+LuluSdNJj4Q+HhH13T7ymDGpx7Arr8w7EjOzVlNEyw/0ZJ22DCJ1CfnviPiwsmE1btiwYTFp0qQ8Nt20zTeH99+HyZPTI6JmZjVG0uMRMax0elmNwUXESxHxUEQ8mdfJv2aNGZPeB6i1xGRm1oKyEoA1Y9So1EaQ3wkwsw7GCaCt+vWDkSPh6qvho4/yjsbMrGzltAW0hyQniuaMGZNaB/1b/TaRZGYdTzkn9lHAC5J+IWm9SgfUIW23Hay2mouBzKxDaTEBRMRBwMbAf4A/SXpI0ti6bhK6VJcucNhhcNddqccwM7MOoNyngN4F/krqFWwQsC/whKTjKxhbx3LooRDhTuPNrMMopw5gT0njgbuB7sAXImI34HPAtyscX8ex+uqw/fbuNN7MOoxy7gC+AvwmIjaKiLMi4nWA7H2A+m0KojFjxsC0aXD//XlHYmbWonLqAA6JiPua+O6f7R9SBzZiBPTtmxqIMzOrceUUAb0n6d2Sz8uSxktaoxpBdhh9+qQXw66/Ht57L+9ozMyaVU4R0K+Bk4GVSf0Bf5vUNPQ1gC91S40ZAx9+CNddl3ckZmbNKicB7BoRF0TEexHxbkRcCOweEdcCy1U4vo7ni1+Eddd1MZCZ1bxyEsAiSftL6pJ99i/6rk19A0vaVdLzkl6U9J22rKtmFDqNf/BBeP75vKMxM2tSOQngQOBg4HVgdjZ8kKTewHFLumFJXYFzgd2AocBoSUOXdH015eCDoWvX9EiomVmNajYBZCfpvSNiz4joHxEDsuEXI+KjiJjYhm1/AXgxIqZFxDxSncLebVhf7Rg0KHUZedllqcMYM7Ma1GwCiIiFVO6kvDLwctH4K9m0xWTNTkySNGnOnDkVCqUCxoyBmTPhmmtg+HCYNSvviMzMFlNOEdADks6RtLWkTQqfdth2Y91nfapOISIujIhhETFswIAB7bDZKtljD+jfH374Q5g4EcaNyzsiM7PFdCtjni2yv8VnsAC2b+O2XwEGF42vArzWxnXWjn79YO5ceOONNH7++enTq5f7DTCzmlDOm8DbNfJp68kf4DFgbUmrS+pBanb6pnZYb22YNg12371hvE8fOPBAmD49v5jMzIqU8ybwQEl/lHRbNj5U0hFt3XBELCA9RXQH8CxwXURMbet6a8agQbDqqg3jc+emZiJWXDG/mMzMipRTB3Ap6SS9Ujb+b+DE9th4RPw9Ij4bEWtGxJntsc6aMns27LNPGv7Sl1wRbGY1pZw6gP4RcZ2k0yBduUtaWOG4Oocbbkh/t9kmFf3cfXe+8ZiZFSnnDuADScuTPaEjaXPgnYpG1dmcfjq88gpcfnnekZiZfaKcBHASqXJ2TUkPAJcD7gmsNXbeGYYNg5/9zC+GmVnNKOcpoCeA4aTHQY8E1o+IpysdWKcipbuAadPg2mvzjsbMDABFtNyem6QtgCEU1RlERNXLM4YNGxaTJk2q9mbbx6JFsNFGqd/gyZNTR/JmZlUg6fGIGFY6vZzHQK8AfglsBWyWfT61ImtBly7pLuCZZ+DGG/OOxsys5TsASc8CQ6OcW4UK69B3AAALF6a+Avr2hUmTUtGQmVmFLfEdADAF8NtL7aFrV/jOd+CJJ+COO/KOxszqXDkJoD/wjKQ7JN1U+FQ6sE7r4INh8GD4yU9SfYCZWU7KeRHsR5UOoq706AGnnALHHw/33ZeaijYzy0GTdwCS1gWIiHuBhyPi3sIH+LhaAXZKRxwBAwfCmZ2v9Qsz6ziaKwK6qmj4oZLvzqtALPWjd2846SS480549NG8ozGzOtVcAlATw42NW2sdfTQst5zvAswsN80lgGhiuLFxa61lloETToCbbkovhpmZVVlzCWAVSWdL+n3RcGH8U3332hI4/nhYemn46U/zjsTM6lBzTwGdXDRc+vZVB34bq4Z85jNwzDHwy1/CGWfAZz+bd0RmVkfKaguoVnT4N4EbM3s2DBkCo0fDJZfkHY2ZdUJteRPYKmngQPj61+GKK2DGjLyjMbM64gRQC04+ObULdNZZeUdiZnXECaAWDB4Mhx4KF18MM2fmHY2Z1YlymoP+rKR/SpqSjW8k6XuVD63OnHoqzJ8Pv/513pGYWZ0o5w7gIuA0YD5A1hvYqEoGVZfWWgtGjYLzz4c338w7GjOrA+UkgD4RUdpegTu2rYTTToMPPoCzz847EjOrA+UkgDckrUn29q+k/QAXVFfCBhvAvvumBPDuu3lHY2adXDkJ4FjgAmBdSa8CJwJHVTSqenb66fD223Ce29szs8oqJwHMiIgdgQHAuhGxVUS06YF1SV+RNFXSIknuX7jYppvCLrukyuAPP8w7GjPrxMpJANMlXQhsDrzfTtudAowA7mun9XUu3/sezJmTmogYPhxmzco7IjPrhMpJAOsAd5GKgqZLOkfSVm3ZaEQ8GxHPt2UdndpWW8E228AvfgETJ8K4cXlHZGadUIsJICI+iojrImIEsDHQF7i34pFlJI2VNEnSpDlz5lRrs/nq3Tt1F/nBB7BoUXo0VErTzczaSVlvAksaLuk84AmgF7B/GcvcJWlKI5+9WxNgRFwYEcMiYtiAAQNas2jHNW1aahyuS/bP07s3HHggTJ+eb1xm1qm02Cm8pOnAk8B1wMkR8UE5K84qjm1JDBoE/fpBoaXWjz6Cvn1hxRXzjcvMOpUWEwDwuYjwQ+nVNnt26jbyrbfg6qvh2WfzjsjMOpkmE4CkUyLiF8BPpE93ARwR31jSjUraF/g96dHSWyU9GRG7LOn6OqUbbkh/33kHJkyA99+HhQuha9d84zKzTqO5O4DCJefj7b3RiBgPjG/v9XZK/frBr36V6gAuugiO8jt4ZtY+WtUjmKQuwNJ5FQl1yh7ByhEBO+wA//oXPP88rLBC3hGZWQeyxD2CSbpKUl9JSwHPAM9LOrml5awdSXDuuemx0FNPzTsaM+skynkMdGh2xb8P8HdgVeDgikZln7beenDSSXDppenlMDOzNionAXSX1J2UAG6MiPlkLYNalX3/+6n3sKOPTp3HmJm1QTkJ4ALgJWAp4D5JqwF+LDQPSy0Fv/sdTJkCv/993tGYWQfXqkrgTxaSukVE1TuFqdtK4GIRsMceqamI556DlVfOOyIzq3FtqQQeKOmPkm7LxocCh1YgRiuHlDqMmT8fvvWtvKMxsw6snCKgS4E7gJWy8X+TOoWxvKy5Zuo+8tpr4a678o7GzDqochJA/4i4DlgEkBX9LKxoVNayU09NieDYY+Hjj/OOxsw6oHISwAeSlqehT+DNgXcqGpW1rFcvOOcc+Pe/U8cxZmatVE4COAm4CVhT0gPA5cDxFY3KyrPrrjByJPzkJ24q2sxarZwOYZ4AhgNbAEcC60fE05UOzMr0m9+kBuJOOCHvSMysgymrQxjgC8DngE2A0ZIOqVxI1iqDB8MPfwg335w+ZmZlavE9AElXAGuSOoUpVP5GW5qDXlJ+D6AJ8+fD5z8PH34IU6dCnz55R2RmNaSp9wDK6RBmGKk9IDf/UKu6d0+NxW23Hfz0p6lOwMysBeUUAU0B3Bdhrdt2WzjoIDjrrNRktJlZC5pMAJJulnQT0B94RtIdkm4qfKoXopXtrLPS46HHHdfQn7CZWROaKwLyw+UdzYorwplnwvHHw4UXwlVXpbeF3Zm8mTWirMbgJA0ENstGH42I1ysaVRNcCVyGhQths81SMdDcuXDkkXDeeXlHZWY5aktjcPsDjwJfAfYHHpG0X/uHaO1i6aVT15EffgiLFsH556cG5Hr3zjsyM6sx5VQCnw5sFhGHRsQhpHcCvl/ZsGyJTZsGBxwA3bLSvR49UofyflPYzEqUkwC6lBT5vFnmcpaHQYOgb9909S/BvHmpsTjXA5hZiXJO5LdnTwAdJukw4FbgtsqGZW0yezYcdVRqKnqppeCWW9I0M7MiLb4IFhEnSxoBbAUIuDAixlc8MltyN9zQMHz//bDlljBiBNx9N/TsmV9cZlZTmnsPYC1JWwJExA0RcVJEfBN4U9KaVYvQ2mbjjeHSS+HBB+GYY/x+gJl9orkioN8C7zUy/cPsuyUm6SxJz0l6WtJ4Scu2ZX3Wgv33h+99Dy65xJ3Jm9knmksAQxpr9jkiJgFD2rjdO4ENImIjUheTp7VxfdaSM86AvfeGk05yN5JmBjSfAHo1812bHiqPiH9kXUsCPAys0pb1WRm6dIErroB11013BC++mHdEZpaz5hLAY5K+XjpR0hHA4+0Yw+E081SRpLGSJkmaNGfOnHbcbB1aZhm46ab0eOhee8G77+YdkZnlqMmmILLmH8YD82g44Q8DegD7RsSsZlcs3UXjrYieHhE3ZvOcnq1zRDnNTbspiHYyYQLstBPsthv87W+pRzEz67Ra3R9ARMwGtpC0HbBBNvnWiLi7nA1GxI4tBHQosAewg/saqLLttoOzz4Zjj4Xvfz/1IWBmdaec9wAmABPac6OSdgVOBYZHxIftuW4r09FHw1NPwc9+BhtuCKNH5x2RmVVZXk06nAMsA9wp6UlJf8gpjvolpUdCt94aDj8cXLRmVnfK6RKy3UXEWnls10r06AF//WtqPnqffVIScJtBZnXDjbrVuwED4MYb4a23YN99U8NxZlYXnAAMPvc5uPxyePjh1Iic6+TN6oITgCUjR8IPf5jaDRo3DoYPh1nNPulrZh2cE4A1+MEPUquhP/pRakV03Li8IzKzCnICsAZLLdXQlHSEu5M0qxUzZ1bkrtwJwBoUupPsVdQM1FZbuTtJszxFwDe+ARMntvtduROANSh0JzlvXkPHMRMnwjXX5BuXWb2JgCefTH17d+kC11+funlt57tyJwBbXKE7yUcegbFjYaWV4JvfhG9/Ox2AZlYZETB5cmqeZd11U2dOEendnB490jx9+sCBB7bbXXkuL4JZDSvuTvKCC2DhQjjxRPjVr+DVV9NTQu5W0qz9TJ0K112XPs89l674t9suXXTtu29KCBdemIpm585Nd+nt9MKmE4A1r2vX1HDc4MFw6qnpDmH8eOjXL+/IzDqOmTNh1Ci49tp08n7uuYaT/tSp6aQ/fDiccEJ6Em+FFRqWLdyVjx2bEsHMme0WVpPNQdciNwedsyuvhDFjYL314LbbYOWV847IrGM45hj4wx9g2LB0FT95cirL33pr+OpX00m/gs2wtLo5aLNPOeggGDgwHaxf+hLcfjsMHZp3VGa1KSIV28yb1zDtscfS327dYMaMVMeWI1cCW+vstBPcdx/Mnw9bbpleGDOzZOHC9OTcSSfBGms0nPy7ZKfaXr1SJe7LL+d+8gcnAFsSG28MDz2U7gZ22im1KGpWr+bNS3fDhafmtt4azj033R1fdBEcemiar3A30I6VuG3lBGBLZsgQeOAB2GQT+MpXUt8CZp1Z8du477+fns0/4IDUou5uu8HVV8O226a/c+bArbfC176W+t4+6qiGxhZrqI0tVwJb23z0UepN7MYb4ZRT0huLBxzQ8LSDWWcxZgxcdhmsumo6iX/8MfTvD3vtlerFdthh8bfoa0hTlcBOANZ2CxfC8centxTXXhv+8x848kg477y8IzNbcgsXwqOPwjbbwIIFn/6+Rw/44INUoVvjmkoALgKytuvaFf70pzT8wgsVeWXdrCpmzUpX+aNGpaKdLbZIiaB/f+jePc1TeBt3xowOcfJvjhOAtY9CQ3KFV9YhPSU0bVp+MZmVKm1Vc/789CTbd7+b6rMGDQ99QAQAABEbSURBVILDDoN774W9907tYL3xBuy3X0oEFXgbN08dO31Z7Sg0JLdgQUoC8+alSuJjj01FQZ3gP4t1Aj/+cTrhjxqVrurvvDNV0nbtmi5YfvrTVKG70UYNj25CRd/GzZPrAKz9jBiREsHYsemtx/vvhxdfTLfMv/0tHHxwKhYyq6Z33knFOfPnf/q7rl1Tcww77NCpmzdxJbDl4/nn4fDD4cEHYffdUwNzq6ySd1TWmc2dm95T+ec/0+exx1LxTdeu6e3cRYtSUc6IEamRwzq4O3UlsOVjnXXSm8O/+x3ccw+sv356OaYDXXhYjSktx1+4MJ3kf/7z9GLicsvB9tuncQlOOw0mTEgXItDwQla/fnVx8m+O7wCseqZNSy/GTJiQbrkvughWXz3vqKyjOfrodCe55ZapaGfCBHj77fTdBhukY2uHHVKS6Nu3YbniIspCOX5x8+edmIuArDYsWpRO/CefnIZ//vPUUmIX34xaM6ZPT3eTTZXjX3FFuuofOLD6sXUANVUEJOnHkp6W9KSkf0jKv1Ukq44uXdJLYlOmpP6Gjz8+vT7/wgvp+wp1fm0dzCuvpJP64YenZkfWWCOd/Hv2TCd8SEU5BxyQ5h092if/JZDXZddZEbFRRHweuAX4QU5xWF5WXTX1KfCnP6W20TfaKFXInXFGRTq/thrSWJKfPTs9c3/kkelt8sGD4ZBDUhMjm26a2pqaMiU9o1/czLLL8dsk9yIgSacBq0bE0S3N6yKgTuq119J/+Mb6HO7VK7U3ZJ3HMcekMvwdd0wn+7vvhmefTd/17ZuSw3bbpU/p8/h1XI7fFjVXByDpTOAQ4B1gu4iY08R8Y4GxAKuuuuqmM2bMqF6QVj2vvZZaFX3wwTQupU5nrr463S1YxzZrVvp3bKwMv0sX+NnP0gl/4407fPMKtajqdQCS7pI0pZHP3gARcXpEDAb+DBzX1Hoi4sKIGBYRwwYMGFCpcC1vK63UcLXXvXu6zX/wwZQEfvrT9Dq+1Y6W6mpefjl1ITp2LKy7brpqnz8/ld8Xruh79kxl96++mlqS3Wwzn/yrrGIJICJ2jIgNGvncWDLrVcDISsVhHUjhdfvHHkuP+m2xRXqs7/TTUxHR2LHwzDN5R2mQmlQo1NVEpBZgL7kkldGvvnq62j/44PSW7dprw1lnwSOPwBFHpOV79UoJYdllXYafo1yKgCStHREvZMPHA8MjYr+WlnMdQJ2aOjW9SHbFFektz513hhNPhF128eOj1da7d/o3aEr//qn55OHD098NN2x4agdchp+TmqoDkPRXYB1gETADOCoiXm1pOSeAOvfGG+mkcc456cSx7rpwwgnpSnOppdK0UaPcGU17e+ONdPX+0EOplcyHH168ffzVVkt3bHvuCeut5/aealBNJYAl5QRgQHr87y9/gd/8Bh5/PL36P3Zsqkj+85/dGU25GkuY8+enx3Iffjid8B9+ODXoB+lK/nOfSz1hPfNMavV1/nzv7w7ACcA6n0JF8TbbNP4Iac+ezRdX1Ltjjkmttu68czqxP/QQTJrU8NjtiiumSvjNN0+fTTdNd1ouxulwnACs85o5M1Ue//3vixdNSOnEtuWW6bPFFqlysrEiio5afFRu3HPnpmftJ09Ob9cuXPjpeaTUp3PhpN/UvrIOp6kE4GeurOMbNCg9Rlpo5nfevNT09CabpE5pLr0Uzj03zbvyyg0JYcstU4Lo1m3xp1o6UnFGadyLFqV2cyZPXvzzwgsNJ/3u3WGZZeC999K0nj1h331TkVpHSn7WZr4DsM6huWKJBQvSSfCBBxo+L7/c/Ppa8wZyW+8eWrv8okWpk52PP2553jXWSE/iFH/WXju1wXThhQ29t7kcv1NzEZBZsZdfTvUH//gHjB8Pb7316Xn692+4uxg0aPHh4mnf/GZq2mBJT6KFphG+9rXUN+3MmekFq+K/xcOzZy9e1FUwcGB6EueLX0wn+vXXh6WXbnybLsevK04AZk05+uh0EuzWLT3Vss02qWnhmTPTk0WFv7NmNV52XkpKj6i25LnnWu4YR4IVVkh3BoMGNfwdNAhuvhnuustP41iLXAdg1pTGOvz+QSMN1C5alJ6JLySFZ5+Fyy5Lj0QuWJAek1x55VSv0KtXy9tdc0146qm0voULU9n8JpvAccfB0KHpZL/CCk03jzBhQkpenayjcqse3wGYtUXh7mFJy9LburxZGWqqQxizTqNw9/Dww+lvazuyaevyZm3gOwAzs07OdwBmZrYYJwAzszrlBGBmVqecAMzM6pQTgJlZnXICMDOrUx3qMVBJc0g9iC2J/kAt9izuuFrHcbWO42qdWo0L2hbbahExoHRih0oAbSFpUmPPwebNcbWO42odx9U6tRoXVCY2FwGZmdUpJwAzszpVTwngwrwDaILjah3H1TqOq3VqNS6oQGx1UwdgZmaLq6c7ADMzK+IEYGZWpzp8ApDUS9Kjkp6SNFXSGY3M01PStZJelPSIpCFF352WTX9e0i5VjuskSc9IelrSPyWtVvTdQklPZp+bqhzXYZLmFG3/a0XfHSrphexzaJXj+k1RTP+W9HbRdxXZX0Xr7yrpX5JuaeS7qh9fZcZV9eOrzLiqfnyVGVcux5eklyRNztb9qfbulZydHUdPS9qk6Lu27a+I6NAfQMDS2XB34BFg85J5jgH+kA2PAq7NhocCTwE9gdWB/wBdqxjXdkCfbPjoQlzZ+Ps57q/DgHMaWfYzwLTs73LZ8HLViqtk/uOBSyq9v4rWfxJwFXBLI99V/fgqM66qH19lxlX146ucuPI6voCXgP7NfL87cFv2f2Rz4JH22l8d/g4gkvez0e7Zp7Rme2/gsmz4emAHScqmXxMRH0fEdOBF4AvViisiJkTEh9now8Aq7bHttsbVjF2AOyPifxHxFnAnsGtOcY0Grm6PbbdE0irAl4GLm5il6sdXOXHlcXyVE1czKnZ8LUFcVTu+yrA3cHn2f+RhYFlJg2iH/dXhEwB8clv3JPA6aYc8UjLLysDLABGxAHgHWL54euaVbFq14ip2BCnLF/SSNEnSw5L2aa+YWhHXyOx283pJg7NpNbG/sqKM1YG7iyZXbH8BvwVOARY18X0ux1cZcRWr2vFVZlxVP77KjCuP4yuAf0h6XNLYRr5var+0eX91igQQEQsj4vOkK5wvSNqgZBY1tlgz06sVVwpOOggYBpxVNHnVSK99HwD8VtKaVYzrZmBIRGwE3EXD1W1N7C9SMcv1EbGwaFpF9pekPYDXI+Lx5mZrZFpFj68y4yrMW7Xjq8y4qn58tWZ/UcXjK7NlRGwC7AYcK2mbku8rdnx1igRQEBFvA/fw6dugV4DBAJK6Af2A/xVPz6wCvFbFuJC0I3A6sFdEfFy0zGvZ32nZshtXK66IeLMolouATbPh3PdXZhQlt+cV3F9bAntJegm4Bthe0pUl8+RxfJUTVx7HV4tx5XR8lbW/MtU8vorX/Townk8XEza1X9q+vypVsVGtDzAAWDYb7g3cD+xRMs+xLF5Jd102vD6LV9JNo/0qgcuJa2NSxeDaJdOXA3pmw/2BF4ChVYxrUNHwvsDD0VDpND2Lb7ls+DPViiv7bh1SpZmqsb9Ktr0tjVdqVv34KjOuqh9fZcZV9eOrnLjyOL6ApYBlioYfBHYtmefLLF4J/Gh77a9udHyDgMskdSXd0VwXEbdIGgdMioibgD8CV0h6kXRlNgogIqZKug54BlgAHBuL3/ZVOq6zgKWBv6Q6Q/4bEXsB6wEXSFqULfvziHiminF9Q9JepH3yP9JTG0TE/yT9GHgsW9e4iPhfFeOCVDl3TWT/AzKV3F+NqoHjq5y48ji+yokrj+OrnLig+sfXQGB89u/TDbgqIm6XdBRARPwB+DvpSaAXgQ+BMdl3bd5fbgrCzKxOdao6ADMzK58TgJlZnXICMDOrU04AZmZ1ygnAzKxOOQFYq0l6v+W5Ppl3W0lbFI3vI2lo0fi47GWlmiXpu21cfrHf3Mj3J0o6pGi8m6Q3JP2sLdttYXt9KrTuAZJur8S6rf05AVilbQtsUTS+D6mVTAAi4gcRcVe1g2qlNiUASn5zsezN4cNJLVQW7Aw8D+yfNSrX3k4EGk0A2XsYSywi5gAzJW3ZlvVYdTgBWLuQtKdSW/j/knSXpIFK7eIfBXxTqa3z4cBewFnZ+JqSLpW0X7aOzSQ9qNQnwKOSlskaiDtL0mNZ42FHNrH9Q7Lvn5J0RTZtNaV28Avt4a+aTb9UqX31ByVNK9r+IEn3ZbFNkbS1pJ8DvbNpf87m+5tSw11TVdR4l6T3JZ2ZxfBwtg+2KP3NJaFvDzwRqRG5gtHA74D/kt78LKz/JUlnSHpCqf34dbPpAyTdmU2/QNIMSf0lLSXp1iyeKZK+KukbwErABEkTiuIeJ+kR4EuSdsj+HSdLukRSz6Lt/1TSQ0oNo20i6Q5J/ym8uJT5G3BgOceN5ay9X7P2p/N/aKRtdNKr6IUXC78G/Cob/hHw7aL5LgX2Kx0HepCaStgsm96X9GbkWOB72bSewCRg9ZJtr0+6Yu6fjX8m+3szcGg2fDjwt6Jt/oV0ATQUeDGb/i3g9Gy4Kw2v6L9fsr3C+nsDU4Dls/EA9syGf1EU92K/uWRdZwDHF433JrXn0if77WcXffdSYV5SHwQXZ8PnAKdlw7tmcfQHRgIXFS3fr2g9/YumB7B/NtyL1MLkZ7Pxy4ETi5Y7Ohv+DfA0sAypGY/Xi9a3MjA57+PUn5Y/vgOw9rIKcIekycDJpJNya6wDzIyIxwAi4t1IV8U7A4coNRP9CKmZ5bVLlt2e1HrjG9myhdfhv0RD0coVwFZFy/wtIhZFeqV/YDbtMWCMpB8BG0bEe03E+g1JT5Ha2B9cFM88oNDT1OPAkDJ+9yBgTtH4HsCESO34/xXYt6RY5oZG1r8VqYEzIuJ24K1s+mRgR0n/J2nriHiniRgWZtuC9O8wPSL+nY1fBhS3TlloMmEyqWOS9yIV+8yVtGz23eukuwyrcU4A1l5+T+rlaUPgSNKVZGuIxpuyFemq9/PZZ/WI+EeZy5YqnufjomEBRMR9pJPdq6S2fQ6hhKRtgR2BL0XE54B/0fBb50dEYRsLoay2tj5i8X01mnTSfol0kl+e1LNXadzF62+0niA7iW9KOln/TNIPmohhbjS0UdRSnUNh+4tYfB8uKoqnF+l3WY1zArD20o904gQ4tGj6e6RigqbGC54DVpK0GUBW/t8NuAM4WlL3bPpnJS1Vsuw/SRWmy2fzfCab/iBZw2ykMumJzf0ApY5AXo+Ii0gNvBX6Xp1f2H72O9+KiA+zMvjNG1lVqaZ+M8CzwFrZ9vuSruZXjYghETGE1NLo6BbWPxHYP1vHzqTiOCStBHwYEVcCvyz6Pc3F8xwwRNJa2fjBwL0tbL/UZ0lFY1bjnABsSfSR9ErR5yRSWf9fJN0PvFE0782kYownJW1NKqo4Oatk/KRCNCLmAV8Ffp8Vr9xJupK8mNSa5hOSpgAXUHJlHRFTgTOBe7Nlf5199Q1Skc7TpBPZCS38rm2BJyX9i1R+/rts+oXA01kl8O1At2ydPyYVA7Wk0d+cuY2GIpYRwN1R1G4/cCOpHfuezaz/DGBnSU+QOhWZSTrJbwg8mhWfnQ78pOj33FaoBC4WEXNJrU3+JSvOWwT8oYzfWGw74NZWLmM5cGugZjmTNB44JSJeWMLlewILI2KBpC8B50fqWS0Xku4D9o7UT63VsM7QH4BZR/cdUmXwEiUAYFXgOkldSBXRX2+vwFpL0gDg1z75dwy+AzAzq1OuAzAzq1NOAGZmdcoJwMysTjkBmJnVKScAM7M69f9N2/28PQmqQgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import subprocess \n",
    "import shlex\n",
    "import glob, os\n",
    "import time\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "j = 0\n",
    "x = []\n",
    "y = []\n",
    "target = \"%%\"\n",
    "for i in list(np.arange(3.00, 5.05, 0.10)):\n",
    "    command_line = f'\\\"C:/Program Files/LAMMPS 64-bit 24Jan2020/bin/lmp_serial\\\" \\\n",
    "-in calc_fcc_ver2.in -var latconst {i:0.02f} -log log_{j:03d}.lammps'\n",
    "    x.append(round(i,3))\n",
    "    file = f\"log_{j:03d}.lammps\"\n",
    "    j += 1\n",
    "    print(command_line)\n",
    "    args = shlex.split(command_line)\n",
    "    p = subprocess.Popen(args) # Success!\n",
    "    p.wait()\n",
    "    with open(file) as f:\n",
    "        for line in f:\n",
    "            if target in line and \"print\" not in line:\n",
    "                exec(line.replace(\"%% \",\"\").replace(\";\",\"\"))\n",
    "                y.append(ecoh)\n",
    "\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "                \n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import *\n",
    "figure()\n",
    "plot(x, y, 'r*-')\n",
    "xlabel('Lattice constant (Angstrom)')\n",
    "ylabel('Cohesive Energy (eV)')\n",
    "title('aluminum cohesive energy evolution')\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm.  It's quite a bit slower when you have to wait within a loop for each subprocess to finish with the `.wait()` command.  I had to put that in the script so it didn't start a subprocess and while it is still running, the Python script would try to open up the log file.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## LAMMPS Input Script Explained\n",
    "\n",
    "In this section, we will break down what LAMMPS is doing. The easy way to do this on your own is to consult the LAMMPS manual for each command or go to the Internet LAMMPS manual, *i.e.*, at https://lammps.sandia.gov\n",
    "\n",
    "So, first, what is different in this example?  Many of the LAMMPS commands have been explained in Tutorial 1.  Here is what I think is important in this tutorial.\n",
    "\n",
    "### `Variable` command\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "lattice fcc \\${latconst} <br>\n",
    "region box block 0 1 0 1 0 1 units lattice <br>\n",
    "create_box 1 box <br><br>\n",
    "lattice fcc \\${latconst} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 <br>\n",
    "create_atoms 1 box <br>\n",
    "replicate 1 1 1\n",
    "</div>\n",
    "\n",
    "OK. Do you remember those variables from the first LAMMPS input script? They can be fed in from the command line as well.  This is a great example of passing in the constant `latconst` from the command line and then using it within other LAMMPS commands.  If you combed through the script, you won't find any `variable` command that sets `latconst` to anything, i.e., because it is meant to be set from the command line.  It's usually a good idea to add a comment about the command line syntax or what variables need to be set outside of the script.\n",
    "\n",
    "### Print `%%` in log file.  Print and retrieve!\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "print \"%% ecoh = ${ecoh};\"\n",
    "</div>\n",
    "\n",
    "You are likely to miss this line if I didn't specifically call it out in the tutorial. This was just an easy way to flag certain lines within the log file, which is filled with all sorts of stuff.  There are a number of different iterations of this and you'll see in the tutorials how this is used in different ways to extract data from the LAMMPS simulations.  I like having all the information in a file, and then being able to extract it when I need to.\n",
    "\n",
    "But Mark, why didn't you just compile LAMMPS with Python and seemlessly pass information back-and-forth?  Well, because this is really easy to set up. Plus, my days of compiling Fortran versions of codes like LAMMPS on various supercomputers has left me permanently scarred.  It's gotten better, I know, but if I can not compile LAMMPS on a **WINDOWS** machine to people how to do a few things... I'll take it! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## FAQs \n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-danger\">\n",
    "Q1: So, what does the `run 0` command actually do? Is it minimization? Is it molecular dynamics?\n",
    "</div>\n",
    "\n",
    "A: The `run 0` command is not minimization or molecular dynamics. The atoms are not displaced at all by this command. This is merely a way for LAMMPS to know that it needs to calculate additionally information for each atom or for the entire system. For instance, before the `run 0` command, LAMMPS may not know the forces on each atom or the energy for the entire system. By executing the `run 0` command, LAMMPS will calculate all the forces and attributes necessary to run molecular dynamics, but will not actually perform the timestep at the end, moving the atoms. Now, the user can print these attributes to screen or to a dumpfile, etc. \n",
    "\n",
    "In some cases, the user will need to tell LAMMPS what values they want LAMMPS to calculate prior to a `run 0` command. For example, if you want the stresses on every atom, you will need to define a compute or a variable command prior to `run 0` and then let LAMMPS know that you need this value by inserting it in to the `thermo_style` command (like in Tutorial 3!). For more information, see the LAMMPS documentation on the run command. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Tutorial Links\n",
    "\n",
    "* [Click here to open Tutorial 1](LAMMPS-Tutorials-01.ipynb)\n",
    "* [Click here to open the next tutorial](LAMMPS-Tutorials-03.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## References \n",
    "\n",
    "1. S. Plimpton, \"Fast Parallel Algorithms for Short-Range Molecular Dynamics,\" J. Comp. Phys., 117, 1-19 (1995)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}