diff --git a/phantomwalk/lib/create_system_dpd.py b/phantomwalk/lib/create_system_dpd.py index 06174e3..6537227 100644 --- a/phantomwalk/lib/create_system_dpd.py +++ b/phantomwalk/lib/create_system_dpd.py @@ -4,8 +4,14 @@ import hoomd import time -from dpd_utils import initialize_snapshot_rand_walk,check_bond_length_equilibration,check_inter_particle_distance,add_hoomd_writers,check_pair_energy,simulation_energy_end - +from phantomwalk.lib.dpd_utils import ( + initialize_snapshot_rand_walk, + check_bond_length_equilibration, + check_inter_particle_distance, + add_hoomd_writers, + simulation_energy_end, + generate_rdf +) def create_polymer_system_dpd( num_pol, @@ -28,7 +34,15 @@ def create_polymer_system_dpd( gsd_file_name='trajectory.gsd', gsd_write_freq=10, log_file_name='log.txt', - log_write_freq=10 + log_write_freq=10, + + # rdf parameters + rdf_start_idx=-1, + rdf_end_idx=None, + rdf_file_name='rdf.txt', + rdf_bins=300, + rdf_r_max=3.0, + rdf_r_min=0.0 ): ''' @@ -79,6 +93,22 @@ def create_polymer_system_dpd( the file that the .txt log file will be saved to log_write_freq : int, default 10 Period to write simulation data to the log file. + rdf_start_idx : int, default -1 + Which frame to use for the start of the RDF average. Supplying None or 0 + will use the first frame. Supports negative indices. + rdf_end_idx : int, default None + Which frame to use for the last frame of the RDF average. Supplying None + will use the last frame. Supports negative indices. + rdf_bins : int, default 300 + How many bins to use for the RDF histogram + rdf_r_max : float, default 3.0 + Maximum interparticle distance to include in the RDF calculation. If + this value is greater than half the box length, then that will be used + instead. + rdf_r_min : float, default 0.0 + Minimum interparticle distance to include in the RDF calculation. + rdf_file_name : str, default 'rdf.txt' + The file to output the RDF data into. ------- Returns @@ -146,6 +176,15 @@ def create_polymer_system_dpd( check_time = time.perf_counter() if (check_time-start_time) > loop_timeout: print("Simulation timed out") + generate_rdf( + start_idx=rdf_start_idx, + end_idx=rdf_end_idx, + bins=rdf_bins, + r_max=rdf_r_max, + r_min=rdf_r_min, + gsd_file_name=gsd_file_name, + output_file_name=rdf_file_name + ) return snap, 0 simulation.run(sim_steps_incr) for writer in simulation.operations.writers: @@ -157,6 +196,15 @@ def create_polymer_system_dpd( check_time = time.perf_counter() if (check_time-start_time) > loop_timeout: print("Simulation timed out") + generate_rdf( + start_idx=rdf_start_idx, + end_idx=rdf_end_idx, + bins=rdf_bins, + r_max=rdf_r_max, + r_min=rdf_r_min, + gsd_file_name=gsd_file_name, + output_file_name=rdf_file_name + ) return snap,0 simulation.run(sim_steps_incr) for writer in simulation.operations.writers: @@ -167,4 +215,13 @@ def create_polymer_system_dpd( end_time = time.perf_counter() total_time = end_time - start_time print("Total build and simulation time:", end_time - start_time) + generate_rdf( + start_idx=rdf_start_idx, + end_idx=rdf_end_idx, + bins=rdf_bins, + r_max=rdf_r_max, + r_min=rdf_r_min, + gsd_file_name=gsd_file_name, + output_file_name=rdf_file_name + ) return snap, total_time diff --git a/phantomwalk/lib/dpd_utils.py b/phantomwalk/lib/dpd_utils.py index f96d09f..3137db2 100644 --- a/phantomwalk/lib/dpd_utils.py +++ b/phantomwalk/lib/dpd_utils.py @@ -4,6 +4,7 @@ import hoomd import time from cmeutils.sampling import is_equilibrated +import csv def initialize_snapshot_rand_walk(num_pol, num_mon, density, bond_length=1.0, seed=1234): ''' @@ -225,3 +226,48 @@ def simulation_energy_end(A,r,r_cut,num_pol,num_mon,density,energy_idx=-1): else: return False +def generate_rdf( + start_idx = -1, + end_idx = None, + bins=300, + r_max=3.0, + r_min=0.0, + gsd_file_name="trajectory.gsd", + output_file_name="rdf.txt" +): + ''' + Generate a radial distribution function for the provided GSD file and save + the data to the provided output file. + + ---------- + Parameters + ---------- + start_idx : int, default -1 + Which frame to use for the start of the RDF average. Supplying None or 0 + will use the first frame. Supports negative indices. + end_idx : int, default None + Which frame to use for the last frame of the RDF average. Supplying None + will use the last frame. Supports negative indices. + bins : int, default 300 + How many bins to use for the histogram + r_max : float, default 3.0 + Maximum interparticle distance to include in the calculation. If this + value is greater than half the box length, then that will be used + instead. + r_min : float, default 0.0 + Minimum interparticle distance to include in the calculation. + gsd_file_name : str, default 'trajectory.gsd' + The file to read the trajectory data used for calculating the RDF. + output_file_name : str, default 'rdf.txt' + The file to output the RDF data into. + ''' + rdf = freud.density.RDF(bins=bins, r_max=r_max, r_min=r_min) + traj = gsd.hoomd.open(gsd_file_name, 'r') + for frame in traj[start_idx:end_idx]: + rdf.compute(system=frame, reset=False) + + np.savetxt( + fname=output_file_name, + X=np.column_stack((rdf.bin_centers, rdf.rdf)), + header="r rdf" + )