From fecf8de08f12e3c14e7fa4a944c2212372b24abf Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Tue, 12 Mar 2024 16:16:43 -0600 Subject: [PATCH 01/10] start tensile test method --- cmeutils/dynamics.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/cmeutils/dynamics.py b/cmeutils/dynamics.py index 32b6a4a..d6f15e0 100644 --- a/cmeutils/dynamics.py +++ b/cmeutils/dynamics.py @@ -8,6 +8,24 @@ from cmeutils import gsd_utils +def tensile_test(gsd_file, period, tensile_axis, ref_energy=1, ref_distance=1): + # Get initial box info and initial stress average + tensor_index_map = {0: 0, 1: 3, 2: 5} + with gsd.hoomd.open(gsd_file) as traj: + n_frames = len(traj) + # init_snap = traj[0] + # init_length = init_snap.configuration.box[tensile_axis] + # Store relevant stress tensor value for each frame + frame_stress_data = np.zeros(n_frames) + for idx, snap in enumerate(traj): + frame_stress_data[idx] = snap.log[ + "md/compute/ThermodynamicQuantities/pressure_tensor" + ][tensor_index_map[tensile_axis]] + + # window_means = np.zeros(n_frames // period) + # window_stds = np.zeros(n_frames // period) + + def msd_from_gsd( gsdfile, atom_types="all", start=0, stop=-1, msd_mode="window" ): From 793e69e533d7b922a3e4dd6e772a52ab8e0e1b36 Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Tue, 12 Mar 2024 20:55:43 -0600 Subject: [PATCH 02/10] use box values to find windows --- cmeutils/dynamics.py | 38 +++++++++++++++++++++++++++++++++----- 1 file changed, 33 insertions(+), 5 deletions(-) diff --git a/cmeutils/dynamics.py b/cmeutils/dynamics.py index d6f15e0..1064829 100644 --- a/cmeutils/dynamics.py +++ b/cmeutils/dynamics.py @@ -4,26 +4,54 @@ import gsd import gsd.hoomd import numpy as np +import scipy from cmeutils import gsd_utils -def tensile_test(gsd_file, period, tensile_axis, ref_energy=1, ref_distance=1): +def tensile_test( + gsd_file, + period, + tensile_axis, + ref_energy=1, + ref_distance=1, + enforce_sampling=False, +): # Get initial box info and initial stress average tensor_index_map = {0: 0, 1: 3, 2: 5} with gsd.hoomd.open(gsd_file) as traj: n_frames = len(traj) - # init_snap = traj[0] - # init_length = init_snap.configuration.box[tensile_axis] + init_snap = traj[0] + init_length = init_snap.configuration.box[tensile_axis] + print(init_length) # Store relevant stress tensor value for each frame frame_stress_data = np.zeros(n_frames) + frame_box_data = np.zeros(n_frames) for idx, snap in enumerate(traj): frame_stress_data[idx] = snap.log[ "md/compute/ThermodynamicQuantities/pressure_tensor" ][tensor_index_map[tensile_axis]] + frame_box_data[idx] = snap.configuration.box[tensile_axis] + + # Perform stress sampling + window_means = np.zeros(n_frames // period) + window_stds = np.zeros(n_frames // period) + window_sems = np.zeros(n_frames // period) + box_lengths = set(frame_box_data) + for idx, box_length in box_lengths: + indices = np.where(frame_box_data == box_length)[0] + stress = frame_stress_data[indices] + if enforce_sampling: # Use cmeutils.sampling, throw error? + pass + else: # Use use the last half of the stress values + cut = -len(stress) // 2 + avg_stress = np.mean(stress[cut:]) + std_stress = np.std(stress[cut:]) + sem_stress = scipy.stats.sem(stress[cut:]) - # window_means = np.zeros(n_frames // period) - # window_stds = np.zeros(n_frames // period) + window_means[idx] = avg_stress + window_stds[idx] = std_stress + window_sems[idx] = sem_stress def msd_from_gsd( From 056061159e42868a0b6d5b8d816ecb3d59478fe9 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Tue, 12 Mar 2024 22:45:33 -0600 Subject: [PATCH 03/10] update workflow, return values --- cmeutils/dynamics.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/cmeutils/dynamics.py b/cmeutils/dynamics.py index 1064829..9b31f58 100644 --- a/cmeutils/dynamics.py +++ b/cmeutils/dynamics.py @@ -11,7 +11,6 @@ def tensile_test( gsd_file, - period, tensile_axis, ref_energy=1, ref_distance=1, @@ -23,7 +22,6 @@ def tensile_test( n_frames = len(traj) init_snap = traj[0] init_length = init_snap.configuration.box[tensile_axis] - print(init_length) # Store relevant stress tensor value for each frame frame_stress_data = np.zeros(n_frames) frame_box_data = np.zeros(n_frames) @@ -34,11 +32,13 @@ def tensile_test( frame_box_data[idx] = snap.configuration.box[tensile_axis] # Perform stress sampling - window_means = np.zeros(n_frames // period) - window_stds = np.zeros(n_frames // period) - window_sems = np.zeros(n_frames // period) box_lengths = set(frame_box_data) - for idx, box_length in box_lengths: + strain = np.zeros(len(box_lengths)) + window_means = np.zeros(len(box_lengths)) + window_stds = np.zeros(len(box_lengths)) + window_sems = np.zeros(len(box_lengths)) + for idx, box_length in enumerate(box_lengths): + strain[idx] = (box_length - init_length) / init_length indices = np.where(frame_box_data == box_length)[0] stress = frame_stress_data[indices] if enforce_sampling: # Use cmeutils.sampling, throw error? @@ -52,6 +52,7 @@ def tensile_test( window_means[idx] = avg_stress window_stds[idx] = std_stress window_sems[idx] = sem_stress + return strain, -window_means, window_stds, window_sems def msd_from_gsd( From 2a1d075fabaede4bab56a3b8a45864fb0ce6fb69 Mon Sep 17 00:00:00 2001 From: chrisjones4 Date: Wed, 13 Mar 2024 11:37:34 -0600 Subject: [PATCH 04/10] replace set with numpy.unique --- cmeutils/dynamics.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/cmeutils/dynamics.py b/cmeutils/dynamics.py index 9b31f58..7773d04 100644 --- a/cmeutils/dynamics.py +++ b/cmeutils/dynamics.py @@ -4,6 +4,7 @@ import gsd import gsd.hoomd import numpy as np +import unyt as u import scipy from cmeutils import gsd_utils @@ -12,8 +13,8 @@ def tensile_test( gsd_file, tensile_axis, - ref_energy=1, - ref_distance=1, + ref_energy=None, + ref_distance=None, enforce_sampling=False, ): # Get initial box info and initial stress average @@ -32,15 +33,20 @@ def tensile_test( frame_box_data[idx] = snap.configuration.box[tensile_axis] # Perform stress sampling - box_lengths = set(frame_box_data) + box_lengths = np.unique(frame_box_data) strain = np.zeros(len(box_lengths)) window_means = np.zeros(len(box_lengths)) window_stds = np.zeros(len(box_lengths)) window_sems = np.zeros(len(box_lengths)) + if ref_energy and ref_distance: + conv_factor = ref_energy.to("J/mol") / (ref_distance.to("m")**3 * u.Avogadros_number_mks) + conv_factor *= 1e-6 + else: + conv_factor = 1 for idx, box_length in enumerate(box_lengths): strain[idx] = (box_length - init_length) / init_length indices = np.where(frame_box_data == box_length)[0] - stress = frame_stress_data[indices] + stress = frame_stress_data[indices] * conv_factor if enforce_sampling: # Use cmeutils.sampling, throw error? pass else: # Use use the last half of the stress values @@ -52,6 +58,8 @@ def tensile_test( window_means[idx] = avg_stress window_stds[idx] = std_stress window_sems[idx] = sem_stress + + return strain, -window_means, window_stds, window_sems From 626bb7e01a428fea78a228ff0aee2d6afad0cc45 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Fri, 15 Mar 2024 22:43:49 -0600 Subject: [PATCH 05/10] update np methods --- cmeutils/dynamics.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/cmeutils/dynamics.py b/cmeutils/dynamics.py index 7773d04..3dfb830 100644 --- a/cmeutils/dynamics.py +++ b/cmeutils/dynamics.py @@ -34,10 +34,10 @@ def tensile_test( # Perform stress sampling box_lengths = np.unique(frame_box_data) - strain = np.zeros(len(box_lengths)) - window_means = np.zeros(len(box_lengths)) - window_stds = np.zeros(len(box_lengths)) - window_sems = np.zeros(len(box_lengths)) + strain = np.zeros_like(box_lengths, dtype=float) + window_means = np.zeros_like(box_lengths, dtype=float) + window_stds = np.zeros_like(box_lengths, dtype=float) + window_sems = np.zeros_like(box_lengths, dtype=float) if ref_energy and ref_distance: conv_factor = ref_energy.to("J/mol") / (ref_distance.to("m")**3 * u.Avogadros_number_mks) conv_factor *= 1e-6 @@ -59,7 +59,6 @@ def tensile_test( window_stds[idx] = std_stress window_sems[idx] = sem_stress - return strain, -window_means, window_stds, window_sems From 9f047dbbcdcf905e1ba0237cbe251603d6b8aeac Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 27 Mar 2024 11:52:40 -0600 Subject: [PATCH 06/10] add some exceptions for ref values --- cmeutils/dynamics.py | 29 +++++++++++++++++++++++------ 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/cmeutils/dynamics.py b/cmeutils/dynamics.py index 3dfb830..5ad06fe 100644 --- a/cmeutils/dynamics.py +++ b/cmeutils/dynamics.py @@ -4,8 +4,8 @@ import gsd import gsd.hoomd import numpy as np -import unyt as u import scipy +import unyt as u from cmeutils import gsd_utils @@ -17,6 +17,28 @@ def tensile_test( ref_distance=None, enforce_sampling=False, ): + if ref_energy or ref_distance: + if not all([ref_energy, ref_distance]): + raise RuntimeError( + "Both ref_energy and ref_distnace must be defined." + ) + if not ( + isinstance(ref_energy, u.array.unyt_quantity) + and isinstance(ref_distance, u.array.unyt_quantity) + ): + raise ValueError( + "ref_energy and ref_distance should be given as " + "unyt.array.unyt_quantity." + ) + # Units of Pa + conv_factor = ref_energy.to("J/mol") / ( + ref_distance.to("m") ** 3 * u.Avogadros_number_mks + ) + # Units of MPa + conv_factor *= 1e-6 + else: + conv_factor = 1 + # Get initial box info and initial stress average tensor_index_map = {0: 0, 1: 3, 2: 5} with gsd.hoomd.open(gsd_file) as traj: @@ -38,11 +60,6 @@ def tensile_test( window_means = np.zeros_like(box_lengths, dtype=float) window_stds = np.zeros_like(box_lengths, dtype=float) window_sems = np.zeros_like(box_lengths, dtype=float) - if ref_energy and ref_distance: - conv_factor = ref_energy.to("J/mol") / (ref_distance.to("m")**3 * u.Avogadros_number_mks) - conv_factor *= 1e-6 - else: - conv_factor = 1 for idx, box_length in enumerate(box_lengths): strain[idx] = (box_length - init_length) / init_length indices = np.where(frame_box_data == box_length)[0] From 2c8f8beaf6cf5d55f87f149f88c2fdc4e2909cc7 Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Thu, 28 Mar 2024 13:51:13 -0600 Subject: [PATCH 07/10] add bootstrap sampling option --- cmeutils/dynamics.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/cmeutils/dynamics.py b/cmeutils/dynamics.py index 5ad06fe..d1806a0 100644 --- a/cmeutils/dynamics.py +++ b/cmeutils/dynamics.py @@ -15,7 +15,7 @@ def tensile_test( tensile_axis, ref_energy=None, ref_distance=None, - enforce_sampling=False, + bootstrap_sampling=False, ): if ref_energy or ref_distance: if not all([ref_energy, ref_distance]): @@ -64,8 +64,22 @@ def tensile_test( strain[idx] = (box_length - init_length) / init_length indices = np.where(frame_box_data == box_length)[0] stress = frame_stress_data[indices] * conv_factor - if enforce_sampling: # Use cmeutils.sampling, throw error? - pass + if bootstrap_sampling: + bootstrap_means = [] + n_data_points = len(strain) + n_samples = 5 + window_size = n_data_points // 5 + for i in range(n_samples): + start = np.random.randint( + low=0, high=(n_data_points - window_size) + ) + window_sample = strain[ + start : start + window_size # noqa: E203 + ] + bootstrap_means.append(window_sample) + avg_stress = np.mean(bootstrap_means) + std_stress = np.std(bootstrap_means) + sem_stress = scipy.stats.sem(bootstrap_means) else: # Use use the last half of the stress values cut = -len(stress) // 2 avg_stress = np.mean(stress[cut:]) From 16ec68d98ad5501740217ec628529f93ef195584 Mon Sep 17 00:00:00 2001 From: chrisjones4 Date: Thu, 28 Mar 2024 14:17:58 -0600 Subject: [PATCH 08/10] fix variable used in averaging --- cmeutils/dynamics.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cmeutils/dynamics.py b/cmeutils/dynamics.py index d1806a0..a137ec8 100644 --- a/cmeutils/dynamics.py +++ b/cmeutils/dynamics.py @@ -65,18 +65,18 @@ def tensile_test( indices = np.where(frame_box_data == box_length)[0] stress = frame_stress_data[indices] * conv_factor if bootstrap_sampling: - bootstrap_means = [] - n_data_points = len(strain) + n_data_points = len(stress) n_samples = 5 window_size = n_data_points // 5 + bootstrap_means = [] for i in range(n_samples): start = np.random.randint( low=0, high=(n_data_points - window_size) ) - window_sample = strain[ + window_sample = stress[ start : start + window_size # noqa: E203 ] - bootstrap_means.append(window_sample) + bootstrap_means.append(np.mean(window_sample)) avg_stress = np.mean(bootstrap_means) std_stress = np.std(bootstrap_means) sem_stress = scipy.stats.sem(bootstrap_means) From 5a929374e0a5030bbd2eb7859788f188b6e6d6bd Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 4 Sep 2024 11:02:58 -0600 Subject: [PATCH 09/10] update frame to freud func --- cmeutils/gsd_utils.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/cmeutils/gsd_utils.py b/cmeutils/gsd_utils.py index 88893d3..5d3233f 100644 --- a/cmeutils/gsd_utils.py +++ b/cmeutils/gsd_utils.py @@ -24,9 +24,17 @@ def frame_to_freud_system(frame, ref_length=None): if ref_length is None: ref_length = 1 box = frame.configuration.box + freud_box = freud.box.Box( + Lx=box[0] * ref_length, + Ly=box[1] * ref_length, + Lz=box[2] * ref_length, + xy=box[3], + xz=box[4], + yz=box[5] + ) box[0:3] *= ref_length xyz = frame.particles.position * ref_length - return freud.locality.NeighborQuery.from_system(system=(box, xyz)) + return freud.locality.AABBQuery(box=freud_box, points=xyz) def get_type_position( From 69a63d457b0acb215c42a8ae31eca0ae5c4aff10 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 4 Sep 2024 17:03:29 +0000 Subject: [PATCH 10/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- cmeutils/gsd_utils.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/cmeutils/gsd_utils.py b/cmeutils/gsd_utils.py index 7575645..cc999a5 100644 --- a/cmeutils/gsd_utils.py +++ b/cmeutils/gsd_utils.py @@ -25,12 +25,12 @@ def frame_to_freud_system(frame, ref_length=None): ref_length = 1 box = frame.configuration.box freud_box = freud.box.Box( - Lx=box[0] * ref_length, - Ly=box[1] * ref_length, - Lz=box[2] * ref_length, - xy=box[3], - xz=box[4], - yz=box[5] + Lx=box[0] * ref_length, + Ly=box[1] * ref_length, + Lz=box[2] * ref_length, + xy=box[3], + xz=box[4], + yz=box[5], ) box[0:3] *= ref_length xyz = frame.particles.position * ref_length