{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "agricultural-poultry",
   "metadata": {},
   "source": [
    "# Testing angles and branches\n",
    "\n",
    "Following https://doi.org/10.1021/acs.jctc.7b00173 , compute the angle distributions.\n",
    "\n",
    "First, consider propane and 2-methylpropane. Also consider very flexible molecules as follows:\n",
    "\n",
    "floppy_angle: 0 - 1 - 2\n",
    "\n",
    "floppy_dihedral: 0 - 1 - 2 - 3\n",
    "\n",
    "floppy_branch:\n",
    "\n",
    "      1\n",
    "      |\n",
    "      0\n",
    "     /  \\\n",
    "    2    3\n",
    "\n",
    "Where the equilibrium angles are all 60 degrees but the harmonic angle potential is not strong, U=k/2(theta-theta0)^2. The bonds are rigid with a length of 1 and the diameters of the beads are 1.\n",
    "\n",
    "The floppy particles with no interactions should have sin(theta) angle distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ba226ff7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import subprocess\n",
    "import argparse\n",
    "import random\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from pyfeasst import physical_constants\n",
    "\n",
    "params = {\n",
    "    \"cubic_side_length\": 90,\n",
    "    \"temperature_K\": 300,\n",
    "    \"trials_per\": 1e5,\n",
    "    \"production\": 1e5,\n",
    "    \"k\": 1,\n",
    "    'angle_bin_width': 0.01,\n",
    "    'angle_bin_center': 100/180*np.pi,\n",
    "    'molecule': 'propane',\n",
    "    'num_jacobian_gaussian': 0,\n",
    "    #'molecule': '2-methylpropane',\n",
    "    #'molecule': 'floppy_branch',\n",
    "    'intra': 1,\n",
    "    'show_plot': False,\n",
    "}\n",
    "\n",
    "def update_params(params):\n",
    "    params['seed'] = random.randrange(int(1e9))\n",
    "    params['beta'] = 1./(params['temperature_K']*physical_constants.MolarGasConstant().value()/1e3) # mol/kJ\n",
    "    if params['molecule'] == 'propane':\n",
    "        params['config'] = \"Configuration cubic_side_length=30 particle_type=pt1:/feasst/particle/propane.txt add_num_pt1_particles=1\"\n",
    "        params['potential'] = \"Potential Model=LennardJones\"\n",
    "        params['branch'] = 'particle_type=pt1 angle=true mobile_site=2 anchor_site=1 anchor_site2=0'\n",
    "        params['grow'] = \"\"\"particle_type=pt1 transfer=true mobile_site=0\n",
    "bond=true mobile_site=1 anchor_site=0\n",
    "angle=true mobile_site=2 anchor_site=1 anchor_site2=0\"\"\"\n",
    "    elif params['molecule'] == 'floppy_angle':\n",
    "        params['branch'] = 'particle_type=pt1 angle=true mobile_site=2 anchor_site=1 anchor_site2=0'\n",
    "        params['grow'] = \"\"\"particle_type=pt1 transfer=true site=0\n",
    "bond=true mobile_site=1 anchor_site=0\n",
    "angle=true mobile_site=2 anchor_site=1 anchor_site2=0\"\"\"\n",
    "    elif params['molecule'] == 'floppy_dihedral':\n",
    "        params['branch'] = \"\"\"particle_type=pt1 angle=true mobile_site=2 anchor_site=1 anchor_site2=0\n",
    "dihedral=true mobile_site=3 anchor_site=2 anchor_site2=1 anchor_site3=0\"\"\"\n",
    "        params['grow'] = \"\"\"particle_type=pt1 transfer=true site=0\n",
    "bond=true mobile_site=1 anchor_site=0\n",
    "angle=true mobile_site=2 anchor_site=1 anchor_site2=0\n",
    "dihedral=true mobile_site=3 anchor_site=2 anchor_site2=1 anchor_site3=0\"\"\"\n",
    "    else:\n",
    "        params['branch'] = 'particle_type=pt1 branch=true mobile_site=2 mobile_site2=3 anchor_site=0 anchor_site2=1'\n",
    "        params['grow'] = \"\"\"particle_type=pt1 transfer=true site=0\n",
    "bond=true mobile_site=1 anchor_site=0\n",
    "branch=true mobile_site=2 mobile_site2=3 anchor_site=0 anchor_site2=1\"\"\"\n",
    "    if params['molecule'] == '2-methylpropane':\n",
    "        params['config'] = \"Configuration cubic_side_length=30 particle_type=pt1:/feasst/particle/2-methylpropane.txt add_num_pt1_particles=1\"\n",
    "        params['potential'] = \"Potential Model=LennardJones\"\n",
    "    elif params['molecule'] == 'floppy_branch':\n",
    "        params['config'] = \"\"\"Configuration cubic_side_length=30 particle_type=pt1:{prefix}.txt add_num_pt1_particles=1 \"\"\".format(**params)\n",
    "        params['potential'] = \"Potential Model=HardSphere\"\n",
    "        if 'intra' in params:\n",
    "            if params['intra'] == 1:\n",
    "                params['potential'] += \"\\nPotential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true\"\n",
    "        with open(params['prefix']+'.txt', 'w') as file1:\n",
    "            file1.write(\"\"\"#FEASST particle file\n",
    "            \n",
    "Site Properties\n",
    "\n",
    "0 epsilon 1 sigma 1 cutoff 1\n",
    "\n",
    "Sites\n",
    "\n",
    "0 0 0 0 0\n",
    "1 0 1 0 0\n",
    "2 0 0 1 0\n",
    "3 0 0 0 1\n",
    "\n",
    "Bonds\n",
    "\n",
    "0 0 0 1\n",
    "0 0 0 2\n",
    "0 0 0 3\n",
    "\n",
    "Bond Properties\n",
    "\n",
    "0 RigidBond length=1 delta=0.0001\n",
    "\n",
    "Angles\n",
    "\n",
    "0 0 1 0 2\n",
    "0 0 1 0 3\n",
    "0 0 2 0 3\n",
    "\n",
    "Angle Properties\n",
    "\n",
    "0 AngleHarmonic equilibrium_degrees=60 k_energy_per_radian_sq={k} num_jacobian_gaussian={num_jacobian_gaussian}\n",
    "\"\"\".format(**params))\n",
    "    elif params['molecule'] == 'floppy_angle':\n",
    "        params['config'] = \"\"\"Configuration cubic_side_length=30 particle_type=pt1:{prefix}.txt add_num_pt1_particles 1 \"\"\".format(**params)\n",
    "        params['potential'] = \"Potential Model=HardSphere\" \n",
    "        with open(params['prefix']+'.txt', 'w') as file1:\n",
    "            file1.write(\"\"\"#FEASST particle file\n",
    "\n",
    "Site Properties\n",
    "\n",
    "0 epsilon 1 sigma 1 cutoff 1\n",
    "\n",
    "Sites\n",
    "\n",
    "0 0 0 0 0\n",
    "1 0 1 0 0\n",
    "2 0 2 0 0\n",
    "\n",
    "Bond Properties\n",
    "\n",
    "0 RigidBond length=1 delta=0.0001\n",
    "\n",
    "Bonds\n",
    "\n",
    "0 0 0 1\n",
    "1 0 1 2\n",
    "\n",
    "Angle Properties\n",
    "\n",
    "0 AngleHarmonic equilibrium_degrees=60 k_energy_per_radian_sq={k}\n",
    "\n",
    "Angles\n",
    "\n",
    "0 0 0 1 2\n",
    "\"\"\".format(**params))\n",
    "    elif params['molecule'] == 'floppy_dihedral':\n",
    "        params['config'] = \"\"\"Configuration cubic_side_length=30 particle_type=pt1:{prefix}.txt add_num_pt1_particles=1 \"\"\".format(**params)\n",
    "        params['potential'] = \"Potential Model=HardSphere\" \n",
    "        with open(params['prefix']+'.txt', 'w') as file1:\n",
    "            file1.write(\"\"\"#FEASST particle file\n",
    "\n",
    "Site Properties\n",
    "\n",
    "0 epsilon=1 sigma=1 cutoff=1\n",
    "\n",
    "Sites\n",
    "\n",
    "0 0 0 0 0\n",
    "1 0 1 0 0\n",
    "2 0 1 1 0\n",
    "3 0 1 1 1\n",
    "\n",
    "Bond Properties\n",
    "\n",
    "0 RigidBond length=1 delta=0.0001\n",
    "\n",
    "Bonds\n",
    "\n",
    "0 0 0 1\n",
    "1 0 1 2\n",
    "2 0 2 3\n",
    "\n",
    "Angle Properties\n",
    "\n",
    "0 AngleHarmonic equilibrium_degrees=60 k_energy_per_radian_sq=0\n",
    "\n",
    "Angles\n",
    "\n",
    "0 0 0 1 2\n",
    "1 0 1 2 3\n",
    "\n",
    "Dihedral Properties\n",
    "\n",
    "0 DihedralTraPPE c0=0 c1=0 c2=0 c3=0\n",
    "\n",
    "Dihedrals\n",
    "\n",
    "0 0 0 1 2 3\n",
    "\n",
    "\"\"\".format(**params))\n",
    "    else:\n",
    "        if params['molecule'] != 'propane':\n",
    "            assert False, 'unrecognized molecule'\n",
    "        \n",
    "def run_sim(params):\n",
    "    print('running')\n",
    "    update_params(params)\n",
    "    with open(params['prefix']+'_grow.txt', 'w') as file1:\n",
    "        file1.write(\"\"\"TrialGrowFile\\n\\n{branch}\"\"\".format(**params))\n",
    "\n",
    "    # write fst script to run a single simulation\n",
    "    with open('script5.txt', \"w\") as myfile: myfile.write(\"\"\"\n",
    "MonteCarlo\n",
    "RandomMT19937 seed={seed}\n",
    "{config}\n",
    "{potential}\n",
    "ThermoParams beta={beta} chemical_potential=1\n",
    "Metropolis\n",
    "TrialGrowFile grow_file={prefix}_grow.txt\n",
    "Let [write]=trials_per_write={trials_per} output_file={prefix}\n",
    "Movie [write].xyz\n",
    "Log [write].txt include_bonds=true\n",
    "Energy [write]_en.txt\n",
    "AnalyzeBonds [write]_bonds{num_jacobian_gaussian}.txt angle_bin_width={angle_bin_width} angle_bin_center={angle_bin_center} dihedral_bin_width={angle_bin_width} dihedral_bin_center={angle_bin_center}\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(\"../../../build/bin/fst < script5.txt > script5.log\", shell=True, executable='/bin/bash')\n",
    "    with open('script5.log', 'r') as file: print(file.read(), '\\n', 'exit:', syscode)\n",
    "\n",
    "def run_gce_ideal_gas_sim(params):\n",
    "    print('running')\n",
    "    update_params(params)\n",
    "    with open(params['prefix']+'_gce_grow.txt', 'w') as file1:\n",
    "        file1.write(\"\"\"TrialGrowFile\\n\\n{grow}\"\"\".format(**params))\n",
    "\n",
    "    # write fst script to run a single simulation\n",
    "    with open('script5.txt', \"w\") as myfile: myfile.write(\"\"\"\n",
    "MonteCarlo\n",
    "RandomMT19937 seed={seed}\n",
    "{config}\n",
    "Potential Model=IdealGas\n",
    "ThermoParams beta={beta} chemical_potential=-25\n",
    "Metropolis\n",
    "#FlatHistogram Macrostate=MacrostateNumParticles width=1 min=0 max=1 Bias=TransitionMatrix min_sweeps=1e3\n",
    "TrialGrowFile grow_file={prefix}_gce_grow.txt\n",
    "Let [write]=trials_per_write {trials_per} output_file={prefix}\n",
    "Movie [write].xyz\n",
    "Log [write].txt include_bonds true\n",
    "Energy [write]_en.txt\n",
    "AnalyzeBonds [write]_bonds{num_jacobian_gaussian}.txt angle_bin_width={angle_bin_width} angle_bin_center={angle_bin_center}\n",
    "CheckEnergy trials_per_update={trials_per} decimal_places=6\n",
    "NumParticles [write]_bonds{num_jacobian_gaussian}_num.txt\n",
    "#CriteriaUpdater trials_per_update=1e4\n",
    "#CriteriaWriter [write]_bonds{num_jacobian_gaussian}_crit.txt\n",
    "#Run until=complete\n",
    "Run num_trials={production}\n",
    "\"\"\".format(**params))\n",
    "\n",
    "    import subprocess\n",
    "    syscode = subprocess.call(\"../../../build/bin/fst < script5.txt > script5.log\", shell=True, executable='/bin/bash')\n",
    "    with open('script5.log', 'r') as file: print(file.read(), '\\n', 'exit:', syscode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c4dc5143",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=51376018\n",
      "# Initializing random number generator with seed: 51376018\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:floppy_angle.txt\n",
      "Potential Model=HardSphere\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=floppy_angle_grow.txt\n",
      "Movie output_file=floppy_angle.xyz trials_per_write=1000000.0\n",
      "Log include_bonds=true output_file=floppy_angle.txt trials_per_write=1000000.0\n",
      "Energy output_file=floppy_angle_en.txt trials_per_write=1000000.0\n",
      "AnalyzeBonds angle_bin_center=1.5707963267948966 angle_bin_width=0.01 dihedral_bin_center=1.5707963267948966 dihedral_bin_width=0.01 output_file=floppy_angle_bonds0.txt trials_per_write=1000000.0\n",
      "CheckEnergy decimal_places=6 trials_per_update=1000000.0\n",
      "Run num_trials=1000000.0\n",
      " \n",
      " exit: 0\n"
     ]
    }
   ],
   "source": [
    "params['k']=0\n",
    "params['prefix'] = params['molecule'] = 'floppy_angle'\n",
    "params['trials_per'] = 1e6\n",
    "params['production'] = 1e6\n",
    "params['angle_bin_center'] = 0.5*np.pi\n",
    "#%timeit \n",
    "run_sim(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6b545813",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('floppy_angle_bonds0.txt', skiprows=10, header=None)\n",
    "if params['show_plot']:\n",
    "    mid=int(len(df[0])/2) # ensure k=0, i=0 is last jg\n",
    "    plt.plot(df[0]/np.pi*180, np.sin(df[0])*np.average(df[1][mid-4:mid+4]), label='sin', color='black', linestyle='solid')\n",
    "    plt.step(df[0]/np.pi*180, df[1])\n",
    "    plt.axvline(90, color='black', linestyle='dotted')\n",
    "    plt.ylabel('degrees')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "256d984b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=863675017\n",
      "# Initializing random number generator with seed: 863675017\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:floppy_dihedral.txt\n",
      "Potential Model=HardSphere\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=floppy_dihedral_grow.txt\n",
      "Movie output_file=floppy_dihedral.xyz trials_per_write=1000000.0\n",
      "Log include_bonds=true output_file=floppy_dihedral.txt trials_per_write=1000000.0\n",
      "Energy output_file=floppy_dihedral_en.txt trials_per_write=1000000.0\n",
      "AnalyzeBonds angle_bin_center=1.5707963267948966 angle_bin_width=0.01 dihedral_bin_center=1.5707963267948966 dihedral_bin_width=0.01 output_file=floppy_dihedral_bonds0.txt trials_per_write=1000000.0\n",
      "CheckEnergy decimal_places=6 trials_per_update=1000000.0\n",
      "Run num_trials=1000000.0\n",
      " \n",
      " exit: 0\n"
     ]
    }
   ],
   "source": [
    "params['k']=0\n",
    "params['prefix'] = params['molecule'] = 'floppy_dihedral'\n",
    "params['trials_per'] = 1e6\n",
    "params['production'] = 1e6\n",
    "params['angle_bin_center'] = 0.5*np.pi\n",
    "#%timeit \n",
    "run_sim(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5690cebc",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('floppy_dihedral_bonds0.txt', skiprows=10, nrows=310, header=None)\n",
    "df2 = pd.read_csv('floppy_dihedral_bonds0.txt', skiprows=328, nrows=310, header=None)\n",
    "if params['show_plot']:\n",
    "    mid=int(len(df[0])/2) # ensure k=0, i=0 is last jg\n",
    "    plt.plot(df[0]/np.pi*180, np.sin(df[0])*np.average(df[1][mid-4:mid+4]), label='sin', color='black', linestyle='solid')\n",
    "    plt.step(df[0]/np.pi*180, df[1], label='angle')\n",
    "    plt.step(df2[0][1:]/np.pi*180, df2[1][1:], label='dihedral')\n",
    "    plt.axvline(90, color='black', linestyle='dotted')\n",
    "    plt.xlabel('degrees')\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a10eb899",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=627491116\n",
      "# Initializing random number generator with seed: 627491116\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:/feasst/particle/propane.txt\n",
      "Potential Model=LennardJones\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=propane_grow.txt\n",
      "Movie output_file=propane.xyz trials_per_write=1000000.0\n",
      "Log include_bonds=true output_file=propane.txt trials_per_write=1000000.0\n",
      "Energy output_file=propane_en.txt trials_per_write=1000000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=propane_bonds0.txt trials_per_write=1000000.0\n",
      "CheckEnergy decimal_places=6 trials_per_update=1000000.0\n",
      "Run num_trials=1000000.0\n",
      " \n",
      " exit: 0\n"
     ]
    }
   ],
   "source": [
    "params['prefix'] = params['molecule'] = 'propane'\n",
    "params['production'] = params['trials_per'] = 1e6\n",
    "params['angle_bin_center'] = 100/180*np.pi\n",
    "#%timeit \n",
    "run_sim(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f2e0ca16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# digitized, poorly, from https://doi.org/10.1021/acs.jctc.7b00173\n",
    "chen = pd.read_csv('../test/data/propane_bonds.txt', comment='#', header=None)\n",
    "df = pd.read_csv('propane_bonds0.txt', skiprows=10, header=None)\n",
    "for d in [df, chen]:\n",
    "    d[1] /= d[1].max()\n",
    "assert np.abs(chen[1].mean() - df[1].mean()) < 0.2\n",
    "assert np.abs(chen[1].std() - df[1].std()) < 0.1\n",
    "if params['show_plot']:\n",
    "    plt.step(chen[0], chen[1], color='black')\n",
    "    plt.step(df[0]/np.pi*180, df[1])\n",
    "    plt.axvline(114, color='black', linestyle='dotted')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "photographic-small",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=592833805\n",
      "# Initializing random number generator with seed: 592833805\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:/feasst/particle/2-methylpropane.txt\n",
      "Potential Model=LennardJones\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=2-methylpropane_grow.txt\n",
      "Movie output_file=2-methylpropane.xyz trials_per_write=1000000.0\n",
      "Log include_bonds=true output_file=2-methylpropane.txt trials_per_write=1000000.0\n",
      "Energy output_file=2-methylpropane_en.txt trials_per_write=1000000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=2-methylpropane_bonds0.txt trials_per_write=1000000.0\n",
      "CheckEnergy decimal_places=6 trials_per_update=1000000.0\n",
      "Run num_trials=1000000.0\n",
      " \n",
      " exit: 0\n"
     ]
    }
   ],
   "source": [
    "params['prefix'] = params['molecule'] = '2-methylpropane'\n",
    "params['production'] = params['trials_per'] = 1e6\n",
    "#%timeit run_sim(params)\n",
    "run_sim(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "293c00fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# digitized, poorly, from https://doi.org/10.1021/acs.jctc.7b00173\n",
    "chen = pd.read_csv('../test/data/2-methylpropane_bonds.txt', comment='#', header=None)\n",
    "df = pd.read_csv('2-methylpropane_bonds0.txt', skiprows=10, header=None)\n",
    "for d in [df, chen]:\n",
    "    d[1] /= d[1].max()\n",
    "assert np.abs(chen[1].mean() - df[1].mean()) < 0.4\n",
    "assert np.abs(chen[1].std() - df[1].std()) < 0.2\n",
    "if params['show_plot']:\n",
    "    plt.step(chen[0], chen[1], color='black')\n",
    "    plt.step(df[0]/np.pi*180, df[1])\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c9c0023a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=656680569\n",
      "# Initializing random number generator with seed: 656680569\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i0.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk0i0.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk0i0.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk0i0_grow.txt\n",
      "Movie output_file=branchk0i0.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk0i0.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk0i0_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk0i0_bonds0.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=791754052\n",
      "# Initializing random number generator with seed: 791754052\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i1.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk0i1.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk0i1.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk0i1_grow.txt\n",
      "Movie output_file=branchk0i1.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk0i1.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk0i1_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk0i1_bonds0.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=570428588\n",
      "# Initializing random number generator with seed: 570428588\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i0.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk0i0.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk0i0.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk0i0_grow.txt\n",
      "Movie output_file=branchk0i0.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk0i0.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk0i0_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk0i0_bonds1.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=633455370\n",
      "# Initializing random number generator with seed: 633455370\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i1.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk0i1.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk0i1.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk0i1_grow.txt\n",
      "Movie output_file=branchk0i1.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk0i1.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk0i1_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk0i1_bonds1.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=192639018\n",
      "# Initializing random number generator with seed: 192639018\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i0.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk0i0.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk0i0.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk0i0_grow.txt\n",
      "Movie output_file=branchk0i0.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk0i0.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk0i0_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk0i0_bonds4.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=267590325\n",
      "# Initializing random number generator with seed: 267590325\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i1.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk0i1.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk0i1.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk0i1_grow.txt\n",
      "Movie output_file=branchk0i1.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk0i1.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk0i1_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk0i1_bonds4.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=385578983\n",
      "# Initializing random number generator with seed: 385578983\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk10i0_grow.txt\n",
      "Movie output_file=branchk10i0.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk10i0_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk10i0_bonds0.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=979224303\n",
      "# Initializing random number generator with seed: 979224303\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i1.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk10i1.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk10i1.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk10i1_grow.txt\n",
      "Movie output_file=branchk10i1.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk10i1.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk10i1_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk10i1_bonds0.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=994542371\n",
      "# Initializing random number generator with seed: 994542371\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk10i0_grow.txt\n",
      "Movie output_file=branchk10i0.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk10i0_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk10i0_bonds1.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=464190821\n",
      "# Initializing random number generator with seed: 464190821\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i1.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk10i1.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk10i1.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk10i1_grow.txt\n",
      "Movie output_file=branchk10i1.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk10i1.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk10i1_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk10i1_bonds1.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=903843675\n",
      "# Initializing random number generator with seed: 903843675\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk10i0_grow.txt\n",
      "Movie output_file=branchk10i0.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk10i0_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk10i0_bonds4.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",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=881915403\n",
      "# Initializing random number generator with seed: 881915403\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i1.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk10i1.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk10i1.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=HardSphere\n",
      "Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=1\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk10i1_grow.txt\n",
      "Movie output_file=branchk10i1.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk10i1.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk10i1_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 dihedral_bin_center=1.7453292519943295 dihedral_bin_width=0.01 output_file=branchk10i1_bonds4.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": [
    "ks=[0, 10]\n",
    "njgs=[0, 1, 4]\n",
    "for k in ks:\n",
    "    params['k'] = k\n",
    "    for njg in njgs:\n",
    "        for intra in [0, 1]:\n",
    "            params['intra'] = intra\n",
    "            params['prefix'] = \"branchk\" + str(params['k'])+'i'+str(params['intra'])\n",
    "            params['molecule'] = 'floppy_branch'\n",
    "            #params['production'] = params['trials_per'] = 1e5\n",
    "            params['trials_per'] = 1e5\n",
    "            params['production'] = 1e5\n",
    "            params['num_jacobian_gaussian']=njg\n",
    "            run_sim(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "sublime-extent",
   "metadata": {},
   "outputs": [],
   "source": [
    "if params['show_plot']:\n",
    "    for k in reversed(ks):\n",
    "        linestyle='solid'\n",
    "        if k == 10: linestyle='dashed'\n",
    "        for intra in [1, 0]:\n",
    "            for njg in njgs:\n",
    "                if njg == 0:\n",
    "                    color='black'\n",
    "                elif njg == 1:\n",
    "                    color='blue'\n",
    "                    if intra == 1: color='green'\n",
    "                elif njg == 4:\n",
    "                    color='red'\n",
    "                    if intra == 1: color='orange'\n",
    "                jg = pd.read_csv(\"branchk\"+str(k)+'i'+str(intra)+'_bonds'+str(njg)+'.txt', skiprows=10, header=None)\n",
    "                plt.step(jg[0]/np.pi*180, jg[1], label='jg '+str(njg)+' intra '+str(intra)+' k='+str(k), linestyle=linestyle, color=color)\n",
    "    mid=int(len(jg[0])/2) # ensure k=0, i=0 is last jg\n",
    "    plt.plot(jg[0]/np.pi*180, np.sin(jg[0])*np.average(jg[1][mid-4:mid+4]), label='sin', color='white', linestyle='dotted')\n",
    "    plt.legend()\n",
    "    plt.xlabel(r'$\\theta$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f53e0c11",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=338534626\n",
      "# Initializing random number generator with seed: 338534626\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=IdealGas\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=-25\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk10i0_gce_grow.txt\n",
      "Movie output_file=branchk10i0.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk10i0_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 output_file=branchk10i0_bonds0.txt trials_per_write=100000.0\n",
      "CheckEnergy decimal_places=6 trials_per_update=100000.0\n",
      "NumParticles output_file=branchk10i0_bonds0_num.txt trials_per_write=100000.0\n",
      "Run num_trials=1000000.0\n",
      " \n",
      " exit: 0\n",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=327568917\n",
      "# Initializing random number generator with seed: 327568917\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=IdealGas\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=-25\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk10i0_gce_grow.txt\n",
      "Movie output_file=branchk10i0.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk10i0_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 output_file=branchk10i0_bonds1.txt trials_per_write=100000.0\n",
      "CheckEnergy decimal_places=6 trials_per_update=100000.0\n",
      "NumParticles output_file=branchk10i0_bonds1_num.txt trials_per_write=100000.0\n",
      "Run num_trials=1000000.0\n",
      " \n",
      " exit: 0\n",
      "running\n",
      "# Usage: /home/hwh/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.16\n",
      "MonteCarlo\n",
      "RandomMT19937 seed=875465004\n",
      "# Initializing random number generator with seed: 875465004\n",
      "Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt\n",
      "# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt\n",
      "# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt\n",
      "# ... and so on (suppressed following rename warnings)\n",
      "Potential Model=IdealGas\n",
      "ThermoParams beta=0.4009078501424201 chemical_potential=-25\n",
      "Metropolis\n",
      "TrialGrowFile grow_file=branchk10i0_gce_grow.txt\n",
      "Movie output_file=branchk10i0.xyz trials_per_write=100000.0\n",
      "Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0\n",
      "Energy output_file=branchk10i0_en.txt trials_per_write=100000.0\n",
      "AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 output_file=branchk10i0_bonds4.txt trials_per_write=100000.0\n",
      "CheckEnergy decimal_places=6 trials_per_update=100000.0\n",
      "NumParticles output_file=branchk10i0_bonds4_num.txt trials_per_write=100000.0\n",
      "Run num_trials=1000000.0\n",
      " \n",
      " exit: 0\n"
     ]
    }
   ],
   "source": [
    "# compare ideal gas GCE\n",
    "ks=[10]#[0, 10]\n",
    "njgs=[0, 1, 4]\n",
    "for k in ks:\n",
    "    params['k'] = k\n",
    "    for njg in njgs:\n",
    "        for intra in [0]:#[0, 1]:\n",
    "            params['intra'] = intra\n",
    "            params['prefix'] = \"branchk\" + str(params['k'])+'i'+str(params['intra'])\n",
    "            params['molecule'] = 'floppy_branch'\n",
    "            #params['production'] = params['trials_per'] = 1e5\n",
    "            params['trials_per'] = 1e5\n",
    "            params['production'] = 1e6\n",
    "            params['num_jacobian_gaussian']=njg\n",
    "            run_gce_ideal_gas_sim(params)"
   ]
  },
  {
   "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
}
