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 unittest
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 add_particles_of_type0 1 particle_type0 /feasst/particle/propane.fstprt"
params['potential'] = "Potential Model LennardJones"
params['branch'] = 'particle_type 0 angle true mobile_site 2 anchor_site 1 anchor_site2 0'
params['grow'] = """particle_type 0 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 0 angle true mobile_site 2 anchor_site 1 anchor_site2 0'
params['grow'] = """particle_type 0 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 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"""
params['grow'] = """particle_type 0 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 0 branch true mobile_site 2 mobile_site2 3 anchor_site 0 anchor_site2 1'
params['grow'] = """particle_type 0 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 add_particles_of_type0 1 particle_type0 /feasst/particle/2-methylpropane.fstprt"
params['potential'] = "Potential Model LennardJones"
elif params['molecule'] == 'floppy_branch':
params['config'] = """Configuration cubic_side_length 30 add_particles_of_type0 1 particle_type0 {prefix}.fstprt""".format(**params)
params['potential'] = "Potential Model HardSphere"
if 'intra' in params:
if params['intra'] == 1:
params['potential'] += "Potential Model HardSphere VisitModel VisitModelIntraMap exclude_bonds true"
with open(params['prefix']+'.fstprt', 'w') as file1:
file1.write("""
4 sites
3 bonds
3 angles
1 site types
1 bond types
1 angle types
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 add_particles_of_type0 1 particle_type0 {prefix}.fstprt""".format(**params)
params['potential'] = "Potential Model HardSphere"
with open(params['prefix']+'.fstprt', 'w') as file1:
file1.write("""
3 sites
2 bonds
1 angles
1 site types
1 bond types
1 angle types
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 add_particles_of_type0 1 particle_type0 {prefix}.fstprt""".format(**params)
params['potential'] = "Potential Model HardSphere"
with open(params['prefix']+'.fstprt', 'w') as file1:
file1.write("""
4 sites
3 bonds
2 angles
1 dihedrals
1 site types
1 bond types
1 angle types
1 dihedral types
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('script.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
Movie trials_per_write {trials_per} output_file {prefix}.xyz
Log trials_per_write {trials_per} output_file {prefix}.txt include_bonds true
Energy trials_per_write {trials_per} output_file {prefix}_en.txt
AnalyzeBonds trials_per_write {trials_per} output_file {prefix}_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} 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)
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('script.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
Movie trials_per_write {trials_per} output_file {prefix}.xyz
Log trials_per_write {trials_per} output_file {prefix}.txt include_bonds true
Energy trials_per_write {trials_per} output_file {prefix}_en.txt
AnalyzeBonds trials_per_write {trials_per} output_file {prefix}_bonds{num_jacobian_gaussian}.txt angle_bin_width {angle_bin_width} angle_bin_center {angle_bin_center}
CheckEnergy trials_per_update {trials_per} tolerance 1e-8
NumParticles trials_per_write {trials_per} output_file {prefix}_bonds{num_jacobian_gaussian}_num.txt
#CriteriaUpdater trials_per_update 1e4
#CriteriaWriter trials_per_write {trials_per} output_file {prefix}_bonds{num_jacobian_gaussian}_crit.txt
#Run until_criteria_complete true
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)
[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
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 878356495
# initializing random number generator with seed: 878356495
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 floppy_angle.fstprt
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 tolerance 1e-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
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 996885932
# initializing random number generator with seed: 996885932
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 floppy_dihedral.fstprt
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 tolerance 1e-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
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 890422165
# initializing random number generator with seed: 890422165
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 /home/hwh/feasst/particle/propane.fstprt
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 tolerance 1e-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_bonds.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
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 775102168
# initializing random number generator with seed: 775102168
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 /home/hwh/feasst/particle/2-methylpropane.fstprt
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 tolerance 1e-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
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 79821632
# initializing random number generator with seed: 79821632
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk0i0.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 733717029
# initializing random number generator with seed: 733717029
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk0i1.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 273331985
# initializing random number generator with seed: 273331985
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk0i0.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 99495162
# initializing random number generator with seed: 99495162
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk0i1.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 610952794
# initializing random number generator with seed: 610952794
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk0i0.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 614753912
# initializing random number generator with seed: 614753912
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk0i1.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 830208665
# initializing random number generator with seed: 830208665
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk10i0.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 59554301
# initializing random number generator with seed: 59554301
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk10i1.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 168818909
# initializing random number generator with seed: 168818909
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk10i0.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 885021345
# initializing random number generator with seed: 885021345
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk10i1.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 575098948
# initializing random number generator with seed: 575098948
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk10i0.fstprt
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 tolerance 1e-8 trials_per_update 100000.0
Run num_trials 100000.0
exit: 0
running
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 786417294
# initializing random number generator with seed: 786417294
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk10i1.fstprt
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 tolerance 1e-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
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 475701026
# initializing random number generator with seed: 475701026
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk10i0.fstprt
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 tolerance 1e-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
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 210635937
# initializing random number generator with seed: 210635937
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk10i0.fstprt
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 tolerance 1e-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
# FEASST version: v0.24.5-13-g5252255c16-dirty-hwh/pybind11
MonteCarlo
RandomMT19937 seed 17488583
# initializing random number generator with seed: 17488583
Configuration add_particles_of_type0 1 cubic_side_length 30 particle_type0 branchk10i0.fstprt
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 tolerance 1e-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!