{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "agricultural-poultry",
   "metadata": {},
   "source": [
    "# Temperature extrapolation of the radius of gyration\n",
    "\n",
    "Simulate a 5-mer simple LJ chain and compute the radius of gyration.\n",
    "\n",
    "See Section III.A of https://doi.org/10.1063/1.5026493 and https://doi.org/10.1063/1.1350578\n",
    "\n",
    "$$\\langle R_g \\rangle = \\frac{\\sum_i R_g^i \\exp(-\\beta U_i)}{\\sum_i \\exp(-\\beta U_i)}$$\n",
    "\n",
    "$$\\frac{d\\langle R_g \\rangle}{d\\beta} = -\\langle R_g U\\rangle + \\langle R_g \\rangle \\langle U \\rangle$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "photographic-small",
   "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=778555011\n",
      "# Initializing random number generator with seed: 778555011\n",
      "Configuration add_num_chain_particles=1 cubic_side_length=30 particle_type=chain:/feasst/plugin/chain/particle/chain5.txt\n",
      "Potential Model=LennardJones VisitModel=VisitModelIntra intra_cut=1\n",
      "ThermoParams beta=1 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=chain5_grow.txt\n",
      "Log output_file=chain5.txt trials_per_write=100000.0\n",
      "Movie output_file=chain5.xyz trials_per_write=100000.0\n",
      "Energy output_file=chain5en.txt trials_per_write=100000.0\n",
      "RadiusOfGyration output_file=chain5rg.txt trials_per_write=100000.0\n",
      "CheckEnergy decimal_places=6 trials_per_update=100000.0\n",
      "Run num_trials=100000.0\n",
      " \n",
      " exit: 0\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import subprocess\n",
    "import argparse\n",
    "import random\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "params = {\n",
    "    \"fstprt\": \"/feasst/plugin/chain/particle/chain5.txt\",\n",
    "    \"cubic_side_length\": 90, \"beta\": 1,\n",
    "    \"trials_per\": 1e5, \"seed\": random.randrange(int(1e9)),\n",
    "    \"production\": 1e5}\n",
    "\n",
    "with open('chain5_grow.txt', 'w') as file1:\n",
    "    file1.write(\"\"\"\n",
    "TrialGrowFile\n",
    "\n",
    "bond=true mobile_site=1 anchor_site=0 particle_type=chain weight=1\n",
    "bond=true mobile_site=2 anchor_site=1\n",
    "bond=true mobile_site=3 anchor_site=2\n",
    "bond=true mobile_site=4 anchor_site=3\"\"\")\n",
    "\n",
    "# write fst script to run a single simulation\n",
    "with open('script4.txt', \"w\") as myfile: myfile.write(\"\"\"\n",
    "MonteCarlo\n",
    "RandomMT19937 seed={seed}\n",
    "Configuration cubic_side_length=30 particle_type=chain:{fstprt} add_num_chain_particles=1 \n",
    "Potential Model=LennardJones VisitModel=VisitModelIntra intra_cut=1\n",
    "ThermoParams beta={beta} chemical_potential=1\n",
    "Metropolis\n",
    "TrialGrowFile grow_file=chain5_grow.txt\n",
    "Let [write]=trials_per_write={trials_per} output_file=chain5\n",
    "Log [write].txt\n",
    "Movie [write].xyz\n",
    "Energy [write]en.txt\n",
    "RadiusOfGyration [write]rg.txt\n",
    "CheckEnergy trials_per_update={trials_per} decimal_places=6\n",
    "Run num_trials={production}\n",
    "\"\"\".format(**params))\n",
    "\n",
    "import subprocess\n",
    "syscode = subprocess.call(\"feasst < script4.txt > script4.log\", shell=True, executable='/bin/bash')\n",
    "with open('script4.log', 'r') as file: print(file.read(), '\\n', 'exit:', syscode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "sublime-extent",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "drg_dbeta -0.09516386530560306\n"
     ]
    }
   ],
   "source": [
    "def rg_fit(beta, beta1, rg1, der1):\n",
    "    return rg1 + (beta - beta1)*der1\n",
    "\n",
    "taylor2001 = pd.read_csv('../test/data/taylor2001rg.csv', header=None)\n",
    "#print(taylor2001)\n",
    "rgdf = pd.read_csv('chain5rg.txt')\n",
    "endf = pd.read_csv('chain5en.txt')\n",
    "rg = rgdf['average'].values[0]\n",
    "en = endf['average'].values[0]\n",
    "assert np.abs(rg - np.sqrt(0.88)) < 0.01\n",
    "assert np.abs(en + 2.04) < 0.06\n",
    "rg_en = rgdf['rgu'].values[0]\n",
    "drg_dbeta = rg*en - rg_en\n",
    "assert np.abs(drg_dbeta + 0.186/2) < 0.01\n",
    "print('drg_dbeta', drg_dbeta)\n",
    "\n",
    "#         import matplotlib.pyplot as plt\n",
    "#         plt.scatter(params['beta'], rg, label='explicit sim')\n",
    "#         plt.plot(1./taylor2001[0], np.sqrt(taylor2001[1]), label='taylor2001')\n",
    "#         plt.plot(1./taylor2001[0],\n",
    "#                  rg_fit(beta=1./taylor2001[0],\n",
    "#                         beta1=params['beta'],\n",
    "#                         rg1=rg,\n",
    "#                         der1=drg_dbeta),\n",
    "#                  label='1st order extrapolation')\n",
    "#         plt.xscale('log')\n",
    "#         plt.ylim([0.2, 1.3])\n",
    "#         plt.legend()\n",
    "#         plt.xlabel(r'$\\beta$', fontsize=16)\n",
    "#         plt.ylabel(r'$R_g$', fontsize=16)\n",
    "#         plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "postal-valentine",
   "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": 5
}
