LAMMPS-Tutorials-08.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LAMMPS Tutorials 08. Simulate a single polymer chain!\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",
    "Thanks to Dmitry Zhuk for writing the first instantiation of this tutorial!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Abstract:\n",
    "<a id=\"Sec1\"></a>\n",
    "\n",
    "In this tutorial, a molecular dynamics simulation in LAMMPS is used to show what happens to a polymer chain at a certain temperature after some time. The chain's movement is caused by a molecular forces between atoms in the chain and by temperature of the chain. These factors are modeled in LAMMPS to show the behavior of this polymer chain. The chain was previously prepared in MATLAB. It contains 100 atoms with bound between neighbor atoms. During the simulation, the chain was equilibrated and fluctuates at temperature. After that the chain was minimized to find it's minimal energy condition. The tutorial explains basic LAMMPS commands and shows how to make a movie using AtomEye and ImageJ. \n",
    "\n",
    "<table><tr><td><img src='https://icme.hpc.msstate.edu/mediawiki/images/8/89/Equ_plus_Min.gif' width=\"300\"><center><br>Equilibration process followed by minimization</center></td></tr></table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Complete the First Tutorials\n",
    "\n",
    "If you have not done so already, you may want to complete the first tutorials available [here](https://github.com/mrkllntschpp/lammps-tutorials). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Step 1: Download Polymer Configuration file\n",
    "\n",
    "The polymer chain of 100 atoms was initially written in MATLAB. The atom's Z coordinate does not vary much, all of them are within 2 Å. The distance between atoms is about 1.5 Å. Essentially, the chain goes from left upper to right lower corner of the box. \n",
    "\n",
    "The data file is shown below and available for download here: https://icme.hpc.msstate.edu/mediawiki/images/e/e1/PE_cl100.txt\n",
    "\n",
    "Just in case you are interested in what it looks like:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "# Model for PE\n",
    "\n",
    "     100     atoms\n",
    "      99     bonds\n",
    "      98     angles\n",
    "      97     dihedrals\n",
    "\n",
    "         1     atom types\n",
    "         1     bond types\n",
    "         1     angle types\n",
    "         1     dihedral types\n",
    "\n",
    "    0.0000   158.5000 xlo xhi\n",
    "    0.0000   158.5000 ylo yhi\n",
    "    0.0000   100.0000 zlo zhi\n",
    "\n",
    "Masses\n",
    "\n",
    "         1          14.02\n",
    "\n",
    "Atoms\n",
    "\n",
    "         1         1         1    5.6240    5.3279   51.6059\n",
    "         2         1         1    7.4995    7.4810   50.2541\n",
    "         3         1         1    8.2322    8.0236   51.2149\n",
    "         4         1         1    9.6108    9.9075   51.7682\n",
    "         5         1         1   11.5481   11.3690   50.4167\n",
    "         6         1         1   12.9409   13.4562   50.2481\n",
    "         7         1         1   14.4708   14.8569   50.0868\n",
    "         8         1         1   16.1916   16.4790   50.5665\n",
    "         9         1         1   17.1338   17.6853   51.8189\n",
    "        10         1         1   19.1109   19.4000   50.3869\n",
    "\n",
    "...\n",
    "\n",
    "Bonds\n",
    "\n",
    "         1         1         1         2\n",
    "         2         1         2         3\n",
    "         3         1         3         4\n",
    "         4         1         4         5\n",
    "         5         1         5         6\n",
    "         6         1         6         7\n",
    "         7         1         7         8\n",
    "         8         1         8         9\n",
    "         9         1         9        10\n",
    "        10         1        10        11\n",
    " \n",
    "...\n",
    "\n",
    "Angles\n",
    "\n",
    "         1         1         1         2         3\n",
    "         2         1         2         3         4\n",
    "         3         1         3         4         5\n",
    "         4         1         4         5         6\n",
    "         5         1         5         6         7\n",
    "         6         1         6         7         8\n",
    "         7         1         7         8         9\n",
    "         8         1         8         9        10\n",
    "\n",
    "...\n",
    "\n",
    "Dihedrals\n",
    "\n",
    "         1         1         1         2         3         4\n",
    "         2         1         2         3         4         5\n",
    "         3         1         3         4         5         6\n",
    "         4         1         4         5         6         7\n",
    "         5         1         5         6         7         8\n",
    "         6         1         6         7         8         9\n",
    "         7         1         7         8         9        10\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the first section defines the numbers of atoms (100), bonds(99), angles (98), and dihedrals (97). Then the number of atom types, bond types, angle types, and dihedral types (each 1 type for simplicity).   Then comes the simulation cell sizes in the x, y, and z directions (denoted by a `lo`(w) and 'hi'(gh) value in each dimension). \n",
    "\n",
    "Then comes the explicit description for the atoms and the connections for the bonds, angles, and dihedrals. In the Atom section, the first column is the atom number, then comes the atom type and a column with information not used in this tutorial (read the description in LAMMPS help). The last three columns are each atom's x, y and z coordinates, correspondingly. \n",
    "\n",
    "The Bond section is intuitive.  Each bond has a number referring to it (column 1 - there are 99 of them), the bond type is used to relate it to a potential function (column 2 - only type 1), and they each connect two atoms, which were labeled in the section just above that (columns 3 and 4).  \n",
    "\n",
    "The Angle section lists the angles between atoms. The first column is the angle number, the second column is the angle ID, then comes the atom numbers between which the angle is defined. Example of the distribution of the atoms can be seen on the picture below.  Note that this picture is defined by the following text above:\n",
    "\n",
    "`         7         1         7         8         9` (angle ID 7, angle type 1, between atoms 7, 8, and 9)\n",
    "\n",
    "The Dihedral angle section is just a little bit more complicated. To define a dihedral angle, four atoms are needed. The syntax is similar to the angles section. The only difference is that one more atom is needed to define a dihedral. Basically, the dihedral angle is the angle between the planes formed by 2 groups of 3 neighbor atoms. An example of the dihedral angle is shown below. Note that this picture is defined by the following text above:\n",
    "\n",
    "`         4         1         4         5         6         7` (dihedral ID 4, dihedral type 1, between atoms 4, 5, 6, and 7)\n",
    "\n",
    "<table><tr><td><img src='https://icme.hpc.msstate.edu/mediawiki/images/f/ff/Atom_angle.jpg' width=\"300\"><center><br>Angles between atoms defined in LAMMPS</center></td><td><img src='https://icme.hpc.msstate.edu/mediawiki/images/a/a8/Dihedral.jpg' width=\"300\"><center><br>Dihedral angles between atoms defined in LAMMPS</center></td></tr></table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Create the LAMMPS input script \n",
    "\n",
    "This input script was run using the January 2020 version of LAMMPS. Changes in some commands in more recent versions may require revision of the input script. This script runs the simulation with the previously generated polymer chain file, which is fed to LAMMPS through the `read_data` command.  To run this script, store it in `in.deform_polymer_chain.txt`and use \n",
    "\n",
    "`lmp_exe < in.deform_polymer_chain.txt` \n",
    "\n",
    "where `lmp_exe` refers to the LAMMPS executable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing in.deform_polymer_chain.txt\n"
     ]
    }
   ],
   "source": [
    "%%writefile in.deform_polymer_chain.txt\n",
    "######################################\n",
    "# LAMMPS INPUT SCRIPT\n",
    "# Polymer Chain Tutorial\n",
    "# Mark Tschopp\n",
    "# The methodology outlined here follows that from Hossain, Tschopp, et al. 2010, Polymer.\n",
    "# The following script requires a LAMMPS data file containing the coordinates and \n",
    "# appropriate bond/angle/dihedral lists for each united atom.                                      \n",
    "# syntax: lmp_exe -in in.deform_polymer_chain.txt\n",
    "\n",
    "######################################\n",
    "# VARIABLES\n",
    "variable fname index PE_cl100.txt\n",
    "variable simname index PE_cl100\n",
    "\n",
    "######################################\n",
    "# INITIALIZATION\n",
    "units real\n",
    "boundary f f f\n",
    "atom_style molecular\n",
    "log log.${simname}.txt\n",
    "read_data ${fname}\n",
    "\n",
    "######################################\n",
    "# DEFINE INTERATOMIC POTENTIAL\n",
    "# Dreiding potential\n",
    "neighbor 0.4 bin\n",
    "neigh_modify every 10 one 10000\n",
    "bond_style      harmonic\n",
    "bond_coeff 1 350 1.53\n",
    "angle_style     harmonic\n",
    "angle_coeff 1 60 109.5\n",
    "dihedral_style multi/harmonic\n",
    "dihedral_coeff 1 1.73 -4.49 0.776 6.99 0.0\n",
    "pair_style lj/cut 10.5\n",
    "pair_coeff 1 1 0.112 4.01 10.5\n",
    "\n",
    "######################################\n",
    "# DEFINE COMPUTES \n",
    "compute csym all centro/atom fcc\n",
    "compute peratom all pe/atom \n",
    "compute eng all pe/atom \n",
    "compute eatoms all reduce sum c_eng \n",
    "\n",
    "#####################################################\n",
    "# EQUILIBRATION\n",
    "# Langevin dynamics at 500 K\n",
    "\n",
    "velocity all create 500.0 1231\n",
    "fix 1 all nve/limit 0.05\n",
    "fix 2 all langevin 500.0 500.0 10.0 904297\n",
    "thermo_style custom step temp \n",
    "thermo          100\n",
    "timestep 1\n",
    "run 10000\n",
    "unfix 1\n",
    "unfix 2\n",
    "write_restart restart.${simname}.dreiding1\n",
    "\n",
    "#####################################################\n",
    "# MINIMIZATION\n",
    "\n",
    "dump 1 all cfg 6 dump.PE_cl100_*.cfg mass type xs ys zs c_csym c_peratom fx fy fz\n",
    "reset_timestep 0 \n",
    "fix 1 all nvt temp 500.0 500.0 100.0\n",
    "thermo 20 \n",
    "thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms \n",
    "min_style cg\n",
    "minimize 1e-15 1e-15 1000 1000 \n",
    "\n",
    "#####################################################\n",
    "print \"All done\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## United atom potential\n",
    "\n",
    "Notice that we didn't have to download a potential file for this simulation.  All of the potential parameters necessary are programmed into those few lines within the LAMMPS script (see below).  The concept of a united atom is that each individual \"united\" atom actually represents multiple atoms, in this case the carbon and two hydrogen atoms from a polyethylene chain (i.e., cl100 represents C$_{100}$H$_{202}$). \n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "neighbor 0.4 bin <br>\n",
    "neigh_modify every 10 one 10000 <br>\n",
    "bond_style      harmonic <br>\n",
    "bond_coeff 1 350 1.53 <br>\n",
    "angle_style     harmonic <br>\n",
    "angle_coeff 1 60 109.5 <br>\n",
    "dihedral_style multi/harmonic <br>\n",
    "dihedral_coeff 1 1.73 -4.49 0.776 6.99 0.0 <br>\n",
    "pair_style lj/cut 10.5 <br>\n",
    "pair_coeff 1 1 0.112 4.01 10.5\n",
    "</div>\n",
    "\n",
    "where the first number after the `pair_coeff`, `bond_coeff`, `angle_coeff`, and `dihedral_coeff` are the \"IDs\" from column 2 of the data file in Step 1, thereby identifying the potential for each of these components."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Run the LAMMPS Simulation\n",
    "\n",
    "Now run the simulation as we have done before using the syntax \n",
    "\n",
    "`lmp_exe < in.deform_polymer_chain.txt`.  \n",
    "\n",
    "On my computer, the 24Jan2020 LAMMPS executable is stored in the `C:\\Program Files\\LAMMPS 64-bit 24Jan2020\\bin\\` folder and is named `lmp_serial.exe`.  \n",
    "\n",
    "The `log.PE_cl100.txt` file should look like the output below. Whoa! 10,000 timesteps...  that must take a while. Nope. The simulation took less than a second on my computer (100,000 timesteps only took about 3 seconds).  Remember that this is only 100 atoms with a very simple potential."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LAMMPS (3 Mar 2020)\n",
      "Reading data file ...\n",
      "  orthogonal box = (0 0 0) to (158.5 158.5 100)\n",
      "  1 by 1 by 1 MPI processor grid\n",
      "  reading atoms ...\n",
      "  100 atoms\n",
      "  scanning bonds ...\n",
      "  1 = max bonds/atom\n",
      "  scanning angles ...\n",
      "  1 = max angles/atom\n",
      "  scanning dihedrals ...\n",
      "  1 = max dihedrals/atom\n",
      "  reading bonds ...\n",
      "  99 bonds\n",
      "  reading angles ...\n",
      "  98 angles\n",
      "  reading dihedrals ...\n",
      "  97 dihedrals\n",
      "Finding 1-2 1-3 1-4 neighbors ...\n",
      "  special bond factors lj:   0          0          0         \n",
      "  special bond factors coul: 0          0          0         \n",
      "  2 = max # of 1-2 neighbors\n",
      "  2 = max # of 1-3 neighbors\n",
      "  4 = max # of 1-4 neighbors\n",
      "  6 = max # of special neighbors\n",
      "  special bonds CPU = 0.000180006 secs\n",
      "  read_data CPU = 0.00569391 secs\n",
      "Neighbor list info ...\n",
      "  update every 10 steps, delay 10 steps, check yes\n",
      "  max neighbors/atom: 10000, page size: 100000\n",
      "  master list distance cutoff = 10.9\n",
      "  ghost atom cutoff = 10.9\n",
      "  binsize = 5.45, bins = 30 30 19\n",
      "  2 neighbor lists, perpetual/occasional/extra = 1 1 0\n",
      "  (1) pair lj/cut, perpetual\n",
      "      attributes: half, newton on\n",
      "      pair build: half/bin/newton\n",
      "      stencil: half/bin/3d/newton\n",
      "      bin: standard\n",
      "  (2) compute centro/atom, occasional\n",
      "      attributes: full, newton on\n",
      "      pair build: full/bin\n",
      "      stencil: full/bin/3d\n",
      "      bin: standard\n",
      "Setting up Verlet run ...\n",
      "  Unit style    : real\n",
      "  Current step  : 0\n",
      "  Time step     : 1\n",
      "Per MPI rank memory allocation (min/avg/max) = 8.099 | 8.099 | 8.099 Mbytes\n",
      "Step Temp \n",
      "       0          500 \n",
      "     100    679.10043 \n",
      "     200    611.54642 \n",
      "     300    641.48308 \n",
      "     400    550.28444 \n",
      "     500    618.64162 \n",
      "     600    574.90254 \n",
      "     700     579.6052 \n",
      "     800    560.73043 \n",
      "     900    589.98744 \n",
      "    1000    429.95347 \n",
      "    1100    540.84598 \n",
      "    1200    571.10339 \n",
      "    1300    515.45725 \n",
      "    1400    590.32144 \n",
      "    1500    545.06421 \n",
      "    1600    525.49468 \n",
      "    1700    641.66642 \n",
      "    1800    629.51826 \n",
      "    1900    529.13537 \n",
      "    2000    587.15636 \n",
      "    2100    579.09598 \n",
      "    2200    516.16292 \n",
      "    2300    568.44108 \n",
      "    2400    510.74808 \n",
      "    2500    568.28071 \n",
      "    2600    566.84316 \n",
      "    2700    424.97451 \n",
      "    2800    509.65558 \n",
      "    2900     660.1342 \n",
      "    3000     518.6186 \n",
      "    3100    581.76475 \n",
      "    3200    585.22442 \n",
      "    3300    480.70202 \n",
      "    3400    567.63538 \n",
      "    3500    596.75494 \n",
      "    3600    541.14352 \n",
      "    3700    478.60105 \n",
      "    3800    531.73688 \n",
      "    3900    523.91918 \n",
      "    4000    505.85932 \n",
      "    4100    618.09178 \n",
      "    4200    517.70664 \n",
      "    4300     615.2756 \n",
      "    4400     508.2206 \n",
      "    4500    526.21919 \n",
      "    4600    511.86477 \n",
      "    4700    566.68262 \n",
      "    4800    521.48111 \n",
      "    4900     558.9255 \n",
      "    5000    589.97877 \n",
      "    5100    526.12526 \n",
      "    5200    518.28703 \n",
      "    5300     533.3145 \n",
      "    5400    476.83486 \n",
      "    5500    548.01645 \n",
      "    5600    498.73744 \n",
      "    5700    595.52205 \n",
      "    5800     575.9437 \n",
      "    5900    545.01307 \n",
      "    6000    548.63026 \n",
      "    6100    521.95715 \n",
      "    6200    477.77285 \n",
      "    6300    438.92221 \n",
      "    6400    473.06068 \n",
      "    6500    515.93384 \n",
      "    6600    488.60385 \n",
      "    6700    445.23567 \n",
      "    6800    552.62501 \n",
      "    6900    566.43421 \n",
      "    7000    509.95734 \n",
      "    7100    531.26926 \n",
      "    7200    577.62881 \n",
      "    7300    479.27758 \n",
      "    7400    506.73024 \n",
      "    7500    537.55261 \n",
      "    7600    517.49232 \n",
      "    7700     506.6552 \n",
      "    7800    509.91288 \n",
      "    7900    550.03833 \n",
      "    8000    510.43975 \n",
      "    8100    490.94551 \n",
      "    8200    549.97739 \n",
      "    8300    477.72934 \n",
      "    8400    444.64217 \n",
      "    8500    527.10749 \n",
      "    8600    548.42555 \n",
      "    8700    489.16576 \n",
      "    8800    512.61258 \n",
      "    8900    526.27395 \n",
      "    9000    469.68818 \n",
      "    9100     493.7436 \n",
      "    9200    519.58576 \n",
      "    9300     608.3647 \n",
      "    9400     515.1291 \n",
      "    9500    538.39326 \n",
      "    9600    512.47775 \n",
      "    9700      557.559 \n",
      "    9800    516.10086 \n",
      "    9900    529.32883 \n",
      "   10000    490.20334 \n",
      "Loop time of 0.382349 on 1 procs for 10000 steps with 100 atoms\n",
      "\n",
      "Performance: 2259.716 ns/day, 0.011 hours/ns, 26154.115 timesteps/s\n",
      "98.5% CPU use with 1 MPI tasks x no OpenMP threads\n",
      "\n",
      "MPI task timing breakdown:\n",
      "Section |  min time  |  avg time  |  max time  |%varavg| %total\n",
      "---------------------------------------------------------------\n",
      "Pair    | 0.049672   | 0.049672   | 0.049672   |   0.0 | 12.99\n",
      "Bond    | 0.21129    | 0.21129    | 0.21129    |   0.0 | 55.26\n",
      "Neigh   | 0.037267   | 0.037267   | 0.037267   |   0.0 |  9.75\n",
      "Comm    | 0.0017045  | 0.0017045  | 0.0017045  |   0.0 |  0.45\n",
      "Output  | 0.00077415 | 0.00077415 | 0.00077415 |   0.0 |  0.20\n",
      "Modify  | 0.078396   | 0.078396   | 0.078396   |   0.0 | 20.50\n",
      "Other   |            | 0.003247   |            |       |  0.85\n",
      "\n",
      "Nlocal:    100 ave 100 max 100 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "Nghost:    0 ave 0 max 0 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "Neighs:    399 ave 399 max 399 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "FullNghs:  0 ave 0 max 0 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "\n",
      "Total # of neighbors = 399\n",
      "Ave neighs/atom = 3.99\n",
      "Ave special neighs/atom = 5.88\n",
      "Neighbor list builds = 518\n",
      "Dangerous builds = 40\n",
      "System init for write_restart ...\n",
      "WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (src/min.cpp:190)\n",
      "Setting up cg style minimization ...\n",
      "  Unit style    : real\n",
      "  Current step  : 0\n",
      "Per MPI rank memory allocation (min/avg/max) = 10.73 | 10.73 | 10.73 Mbytes\n",
      "Step PotEng Lx Ly Lz Press Pxx Pyy Pzz c_eatoms \n",
      "       0     939.0478        158.5        158.5          100   -103.59827   -153.62204   -155.64402   -1.5287482     939.0478 \n",
      "      20    690.94736        158.5        158.5          100   -87.999784   -133.21524   -133.84363    3.0595213    690.94736 \n",
      "      40    526.91489        158.5        158.5          100   -74.970077   -114.48541    -112.5078    2.0829765    526.91489 \n",
      "      60    356.03461        158.5        158.5          100   -59.462004   -91.044182   -89.488953    2.1471239    356.03461 \n",
      "      80    207.25645        158.5        158.5          100    -43.58014   -66.605544   -66.129841     1.994964    207.25645 \n",
      "     100    109.13827        158.5        158.5          100   -30.415924   -46.657818   -46.945286    2.3553312    109.13827 \n",
      "     120    68.819378        158.5        158.5          100   -23.379241   -36.472647   -36.695413    3.0303379    68.819378 \n",
      "     140    38.506783        158.5        158.5          100   -17.538345   -27.573219   -27.784913    2.7430971    38.506783 \n",
      "     160     17.27959        158.5        158.5          100   -11.671201   -18.719908   -18.842429    2.5487348     17.27959 \n",
      "     180    2.8675512        158.5        158.5          100   -7.1491533   -12.263057   -11.770757    2.5863542    2.8675512 \n",
      "     200   -2.1848489        158.5        158.5          100   -4.9654302   -8.9759364   -8.6858968    2.7655427   -2.1848489 \n",
      "     220   -7.0009826        158.5        158.5          100   -2.1080084   -4.8272342   -4.3216139    2.8248228   -7.0009826 \n",
      "     240   -8.9832643        158.5        158.5          100  -0.32751718   -2.2244082   -1.6617359    2.9035926   -8.9832643 \n",
      "     260   -9.7033977        158.5        158.5          100   0.35491076   -1.2217998   -0.5872161    2.8737482   -9.7033977 \n",
      "     280   -10.128971        158.5        158.5          100   0.72642708   -0.6765529  -0.02354524    2.8793794   -10.128971 \n",
      "     300   -10.449359        158.5        158.5          100    1.1484347 -0.071303239    0.6107478    2.9058594   -10.449359 \n",
      "     320   -10.581447        158.5        158.5          100     1.292233   0.16741435   0.81522263     2.894062   -10.581447 \n",
      "     340   -10.740691        158.5        158.5          100    1.2370894   0.17330662   0.83426616    2.7036955   -10.740691 \n",
      "     360   -10.808291        158.5        158.5          100    1.5689286   0.57485581    1.2351347    2.8967954   -10.808291 \n",
      "     380   -11.132365        158.5        158.5          100    1.9304894    1.1369933    1.8133551    2.8411199   -11.132365 \n",
      "     400   -11.168527        158.5        158.5          100    1.9277665    1.1519308    1.8280985    2.8032701   -11.168527 \n",
      "     407   -11.181251        158.5        158.5          100    1.9839572    1.2351284    1.8856273    2.8311159   -11.181251 \n",
      "Loop time of 0.080941 on 1 procs for 407 steps with 100 atoms\n",
      "\n",
      "75.5% CPU use with 1 MPI tasks x no OpenMP threads\n",
      "\n",
      "Minimization stats:\n",
      "  Stopping criterion = max force evaluations\n",
      "  Energy initial, next-to-last, final = \n",
      "          939.04780435     -11.1764729995     -11.1812508351\n",
      "  Force two-norm initial, final = 401.034 2.70442\n",
      "  Force max component initial, final = 70.02 0.542305\n",
      "  Final line search alpha, max atom move = 0.0103102 0.00559128\n",
      "  Iterations, force evaluations = 407 1004\n",
      "\n",
      "MPI task timing breakdown:\n",
      "Section |  min time  |  avg time  |  max time  |%varavg| %total\n",
      "---------------------------------------------------------------\n",
      "Pair    | 0.014533   | 0.014533   | 0.014533   |   0.0 | 17.96\n",
      "Bond    | 0.023352   | 0.023352   | 0.023352   |   0.0 | 28.85\n",
      "Neigh   | 0.0042591  | 0.0042591  | 0.0042591  |   0.0 |  5.26\n",
      "Comm    | 0.00019002 | 0.00019002 | 0.00019002 |   0.0 |  0.23\n",
      "Output  | 0.035332   | 0.035332   | 0.035332   |   0.0 | 43.65\n",
      "Modify  | 0          | 0          | 0          |   0.0 |  0.00\n",
      "Other   |            | 0.003274   |            |       |  4.05\n",
      "\n",
      "Nlocal:    100 ave 100 max 100 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "Nghost:    0 ave 0 max 0 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "Neighs:    470 ave 470 max 470 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "FullNghs:  940 ave 940 max 940 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "\n",
      "Total # of neighbors = 940\n",
      "Ave neighs/atom = 9.4\n",
      "Ave special neighs/atom = 5.88\n",
      "Neighbor list builds = 55\n",
      "Dangerous builds = 0\n",
      "All done\n",
      "Total wall time: 0:00:00\n",
      "--- 1.1596903800964355 seconds ---\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "start_time = time.time()\n",
    "!lmp_serial < in.deform_polymer_chain.txt\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, the dump files should be able to be used to examine what is happening during the simulation.  Notice that the temperature jumps around quite a bit, a function of the small number of (unconstrained) atoms and the high temperature.  In the minimization stage, the energy starts very high and is rapidly reduced to a minimum energy, i.e., the minimum energy state at 0 K.\n",
    "\n",
    "<table><tr><td><img src='https://icme.hpc.msstate.edu/mediawiki/images/8/89/Equ_plus_Min.gif' width=\"300\"><center><br>Equilibration process followed by minimization</center></td></tr></table>\n"
   ]
  },
  {
   "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",
    "In general, this script does equilibration and minimization to the polymer chain. Polymer chain data file named `PE_cl100.txt` should be in the present working directory. The line `dump 1 all cfg 6 dump.PE_cl100_*.cfg mass type xs ys zs c_csym c_peratom fx fy fz` is used to output information during the simulation; this can be moved before the equilibration part of the script to output the process of equilibration. Please note that you need to change the time step for the dump command so it won't give too many or too few dump files. It may be reasonable to choose an `every n timesteps` number to correspond to roughly 100 snapshots. \n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "###################################### <br>\n",
    "# VARIABLES <br>\n",
    "variable fname index PE_cl100.txt <br>\n",
    "variable simname index PE_cl100 <br> <br>\n",
    "###################################### <br>\n",
    "# INITIALIZATION <br>\n",
    "units real <br>\n",
    "boundary f f f <br>\n",
    "atom_style molecular <br>\n",
    "log log.\\${simname}.txt <br>\n",
    "read_data \\${fname}\n",
    "</div>\n",
    "\n",
    "This part is just opens the data file, defines the boundary conditions, units, logfile's name, etc. \n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "###################################### <br>\n",
    "# DEFINE INTERATOMIC POTENTIAL <br>\n",
    "# Dreiding potential <br>\n",
    "neighbor 0.4 bin <br>\n",
    "neigh_modify every 10 one 10000 <br>\n",
    "bond_style      harmonic <br>\n",
    "bond_coeff 1 350 1.53 <br>\n",
    "angle_style     harmonic <br>\n",
    "angle_coeff 1 60 109.5 <br>\n",
    "dihedral_style multi/harmonic <br>\n",
    "dihedral_coeff 1 1.73 -4.49 0.776 6.99 0.0 <br>\n",
    "pair_style lj/cut 10.5 <br>\n",
    "pair_coeff 1 1 0.112 4.01 10.5\n",
    "</div>\n",
    "\n",
    "As mentioned earlier, this section gives the information about the potential functions/parameters to be used for the bonds, angles, dihedrals in a chain. Bond_style and bond_coeff defines the type on the force field between atoms and a magnitude of this fields. \"1\" here corresponds to the second column of the \"Bonds\" section of the data file. Thus, every atom pair with \"1\" in the second column has properties that correspond to that potential function during the simulation. The angles and dihedrals are similar. Angle_* and dihedral_* lines define the angles and dihedral angles between atoms in the polymer chain. Pair_* commands used to define pair potentials between atoms that are within a cutoff distance.\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "###################################### <br>\n",
    "# DEFINE COMPUTES  <br>\n",
    "compute csym all centro/atom fcc <br>\n",
    "compute peratom all pe/atom  <br>\n",
    "compute eng all pe/atom  <br>\n",
    "compute eatoms all reduce sum c_eng\n",
    "</div>\n",
    "\n",
    "This commands are used to define which properties are to be calculated during the simulation. For example, `pe/atom` simply means to calculate the potential energy for each atom. Information about other possible properties to calculate can be found in the LAMMPS documentation.\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "##################################################### <br>\n",
    "# EQUILIBRATION <br>\n",
    "# Langevin dynamics at 500 K <br>\n",
    "velocity all create 500.0 1231 <br>\n",
    "fix 1 all nve/limit 0.05 <br>\n",
    "fix 2 all langevin 500.0 500.0 10.0 904297 <br>\n",
    "thermo_style custom step temp  <br>\n",
    "thermo          100 <br>\n",
    "timestep 1 <br>\n",
    "run 10000 <br>\n",
    "unfix 1 <br>\n",
    "unfix 2 <br>\n",
    "write_restart restart.${simname}.dreiding1 <br>\n",
    "</div>\n",
    "\n",
    "The equilibration step allow atoms in the polymer chain to move freely at a certain temperature. The `velocity` command instantiates the temperature for the various atoms within the chain. The `fix` command performs the check of the system every time step and updates the velocities and positions of the atoms. The `thermo` commands define the thermodynamic output during the simulation. The `run` command starts the simulation.  After the 10,000 timesteps are completed, the `unfix` commands remove the fixes; this is just a good book-keeping practice when the LAMMPS input scripts start to get a little more complicated.\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "##################################################### <br>\n",
    "# MINIMIZATION <br>\n",
    "dump 1 all cfg 6 dump.PE_cl100_*.cfg mass type xs ys zs c_csym c_peratom fx fy fz <br>\n",
    "reset_timestep 0  <br>\n",
    "fix 1 all nvt temp 500.0 500.0 100.0 <br>\n",
    "thermo 20  <br>\n",
    "thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms  <br>\n",
    "min_style cg <br>\n",
    "minimize 1e-15 1e-15 1000 1000\n",
    "</div>\n",
    "\n",
    "The other step in the program is the minimization. This process finds the minimum energy condition for this configuration. The parameters for minimize command include: (left to right) stopping tolerance for energy, stopping tolerance for force, max iterations of minimization routine and max number of force/energy evaluations, where stopping tolerance for energy is the first parameter and the max number of force/energy evaluations is the last one. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## FAQs \n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-danger\">\n",
    "    <strong>Question 1</strong>: None yet.  \n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Links\n",
    "\n",
    "* [Click here to open Tutorial 1](LAMMPS-Tutorials-01.ipynb)\n",
    "* [Click here to open Tutorial 2](LAMMPS-Tutorials-02.ipynb)\n",
    "* [Click here to open Tutorial 3](LAMMPS-Tutorials-03.ipynb)\n",
    "* [Click here to open Tutorial 4](LAMMPS-Tutorials-04.ipynb)\n",
    "* [Click here to open Tutorial 5](LAMMPS-Tutorials-05.ipynb)\n",
    "* [Click here to open Tutorial 6](LAMMPS-Tutorials-06.ipynb)\n",
    "* [Click here to open Tutorial 7](LAMMPS-Tutorials-07.ipynb)\n",
    "* [Click here to open the next tutorial](LAMMPS-Tutorials-09.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). \n",
    "1. Hossain, D., Tschopp, M.A. (corresponding author), Ward, D.K., Bouvard, J.L., Wang, P., Horstemeyer, M.F., \"[Molecular dynamics simulations of deformation mechanisms of amorphous polyethylene](https://doi.org/10.1016/j.polymer.2010.10.009),\" Polymer, 51 (2010) 6071-6083."
   ]
  }
 ],
 "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": 4
}