BUG: ClusterGrowth: LAMMPS handle coordinates slightly outside box

[kay-ro]

Description:

Since lammps may be storing particle positions that are slightly outside the simulation box, we have to handle coordinates that are slightly outside box_boundary.

example data:
gyrationstuffex.zip

Code for reproduction:

test2clg = amep.evaluate.ClusterGrowth(traj1, mode="largest",skip=0.2,min_size=200)

Error message:

---

ValueError Traceback (most recent call last)
Cell In[20], line 14
10 #test1clg = amep.evaluate.ClusterGrowth(traj1, mode="mean",skip=0.2,min_size=200)
11 #test3clg = amep.evaluate.ClusterGrowth(traj1, mode="weighted mean",skip=0.2,min_size=200)
13 print("LÜLÜLÜ")
---> 14 test2clg = amep.evaluate.ClusterGrowth(traj1, mode="largest",skip=0.2,min_size=200)
15 #axstest.plot(test1clg.times, test1clg.frames, color = "blue", label="mean")
16 #axstest.plot(test3clg.times, test1clg.frames, color = "green", label="weighted mean")
17 axstest.plot(test2clg.times, test2clg.frames, color = "orange", label="largest")

File ~/anaconda3/envs/new-amep/lib/python3.11/site-packages/amep/evaluate.py:3523, in ClusterGrowth.init(self, traj, skip, nav, min_size, ptype, ftype, mode, use_density, **kwargs)
3515 raise ValueError(
3516 "amep.evaluate.ClusterGrowth: Invalid value for mode. Got "
3517 f"{self.__mode}. Please choose one of ['average', 'mean', "
3518 "'weighted mean']."
3519 )
3521 if isinstance(self.__traj, ParticleTrajectory):
3522 # calculation for particles
-> 3523 self.__frames, self.__avg, self.__indices = average_func(
3524 self.__compute_particles,
3525 traj,
3526 skip=self.__skip,
3527 nr=self.__nav,
3528 indices=True
3529 )
3530 print(self.__frames, self.__avg, self.__indices)
3531 elif isinstance(self.__traj, FieldTrajectory):
3532 # calculation for fields

File ~/anaconda3/envs/new-amep/lib/python3.11/site-packages/amep/utils.py:142, in average_func(func, data, skip, nr, indices, **kwargs)
140 evaluated_indices = np.array(np.ceil(np.linspace(skip*N, N-1, nr)), dtype=int)
141 print(evaluated_indices)
--> 142 func_result = [func(x, **kwargs) for x in tqdm(data[evaluated_indices])]
143 evaluated = np.array(func_result)
144 if indices:

File ~/anaconda3/envs/new-amep/lib/python3.11/site-packages/amep/utils.py:142, in (.0)
140 evaluated_indices = np.array(np.ceil(np.linspace(skip*N, N-1, nr)), dtype=int)
141 print(evaluated_indices)
--> 142 func_result = [func(x, **kwargs) for x in tqdm(data[evaluated_indices])]
143 evaluated = np.array(func_result)
144 if indices:

File ~/anaconda3/envs/new-amep/lib/python3.11/site-packages/amep/evaluate.py:3570, in ClusterGrowth.__compute_particles(self, frame)
3556 r'''
3557 Calculation for a single frame.
3558
(...)
3567 Largest, mean, or weighted mean cluster size.
3568 '''
3569 # determine clusters for particles
-> 3570 sorted_clusters, idx = identify(
3571 frame.coords(ptype=self.__ptype),
3572 frame.box,
3573 **self.__kwargs
3574 )
3575 sizes = cluster_sizes(sorted_clusters)
3577 # get relevant values

File ~/anaconda3/envs/new-amep/lib/python3.11/site-packages/amep/cluster.py:977, in identify(coords, box_boundary, sizes, pbc, rmax)
912 r"""
913 Identify clusters from particle coordinates, and respective sizes.
914
(...)
969
970 """
971 # extracting relevant information, i.e, particles belong to the same
972 # cluster, if their distance is smaller than rmax; here we generate an
973 # array of shape (N,2) that contains particle pairs, where N is the number
974 # of pairs, the first column contains the index of the first particle i,
975 # the second column of the second particle j, where particle i and j are
976 # closer than rmax to each other
--> 977 pairs = find_pairs(
978 coords, box_boundary, pbc=pbc, rmax=rmax, sizes=sizes
979 )
980 # results list
981 clusters = []

File ~/anaconda3/envs/new-amep/lib/python3.11/site-packages/amep/pbc.py:750, in find_pairs(coords, box_boundary, ids, sizes, other_coords, other_ids, other_sizes, pbc, rmax)
747 if sizes is None:
748 if other_coords is None:
749 # create KDTree object
--> 750 tree = kdtree(coords, box_boundary, pbc=pbc)
752 # get list of pairs (particle indices)
753 pairs = np.asarray(tree.query_pairs(rmax, output_type='ndarray'))

File ~/anaconda3/envs/new-amep/lib/python3.11/site-packages/amep/pbc.py:651, in kdtree(coords, box_boundary, pbc)
648 coords[coords[:,1]==box[1],1]=0
649 coords[coords[:,2]==box[2],2]=0
--> 651 return KDTree(coords, boxsize=box)

File ~/anaconda3/envs/new-amep/lib/python3.11/site-packages/scipy/spatial/_kdtree.py:360, in KDTree.init(self, data, leafsize, compact_nodes, copy_data, balanced_tree, boxsize)
357 raise TypeError("KDTree does not work with complex data")
359 # Note KDTree has different default leafsize from cKDTree
--> 360 super().init(data, leafsize, compact_nodes, copy_data,
361 balanced_tree, boxsize)

File _ckdtree.pyx:587, in scipy.spatial._ckdtree.cKDTree.init()

ValueError: Some input data are greater than the size of the periodic box.

Python and AMEP versions:

1.0.2

Additional information:

Lammps:
https://docs.lammps.org/dump_modify.html
The pbc keyword applies to all the dump styles. As explained on the dump doc page, atom coordinates in a dump file may be slightly outside the simulation box. This is because periodic boundary conditions are enforced only on timesteps when neighbor lists are rebuilt, which will not typically coincide with the timesteps dump snapshots are written. If the setting of this keyword is set to yes, then all atoms will be remapped to the periodic box before the snapshot is written, then restored to their original position. If it is set to no they will not be. The no setting is the default because it requires no extra computation.

this may not be specific to evaluate.ClusterGrowth! keep an eye out for this!

How did you install AMEP?

None

[kay-ro]

a quick investigation shows that this might be caused by floating point operation errors. pbc.fold() returns a np.float32 array but if a function were to use a 64 bit float, this might cause issues.

[kay-ro]

Order of execution of moving coordinates and folding them in the box is important. issue fixed in branch gromacs_reader. Will be merged soon.