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 numpy as np
import pandas as pd

params = {
    "fstprt": "/feasst/plugin/chain/particle/chain5.txt",
    "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 chain 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('script4.txt', "w") as myfile: myfile.write("""
MonteCarlo
RandomMT19937 seed={seed}
Configuration cubic_side_length=30 particle_type=chain:{fstprt} add_num_chain_particles=1
Potential Model=LennardJones VisitModel=VisitModelIntra intra_cut=1
ThermoParams beta={beta} chemical_potential=1
Metropolis
TrialGrowFile file_name=chain5_grow.txt
Let [write]=trials_per_write={trials_per} output_file=chain5
Log [write].txt
Movie [write].xyz
Energy [write]en.txt
RadiusOfGyration [write]rg.txt
CheckEnergy trials_per_update={trials_per} decimal_places=8
Run num_trials={production}
""".format(**params))

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

# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=199934304
# Initializing random number generator with seed: 199934304
Configuration add_num_chain_particles=1 cubic_side_length=30 particle_type=chain:/feasst/plugin/chain/particle/chain5.txt
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:416] TrialGrowFile::file_name renamed to grow_file.
Log output_file=chain5.txt trials_per_write=100000.0
Movie output_file=chain5.xyz trials_per_write=100000.0
Energy output_file=chain5en.txt trials_per_write=100000.0
RadiusOfGyration output_file=chain5rg.txt trials_per_write=100000.0
CheckEnergy decimal_places=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

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]
assert np.abs(rg - np.sqrt(0.88)) < 0.01
assert np.abs(en + 2.04) < 0.06
rg_en = rgdf['rgu'].values[0]
drg_dbeta = rg*en - rg_en
assert np.abs(drg_dbeta + 0.186/2) < 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()

drg_dbeta -0.09754849439745095

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