{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reference configuration of SPC/E water in non-cuboid domain\n",
    "\n",
    "Compare the potential energy against the [NIST SRSW](https://www.nist.gov/mml/csd/chemical-informatics-group/spce-water-reference-calculations-non-cuboid-cell-10a-cutoff) values. This tutorial requires the use of [Ewald](../doc/Ewald_arguments.rst)."
   ]
  },
  {
   "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",
      "Configuration particle_type=spce:/feasst/particle/spce_new.txt side_length=30.0,28.97777478867205,29.51512917398008 xy=7.764571353075622 xyz_file=../../../plugin/charge/test/data/spce_triclinic_sample_periodic1.xyz xz=-2.6146722824297473 yz=-4.692615336756641\n",
      "Potential Model=LennardJones\n",
      "Potential VisitModel=LongRangeCorrections\n",
      "Potential VisitModel=Ewald alpha=0.2850 kxmax=7 kymax=7 kzmax=7\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 2.3037\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 2.3037\n",
      "Potential Model=ChargeScreened\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 2.3037\n",
      "Potential Model=ChargeScreenedIntra VisitModel=VisitModelBond\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 2.3037\n",
      "Potential Model=ChargeSelf\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 2.3037\n",
      "ThermoParams beta=1000000\n",
      "Metropolis\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 2.3037\n",
      "Log clear_file=true max_precision=true output_file=en.csv\n",
      "Run num_trials=1\n",
      "#Warn 0 [:802] No Trials to attempt.\n",
      " \n",
      " exit: 0\n"
     ]
    }
   ],
   "source": [
    "script=\"\"\"\n",
    "MonteCarlo\n",
    "Configuration side_length=30.0,28.97777478867205,29.51512917398008 \\\n",
    "    xy=7.764571353075622 xz=-2.6146722824297473 yz=-4.692615336756641 \\\n",
    "    xyz_file=../../../plugin/charge/test/data/spce_triclinic_sample_periodic1.xyz \\\n",
    "    particle_type=spce:/feasst/particle/spce_new.txt\n",
    "Potential Model=LennardJones\n",
    "Potential VisitModel=LongRangeCorrections\n",
    "Potential VisitModel=Ewald alpha=0.2850 kxmax=7 kymax=7 kzmax=7\n",
    "Potential Model=ChargeScreened\n",
    "Potential Model=ChargeScreenedIntra VisitModel=VisitModelBond\n",
    "Potential Model=ChargeSelf\n",
    "ThermoParams beta=1000000\n",
    "Metropolis\n",
    "Log output_file=en.csv max_precision=true clear_file=true\n",
    "Run num_trials=1\n",
    "\"\"\"\n",
    "\n",
    "with open('script.txt', 'w') as file: file.write(script)\n",
    "import subprocess\n",
    "syscode = subprocess.call(\"feasst < script.txt > script.log\", shell=True, executable='/bin/bash')\n",
    "with open('script.log', 'r') as file: print(file.read(), '\\n', 'exit:', syscode)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the test passes, the energy is within the tolerance of the SRSW value and the two ensemble average methods agreed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "df=pd.read_csv('en.csv')\n",
    "assert np.abs(df['LennardJones'][0] - 931.15451) < 1e-4\n",
    "assert np.abs(df['LongRangeCorrections'][0] + 34.16569) < 1e-4\n",
    "assert np.abs(df['Ewald'][0] - 371.46525) < 1e-4\n",
    "assert np.abs(df['ChargeScreened'][0] + 6046.43627) < 1e-4\n",
    "assert np.abs(df['ChargeScreenedIntra'][0] - 95078.89447) < 1e-4\n",
    "assert np.abs(df['ChargeSelf'][0] + 96297.75579) < 1e-4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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",
      "Configuration particle_type=spce:/feasst/particle/spce_new.txt side_length=36.0,36.0,31.17691453623979 xy=0.0 xyz_file=../../../plugin/charge/test/data/spce_monoclinic_sample_periodic4.xyz xz=18.0 yz=0.0\n",
      "Potential Model=LennardJones\n",
      "Potential VisitModel=LongRangeCorrections\n",
      "Potential VisitModel=Ewald alpha=0.2850 kxmax=7 kymax=7 kzmax=7\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 1.99017\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 1.99017\n",
      "Potential Model=ChargeScreened\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 1.99017\n",
      "Potential Model=ChargeScreenedIntra VisitModel=VisitModelBond\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 1.99017\n",
      "Potential Model=ChargeSelf\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 1.99017\n",
      "ThermoParams beta=1000000\n",
      "Metropolis\n",
      "# Ewald alpha: 0.285\n",
      "# Ewald kmax_squared 1.99017\n",
      "Log clear_file=true max_precision=true output_file=en.csv\n",
      "Run num_trials=1\n",
      "#Warn 0 [:802] No Trials to attempt.\n",
      " \n",
      " exit: 0\n"
     ]
    }
   ],
   "source": [
    "script=\"\"\"\n",
    "MonteCarlo\n",
    "Configuration side_length=36.0,36.0,31.17691453623979 \\\n",
    "    xy=0.0 xz=18.0 yz=0.0 \\\n",
    "    xyz_file=../../../plugin/charge/test/data/spce_monoclinic_sample_periodic4.xyz \\\n",
    "    particle_type=spce:/feasst/particle/spce_new.txt\n",
    "Potential Model=LennardJones\n",
    "Potential VisitModel=LongRangeCorrections\n",
    "Potential VisitModel=Ewald alpha=0.2850 kxmax=7 kymax=7 kzmax=7\n",
    "Potential Model=ChargeScreened\n",
    "Potential Model=ChargeScreenedIntra VisitModel VisitModelBond\n",
    "Potential Model=ChargeSelf\n",
    "ThermoParams beta=1000000\n",
    "Metropolis\n",
    "Log output_file=en.csv max_precision=true clear_file=true\n",
    "Run num_trials=1\n",
    "\"\"\"\n",
    "\n",
    "with open('script.txt', 'w') as file: file.write(script)\n",
    "import subprocess\n",
    "syscode = subprocess.call(\"feasst < script.txt > script.log\", shell=True, executable='/bin/bash')\n",
    "with open('script.log', 'r') as file: print(file.read(), '\\n', 'exit:', syscode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "df=pd.read_csv('en.csv')\n",
    "assert np.abs(df['LennardJones'][0] - 208.07025846554663) < 1e-4\n",
    "assert np.abs(df['LongRangeCorrections'][0] + 1.35601402285723) < 1e-4\n",
    "# Note that the Ewald value is slightly different than SRSW due to change in kmax\n",
    "assert np.abs(df['Ewald'][0] - 185.60955025935957) < 1e-4\n",
    "assert np.abs(df['ChargeScreened'][0] + 1425.6143894337906) < 1e-2\n",
    "assert np.abs(df['ChargeScreenedIntra'][0] - 23769.71831128121) < 1e-2\n",
    "assert np.abs(df['ChargeSelf'][0] + 24074.443366236523) < 1e-2"
   ]
  },
  {
   "cell_type": "markdown",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
