IPKM/mdevaluate
6
1
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Changed input of hbonds to atom_indices

Browse Source
This commit is contained in:
Sebastian Kloth 2024-01-16 13:38:13 +01:00
parent 9200d31af6
commit 95b46c43be
1 changed files with 27 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
src/mdevaluate/distribution.py
View File

@ -12,6 +12,7 @@ from .coordinates import (
Coordinates, Coordinates,
CoordinateFrame, CoordinateFrame,
next_neighbors, next_neighbors,
number_of_neighbors,
) )
from .autosave import autosave_data from .autosave import autosave_data
from .pbc import pbc_diff, pbc_points from .pbc import pbc_diff, pbc_points
@ -57,7 +58,7 @@ def time_distribution(
function: Callable, function: Callable,
coordinates: Coordinates, coordinates: Coordinates,
coordinates_b: Optional[Coordinates] = None, coordinates_b: Optional[Coordinates] = None,
skip: float = 0.1, skip: float = 0,
segments: int = 100, segments: int = 100,
) -> Tuple[NDArray, List]: ) -> Tuple[NDArray, List]:
""" """
@ -88,7 +89,7 @@ def time_distribution(
return times, result return times, result
def gr( def rdf(
atoms_a: CoordinateFrame, atoms_a: CoordinateFrame,
atoms_b: Optional[CoordinateFrame] = None, atoms_b: Optional[CoordinateFrame] = None,
bins: Optional[ArrayLike] = None, bins: Optional[ArrayLike] = None,
@ -122,10 +123,9 @@ def gr(
if bins is None: if bins is None:
bins = np.arange(0, 1, 0.01) bins = np.arange(0, 1, 0.01)
box = atoms_b.box particles_in_volume = int(
n = len(atoms_a) / np.prod(np.diag(box)) np.max(number_of_neighbors(atoms_a, query_atoms=atoms_b, r_max=bins[-1])) * 1.1
V = 4 / 3 * np.pi * bins[-1] ** 3 )
particles_in_volume = int(n * V * 1.1)
distances, indices = next_neighbors( distances, indices = next_neighbors(
atoms_a, atoms_a,
atoms_b, atoms_b,
@ -144,7 +144,7 @@ def gr(
) )
distances = np.concatenate(new_distances) distances = np.concatenate(new_distances)
else: else:
distances = distances.flatten() distances = [d for dist in distances for d in dist]
hist, bins = np.histogram(distances, bins=bins, range=(0, bins[-1]), density=False) hist, bins = np.histogram(distances, bins=bins, range=(0, bins[-1]), density=False)
hist = hist / len(atoms_a) hist = hist / len(atoms_a)
@ -310,26 +310,27 @@ def next_neighbor_distribution(
def hbonds( def hbonds(
D: CoordinateFrame, atoms: CoordinateFrame,
H: CoordinateFrame, donator_indices: ArrayLike,
A: CoordinateFrame, hydrogen_indices: ArrayLike,
acceptor_indices: ArrayLike,
DA_lim: float = 0.35, DA_lim: float = 0.35,
HA_lim: float = 0.35, HA_lim: float = 0.35,
min_cos: float = np.cos(30 * np.pi / 180), max_angle_deg: float = 30,
full_output: bool = False, full_output: bool = False,
) -> Union[NDArray, tuple[NDArray, NDArray, NDArray]]: ) -> Union[NDArray, tuple[NDArray, NDArray, NDArray]]:
""" """
Compute h-bond pairs Compute h-bond pairs
Args: Args:
D: Set of coordinates for donators. atoms: Set of all coordinates for a frame.
H: Set of coordinates for hydrogen atoms. Should have the same donator_indices: Set of indices for donators.
hydrogen_indices: Set of indices for hydrogen atoms. Should have the same
length as D. length as D.
A: Set of coordinates for acceptors. acceptor_indices: Set of indices for acceptors.
DA_lim (opt.): Minimum distance beteen donator and acceptor. DA_lim (opt.): Minimum distance between donator and acceptor.
HA_lim (opt.): Minimum distance beteen hydrogen and acceptor. HA_lim (opt.): Minimum distance between hydrogen and acceptor.
min_cos (opt.): Minimum cosine for the HDA angle. Default is max_angle_deg (opt.): Maximum angle in degree for the HDA angle.
equivalent to a maximum angle of 30 degree.
full_output (opt.): Returns additionally the cosine of the full_output (opt.): Returns additionally the cosine of the
angles and the DA distances angles and the DA distances
@ -362,6 +363,10 @@ def hbonds(
pairs[:, 1] = pind[pairs[:, 1]] pairs[:, 1] = pind[pairs[:, 1]]
return pairs return pairs
D = atoms[donator_indices]
H = atoms[hydrogen_indices]
A = atoms[acceptor_indices]
min_cos = np.cos(max_angle_deg * np.pi / 180)
box = D.box box = D.box
if len(D) <= len(A): if len(D) <= len(A):
pairs = dist_DltA(D, A, DA_lim) pairs = dist_DltA(D, A, DA_lim)