{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3cc0aea6",
   "metadata": {},
   "source": [
    "# Calculation of bulk properties with molecular dynamics\n",
    "\n",
    "In this notebook, we will use a molecular dynamics simulation to compute the thermal expansion coefficient of solid Xe at 77 K and 1 atm pressure. In the isothermal-isobaric ensemble (number of atoms, temperature, and pressure are all held constant - the NPT ensemble), the thermal expansion coefficient can be related to the volume, enthalpy, pressure, and temperature as:\n",
    "\n",
    "$$ \\alpha_V = \\frac{1}{Vk_BT^2}\\Bigl(\\langle HV \\rangle - \\langle H\\rangle \\langle V \\rangle \\Bigr) $$\n",
    "\n",
    "where $k_B$ is Boltzmann's constant and the enthalpy is:\n",
    "\n",
    "$$ H = E + PV $$\n",
    "\n",
    "The quantities inside the angle brackets represent average values. So $\\langle H \\rangle \\langle V \\rangle$ is the average enthalpy times the average volume, while $\\langle HV \\rangle$ is the average of the product of $H$ and $V$ for every frame.\n",
    "\n",
    "To perform this calculation, the `xe-729.pdb` file contains a lattice of 729 Xe atoms, and the `xenon.xml` file contains the Lennard-Jones parameters for the Xe atom. Because we want to simulate bulk Xe, and not the properties of a Xe nanocrystal, we will use a periodic boundary condition with a cutoff of 1.2 nm. Also, to maintain constant temperature, a Langevin thermostat is used with a target temperature of 77 K and a collision rate of 1 ps$^{-1}$. Finally, constant pressure is maintained with a Monte Carlo barostat with an attempt interval set to 25 steps.\n",
    "\n",
    "In the code below, every 1000 steps output is written to a `teajectory.pdb` file, which can be viewed with VMD or another visualization package. In addition, the step number, time, PE, KE, total energy, temperature, volume, and density are written to a `data.csv` file every 1000 steps and to standard output every 10000 steps.\n",
    "\n",
    "The simulation is initialized with random velocities drawn from a 77 K distribution, then allowed to minimize and equilibrate before the production run is started. The production run goes for 1000000 steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1141df35",
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################################################################\n",
    "# this script was generated by openmm-builder. to customize it further,\n",
    "# you can save the file to disk and edit it with your favorite editor.\n",
    "##########################################################################\n",
    "\n",
    "from __future__ import print_function\n",
    "from openmm import app\n",
    "import openmm as mm\n",
    "from simtk import unit\n",
    "from sys import stdout\n",
    "\n",
    "pdb = app.PDBFile('xe-729.pdb')\n",
    "forcefield = app.ForceField('xenon.xml')\n",
    "\n",
    "system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffPeriodic, \n",
    "    nonbondedCutoff=1.2*unit.nanometers, constraints=None, rigidWater=False)\n",
    "integrator = mm.LangevinIntegrator(77*unit.kelvin, 1.0/unit.picoseconds, \n",
    "    5.0*unit.femtoseconds)\n",
    "system.addForce(mm.MonteCarloBarostat(1*unit.atmospheres, 77*unit.kelvin, 25))\n",
    "\n",
    "platform = mm.Platform.getPlatformByName('CPU')\n",
    "simulation = app.Simulation(pdb.topology, system, integrator, platform)\n",
    "simulation.context.setPositions(pdb.positions)\n",
    "\n",
    "print('Minimizing...')\n",
    "simulation.minimizeEnergy(maxIterations=100)\n",
    "\n",
    "simulation.context.setVelocitiesToTemperature(77*unit.kelvin)\n",
    "print('Equilibrating...')\n",
    "simulation.step(100)\n",
    "\n",
    "simulation.reporters.append(app.PDBReporter('trajectory.pdb', 1000))\n",
    "simulation.reporters.append(app.StateDataReporter(stdout, 10000, step=True, \n",
    "    time=True, potentialEnergy=True, kineticEnergy=True, totalEnergy=True, \n",
    "    temperature=True, volume=True, density=True, separator='\\t'))\n",
    "simulation.reporters.append(app.StateDataReporter('data.csv', 1000, step=True, \n",
    "    time=True, potentialEnergy=True, kineticEnergy=True, totalEnergy=True, \n",
    "    temperature=True, volume=True, density=True, separator=','))\n",
    "\n",
    "print('Running Production...')\n",
    "simulation.step(10000)\n",
    "print('Done!')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a760b563",
   "metadata": {},
   "source": [
    "At this point, you can open the `trajectory.pdb` file in VMD to ensure that you see a Xe solid, where the atoms are vibrating about their lattice sites. Under \"Extensions > Analysis\" you can find the menu for calculating the radial distribution function between atoms: setting `name = Xe1` for both selection fields will compute the function.\n",
    "\n",
    "Now we can analyze the results of the simulation and compute the thermal expansion coefficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8e0c018",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_csv('data.csv')\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9897b54",
   "metadata": {},
   "source": [
    "The thermostat does not force the temperature to 77 K at every frame; it works by altering the equation of motion to add random forces to slow particles and friction to fast ones, such that over time, the average temperature approaches 77 K. We can visualize a histogram of the temperature (and indeed the other properties) to see its distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f47285dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "ax.hist(df['Temperature (K)'],bins=100)\n",
    "ax.set_xlabel('Temperature (K)')\n",
    "ax.set_ylabel('Frames')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "699d78d6",
   "metadata": {},
   "source": [
    "Now we are ready to compute the thermal expansion coefficient. Because the units of pressure in the simulation are atm and the units of volume are nm$^3$, some unit conversions are needed for calculating enthalpy: $PV$ must have the same units as $E$ (kJ/mol). This is accomplished by converting atm to Pa (J/m$^3$), converting atoms to moles with Avogadro's number, converting J to kJ, and nm$^3$ to m$^3$.\n",
    "\n",
    "Unit conversions are available in the scipy.constants module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4eff2f69",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.constants as spc\n",
    "\n",
    "atm = spc.value('standard atmosphere') # 101325\n",
    "NA = spc.value('Avogadro constant')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10a0b817",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "V = df['Box Volume (nm^3)']\n",
    "T = df['Temperature (K)']\n",
    "E = df['Total Energy (kJ/mole)']\n",
    "\n",
    "Vavg = np.mean(V)\n",
    "Tavg = np.mean(T)\n",
    "Eavg = np.mean(E)\n",
    "P = 1.0\n",
    "print(Vavg,Tavg,Eavg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b8064f1",
   "metadata": {},
   "source": [
    "Now we can compute the enthalpy in each frame, the total average enthalpy, the product of enthalpy and volume ($HV$) for each frame, and the average $\\langle HV \\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f71b547e",
   "metadata": {},
   "outputs": [],
   "source": [
    "H = E + P*V*atm*NA/spc.kilo*(spc.nano**3) #units kJ/mol\n",
    "Havg = Eavg + P*Vavg*atm*NA/spc.kilo*(spc.nano**3) #units kJ/mol\n",
    "HVavg = np.mean(H*V) #units kJ nm^3 / mol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3739aaeb",
   "metadata": {},
   "source": [
    "And finally, the thermal expansion coefficient $\\alpha$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6874ed49",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'alpha = {1/(Vavg*spc.value(\"Boltzmann constant\")*NA/spc.kilo*Tavg**2)*(HVavg - Havg*Vavg)} 1/K')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76c69aa6",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}