{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lennard Jones in the canonical ensemble\n",
    "\n",
    "In this example, the [average energy reported in the NIST SRSW](https://mmlapps.nist.gov/srs/LJ_PURE/mc.htm) at $T^*=0.9$ and $\\rho^*=0.003$ is computed with Monte Carlo simulation of a bulk Lennard Jones fluid ([Section IIIB of the manuscript](https://doi.org/10.1063/5.0224283)).\n",
    "The simulation initializes a desired number of particles in cubic periodic boundaries, tunes the maximum trial displacement during equilibration, and finally obtains the ensemble average energy.\n",
    "\n",
    "We use Python to generate a string-based input to FEASST for the following conveniences:\n",
    "\n",
    " - Variables required by FEASST arguments may need to be converted from variables that the user prefers to input, and these FEASST arguments can be generated by Python to machine precision. For example, if the user prefers to input temperature and state a temperature of $T^*=0.9$ in their manuscript, but FEASST requires $\\beta=\\frac{1}{k_B T}$, then $T^*=0.9$ is easier to input than $\\beta=1.111111111111111$.\n",
    " - Python generation of input files may create multiple related simulations from a single source. For example, the next tutorial simulates a range of densities with a single Python file.\n",
    " - Python provides formatted strings and argument parsing.\n",
    "\n",
    "In this example, we begin a simulation by defining a string variable `fstin`.\n",
    "Optional comments can be included for future reference.\n",
    "First, a [MonteCarlo](../plugin/monte_carlo/doc/MonteCarlo_arguments.rst) object is created.\n",
    "Then, the random number generator is given a seed.\n",
    "The available arguments of [RandomMT19937](../plugin/math/doc/RandomMT19937_arguments.rst) can be found in the documentation of that object, or its base-class, [Random](../plugin/math/doc/Random_arguments.rst). \n",
    "Arguments of objects are input as pairs of argument names then values separated by an equal sign (e.g., \"name=value\").\n",
    "Each argument pair is separated by a space.\n",
    "In this case, `seed` is the argument name and `1572362164` is the argument value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Usage: feasst < file.txt\n",
      "# For more information, use the command \"feasst-menu\"\n",
      "# Exit with ctrl-c\n",
      "FEASST version 0.25.19\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=1572362164\n",
      "# Initializing random number generator with seed: 1572362164\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# This and all tutorial code blocks are in Python and included in your download of FEASST.\n",
    "# For example, https://pages.nist.gov/feasst/tutorial/tutorial.html is in $HOME/feasst/tutorial/tutorial.ipynb . Usage: jupyter notebook tutorial.ipynb\n",
    "# And https://pages.nist.gov/feasst/tutorial/launch.html is in $HOME/feasst/tutorial/launch.py . Usage: python launch.py\n",
    "\n",
    "import subprocess\n",
    "\n",
    "fstin=\"\"\"\n",
    "# comments begin with the '#' symbol and are not displayed in the output\n",
    "MonteCarlo\n",
    "RandomMT19937 seed=1572362164\n",
    "\"\"\"\n",
    "\n",
    "def run(fstin):\n",
    "    with open('fstin.txt', 'w') as file: file.write(fstin)\n",
    "    syscode = subprocess.call(\"feasst < fstin.txt > fstout.txt\", shell=True, executable='/bin/bash')\n",
    "    with open('fstout.txt', 'r') as file:\n",
    "        print(file.read())\n",
    "run(fstin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output FEASST may be used to reproduce the simulation.\n",
    "\n",
    "If there is a typo in the arguments, this may result in an exception which will print to the terminal. You may test this now by changing the `seed` argument to `sd`. Because the argument `sd` is not recognized by [RandomMT19937](../plugin/math/doc/RandomMT19937_arguments.rst), a nonzero system code should be returned and a verbose error should appear in your terminal.\n",
    "\n",
    "The next step is to add a [Configuration](../plugin/configuration/doc/Configuration_arguments.rst).\n",
    "In this example, a simple cubic periodic box is defined based on the number of particles and density.\n",
    "Because the cubic length needs to be input to high precision, a Python format string is used to input a variable enclosed in curly brackets.\n",
    "\n",
    "In addition, a particle type \"lj\" is defined by the file `lj.txt`.\n",
    "Any argument value beginning with \"/feasst\" will have that beginning replaced with the [installation directory of FEASST](../plugin/utils/doc/definitions.rst)).\n",
    "See [Particle](../particle/README.rst) for more information about the format of the particle file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Usage: feasst < file.txt\n",
      "# For more information, use the command \"feasst-menu\"\n",
      "# Exit with ctrl-c\n",
      "FEASST version 0.25.19\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=1572362164\n",
      "# Initializing random number generator with seed: 1572362164\n",
      "Configuration cubic_side_length=55.03212081491043 particle_type=lj:/feasst/particle/lj_new.txt\n",
      "\n"
     ]
    }
   ],
   "source": [
    "num_particles=500\n",
    "density=0.003\n",
    "length=(num_particles/density)**(1./3.)\n",
    "fstin+=f\"\"\"\n",
    "Configuration cubic_side_length={length} particle_type=lj:/feasst/particle/lj_new.txt\"\"\"\n",
    "run(fstin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two [Potentials](../plugin/system/doc/Potential_arguments.rst) are introduced starting with the pair-wise [Lennard-Jones (LJ) model](../plugin/system/doc/LennardJones_arguments.rst).\n",
    "Because the domain is large compared to the potential cutoff, a [cell list](../plugin/system/doc/VisitModelCell_arguments.rst) speeds up the energy calculation.\n",
    "Then, [long-range corrections](../plugin/system/doc/LongRangeCorrections_arguments.rst) approximately account for the cut off of the LJ potential by assuming a pair-wise radial distribution function of unity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fstin+=\"\"\"\n",
    "Potential Model=LennardJones VisitModel=VisitModelCell min_length=max_cutoff\n",
    "Potential VisitModel=LongRangeCorrections\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[ThermoParams](../plugin/system/doc/ThermoParams_arguments.rst) set temperature and the chemical potential of each particle type.\n",
    "A chemical potential is required when the initial configuration is generated with grand canonical particle additions, but the chemical potential will not contribute to canonical ensemble production simulations. [Metropolis](../plugin/monte_carlo/doc/Metropolis_arguments.rst) acceptance criteria is also utilized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Usage: feasst < file.txt\n",
      "# For more information, use the command \"feasst-menu\"\n",
      "# Exit with ctrl-c\n",
      "FEASST version 0.25.19\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=1572362164\n",
      "# Initializing random number generator with seed: 1572362164\n",
      "Configuration cubic_side_length=55.03212081491043 particle_type=lj:/feasst/particle/lj_new.txt\n",
      "Potential Model=LennardJones VisitModel=VisitModelCell min_length=max_cutoff\n",
      "Potential VisitModel=LongRangeCorrections\n",
      "ThermoParams beta=1.1111111111111112 chemical_potential=-1\n",
      "Metropolis\n",
      "\n"
     ]
    }
   ],
   "source": [
    "fstin+=f\"\"\"\n",
    "ThermoParams beta={1./0.9} chemical_potential=-1\n",
    "Metropolis\"\"\"\n",
    "run(fstin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A [TrialTranslate](../plugin/monte_carlo/doc/TrialTranslate_arguments.rst) attempts translation of a random particle by a random distance bound in each dimension by a `tunable_param`.\n",
    "The desired acceptance ratio, `tunable_target_acceptance`, is obtained by adjusting this parameter with [Tune](../plugin/steppers/doc/Tune_arguments.rst)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "fstin+=\"\"\"\n",
    "TrialTranslate weight=1 tunable_param=2 tunable_target_acceptance=0.2\n",
    "Tune\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[CheckEnergy](../plugin/steppers/doc/CheckEnergy_arguments.rst) asserts that the optimized energy calculations match the unoptimized calculations within a tolerance.\n",
    "This check is only performed every `trials_per_update` because the unoptimized energy calculation is expensive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "fstin+=\"\"\"\n",
    "CheckEnergy trials_per_update=1e4 decimal_places=8\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[TrialAdd](../plugin/monte_carlo/doc/TrialTranslate_arguments.rst), generates an initial configuration with a [Run](../plugin/actions/doc/Run_arguments.rst) with the desired number of particles. Afterward, we [Remove](../plugin/actions/doc/Remove_arguments.rst) [TrialAdd](../plugin/monte_carlo/doc/TrialTranslate_arguments.rst) to simulate the canonical ensemble."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "fstin+=f\"\"\"\n",
    "Run until_num_particles={num_particles} Trial=TrialAdd weight=2 particle_type=lj\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation [Run](../plugin/actions/doc/Run_arguments.rst) for a number of equilibration trials before [Remove](../plugin/actions/doc/Remove_arguments.rst) disables tuning of the maximum displacement ([to satisfy detailed balance](https://github.com/usnistgov/best-practices-mc/))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "fstin+=\"\"\"\n",
    "Run num_trials=1e5\n",
    "Remove name=Tune\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, initialize a production simulation.\n",
    "First, a [Log](../plugin/steppers/doc/Log_arguments.rst) periodically outputs the instantaneous properties of the simulation.\n",
    "[Log](../plugin/steppers/doc/Log_arguments.rst) is useful for monitoring simulation progress, but is not recommended for calculations.\n",
    "Instead, [Energy](../plugin/steppers/doc/Energy_arguments.rst) computes the average energy by accumulating its value after every trial (more often than once every line in the [Log](../plugin/steppers/doc/Log_arguments.rst) file).\n",
    "[Movie](../plugin/steppers/doc/Movie_arguments.rst) produces a trajectory in the [XYZ file format](../plugin/configuration/doc/PrinterXYZ_arguments.rst).\n",
    "Additional [Analyze](../plugin/monte_carlo/doc/Analyze_arguments.rst) or [Modify](../plugin/monte_carlo/doc/Modify_arguments.rst) may be added to perform some task contingent upon the number of attempted trials.\n",
    "\n",
    "Because all three of these [Analyze](../plugin/monte_carlo/doc/Analyze_arguments.rst) share similar arguments, [Let](../plugin/utils/doc/Let_arguments.rst) minimizes redundancy by defining a string of characters as the variable \\[write\\]. Variables must begin and end with square brackets followed by an equal sign to set its value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "write='trials_per_write=1e4 output_file=lj'\n",
    "fstin+=f\"\"\"\n",
    "Log {write}.csv\n",
    "Energy {write}_en.csv\n",
    "Movie {write}.xyz\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next code block is equivalent to the previous, but demonstrates the use of [For](../plugin/utils/doc/For_arguments.rst) loops. Each [For](../plugin/utils/doc/For_arguments.rst) executes the following lines with given string substitutions until [EndFor](../plugin/utils/doc/EndFor_arguments.rst), and then repeats those lines again for each comma-separated list of arguments. Within the same [For](../plugin/utils/doc/For_arguments.rst) loop, multiple string substitutions are defined with colon-separated values. For example, \\[Analyze\\]=Log and \\[extension\\]=.csv in the first iteraction of the loop, as can be seen in the output. [If](../plugin/utils/doc/If_arguments.rst), [Else](../plugin/utils/doc/Else_arguments.rst) and [EndIf](../plugin/utils/doc/EndIf_arguments.rst) may also be used for optional arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Usage: feasst < file.txt\n",
      "# For more information, use the command \"feasst-menu\"\n",
      "# Exit with ctrl-c\n",
      "FEASST version 0.25.19\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=1572362164\n",
      "# Initializing random number generator with seed: 1572362164\n",
      "Configuration cubic_side_length=55.03212081491043 particle_type=lj:/feasst/particle/lj_new.txt\n",
      "Potential Model=LennardJones VisitModel=VisitModelCell min_length=max_cutoff\n",
      "Potential VisitModel=LongRangeCorrections\n",
      "ThermoParams beta=1.1111111111111112 chemical_potential=-1\n",
      "Metropolis\n",
      "TrialTranslate tunable_param=2 tunable_target_acceptance=0.2 weight=1\n",
      "Tune\n",
      "CheckEnergy decimal_places=8 trials_per_update=1e4\n",
      "Run Trial=TrialAdd particle_type=lj until_num_particles=500 weight=2\n",
      "Run num_trials=1e5\n",
      "Remove name=Tune\n",
      "Log output_file=lj.csv trials_per_write=1e4\n",
      "Energy output_file=lj_en.csv trials_per_write=1e4\n",
      "Movie output_file=lj.xyz trials_per_write=1e4\n",
      "Log output_file=lj_with_for.csv trials_per_write=1e4\n",
      "Energy output_file=lj_with_for_en.csv trials_per_write=1e4\n",
      "Movie output_file=lj_with_for.xyz trials_per_write=1e4\n",
      "Run num_trials=1e6\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for pair in [['Log', '.csv'], ['Energy', '_en.csv'], ['Movie', '.xyz']]:\n",
    "    fstin+=f\"\\n{pair[0]} trials_per_write=1e4 output_file=lj_with_for{pair[1]}\"\n",
    "fstin+=\"\\nRun num_trials=1e6\"\n",
    "run(fstin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the [average energy to the NIST SRSW](https://mmlapps.nist.gov/srs/LJ_PURE/mc.htm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U_FEASST/N= -0.030332409560752716 +/- 0.0002319803861581972\n",
      "U_SRSW/N= -0.029787 +/- 3.21e-05\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "df=pd.read_csv('lj_en.csv')\n",
    "print('U_FEASST/N=', df['average'][0]/num_particles, '+/-', df['block_stdev'][0]/num_particles)\n",
    "print('U_SRSW/N=', -2.9787E-02, '+/-', 3.21E-05)\n",
    "import math\n",
    "assert abs(-2.9787E-02-df['average'][0]/num_particles) < 3*math.sqrt(3.21E-05**2+(df['block_stdev'][0]/num_particles)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should also find the `lj.xyz` trajectory file with an automatically-generated `lj.xyz.vmd` file for use with VMD (e.g., `vmd -e lj.xyz.vmd`).\n",
    "\n",
    "Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please [contact](../CONTACT.rst) us. Any feedback is appreciated!"
   ]
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
