Predicting Phase Behavior.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b43c086c-80c0-40ac-9397-0f2ef1c8cc9d",
   "metadata": {},
   "source": [
    "# Predicting phase behavior\n",
    "\n",
    "This notebook will focus on using molecular dynamics (MD) simulation to predict the phase behavior of point particles that interact via the Lennard-Jones potential. \n",
    "\n",
    "To perform these simulations we will use HOOMD-Blue, as we can easily and quickly run this from the jupyter notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6517d0b0-8257-4bbb-925d-b76b31118808",
   "metadata": {},
   "source": [
    "## Single component Lennard-Jones system\n",
    "\n",
    "First, let us consider a single component system composed of Lennard-Jones (LJ) spheres, identical to what we have previously examined in prior notebooks.  Here, the goal is to use simulation to perform simulations of the LJ particles in different phases: solid, liquid and gas. \n",
    "\n",
    "We will do this by modifying the density and temperature of the simulation to move to different points in the phase diagram. \n",
    "\n",
    "The follow paper provides a clear phase diagram that we can use to direct our simulations:\n",
    ">T. E. Karakasidis, A. Fragkou, and A. Liakopoulos, \"System dynamics revealed by recurrence quantification analysis: Application to molecular dynamics simulations\" Phys. Rev. E 76, 021120 – Published 24 August 2007\n",
    "https://journals.aps.org/pre/abstract/10.1103/PhysRevE.76.021120\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7213dc5d-48d6-4267-8c44-9f9eb9feb14c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(url='img/lj_phase_diagram.png')  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c09723e-1d0b-480d-8e66-418b68e322e3",
   "metadata": {},
   "source": [
    "### Radial Distribution Function \n",
    "A common way to help evaluate the phase of a system is to use the radial distribution function (RDF, sometimes referred to as g(r)). In the RDF, we create a histogram of the distance between all particle pairs, normalized by the density of the system multipied by the volume of the bin (i.e., $\\rho *4 \\pi r^2dr$ where $\\rho$ is the number density and dr is the width of the bin). Systems where particle positions are strongly correlated will show larger values in the RDF. An example of the RDF for different phases is shown below.  See [this link](https://chem.libretexts.org/Bookshelves/Biological_Chemistry/Concepts_in_Biophysical_Chemistry_(Tokmakoff)/01%3A_Water_and_Aqueous_Solutions/01%3A_Fluids/1.02%3A_Radial_Distribution_Function) for more info on the RDF.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58a58a90-98ff-442c-9bef-2f67f6608b0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(url='https://upload.wikimedia.org/wikipedia/commons/b/b8/Simulated_Radial_Distribution_Functions_for_Solid%2C_Liquid%2C_and_Gaseous_Argon.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4817454d-bfd4-4db2-87d5-7ba5bb19f532",
   "metadata": {},
   "source": [
    "### Simulation workflow\n",
    "\n",
    "We'll need to modify the prior simulation script to allow us to change the box volume, in order to be able to adjust the system density.  A tutorial for modifying the box that provides step-by-step details is available in the (HOOMD-Blue documentation)[https://hoomd-blue.readthedocs.io/en/v3.2.0/tutorial/01-Introducing-Molecular-Dynamics/03-Compressing-the-System.html]. \n",
    "\n",
    "\n",
    "The structure of the script is also now slightly changed compared to the prior examples. In the simulation runtime section, we will now have 3 distinct stages:\n",
    "1) compression  -- adjust the box to achieve the target density\n",
    "2) relaxation -- allow the system to reach steady state, especially important after the box size perturbation\n",
    "3) production -- final portion of the simulation for data collection\n",
    "\n",
    "It is generally good practice to break up a simulation up in sub steps.  This can help to avoid accidentically utilizing non-equilibrium data when reporting properties. This can also make it more obvious if there is an issue, at what point during the simulation workflow is this occuring.    This may often be required to ensure stability, allowing each substep to use different init parameters (e.g., temperature, timestep, integration algorithm, etc.)  For example, the initial configuration may not be especially favorable, with individual particles experiencing large forces; to relieve the system of these high energy/large forces we often use a much smaller timestep initially after performing energy minimization.  Additionally, it is often common practice to start a simulation at higher temperature than the desired target temperature, to help decorrelation from the initial configuration and speed the process of reaching equilibrium (as temperature increases, so does system mobility).\n",
    "\n",
    "#### Simulation outputs\n",
    "This simulation run will generate several GSD files. \n",
    "* thermo_full.gsd -- thermodynamic quantities logged during the entire duration of the simulation\n",
    "* thermo_prod.gsd -- thermodynamic quantities logged during the final production portion of the simulation\n",
    "* trajectory_full.gsd -- trajectory logged during the entire duration of the simulation\n",
    "* trajectory_prod.gsd -- trajectory logged during the final production portion of the simulation\n",
    "\n",
    "We will also print to the screen the current timestep and timesteps-per-second periodically during the simulation run, with print statements also noting when each of the substeps have been completed.\n",
    "\n",
    "We've also made it easier to modify by placing variables for the variables we will be adjust at the top of the script file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9e0a634-1be2-421b-b4d1-0bb789828185",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hoomd\n",
    "import gsd.hoomd\n",
    "\n",
    "system_temperature = 1.0\n",
    "final_rho = 0.5\n",
    "compression_time = 20000\n",
    "relax_time = 50000\n",
    "production_time = 50000\n",
    "\n",
    "\n",
    "\n",
    "cpu = hoomd.device.CPU()\n",
    "sim = hoomd.Simulation(device=cpu, seed=1)\n",
    "\n",
    "#########################\n",
    "# system initialization\n",
    "#########################\n",
    "\n",
    "sim.create_state_from_gsd(filename='datafiles/lj.gsd')\n",
    "sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=system_temperature)\n",
    "\n",
    "\n",
    "#########################\n",
    "# interaction definition\n",
    "#########################\n",
    "\n",
    "cell = hoomd.md.nlist.Cell(buffer=1.0)\n",
    "lj = hoomd.md.pair.LJ(nlist=cell)\n",
    "lj.params[('LJ', 'LJ')] = dict(epsilon=1, sigma=1)\n",
    "lj.r_cut[('LJ', 'LJ')] = 2.5\n",
    "\n",
    "#########################\n",
    "# integrator setup\n",
    "#########################\n",
    "\n",
    "\n",
    "nvt = hoomd.md.methods.NVT(kT=system_temperature, filter=hoomd.filter.All(), tau=1.0)\n",
    "integrator = hoomd.md.Integrator(dt=0.005)\n",
    "integrator.forces.append(lj)\n",
    "integrator.methods.append(nvt)\n",
    "sim.operations.integrator = integrator\n",
    "\n",
    "#########################\n",
    "# runtime parameters\n",
    "#########################\n",
    "\n",
    "# add a logger to output the current step and timesteps per second (TPS) \n",
    "logger = hoomd.logging.Logger(categories=['scalar', 'string'])\n",
    "logger.add(sim, quantities=['timestep', 'tps'])\n",
    "table = hoomd.write.Table(trigger=hoomd.trigger.Periodic(period=5000),\n",
    "                          logger=logger)\n",
    "sim.operations.writers.append(table)\n",
    "\n",
    "# add a logger to output the thermodynamic data during the entire run \n",
    "thermodynamic_properties = hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All())\n",
    "sim.operations.computes.append(thermodynamic_properties)\n",
    "\n",
    "log_thermo_full = hoomd.logging.Logger()\n",
    "log_thermo_full.add(thermodynamic_properties)\n",
    "gsd_thermo_full = hoomd.write.GSD(filename='thermo_full.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(100),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.Null(),\n",
    "                             log=log_thermo_full)\n",
    "sim.operations.writers.append(gsd_thermo_full)\n",
    "\n",
    "\n",
    "log_traj_full = hoomd.logging.Logger()\n",
    "gsd_traj_full = hoomd.write.GSD(filename='trajectory_full.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(5000),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.All(),\n",
    "                             log=log_traj_full)\n",
    "sim.operations.writers.append(gsd_traj_full)\n",
    "\n",
    "# ================\n",
    "# compress the box\n",
    "# ================\n",
    "\n",
    "\n",
    "# ramp will define a switching factor for adjusting from the initial to the final box size\n",
    "ramp = hoomd.variant.Ramp(A=0, B=1, t_start=sim.timestep, t_ramp=compression_time)\n",
    "\n",
    "initial_box = sim.state.box\n",
    "final_box = hoomd.Box.from_box(initial_box)  # make a copy of initial_box\n",
    "final_box.volume = sim.state.N_particles / final_rho\n",
    "box_resize_trigger = hoomd.trigger.Periodic(10)\n",
    "box_resize = hoomd.update.BoxResize(box1=initial_box,\n",
    "                                    box2=final_box,\n",
    "                                    variant=ramp,\n",
    "                                    trigger=box_resize_trigger)\n",
    "sim.operations.updaters.append(box_resize)\n",
    "sim.run(compression_time+1)\n",
    "if sim.state.box == final_box:\n",
    "    print(\"system box size matches target\")\n",
    "    \n",
    "sim.operations.updaters.remove(box_resize)\n",
    "print(\"compression complete\")\n",
    "\n",
    "# ================\n",
    "# relaxation\n",
    "# ================\n",
    "\n",
    "sim.run(relax_time)\n",
    "print(\"relaxation complete\")\n",
    "\n",
    "# ================\n",
    "# production \n",
    "# ================\n",
    "\n",
    "log_thermo = hoomd.logging.Logger()\n",
    "log_thermo.add(thermodynamic_properties)\n",
    "gsd_thermo = hoomd.write.GSD(filename='thermo_prod.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(100),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.Null(),\n",
    "                             log=log_thermo)\n",
    "sim.operations.writers.append(gsd_thermo)\n",
    "\n",
    "log_traj = hoomd.logging.Logger()\n",
    "gsd_traj = hoomd.write.GSD(filename='trajectory_prod.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(5000),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.All(),\n",
    "                             log=log_traj)\n",
    "sim.operations.writers.append(gsd_traj)\n",
    "\n",
    "sim.run(production_time)\n",
    "\n",
    "hoomd.write.GSD.write(state=sim.state, filename='final.gsd', mode='wb')\n",
    "print(\"simulation complete\")\n",
    "\n",
    "#delete the instances we defined ensure writers are closed\n",
    "del sim, gsd_thermo, gsd_traj, thermodynamic_properties, log_thermo, log_traj\n",
    "del log_thermo_full, gsd_thermo_full, logger, table, log_traj_full, gsd_traj_full\n",
    "del integrator, nvt, lj, cell, cpu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc34c0c3-384d-45cb-a55b-2fd797d98b99",
   "metadata": {},
   "source": [
    "Take note of the value of the timesteps per second (TPS) during the simulation run.  Higher values mean the code is running more efficiently.  For a system where we compress the system (i.e., density is higher than the initial density of 0.125 in the lj.gsd data file), the TPS value drops during the compression stage.  Recall the discussion of neighorlists and cutoffs; as the system density increases, the number of neighbors we need to consider will also increase. \n",
    "\n",
    "# Exercises\n",
    "\n",
    "## Exercise 1: fluid phase simulation\n",
    "\n",
    "Pick a few statepoints from the phase diagram to simulate a fluid phase, adjusting the density and temperature as needed in the simulation code above\n",
    " * Visualize the trajectory using the code below. \n",
    " * Examine the full thermodynamic output (see below) to ensure that the system has reached a steady state before the production run.\n",
    "  * As discussed earlier, this could be indicated by convergence of a quantity such as PE; for the sake of this example, visual inspection of the plot is more than sufficient, although more rigorous, mathematical approaches can be taken (which have the advantage of allowing automated checks for convergence).\n",
    "  * Add in code to read in the production thermo data and calculate the mean and std of the potential energy.\n",
    " * Use the following code to calculate and plot the RDF for the production data. How does this compare the example plots above?  Does the RDF you plotted seem to indicate a fluid phase? Why?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bed149ef-b46c-4bc1-a325-1f8a8f4fcbf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "matplotlib.style.use('default')\n",
    "import mdtraj as md\n",
    "from fresnel_render import render\n",
    "from fresnel_render import render_movie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a233770f-eedb-40d2-b2aa-b1f1306d2ac6",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj_log = gsd.hoomd.open('trajectory_prod.gsd', 'rb')\n",
    "\n",
    "render_movie(traj_log)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7ccd083-df78-4d15-a01f-a0ad6a367e7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "thermo_log = gsd.hoomd.open('thermo_full.gsd', 'rb')\n",
    "\n",
    "time = []\n",
    "pe = []\n",
    "for frame in thermo_log:\n",
    "    time.append(frame.configuration.step)\n",
    "    pe.append(frame.log['md/compute/ThermodynamicQuantities/potential_energy'][0])\n",
    "\n",
    "n = thermo_log[-1].log['md/compute/ThermodynamicQuantities/num_particles'][0]\n",
    "vol = thermo_log[-1].log['md/compute/ThermodynamicQuantities/volume'][0]\n",
    "\n",
    "print(\"density: \", n/vol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2672c427-24de-4e3c-bdcb-4c9313626b3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time, pe)\n",
    "plt.ylabel('potential energy')\n",
    "plt.xlabel('timestep')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e1c3a34-e504-489f-9382-d4a246ca6d83",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = md.load('trajectory_prod.gsd')\n",
    "r_traj_fluid, g_r_traj_fluid = md.compute_rdf(traj,  r_range=[0, traj.unitcell_lengths[0][1]/2.0], \n",
    "                                pairs = traj.top.select_pairs(\"all\", \"all\"), \n",
    "                                  bin_width=0.05, periodic=True, opt=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98a55211-dadc-4cc3-9648-39d6c0f146ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(r_traj_fluid, g_r_traj_fluid)\n",
    "plt.ylabel('g(r)')\n",
    "plt.xlabel('r')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "112c4b74-1b9a-4113-8743-dff9c0c80606",
   "metadata": {},
   "source": [
    "## Exercise 2: Simulation of a solid phase\n",
    "\n",
    "Now, adjust the density to 1.4 and temperature to be 1.0 in the code and rerun the simulation. As you can see in the phase diagram (shown again below), this should be a solid crystalline phase for the LJ particles.  To avoid having to move around in the notebook, the same code is reproduced below. \n",
    " * Render the simulation trajectory of the production step as movie. Does the phase look like a well ordered solid?\n",
    " * Calculate and plot the RDF. How does this compare the example plots above?  Would you conclude from both RDF and visual inspection that you have formed a solid phase? Why or why not?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "112bc4cf-ace6-482e-b2fc-a6cd8b98be9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(url='img/lj_phase_diagram.png')  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10c2ea2e-b582-44db-8ccd-ae78c2be0619",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hoomd\n",
    "import gsd.hoomd\n",
    "\n",
    "system_temperature = 1.0\n",
    "final_rho = 0.5\n",
    "compression_time = 20000\n",
    "relax_time = 50000\n",
    "production_time = 50000\n",
    "\n",
    "\n",
    "\n",
    "cpu = hoomd.device.CPU()\n",
    "sim = hoomd.Simulation(device=cpu, seed=1)\n",
    "\n",
    "#########################\n",
    "# system initialization\n",
    "#########################\n",
    "\n",
    "sim.create_state_from_gsd(filename='datafiles/lj.gsd')\n",
    "sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=system_temperature)\n",
    "\n",
    "\n",
    "#########################\n",
    "# interaction definition\n",
    "#########################\n",
    "\n",
    "cell = hoomd.md.nlist.Cell(buffer=1.0)\n",
    "lj = hoomd.md.pair.LJ(nlist=cell)\n",
    "lj.params[('LJ', 'LJ')] = dict(epsilon=1, sigma=1)\n",
    "lj.r_cut[('LJ', 'LJ')] = 2.5\n",
    "\n",
    "#########################\n",
    "# integrator setup\n",
    "#########################\n",
    "\n",
    "\n",
    "nvt = hoomd.md.methods.NVT(kT=system_temperature, filter=hoomd.filter.All(), tau=1.0)\n",
    "integrator = hoomd.md.Integrator(dt=0.005)\n",
    "integrator.forces.append(lj)\n",
    "integrator.methods.append(nvt)\n",
    "sim.operations.integrator = integrator\n",
    "\n",
    "#########################\n",
    "# runtime parameters\n",
    "#########################\n",
    "\n",
    "# add a logger to output the current step and timesteps per second (TPS) \n",
    "logger = hoomd.logging.Logger(categories=['scalar', 'string'])\n",
    "logger.add(sim, quantities=['timestep', 'tps'])\n",
    "table = hoomd.write.Table(trigger=hoomd.trigger.Periodic(period=5000),\n",
    "                          logger=logger)\n",
    "sim.operations.writers.append(table)\n",
    "\n",
    "# add a logger to output the thermodynamic data during the entire run \n",
    "thermodynamic_properties = hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All())\n",
    "sim.operations.computes.append(thermodynamic_properties)\n",
    "\n",
    "log_thermo_full = hoomd.logging.Logger()\n",
    "log_thermo_full.add(thermodynamic_properties)\n",
    "gsd_thermo_full = hoomd.write.GSD(filename='thermo_full.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(100),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.Null(),\n",
    "                             log=log_thermo_full)\n",
    "sim.operations.writers.append(gsd_thermo_full)\n",
    "\n",
    "\n",
    "log_traj_full = hoomd.logging.Logger()\n",
    "gsd_traj_full = hoomd.write.GSD(filename='trajectory_full.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(5000),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.All(),\n",
    "                             log=log_traj_full)\n",
    "sim.operations.writers.append(gsd_traj_full)\n",
    "\n",
    "# ================\n",
    "# compress the box\n",
    "# ================\n",
    "\n",
    "\n",
    "# ramp will define a switching factor for adjusting from the initial to the final box size\n",
    "ramp = hoomd.variant.Ramp(A=0, B=1, t_start=sim.timestep, t_ramp=compression_time)\n",
    "\n",
    "initial_box = sim.state.box\n",
    "final_box = hoomd.Box.from_box(initial_box)  # make a copy of initial_box\n",
    "final_box.volume = sim.state.N_particles / final_rho\n",
    "box_resize_trigger = hoomd.trigger.Periodic(10)\n",
    "box_resize = hoomd.update.BoxResize(box1=initial_box,\n",
    "                                    box2=final_box,\n",
    "                                    variant=ramp,\n",
    "                                    trigger=box_resize_trigger)\n",
    "sim.operations.updaters.append(box_resize)\n",
    "sim.run(compression_time+1)\n",
    "if sim.state.box == final_box:\n",
    "    print(\"system box size matches target\")\n",
    "    \n",
    "sim.operations.updaters.remove(box_resize)\n",
    "print(\"compression complete\")\n",
    "\n",
    "# ================\n",
    "# relaxation\n",
    "# ================\n",
    "\n",
    "sim.run(relax_time)\n",
    "print(\"relaxation complete\")\n",
    "\n",
    "# ================\n",
    "# production \n",
    "# ================\n",
    "\n",
    "log_thermo = hoomd.logging.Logger()\n",
    "log_thermo.add(thermodynamic_properties)\n",
    "gsd_thermo = hoomd.write.GSD(filename='thermo_prod.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(100),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.Null(),\n",
    "                             log=log_thermo)\n",
    "sim.operations.writers.append(gsd_thermo)\n",
    "\n",
    "log_traj = hoomd.logging.Logger()\n",
    "gsd_traj = hoomd.write.GSD(filename='trajectory_prod.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(5000),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.All(),\n",
    "                             log=log_traj)\n",
    "sim.operations.writers.append(gsd_traj)\n",
    "\n",
    "sim.run(production_time)\n",
    "\n",
    "hoomd.write.GSD.write(state=sim.state, filename='final.gsd', mode='wb')\n",
    "print(\"simulation complete\")\n",
    "\n",
    "#delete the instances we defined ensure writers are closed\n",
    "del sim, gsd_thermo, gsd_traj, thermodynamic_properties, log_thermo, log_traj\n",
    "del log_thermo_full, gsd_thermo_full, logger, table, log_traj_full, gsd_traj_full\n",
    "del integrator, nvt, lj, cell, cpu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b10ca16-2454-4fc1-99f8-a548c27e6c02",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj_log = gsd.hoomd.open('trajectory_prod.gsd', 'rb')\n",
    "\n",
    "#render the last frame\n",
    "render(traj_log[-1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5774176f-bbe0-43f9-a774-848b8c7b784a",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = md.load('trajectory_prod.gsd')\n",
    "r_traj_solid1, g_r_traj_solid1 = md.compute_rdf(traj,  r_range=[0, traj.unitcell_lengths[0][1]/2.0], \n",
    "                                pairs = traj.top.select_pairs(\"all\", \"all\"), \n",
    "                                  bin_width=0.01, periodic=True, opt=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82299eda-9d4d-4000-beaa-4c15e7be0ea3",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(r_traj_solid1, g_r_traj_solid1)\n",
    "plt.ylabel('g(r)')\n",
    "plt.xlabel('r')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55f04d9e-4f82-4871-a50d-e492ee170aa1",
   "metadata": {},
   "source": [
    "## Exercise 3: Equilibrating the system\n",
    "\n",
    "Now, set the statepoint to be T = 0.25 and density = 1.4, and reduce and rerun the simulation in the code below.  Note, this code has the time for the compression and relaxation steps set to be 5000  and 20000 steps, respectively,  intentionally. \n",
    "* Calculate the RDF below; we will save the output into a new set of variables and plot this against the T=1.0 and density = 1.4 data from above. \n",
    "    * What differences do you observe?  While we would not necessarily expect these plots to be identical, if anything we would expect the lower temperature state to demonstrate more well defined ordering; clearly it is not. \n",
    "* Render the last frame of the trajectory.\n",
    "    * How does this visually compare to the higher temperature state and what you'd expect a crystal to look like?\n",
    " \n",
    "The differences are related to factors: we rapidly increased the rate of compression while simultaneously lowering the temperature.  This effectively has trapped the system in a metastable state that it has trouble escaping. \n",
    " * Increase the time over which compression occurs and rerun the simulation.  At what value you see better agreement between the two RDFs? \n",
    " * Take note of the mean and stdev of the PE for different values of compression time. How do they compare and do you see convergence as the time increases?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41a7a046-8b96-493d-bb9c-eb709c7012b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hoomd\n",
    "import gsd.hoomd\n",
    "\n",
    "system_temperature = 1.0\n",
    "final_rho = 0.5\n",
    "compression_time = 5000\n",
    "relax_time = 20000\n",
    "production_time = 50000\n",
    "\n",
    "\n",
    "\n",
    "cpu = hoomd.device.CPU()\n",
    "sim = hoomd.Simulation(device=cpu, seed=1)\n",
    "\n",
    "#########################\n",
    "# system initialization\n",
    "#########################\n",
    "\n",
    "sim.create_state_from_gsd(filename='datafiles/lj.gsd')\n",
    "sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=system_temperature)\n",
    "\n",
    "\n",
    "#########################\n",
    "# interaction definition\n",
    "#########################\n",
    "\n",
    "cell = hoomd.md.nlist.Cell(buffer=1.0)\n",
    "lj = hoomd.md.pair.LJ(nlist=cell)\n",
    "lj.params[('LJ', 'LJ')] = dict(epsilon=1, sigma=1)\n",
    "lj.r_cut[('LJ', 'LJ')] = 2.5\n",
    "\n",
    "#########################\n",
    "# integrator setup\n",
    "#########################\n",
    "\n",
    "\n",
    "nvt = hoomd.md.methods.NVT(kT=system_temperature, filter=hoomd.filter.All(), tau=1.0)\n",
    "integrator = hoomd.md.Integrator(dt=0.005)\n",
    "integrator.forces.append(lj)\n",
    "integrator.methods.append(nvt)\n",
    "sim.operations.integrator = integrator\n",
    "\n",
    "#########################\n",
    "# runtime parameters\n",
    "#########################\n",
    "\n",
    "# add a logger to output the current step and timesteps per second (TPS) \n",
    "logger = hoomd.logging.Logger(categories=['scalar', 'string'])\n",
    "logger.add(sim, quantities=['timestep', 'tps'])\n",
    "table = hoomd.write.Table(trigger=hoomd.trigger.Periodic(period=5000),\n",
    "                          logger=logger)\n",
    "sim.operations.writers.append(table)\n",
    "\n",
    "# add a logger to output the thermodynamic data during the entire run \n",
    "thermodynamic_properties = hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All())\n",
    "sim.operations.computes.append(thermodynamic_properties)\n",
    "\n",
    "log_thermo_full = hoomd.logging.Logger()\n",
    "log_thermo_full.add(thermodynamic_properties)\n",
    "gsd_thermo_full = hoomd.write.GSD(filename='thermo_full.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(100),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.Null(),\n",
    "                             log=log_thermo_full)\n",
    "sim.operations.writers.append(gsd_thermo_full)\n",
    "\n",
    "\n",
    "log_traj_full = hoomd.logging.Logger()\n",
    "gsd_traj_full = hoomd.write.GSD(filename='trajectory_full.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(5000),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.All(),\n",
    "                             log=log_traj_full)\n",
    "sim.operations.writers.append(gsd_traj_full)\n",
    "\n",
    "# ================\n",
    "# compress the box\n",
    "# ================\n",
    "\n",
    "\n",
    "# ramp will define a switching factor for adjusting from the initial to the final box size\n",
    "ramp = hoomd.variant.Ramp(A=0, B=1, t_start=sim.timestep, t_ramp=compression_time)\n",
    "\n",
    "initial_box = sim.state.box\n",
    "final_box = hoomd.Box.from_box(initial_box)  # make a copy of initial_box\n",
    "final_box.volume = sim.state.N_particles / final_rho\n",
    "box_resize_trigger = hoomd.trigger.Periodic(10)\n",
    "box_resize = hoomd.update.BoxResize(box1=initial_box,\n",
    "                                    box2=final_box,\n",
    "                                    variant=ramp,\n",
    "                                    trigger=box_resize_trigger)\n",
    "sim.operations.updaters.append(box_resize)\n",
    "sim.run(compression_time+1)\n",
    "if sim.state.box == final_box:\n",
    "    print(\"system box size matches target\")\n",
    "    \n",
    "sim.operations.updaters.remove(box_resize)\n",
    "print(\"compression complete\")\n",
    "\n",
    "# ================\n",
    "# relaxation\n",
    "# ================\n",
    "\n",
    "sim.run(relax_time)\n",
    "print(\"relaxation complete\")\n",
    "\n",
    "# ================\n",
    "# production \n",
    "# ================\n",
    "\n",
    "log_thermo = hoomd.logging.Logger()\n",
    "log_thermo.add(thermodynamic_properties)\n",
    "gsd_thermo = hoomd.write.GSD(filename='thermo_prod.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(100),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.Null(),\n",
    "                             log=log_thermo)\n",
    "sim.operations.writers.append(gsd_thermo)\n",
    "\n",
    "log_traj = hoomd.logging.Logger()\n",
    "gsd_traj = hoomd.write.GSD(filename='trajectory_prod.gsd',\n",
    "                             trigger=hoomd.trigger.Periodic(5000),\n",
    "                             mode='wb',\n",
    "                             filter=hoomd.filter.All(),\n",
    "                             log=log_traj)\n",
    "sim.operations.writers.append(gsd_traj)\n",
    "\n",
    "sim.run(production_time)\n",
    "\n",
    "hoomd.write.GSD.write(state=sim.state, filename='final.gsd', mode='wb')\n",
    "print(\"simulation complete\")\n",
    "\n",
    "#delete the instances we defined ensure writers are closed\n",
    "del sim, gsd_thermo, gsd_traj, thermodynamic_properties, log_thermo, log_traj\n",
    "del log_thermo_full, gsd_thermo_full, logger, table, log_traj_full, gsd_traj_full\n",
    "del integrator, nvt, lj, cell, cpu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30751d6d-a90c-470a-87cd-a21ba348ef8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = md.load('trajectory_prod.gsd')\n",
    "r_traj_solid2, g_r_traj_solid2 = md.compute_rdf(traj,  r_range=[0, traj.unitcell_lengths[0][1]/2.0], \n",
    "                                pairs = traj.top.select_pairs(\"all\", \"all\"), \n",
    "                                  bin_width=0.01, periodic=True, opt=True)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30c0e86a-2588-43dc-afc8-eb3c3ebd0601",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(r_traj_solid1, g_r_traj_solid1, label='T=1.0, rho=1.4')\n",
    "plt.plot(r_traj_solid2, g_r_traj_solid2, label='T=0.25, rho=1.4')\n",
    "\n",
    "\n",
    "plt.ylabel('g(r)')\n",
    "plt.xlabel('r')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ea46e45-a266-4d92-8c65-99c5cc638842",
   "metadata": {},
   "outputs": [],
   "source": [
    "thermo_log = gsd.hoomd.open('thermo_prod.gsd', 'rb')\n",
    "import numpy as np\n",
    "time = []\n",
    "pe = []\n",
    "for frame in thermo_log:\n",
    "    time.append(frame.configuration.step)\n",
    "    pe.append(frame.log['md/compute/ThermodynamicQuantities/potential_energy'][0])\n",
    "\n",
    "print('mean pe: ', np.mean(np.array(pe)), ' stdev: ', np.std(np.array(pe)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed1ad6a5-a036-45b6-89ce-d725f231241e",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj_log = gsd.hoomd.open('trajectory_prod.gsd', 'rb')\n",
    "\n",
    "render(traj_log[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e60485cb-e628-4017-a3bf-4412b127329f",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}