Testing angles and branches
Following https://doi.org/10.1021/acs.jctc.7b00173 , compute the angle distributions.
First, consider propane and 2-methylpropane. Also consider very flexible molecules as follows:
floppy_angle: 0 - 1 - 2
floppy_dihedral: 0 - 1 - 2 - 3
floppy_branch:
1
|
0
/ \
2 3
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.
The floppy particles with no interactions should have sin(theta) angle distributions.
[1]:
import sys
import subprocess
import argparse
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pyfeasst import physical_constants
params = {
"cubic_side_length": 90,
"temperature_K": 300,
"trials_per": 1e5,
"production": 1e5,
"k": 1,
'angle_bin_width': 0.01,
'angle_bin_center': 100/180*np.pi,
'molecule': 'propane',
'num_jacobian_gaussian': 0,
#'molecule': '2-methylpropane',
#'molecule': 'floppy_branch',
'intra': 1,
'show_plot': False,
}
def update_params(params):
params['seed'] = random.randrange(int(1e9))
params['beta'] = 1./(params['temperature_K']*physical_constants.MolarGasConstant().value()/1e3) # mol/kJ
if params['molecule'] == 'propane':
params['config'] = "Configuration cubic_side_length=30 particle_type=pt1:/feasst/particle/propane.txt add_num_pt1_particles=1"
params['potential'] = "Potential Model=LennardJones"
params['branch'] = 'particle_type pt1 angle true mobile_site 2 anchor_site 1 anchor_site2 0'
params['grow'] = """particle_type pt1 transfer true mobile_site 0
bond true mobile_site 1 anchor_site 0
angle true mobile_site 2 anchor_site 1 anchor_site2 0"""
elif params['molecule'] == 'floppy_angle':
params['branch'] = 'particle_type pt1 angle true mobile_site 2 anchor_site 1 anchor_site2 0'
params['grow'] = """particle_type pt1 transfer true site 0
bond true mobile_site 1 anchor_site 0
angle true mobile_site 2 anchor_site 1 anchor_site2 0"""
elif params['molecule'] == 'floppy_dihedral':
params['branch'] = """particle_type pt1 angle true mobile_site 2 anchor_site 1 anchor_site2 0
dihedral true mobile_site 3 anchor_site 2 anchor_site2 1 anchor_site3 0"""
params['grow'] = """particle_type pt1 transfer true site 0
bond true mobile_site 1 anchor_site 0
angle true mobile_site 2 anchor_site 1 anchor_site2 0
dihedral true mobile_site 3 anchor_site 2 anchor_site2 1 anchor_site3 0"""
else:
params['branch'] = 'particle_type pt1 branch true mobile_site 2 mobile_site2 3 anchor_site 0 anchor_site2 1'
params['grow'] = """particle_type pt1 transfer true site 0
bond true mobile_site 1 anchor_site 0
branch true mobile_site 2 mobile_site2 3 anchor_site 0 anchor_site2 1"""
if params['molecule'] == '2-methylpropane':
params['config'] = "Configuration cubic_side_length=30 particle_type=pt1:/feasst/particle/2-methylpropane.txt add_num_pt1_particles=1"
params['potential'] = "Potential Model=LennardJones"
elif params['molecule'] == 'floppy_branch':
params['config'] = """Configuration cubic_side_length=30 particle_type=pt1:{prefix}.txt add_num_pt1_particles=1 """.format(**params)
params['potential'] = "Potential Model=HardSphere"
if 'intra' in params:
if params['intra'] == 1:
params['potential'] += "\nPotential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true"
with open(params['prefix']+'.txt', 'w') as file1:
file1.write("""#FEASST particle file
Site Properties
0 epsilon 1 sigma 1 cutoff 1
Sites
0 0 0 0 0
1 0 1 0 0
2 0 0 1 0
3 0 0 0 1
Bonds
0 0 0 1
0 0 0 2
0 0 0 3
Bond Properties
0 RigidBond length 1 delta 0.0001
Angles
0 0 1 0 2
0 0 1 0 3
0 0 2 0 3
Angle Properties
0 AngleHarmonic equilibrium_degrees 60 k_energy_per_radian_sq {k} num_jacobian_gaussian {num_jacobian_gaussian}
""".format(**params))
elif params['molecule'] == 'floppy_angle':
params['config'] = """Configuration cubic_side_length=30 particle_type=pt1:{prefix}.txt add_num_pt1_particles 1 """.format(**params)
params['potential'] = "Potential Model=HardSphere"
with open(params['prefix']+'.txt', 'w') as file1:
file1.write("""#FEASST particle file
Site Properties
0 epsilon 1 sigma 1 cutoff 1
Sites
0 0 0 0 0
1 0 1 0 0
2 0 2 0 0
Bond Properties
0 RigidBond length 1 delta 0.0001
Bonds
0 0 0 1
1 0 1 2
Angle Properties
0 AngleHarmonic equilibrium_degrees 60 k_energy_per_radian_sq {k}
Angles
0 0 0 1 2
""".format(**params))
elif params['molecule'] == 'floppy_dihedral':
params['config'] = """Configuration cubic_side_length=30 particle_type=pt1:{prefix}.txt add_num_pt1_particles=1 """.format(**params)
params['potential'] = "Potential Model=HardSphere"
with open(params['prefix']+'.txt', 'w') as file1:
file1.write("""#FEASST particle file
Site Properties
0 epsilon 1 sigma 1 cutoff 1
Sites
0 0 0 0 0
1 0 1 0 0
2 0 1 1 0
3 0 1 1 1
Bond Properties
0 RigidBond length 1 delta 0.0001
Bonds
0 0 0 1
1 0 1 2
2 0 2 3
Angle Properties
0 AngleHarmonic equilibrium_degrees 60 k_energy_per_radian_sq 0
Angles
0 0 0 1 2
1 0 1 2 3
Dihedral Properties
0 DihedralTraPPE c0 0 c1 0 c2 0 c3 0
Dihedrals
0 0 0 1 2 3
""".format(**params))
else:
if params['molecule'] != 'propane':
assert False, 'unrecognized molecule'
def run_sim(params):
print('running')
update_params(params)
with open(params['prefix']+'_grow.txt', 'w') as file1:
file1.write("""TrialGrowFile\n\n{branch}""".format(**params))
# write fst script to run a single simulation
with open('script5.txt', "w") as myfile: myfile.write("""
MonteCarlo
RandomMT19937 seed={seed}
{config}
{potential}
ThermoParams beta={beta} chemical_potential=1
Metropolis
TrialGrowFile grow_file={prefix}_grow.txt
Let [write]=trials_per_write={trials_per} output_file={prefix}
Movie [write].xyz
Log [write].txt include_bonds=true
Energy [write]_en.txt
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}
CheckEnergy trials_per_update={trials_per} decimal_places=8
Run num_trials={production}
""".format(**params))
import subprocess
syscode = subprocess.call("../../../build/bin/fst < script5.txt > script5.log", shell=True, executable='/bin/bash')
with open('script5.log', 'r') as file: print(file.read(), '\n', 'exit:', syscode)
def run_gce_ideal_gas_sim(params):
print('running')
update_params(params)
with open(params['prefix']+'_gce_grow.txt', 'w') as file1:
file1.write("""TrialGrowFile\n\n{grow}""".format(**params))
# write fst script to run a single simulation
with open('script5.txt', "w") as myfile: myfile.write("""
MonteCarlo
RandomMT19937 seed={seed}
{config}
Potential Model=IdealGas
ThermoParams beta={beta} chemical_potential=-25
Metropolis
#FlatHistogram Macrostate=MacrostateNumParticles width=1 min=0 max=1 Bias=TransitionMatrix min_sweeps=1e3
TrialGrowFile grow_file={prefix}_gce_grow.txt
Let [write]=trials_per_write {trials_per} output_file={prefix}
Movie [write].xyz
Log [write].txt include_bonds true
Energy [write]_en.txt
AnalyzeBonds [write]_bonds{num_jacobian_gaussian}.txt angle_bin_width={angle_bin_width} angle_bin_center={angle_bin_center}
CheckEnergy trials_per_update={trials_per} decimal_places=8
NumParticles [write]_bonds{num_jacobian_gaussian}_num.txt
#CriteriaUpdater trials_per_update=1e4
#CriteriaWriter [write]_bonds{num_jacobian_gaussian}_crit.txt
#Run until=complete
Run num_trials={production}
""".format(**params))
import subprocess
syscode = subprocess.call("../../../build/bin/fst < script5.txt > script5.log", shell=True, executable='/bin/bash')
with open('script5.log', 'r') as file: print(file.read(), '\n', 'exit:', syscode)
[2]:
params['k']=0
params['prefix'] = params['molecule'] = 'floppy_angle'
params['trials_per'] = 1e6
params['production'] = 1e6
params['angle_bin_center'] = 0.5*np.pi
#%timeit
run_sim(params)
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=464498080
# Initializing random number generator with seed: 464498080
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:floppy_angle.txt
Potential Model=HardSphere
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=floppy_angle_grow.txt
Movie output_file=floppy_angle.xyz trials_per_write=1000000.0
Log include_bonds=true output_file=floppy_angle.txt trials_per_write=1000000.0
Energy output_file=floppy_angle_en.txt trials_per_write=1000000.0
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
CheckEnergy decimal_places=8 trials_per_update=1000000.0
Run num_trials=1000000.0
exit: 0
[3]:
df = pd.read_csv('floppy_angle_bonds0.txt', skiprows=10, header=None)
if params['show_plot']:
mid=int(len(df[0])/2) # ensure k=0, i=0 is last jg
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')
plt.step(df[0]/np.pi*180, df[1])
plt.axvline(90, color='black', linestyle='dotted')
plt.ylabel('degrees')
[4]:
params['k']=0
params['prefix'] = params['molecule'] = 'floppy_dihedral'
params['trials_per'] = 1e6
params['production'] = 1e6
params['angle_bin_center'] = 0.5*np.pi
#%timeit
run_sim(params)
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=106926204
# Initializing random number generator with seed: 106926204
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:floppy_dihedral.txt
Potential Model=HardSphere
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=floppy_dihedral_grow.txt
Movie output_file=floppy_dihedral.xyz trials_per_write=1000000.0
Log include_bonds=true output_file=floppy_dihedral.txt trials_per_write=1000000.0
Energy output_file=floppy_dihedral_en.txt trials_per_write=1000000.0
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
CheckEnergy decimal_places=8 trials_per_update=1000000.0
Run num_trials=1000000.0
exit: 0
[5]:
df = pd.read_csv('floppy_dihedral_bonds0.txt', skiprows=10, nrows=310, header=None)
df2 = pd.read_csv('floppy_dihedral_bonds0.txt', skiprows=328, nrows=310, header=None)
if params['show_plot']:
mid=int(len(df[0])/2) # ensure k=0, i=0 is last jg
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')
plt.step(df[0]/np.pi*180, df[1], label='angle')
plt.step(df2[0][1:]/np.pi*180, df2[1][1:], label='dihedral')
plt.axvline(90, color='black', linestyle='dotted')
plt.xlabel('degrees')
plt.legend()
[6]:
params['prefix'] = params['molecule'] = 'propane'
params['production'] = params['trials_per'] = 1e6
params['angle_bin_center'] = 100/180*np.pi
#%timeit
run_sim(params)
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=742220871
# Initializing random number generator with seed: 742220871
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:/feasst/particle/propane.txt
Potential Model=LennardJones
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=propane_grow.txt
Movie output_file=propane.xyz trials_per_write=1000000.0
Log include_bonds=true output_file=propane.txt trials_per_write=1000000.0
Energy output_file=propane_en.txt trials_per_write=1000000.0
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
CheckEnergy decimal_places=8 trials_per_update=1000000.0
Run num_trials=1000000.0
exit: 0
[7]:
# digitized, poorly, from https://doi.org/10.1021/acs.jctc.7b00173
chen = pd.read_csv('../test/data/propane_bonds.txt', comment='#', header=None)
df = pd.read_csv('propane_bonds0.txt', skiprows=10, header=None)
for d in [df, chen]:
d[1] /= d[1].max()
assert np.abs(chen[1].mean() - df[1].mean()) < 0.2
assert np.abs(chen[1].std() - df[1].std()) < 0.1
if params['show_plot']:
plt.step(chen[0], chen[1], color='black')
plt.step(df[0]/np.pi*180, df[1])
plt.axvline(114, color='black', linestyle='dotted')
[8]:
params['prefix'] = params['molecule'] = '2-methylpropane'
params['production'] = params['trials_per'] = 1e6
#%timeit run_sim(params)
run_sim(params)
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=56179832
# Initializing random number generator with seed: 56179832
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:/feasst/particle/2-methylpropane.txt
Potential Model=LennardJones
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=2-methylpropane_grow.txt
Movie output_file=2-methylpropane.xyz trials_per_write=1000000.0
Log include_bonds=true output_file=2-methylpropane.txt trials_per_write=1000000.0
Energy output_file=2-methylpropane_en.txt trials_per_write=1000000.0
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
CheckEnergy decimal_places=8 trials_per_update=1000000.0
Run num_trials=1000000.0
exit: 0
[9]:
# digitized, poorly, from https://doi.org/10.1021/acs.jctc.7b00173
chen = pd.read_csv('../test/data/2-methylpropane_bonds.txt', comment='#', header=None)
df = pd.read_csv('2-methylpropane_bonds0.txt', skiprows=10, header=None)
for d in [df, chen]:
d[1] /= d[1].max()
assert np.abs(chen[1].mean() - df[1].mean()) < 0.4
assert np.abs(chen[1].std() - df[1].std()) < 0.2
if params['show_plot']:
plt.step(chen[0], chen[1], color='black')
plt.step(df[0]/np.pi*180, df[1])
plt.legend()
[10]:
ks=[0, 10]
njgs=[0, 1, 4]
for k in ks:
params['k'] = k
for njg in njgs:
for intra in [0, 1]:
params['intra'] = intra
params['prefix'] = "branchk" + str(params['k'])+'i'+str(params['intra'])
params['molecule'] = 'floppy_branch'
#params['production'] = params['trials_per'] = 1e5
params['trials_per'] = 1e5
params['production'] = 1e5
params['num_jacobian_gaussian']=njg
run_sim(params)
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=546964109
# Initializing random number generator with seed: 546964109
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i0.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk0i0.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk0i0.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk0i0_grow.txt
Movie output_file=branchk0i0.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk0i0.txt trials_per_write=100000.0
Energy output_file=branchk0i0_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=965216498
# Initializing random number generator with seed: 965216498
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i1.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk0i1.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk0i1.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk0i1_grow.txt
Movie output_file=branchk0i1.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk0i1.txt trials_per_write=100000.0
Energy output_file=branchk0i1_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=49576794
# Initializing random number generator with seed: 49576794
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i0.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk0i0.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk0i0.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk0i0_grow.txt
Movie output_file=branchk0i0.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk0i0.txt trials_per_write=100000.0
Energy output_file=branchk0i0_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=833521789
# Initializing random number generator with seed: 833521789
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i1.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk0i1.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk0i1.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk0i1_grow.txt
Movie output_file=branchk0i1.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk0i1.txt trials_per_write=100000.0
Energy output_file=branchk0i1_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=172008195
# Initializing random number generator with seed: 172008195
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i0.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk0i0.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk0i0.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk0i0_grow.txt
Movie output_file=branchk0i0.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk0i0.txt trials_per_write=100000.0
Energy output_file=branchk0i0_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=926832678
# Initializing random number generator with seed: 926832678
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk0i1.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk0i1.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk0i1.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk0i1_grow.txt
Movie output_file=branchk0i1.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk0i1.txt trials_per_write=100000.0
Energy output_file=branchk0i1_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=509031663
# Initializing random number generator with seed: 509031663
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk10i0_grow.txt
Movie output_file=branchk10i0.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0
Energy output_file=branchk10i0_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=539809929
# Initializing random number generator with seed: 539809929
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i1.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk10i1.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk10i1.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk10i1_grow.txt
Movie output_file=branchk10i1.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk10i1.txt trials_per_write=100000.0
Energy output_file=branchk10i1_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=275031964
# Initializing random number generator with seed: 275031964
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk10i0_grow.txt
Movie output_file=branchk10i0.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0
Energy output_file=branchk10i0_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=486210854
# Initializing random number generator with seed: 486210854
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i1.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk10i1.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk10i1.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk10i1_grow.txt
Movie output_file=branchk10i1.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk10i1.txt trials_per_write=100000.0
Energy output_file=branchk10i1_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=803034985
# Initializing random number generator with seed: 803034985
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk10i0_grow.txt
Movie output_file=branchk10i0.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0
Energy output_file=branchk10i0_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=274724731
# Initializing random number generator with seed: 274724731
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i1.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk10i1.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk10i1.txt
# ... and so on (suppressed following rename warnings)
Potential Model=HardSphere
Potential Model=HardSphere VisitModel=VisitModelIntraMap exclude_bonds=true
ThermoParams beta=0.4009078501424201 chemical_potential=1
Metropolis
TrialGrowFile grow_file=branchk10i1_grow.txt
Movie output_file=branchk10i1.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk10i1.txt trials_per_write=100000.0
Energy output_file=branchk10i1_en.txt trials_per_write=100000.0
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
CheckEnergy decimal_places=8 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
[11]:
if params['show_plot']:
for k in reversed(ks):
linestyle='solid'
if k == 10: linestyle='dashed'
for intra in [1, 0]:
for njg in njgs:
if njg == 0:
color='black'
elif njg == 1:
color='blue'
if intra == 1: color='green'
elif njg == 4:
color='red'
if intra == 1: color='orange'
jg = pd.read_csv("branchk"+str(k)+'i'+str(intra)+'_bonds'+str(njg)+'.txt', skiprows=10, header=None)
plt.step(jg[0]/np.pi*180, jg[1], label='jg '+str(njg)+' intra '+str(intra)+' k='+str(k), linestyle=linestyle, color=color)
mid=int(len(jg[0])/2) # ensure k=0, i=0 is last jg
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')
plt.legend()
plt.xlabel(r'$\theta$')
[12]:
# compare ideal gas GCE
ks=[10]#[0, 10]
njgs=[0, 1, 4]
for k in ks:
params['k'] = k
for njg in njgs:
for intra in [0]:#[0, 1]:
params['intra'] = intra
params['prefix'] = "branchk" + str(params['k'])+'i'+str(params['intra'])
params['molecule'] = 'floppy_branch'
#params['production'] = params['trials_per'] = 1e5
params['trials_per'] = 1e5
params['production'] = 1e6
params['num_jacobian_gaussian']=njg
run_gce_ideal_gas_sim(params)
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=831621889
# Initializing random number generator with seed: 831621889
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt
# ... and so on (suppressed following rename warnings)
Potential Model=IdealGas
ThermoParams beta=0.4009078501424201 chemical_potential=-25
Metropolis
TrialGrowFile grow_file=branchk10i0_gce_grow.txt
Movie output_file=branchk10i0.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0
Energy output_file=branchk10i0_en.txt trials_per_write=100000.0
AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 output_file=branchk10i0_bonds0.txt trials_per_write=100000.0
CheckEnergy decimal_places=8 trials_per_update=100000.0
NumParticles output_file=branchk10i0_bonds0_num.txt trials_per_write=100000.0
Run num_trials=1000000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=273644159
# Initializing random number generator with seed: 273644159
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt
# ... and so on (suppressed following rename warnings)
Potential Model=IdealGas
ThermoParams beta=0.4009078501424201 chemical_potential=-25
Metropolis
TrialGrowFile grow_file=branchk10i0_gce_grow.txt
Movie output_file=branchk10i0.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0
Energy output_file=branchk10i0_en.txt trials_per_write=100000.0
AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 output_file=branchk10i0_bonds1.txt trials_per_write=100000.0
CheckEnergy decimal_places=8 trials_per_update=100000.0
NumParticles output_file=branchk10i0_bonds1_num.txt trials_per_write=100000.0
Run num_trials=1000000.0
exit: 0
running
# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=392305367
# Initializing random number generator with seed: 392305367
Configuration add_num_pt1_particles=1 cubic_side_length=30 particle_type=pt1:branchk10i0.txt
# Renamed Bond name:0->1 for particle_type:0 in:branchk10i0.txt
# Renamed Bond name:0->2 for particle_type:0 in:branchk10i0.txt
# ... and so on (suppressed following rename warnings)
Potential Model=IdealGas
ThermoParams beta=0.4009078501424201 chemical_potential=-25
Metropolis
TrialGrowFile grow_file=branchk10i0_gce_grow.txt
Movie output_file=branchk10i0.xyz trials_per_write=100000.0
Log include_bonds=true output_file=branchk10i0.txt trials_per_write=100000.0
Energy output_file=branchk10i0_en.txt trials_per_write=100000.0
AnalyzeBonds angle_bin_center=1.7453292519943295 angle_bin_width=0.01 output_file=branchk10i0_bonds4.txt trials_per_write=100000.0
CheckEnergy decimal_places=8 trials_per_update=100000.0
NumParticles output_file=branchk10i0_bonds4_num.txt trials_per_write=100000.0
Run num_trials=1000000.0
exit: 0
Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!