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=6
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=6
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/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=51376018
# Initializing random number generator with seed: 51376018
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=6 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/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=863675017
# Initializing random number generator with seed: 863675017
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=6 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/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=627491116
# Initializing random number generator with seed: 627491116
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=6 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/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=592833805
# Initializing random number generator with seed: 592833805
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=6 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/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=656680569
# Initializing random number generator with seed: 656680569
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=791754052
# Initializing random number generator with seed: 791754052
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=570428588
# Initializing random number generator with seed: 570428588
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=633455370
# Initializing random number generator with seed: 633455370
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=192639018
# Initializing random number generator with seed: 192639018
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=267590325
# Initializing random number generator with seed: 267590325
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=385578983
# Initializing random number generator with seed: 385578983
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=979224303
# Initializing random number generator with seed: 979224303
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=994542371
# Initializing random number generator with seed: 994542371
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=464190821
# Initializing random number generator with seed: 464190821
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=903843675
# Initializing random number generator with seed: 903843675
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=6 trials_per_update=100000.0
Run num_trials=100000.0
exit: 0
running
# Usage: /home/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=881915403
# Initializing random number generator with seed: 881915403
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=6 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/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=338534626
# Initializing random number generator with seed: 338534626
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=6 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/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=327568917
# Initializing random number generator with seed: 327568917
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=6 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/hwh/feasst/build/bin/fst < file.txt
FEASST version 0.25.16
MonteCarlo
RandomMT19937 seed=875465004
# Initializing random number generator with seed: 875465004
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=6 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!