{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "unknown-specific",
   "metadata": {},
   "source": [
    "# Lennard Jones potential test\n",
    "\n",
    "First, the potential energy between two particles is compared to the analytical value.\n",
    "Next, the energy of an entire configuration is compared against the [NIST SRSW](https://www.nist.gov/programs-projects/nist-standard-reference-simulation-website) reference calculations.\n",
    "\n",
    "The energy of a reference configuration is obtained with many decimal places using the [Log](../../steppers/doc/Log_arguments.rst)::`max_precision` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "split-sculpture",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Usage: /home/user/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.13\n",
      "MonteCarlo\n",
      "Configuration particle_type=lj:/feasst/particle/lj_new.txt xyz_file=two.xyz\n",
      "Potential Model=LennardJones\n",
      "Potential VisitModel=LongRangeCorrections\n",
      "ThermoParams beta=1000000\n",
      "Metropolis\n",
      "Log clear_file=true max_precision=true output_file=lj.csv\n",
      "Run num_trials=1\n",
      "#Warn 0 [plugin/monte_carlo/src/monte_carlo.cpp:786] No Trials to attempt.\n",
      "# Usage: /home/user/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.13\n",
      "MonteCarlo\n",
      "Configuration cubic_side_length=8 particle_type=lj:/feasst/particle/lj_new.txt xyz_file=/home/user/feasst/plugin/configuration/test/data/lj_sample_config_periodic4.xyz\n",
      "Potential Model=LennardJones\n",
      "Potential VisitModel=LongRangeCorrections\n",
      "ThermoParams beta=1000000\n",
      "Metropolis\n",
      "Log clear_file=true max_precision=true output_file=lj.csv\n",
      "Run num_trials=1\n",
      "#Warn 0 [plugin/monte_carlo/src/monte_carlo.cpp:786] No Trials to attempt.\n",
      "# Usage: /home/user/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.13\n",
      "MonteCarlo\n",
      "Configuration particle_type=lj:/feasst/particle/lj_new.txt side_length=10.0,9.84807753012208,9.64974312607518 xy=1.7364817766693041 xyz_file=/home/user/feasst/plugin/configuration/test/data/lj_triclinic_sample_config_periodic3.xyz xz=2.5881904510252074 yz=0.42863479791864567\n",
      "Potential Model=LennardJones\n",
      "Potential VisitModel=LongRangeCorrections\n",
      "ThermoParams beta=1000000\n",
      "Metropolis\n",
      "Log clear_file=true max_precision=true output_file=lj.csv\n",
      "Run num_trials=1\n",
      "#Warn 0 [plugin/monte_carlo/src/monte_carlo.cpp:786] No Trials to attempt.\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import sys\n",
    "import subprocess\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# Run a feasst simulation with the given parameters\n",
    "def run_fst(params):\n",
    "    with open(\"launch0.txt\", \"w\") as myfile: myfile.write(\"\"\"\n",
    "MonteCarlo\n",
    "Configuration {config_params}\n",
    "Potential Model=LennardJones\n",
    "Potential VisitModel=LongRangeCorrections\n",
    "ThermoParams beta=1000000\n",
    "Metropolis\n",
    "Log output_file=lj.csv max_precision=true clear_file=true\n",
    "Run num_trials=1\n",
    "\"\"\".format(**params))\n",
    "    syscode = subprocess.call(\"../../../build/bin/fst < launch0.txt\", shell=True, executable='/bin/bash')\n",
    "    if syscode > 0: sys.exit(1)\n",
    "\n",
    "\"\"\"Test the LJ potential against analytical calculation of two particles\"\"\"\n",
    "params = {\"displacement\": 1.2345}\n",
    "with open(\"two.xyz\", \"w\") as myfile: myfile.write(\n",
    "\"\"\"2\n",
    "-1 8 8 8\n",
    "0 0 0 0\n",
    "1 0 0 {displacement}\"\"\".format(**params))\n",
    "run_fst({\"config_params\": \"particle_type=lj:/feasst/particle/lj_new.txt xyz_file=two.xyz\"})\n",
    "df = pd.read_csv('lj.csv')\n",
    "assert 2 == df['num_particles_lj'][0]\n",
    "\n",
    "# compute the expected analytical LJ and LRC energies\n",
    "enlj = 4*(params[\"displacement\"]**(-12) - params[\"displacement\"]**(-6))\n",
    "rcut = 3 #mc.system().configuration().model_params().select(\"cutoff\").value(0)\n",
    "enlrc = (8./3.)*math.pi*2**2/ 8**3*((1./3.)*rcut**(-9) - rcut**(-3))\n",
    "\n",
    "# Compare the analytical results with the FEASST computed energies.\n",
    "assert np.abs(enlj - df['LennardJones'][0]) < 12\n",
    "assert np.abs(enlrc - df['LongRangeCorrections'][0]) < 12\n",
    "assert np.abs(enlj + enlrc - df['energy'][0]) < 12\n",
    "\n",
    "\"\"\"Test the LJ potential against a configuration of 30 particles.\n",
    "In particular, the 4th configuration of the LJ SRSW reference:\n",
    "https://www.nist.gov/mml/csd/chemical-informatics-research-group/lennard-jones-fluid-reference-calculations\n",
    "\"\"\"\n",
    "run_fst({\"config_params\": \"cubic_side_length=8 particle_type=lj:/feasst/particle/lj_new.txt \\\n",
    "          xyz_file=/feasst/plugin/configuration/test/data/lj_sample_config_periodic4.xyz\"})\n",
    "df = pd.read_csv('lj.csv')\n",
    "assert df['num_particles_lj'][0] == 30\n",
    "enlj = -16.790321304625856\n",
    "enlrc = -0.5451660014945704\n",
    "assert np.abs(enlj - df['LennardJones'][0]) < 12\n",
    "assert np.abs(enlrc - df['LongRangeCorrections'][0]) < 12\n",
    "assert np.abs(enlj + enlrc - df['energy'][0]) < 12\n",
    "\n",
    "\"\"\"Test the LJ potential against a configuration of 300 particles in a trinclinic cell.\n",
    "In particular, the 3th configuration of the triclinic LJ SRSW reference:\n",
    "https://www.nist.gov/mml/csd/chemical-informatics-group/lennard-jones-fluid-reference-calculations-non-cuboid-cell\n",
    "\"\"\"\n",
    "run_fst({\"config_params\": \"side_length=10.0,9.84807753012208,9.64974312607518 \\\n",
    "    xy=1.7364817766693041 xz=2.5881904510252074 yz=0.42863479791864567 \\\n",
    "    particle_type=lj:/feasst/particle/lj_new.txt \\\n",
    "    xyz_file=/feasst/plugin/configuration/test/data/lj_triclinic_sample_config_periodic3.xyz\"})\n",
    "df = pd.read_csv('lj.csv')\n",
    "assert df['num_particles_lj'][0] == 300\n",
    "enlj = -505.78567945268367\n",
    "enlrc = -29.37186430697248\n",
    "assert np.abs(enlj - df['LennardJones'][0]) < 12\n",
    "assert np.abs(enlrc - df['LongRangeCorrections'][0]) < 12\n",
    "assert np.abs(enlj + enlrc - df['energy'][0]) < 12"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "numeric-tattoo",
   "metadata": {},
   "source": [
    "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
