Skip to content

Commit

Permalink
add diffraction_pattern method and unit test
Browse files Browse the repository at this point in the history
  • Loading branch information
chrisjonesBSU committed Oct 13, 2023
1 parent 8124963 commit 91af6e5
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 3 deletions.
53 changes: 50 additions & 3 deletions cmeutils/structure.py
Original file line number Diff line number Diff line change
Expand Up @@ -476,7 +476,7 @@ def structure_factor(
method="direct",
ref_length=None,
):
"""
"""Uses freud's diffraction module to find 1D structure factors.
Parameters
----------
Expand Down Expand Up @@ -508,6 +508,9 @@ def structure_factor(
freud.diffraction.StaticStructureFactorDebye
"""

if not ref_length:
ref_length = 1

if method.lower() == "direct":
sf = freud.diffraction.StaticStructureFactorDirect(
bins=bins, k_max=k_max, k_min=k_min
Expand All @@ -520,15 +523,59 @@ def structure_factor(
raise ValueError(
f"Optional methods are `debye` or `direct`, you chose {method}"
)
if not ref_length:
ref_length = 1
with gsd.hoomd.open(gsdfile, mode="r") as trajectory:
for frame in trajectory[start:stop]:
system = frame_to_freud_system(frame=frame, ref_length=ref_length)
sf.compute(system=system, reset=False)
return sf


def diffraction_pattern(
gsdfile,
views,
start=0,
stop=-1,
ref_length=None,
grid_size=1024,
output_size=None,
):
"""Uses freud's diffraction module to find 2D diffraction patterns.
Parameters
----------
gsdfile : str, required
File path to the GSD trajectory.
views : list, required
List of orientations (quarternions) to average over.
See cmeutils.structure.get_quarternions
start : int, default 0
Starting frame index for accumulating the Sq. Negative numbers index
from the end.
stop : int, optional default None
Final frame index for accumulating the Sq. If None, the last frame
will be used.
ref_length : float, optional, default None
Set a reference length to convert from reduced units to real units.
If None, uses 1 by default.
Returns
-------
freud.diffraction.DiffractionPattern
"""

if not ref_length:
ref_length = 1
dp = freud.diffraction.DiffractionPattern(
grid_size=grid_size, output_size=output_size
)
with gsd.hoomd.open(gsdfile) as trajectory:
for frame in trajectory[start:stop]:
system = frame_to_freud_system(frame=frame, ref_length=ref_length)
for view in views:
dp.compute(system=system, view_orientation=view, reset=False)
return dp


def get_centers(gsdfile, new_gsdfile):
"""Create a gsd file containing the molecule centers from an existing gsd
file.
Expand Down
6 changes: 6 additions & 0 deletions cmeutils/tests/test_structure.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
all_atom_rdf,
angle_distribution,
bond_distribution,
diffraction_pattern,
dihedral_distribution,
get_centers,
get_quaternions,
Expand Down Expand Up @@ -187,6 +188,11 @@ def test_angle_range_outside(self, p3ht_gsd):
theta_max=180,
)

def test_diffraction_pattern(self, gsdfile_bond):
views = get_quaternions(n_views=5)
dp = diffraction_pattern(gsdfile_bond, views=views)
assert isinstance(dp, freud.diffraction.DiffractionPattern)

def test_structure_factor_direct(self, gsdfile_bond):
sf = structure_factor(gsdfile_bond, k_min=0.2, k_max=5)
assert isinstance(sf, freud.diffraction.StaticStructureFactorDirect)
Expand Down

0 comments on commit 91af6e5

Please sign in to comment.