Simulation of a single freely-jointed chain of 10,000 monomers

import math
import numpy as np
import pandas as pd

params = {"num_monomers": 10000}
params["length"] = 3*params["num_monomers"]

script="""
MonteCarlo
Configuration cubic_side_length={length} particle_type=chain:chain.txt.swp
Potential Model=IdealGas
ThermoParams beta=1 chemical_potential=1
Metropolis
TrialAdd particle_type=chain
Run until_num_particles=1
Remove name=TrialAdd
TrialGrowFile grow_file=grow_chain.txt.swp
Let [write]=trials_per_write=1 output_file=chain
Log [write].csv
Movie [write].xyz
EndToEndDistance output_file=end_to_end.txt
RadiusOfGyration output_file=rg0.txt
Run num_trials=1e2
""".format(**params)

def linear_single_site(n, filename):
    # write fstprt file
    with open (filename, "w") as f:
        f.write(
"# FEASST particle file\n\
\n\
Site Properties\n\
\n\
0 epsilon 1 sigma 1.0 cutoff 3.0\n\
\n\
Bond Properties\n\
\n\
0 RigidBond length 1.0 delta 0.00001\n\
\n\
Sites\n\
\n\
")

        for i in range(n):
            f.write(str(i) + " 0 " + str(i) + ".0 0.0 0.0\n")
        f.write("\nBonds\n\n")
        for i in range(n - 1):
            f.write(str(i) + " 0 " + str(i) + " " + str(i + 1) + "\n")
        f.write("\n")

    # write TrialGrowFile
    with open('grow_' + filename, 'w') as f:
        f.write("TrialGrowFile\n\nparticle_type chain regrow true site 0\n")
        for i in range(1, n):
            f.write("bond true mobile_site "+str(i)+" anchor_site "+str(i-1)+"\n")

linear_single_site(params["num_monomers"], "chain.txt.swp")

with open('script0.txt', 'w') as file: file.write(script)
import subprocess
syscode = subprocess.call("../../../build/bin/fst < script0.txt > script0.log", shell=True, executable='/bin/bash')
with open('script0.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
Configuration cubic_side_length=30000 particle_type=chain:chain.txt.swp
Potential Model=IdealGas
ThermoParams beta=1 chemical_potential=1
Metropolis
TrialAdd particle_type=chain
Run until_num_particles=1
# Initializing random number generator with seed: 1748975657
Remove name=TrialAdd
TrialGrowFile grow_file=grow_chain.txt.swp
Log output_file=chain.csv trials_per_write=1
Movie output_file=chain.xyz trials_per_write=1
EndToEndDistance output_file=end_to_end.txt
RadiusOfGyration output_file=rg0.txt
Run num_trials=1e2

 exit: 0

end_to_end=pd.read_csv('end_to_end.txt')
assert np.abs(end_to_end['average'][0] - math.sqrt(params['num_monomers'])) < 20
rg=pd.read_csv('rg0.txt')
assert np.abs(rg['average'][0] - np.sqrt(params['num_monomers']/6)) < 200

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