Glass Transition Temperature Working.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 0,
        "width": 12
       },
       "report_default": {}
      }
     }
    }
   },
   "source": [
    "# Calculate the glass transition temperature, $T_g$, of an atomistic amorphous system using simulation tools at nanoHUB.org"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# System setup\n",
    "import sys, os, tempfile\n",
    "\n",
    "#sys.path.append('/apps/share64/debian7/ipynb_tools')\n",
    "sys.path.append('/Storage/davidelbert/persistent/hackathon2017/ben/nanoHUB_remote')\n",
    "\n",
    "# Run this calculation in a temporary directory\n",
    "workdir = tempfile.mkdtemp(dir=os.getcwd())\n",
    "#%cd $workdir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nanoHUB_remote\n",
    "\n",
    "# In mysecrets.py, set your web app and account secrets.\n",
    "# \n",
    "# Example, with user credentials:\n",
    "# auth_data = {\n",
    "#    'client_id': '',       # Get this when you register a web app\n",
    "#    'client_secret': '',   # Get this when you register a web app\n",
    "#    'grant_type': 'password',\n",
    "#    'username': '',        # Get this when you create a nanoHUB account\n",
    "#    'password': ''         # Get this when you create a nanoHUB account\n",
    "# }\n",
    "#\n",
    "# This design is strictly for convenience so that this notebook can be \n",
    "# shared without storing any secrets.\n",
    "from nanoHUB_remote.mysecrets import auth_data\n",
    "\n",
    "# Authenticate; use headers in all subsequent steps\n",
    "headers = nanoHUB_remote.authenticate(auth_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "No module named hublib",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mImportError\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-9cd930afdfda>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mhublib\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mImportError\u001b[0m: No module named hublib"
     ]
    }
   ],
   "source": [
    "import hublib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hubwf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 4,
        "width": 4
       },
       "report_default": {}
      }
     }
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Run PolymerModeler to pack chains into a box\n",
    "import hubwf\n",
    "\n",
    "# Create an interface to the installed PolymerModeler tool\n",
    "task = hubwf.RapptureTask('polymod')\n",
    "\n",
    "# Select the polymer\n",
    "task.set_input_value('Polymer', 'PMMA (atactic)')\n",
    "\n",
    "# Set the length of the chains (monomers per chain)\n",
    "task.set_input_value('Monomers per chain', 10)\n",
    "\n",
    "# Set the number of chains\n",
    "task.set_input_value('Number of chains', 30)   \n",
    "\n",
    "# Have the tool generate input files for LAMMPS\n",
    "task.set_input_value('Simulation choice', 'LAMMPS input files only')\n",
    "\n",
    "# Generate LAMMPS commands to relax the structure\n",
    "task.set_input_value('Minimization levels', 3)\n",
    "task.set_input_value('Minimization steps', 5000)\n",
    "\n",
    "# Generate LAMMPS commands to thermalize at 800K\n",
    "task.set_input_value('Ensemble', 'npt')\n",
    "task.set_input_value('Number of MD steps', 10000)\n",
    "task.set_input_value('Temperature', 800, exact=True)\n",
    "\n",
    "# Run the tool\n",
    "task.run()\n",
    "print 'Finished running PolymerModeler'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 12,
        "hidden": false,
        "row": 8,
        "width": 11
       },
       "report_default": {}
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "# View the initial structure built by PolymerModeler\n",
    "import nglview as nv\n",
    "import subprocess as sp\n",
    "\n",
    "tmp_pdb = 'temp.pdb'\n",
    "polymod_pdb = 'polymod.pdb'\n",
    "with open(tmp_pdb, 'w') as f:\n",
    "    f.write(task.get_output_pdb('Built structures, unwrapped'))\n",
    "sp.call(['babel', tmp_pdb, polymod_pdb])  # add CONECT records\n",
    "v = nv.show_structure_file(polymod_pdb)\n",
    "v.representations = [{'type': 'ball+stick', 'params': {}}]\n",
    "v.parameters = {'backgroundColor': 'black'}\n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Relax and thermalize the initial structure built by PolymerModeler\n",
    "import hublib.use\n",
    "%use lammps-15May15\n",
    "\n",
    "nodes = 2\n",
    "walltime = '1:00:00'\n",
    "datafile = 'polymer_relax.data'\n",
    "with open(datafile, 'w') as f:\n",
    "    f.write(task.get_output_value('LAMMPS data file'))\n",
    "infile = 'lammps_relax.in'\n",
    "with open(infile, 'w') as f:\n",
    "    f.write(task.get_output_value('LAMMPS input file'))\n",
    "\n",
    "!submit -n $nodes -w $walltime -i $datafile lammps-15May15-parallel -in $infile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Cool and compress the structure with LAMMPS\n",
    "infile = 'lammps_cool.in'\n",
    "T_steps = 10000\n",
    "lammps_input = \"\"\"\n",
    "read_restart      restart.lammps\n",
    "\n",
    "neighbor          2.0 bin\n",
    "thermo_style      custom temp density vol etotal ke pe ebond eangle edihed eimp evdwl ecoul elong press pxx pyy pzz pxy pxz pyz lx ly lz\n",
    "thermo            100\n",
    "thermo_modify     flush yes line multi format float %16.6g\n",
    "\n",
    "timestep          4\n",
    "neigh_modify      every 1 delay 5\n",
    "kspace_style      pppm 1e-4\n",
    "run_style         respa 3 2 2 bond 1 pair 2 kspace 3\n",
    "velocity          all create 800.0 1214556 mom yes rot yes dist gaussian\n",
    "\n",
    "# Loop\n",
    "variable          i loop 26\n",
    "label             loopa\n",
    "reset_timestep    0\n",
    "variable          T equal 800+($i-1)*(-20)\n",
    "\n",
    "fix               1 all npt temp $T $T 100.0 iso 1.0 1.0 1000.\n",
    "dump              2 all custom 500 strcool_$TK.dump id type x y z\n",
    "run               {0}\n",
    "write_data        str_cool$TK.data\n",
    "unfix             1\n",
    "undump            2\n",
    "next              i\n",
    "jump              {1} loopa\n",
    "\"\"\".format(T_steps, infile)\n",
    "\n",
    "with open(infile, 'w' ) as f:\n",
    "    f.write(lammps_input)\n",
    "walltime = '12:00:00'\n",
    "\n",
    "!submit -n $nodes -w $walltime -i restart.lammps lammps-15May15-parallel -in $infile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the density and temperature data from LAMMPS output; find Tg\n",
    "from Tgcalc import getVDT, getVDTdata, getTg_BL\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook\n",
    "\n",
    "getVDT('log.lammps', T_steps)\n",
    "data_file = 'log.rec_VDT'\n",
    "\n",
    "# Plot density vs. T\n",
    "T, rho, err = getVDTdata(data_file)\n",
    "plt.style.use('ggplot')\n",
    "plt.plot(T, rho, 'b-')\n",
    "plt.xlabel('T (K)')\n",
    "plt.ylabel('Density (g/cc)')\n",
    "plt.show()\n",
    "\n",
    "# Calculate Tg\n",
    "print 'Tg from bilinear fit:', getTg_BL(data_file)"
   ]
  }
 ],
 "metadata": {
  "extensions": {
   "jupyter_dashboards": {
    "activeView": "grid_default",
    "version": 1,
    "views": {
     "grid_default": {
      "cellMargin": 10,
      "defaultCellHeight": 20,
      "maxColumns": 12,
      "name": "grid",
      "type": "grid"
     },
     "report_default": {
      "name": "report",
      "type": "report"
     }
    }
   }
  },
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}