Temperature extrapolation of the radius of gyration

Simulate a 5-mer simple LJ chain and compute the radius of gyration.

See Section III.A of https://doi.org/10.1063/1.5026493 and https://doi.org/10.1063/1.1350578

\[\langle R_g \rangle = \frac{\sum_i R_g^i \exp(-\beta U_i)}{\sum_i \exp(-\beta U_i)}\]

\[\frac{d\langle R_g \rangle}{d\beta} = -\langle R_g U\rangle + \langle R_g \rangle \langle U \rangle\]

import sys
import subprocess
import argparse
import random
import unittest
import numpy as np
import pandas as pd

params = {
    "fstprt": "/feasst/plugin/chain/particle/chain5.fstprt",
    "cubic_side_length": 90, "beta": 1,
    "trials_per": 1e5, "seed": random.randrange(int(1e9)),
    "production": 1e5}

with open('chain5_grow.txt', 'w') as file1:
    file1.write("""
TrialGrowFile

bond true mobile_site 1 anchor_site 0 particle_type 0 weight 1
bond true mobile_site 2 anchor_site 1
bond true mobile_site 3 anchor_site 2
bond true mobile_site 4 anchor_site 3""")

# write fst script to run a single simulation
with open('script.txt', "w") as myfile: myfile.write("""
MonteCarlo
RandomMT19937 seed {seed}
Configuration cubic_side_length 30 add_particles_of_type0 1 particle_type0 {fstprt}
Potential Model LennardJones VisitModel VisitModelIntra intra_cut 1
ThermoParams beta {beta} chemical_potential 1
Metropolis
TrialGrowFile file_name chain5_grow.txt
Movie trials_per_write {trials_per} file_name chain5.xyz
Log trials_per_write {trials_per} file_name chain5.txt
Energy trials_per_write {trials_per} file_name chain5en.txt
RadiusOfGyration trials_per_write {trials_per} file_name chain5rg.txt
CheckEnergy trials_per_update {trials_per} tolerance 1e-8
Run num_trials {production}
""".format(**params))

import subprocess
syscode = subprocess.call("../../../build/bin/fst < script.txt > script.log", shell=True, executable='/bin/bash')
with open('script.log', 'r') as file: print(file.read(), '\n', 'exit:', syscode)

# FEASST version: v0.24.2-6-gf3cf33172b-dirty-hwh/whatnow
MonteCarlo
RandomMT19937 seed 101957491
# initializing random number generator with seed: 101957491
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 /home/hwh/feasst/plugin/chain/particle/chain5.fstprt
Potential Model LennardJones VisitModel VisitModelIntra intra_cut 1
ThermoParams beta 1 chemical_potential 1
Metropolis
TrialGrowFile file_name chain5_grow.txt
# Warn 0 [plugin/chain/src/trial_grow.cpp:431] TrialGrowFile::file_name renamed to grow_file.
Movie file_name chain5.xyz trials_per_write 100000.0
# Warn 0 [plugin/monte_carlo/src/stepper.cpp:22] Stepper argument file_name was renamed to output_file.
Log file_name chain5.txt trials_per_write 100000.0
# Warn 0 [plugin/monte_carlo/src/stepper.cpp:22] Stepper argument file_name was renamed to output_file.
Energy file_name chain5en.txt trials_per_write 100000.0
# Warn 0 [plugin/monte_carlo/src/stepper.cpp:22] Stepper argument file_name was renamed to output_file.
RadiusOfGyration file_name chain5rg.txt trials_per_write 100000.0
# Warn 0 [plugin/monte_carlo/src/stepper.cpp:22] Stepper argument file_name was renamed to output_file.
CheckEnergy tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0

 exit: 0

def rg_fit(beta, beta1, rg1, der1):
    return rg1 + (beta - beta1)*der1

class TestFreelyJointedIdealChain(unittest.TestCase):
    def test(self):
        taylor2001 = pd.read_csv('../test/data/taylor2001rg.csv', header=None)
        #print(taylor2001)
        rgdf = pd.read_csv('chain5rg.txt')
        endf = pd.read_csv('chain5en.txt')
        rg = rgdf['average'].values[0]
        en = endf['average'].values[0]
        self.assertAlmostEqual(rg, np.sqrt(0.88), delta=0.01)
        self.assertAlmostEqual(en, -2.04, delta=0.06)
        rg_en = rgdf['rgu'].values[0]
        drg_dbeta = rg*en - rg_en
        self.assertAlmostEqual(drg_dbeta, -0.186/2, delta=0.01)
        print('drg_dbeta', drg_dbeta)

#         import matplotlib.pyplot as plt
#         plt.scatter(params['beta'], rg, label='explicit sim')
#         plt.plot(1./taylor2001[0], np.sqrt(taylor2001[1]), label='taylor2001')
#         plt.plot(1./taylor2001[0],
#                  rg_fit(beta=1./taylor2001[0],
#                         beta1=params['beta'],
#                         rg1=rg,
#                         der1=drg_dbeta),
#                  label='1st order extrapolation')
#         plt.xscale('log')
#         plt.ylim([0.2, 1.3])
#         plt.legend()
#         plt.xlabel(r'$\beta$', fontsize=16)
#         plt.ylabel(r'$R_g$', fontsize=16)
#         plt.show()


unittest.main(argv=[''], verbosity=2, exit=False)

test (__main__.TestFreelyJointedIdealChain) ... ok

----------------------------------------------------------------------
Ran 1 test in 0.003s

OK

drg_dbeta -0.09925819217655563

<unittest.main.TestProgram at 0x7f2aee30abc0>

Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!