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

Fixed rdf

Browse Source
This commit is contained in:
Sebastian Kloth 2024-01-23 11:01:43 +01:00
parent 5c17e04b38
commit 02fed343f0
2 changed files with 11 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/mdevaluate/coordinates.py
View File

@ -646,7 +646,7 @@ def next_neighbors(
return distances_new, indices_new return distances_new, indices_new
else: else:
atoms_pbc, atoms_pbc_index = pbc_points( atoms_pbc, atoms_pbc_index = pbc_points(
query_atoms, box, thickness=distance_upper_bound + 0.1, index=True, **kwargs atoms, box, thickness=distance_upper_bound + 0.1, index=True, **kwargs
) )
tree = KDTree(atoms_pbc) tree = KDTree(atoms_pbc)
distances, indices = tree.query( distances, indices = tree.query(
@ -692,12 +692,16 @@ def number_of_neighbors(
elif not distinct: elif not distinct:
dnn = 1 dnn = 1
print(atoms[:5])
print(query_atoms[:5])
print(" ")
box = atoms.box box = atoms.box
if np.all(np.diag(np.diag(box)) == box): if np.all(np.diag(np.diag(box)) == box):
atoms = atoms % np.diag(box) atoms = atoms % np.diag(box)
tree = KDTree(atoms, boxsize=np.diag(box)) tree = KDTree(atoms, boxsize=np.diag(box))
else: else:
atoms_pbc = pbc_points(query_atoms, box, thickness=r_max + 0.1, **kwargs) atoms_pbc = pbc_points(atoms, box, thickness=r_max + 0.1, **kwargs)
tree = KDTree(atoms_pbc) tree = KDTree(atoms_pbc)
num_of_neighbors = tree.query_ball_point(query_atoms, r_max, return_length=True) num_of_neighbors = tree.query_ball_point(query_atoms, r_max, return_length=True)

7
src/mdevaluate/distribution.py
View File

@ -126,6 +126,9 @@ def rdf(
particles_in_volume = int( particles_in_volume = int(
np.max(number_of_neighbors(atoms_a, query_atoms=atoms_b, r_max=bins[-1])) * 1.1 np.max(number_of_neighbors(atoms_a, query_atoms=atoms_b, r_max=bins[-1])) * 1.1
) )
print(atoms_a[:5])
print(atoms_b[:5])
print(" ")
distances, indices = next_neighbors( distances, indices = next_neighbors(
atoms_a, atoms_a,
atoms_b, atoms_b,
@ -147,9 +150,9 @@ def rdf(
distances = [d for dist in distances for d in dist] 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_b)
hist = hist / (4 / 3 * np.pi * bins[1:] ** 3 - 4 / 3 * np.pi * bins[:-1] ** 3) hist = hist / (4 / 3 * np.pi * bins[1:] ** 3 - 4 / 3 * np.pi * bins[:-1] ** 3)
n = len(atoms_b) / np.prod(np.diag(atoms_b.box)) n = len(atoms_a) / np.prod(np.diag(atoms_a.box))
hist = hist / n hist = hist / n
return hist return hist