Session4-Monte_Carlo/Monte_Carlo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo Simulation\n",
    "\n",
    "For today's session, we will no longer be performing simulations using molecular dynamics. Instead, we will be using the second of the \"Big 2\" simulation approaches: Monte Carlo.\n",
    "\n",
    "Monte Carlo (MC) offers several advantages and disadvantages over molecular dynamics (MD), and it will depend on the system and problem of interest as to which simulation approach is appropriate. For systems where large energy barriers can result in systems that are easily trapped in local minima, MC can provide a significant advantage over MD. However, typically there is no \"time\" in a Monte Carlo simulation, thus properties that require time-averaging can no be assessed via MC.\n",
    "\n",
    "In this session, you will introduced to the basics of MC through simulations of a simple Lennard-Jones system (comparable to the system you performed MD simulations of on Tuesday). While several open-source, highly-optimized MD codes exist (you have already been introduced to HOOMD and GROMACS), relatively few MC codes exist. As a result, for this session we will be using an ad hoc MC code written in C++, that has been wrapped in Python for clarity. We will explore some of the underlying code and observe how changing various run parameters effect the result of the simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing the MC Code\n",
    "\n",
    "The MC code is contained in several files located in the `src` directory. Modules exist for the LJ force field, the LJ system, and the Monte Carlo algorithm. We will import these three classes here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src import ForceField, LJBox, MonteCarlo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the force field\n",
    "\n",
    "Let's define the Lennard-Jones parameters that will dictate how our particles will interact with one another. In the box below we will define values for `sigma`, `epsilon`, and `cutoff`. Recall these parameters from the Lennard-Jones equation:\n",
    "\n",
    "\\begin{align}\n",
    "U_{LJ} = 4 \\epsilon \\left [ {\\left ( \\frac{\\sigma}{r} \\right )} ^{12} - {\\left ( \\frac{\\sigma}{r} \\right )} ^{6} \\right ],\\: r \\le r_{cut}\n",
    "\\end{align}\n",
    "\n",
    "We will begin by using a sigma of 1.0, an epsilon of 1.0, and a cutoff of 2.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 1.0 \n",
    "epsilon = 1.0 \n",
    "cutoff = 2.5 \n",
    "\n",
    "forcefield = ForceField(sigma=sigma, epsilon=epsilon, cutoff=cutoff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this code, force field parameters are stored within a `ForceField` object. Let's take a quick look at this (simple) class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat src/forcefield.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see above, the `ForceField` class simply features a constructor to load force field information into several class attributes and a `show` method to plot the force field.\n",
    "\n",
    "Let's use the `show` method here to obtain a better appreciation for the force field we've defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forcefield.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define system parameters\n",
    "\n",
    "Now we want to define our system parameters. For this system these parameters will be the total number of particles and the number density of particles to place in the box.\n",
    "\n",
    "We will begin by defining a system of 216 particles and a number density of 0.7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_particles = 216 \n",
    "number_density = 0.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `LJBox` class will take these system parameters and use mBuild to create the system. Let's take a quick look at that class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat src/ljbox.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = LJBox(n_particles, number_density)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define Monte Carlo parameters\n",
    "\n",
    "Now let's define the parameters for our Monte Carlo simulation. There are several parameters that we need to define:\n",
    "  - Temperature &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; Recall $\\beta = 1 / {k_B T}$\n",
    "  - $\\Delta_{x, init}$ &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp;&nbsp;&nbsp; Initial maximum particle displacement\n",
    "  - Acceptance probability &nbsp; &nbsp; &nbsp; At what probability do we want move to be accepted?\n",
    "  - $n_{relax}$ &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Number of steps to relax from the initial configuration\n",
    "  - $n_{MC}$ &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp; Number of MC steps to perform after the relaxation period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = 1.2\n",
    "dx = 0.1            # Initial maximum displacement\n",
    "target = 0.5        # Target acceptance probabality\n",
    "n_relax = 10000     # Number of steps to relax from initial configuration\n",
    "n_mc = 100000       # Total number of MC steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will feed these parameters into the constructor for a `MonteCarlo` object. This class provides a wrapper around the C++ code that performs the Monte Carlo simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mc = MonteCarlo(forcefield=forcefield,\n",
    "                system=system,\n",
    "                dx=dx,\n",
    "                temperature=temperature,\n",
    "                target=target)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Relax the system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will perform a relaxation, both to move our system out of a high energy initial state, and to determine the optimum value of `dx` to yield the target acceptance probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "mc.relax(n_relax, adjust_freq=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's view the trajectory of the relaxation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview as nv\n",
    "import mdtraj as md\n",
    "traj = md.load('relax.xyz', top='system_init.gro')\n",
    "nv.show_mdtraj(traj, representations=[{'type': 'spacefill', 'params': {'radius': 0.5}}])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Perform production Monte Carlo run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have optimized `dx` we can perform a production run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mc.run(n_mc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### View the trajectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview as nv\n",
    "import mdtraj as md\n",
    "traj = md.load('run.xyz', top='system_init.gro')\n",
    "nv.show_mdtraj(traj, representations=[{'type': 'spacefill', 'params': {'radius': 0.5}}])"
   ]
  }
 ],
 "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}