Clear Filters
Clear Filters

Efficient management of interacting atoms based on their mutual distance

5 views (last 30 days)
Dear all
Currently, I am dealing with a dataset, which I attach here in the "File.txt", file, where the first column is an atomic label ranging from 0 to 95900 (that it is, to the number of atoms-1). Each row represents the spatial coordinates of each i-th atom of the system, denoting the second, third, and fourth columns the coordinate values, in Angstroms, along x-, y-, and z-th directions. It is important to note that the fourth column, that it is, the z-th spatial coordinate, only can have two different values, 3.27645 and 9.82515. So basically, I have two monolayers at two different heights.
My goal is, for each atom, to loop over the first ten nearest interlayer atoms. That it is, for the i-th atom, to loop over the atoms with the lowest deviations in the xy plane from (xi,yi) and with different value of the z-th component. Importantly, if, for example, during the looping one founds that atom with id 4 and 560 have in common the atom with 1100 as one of the ten nearest neighbors to them, it would be very desirable that when looping in for 1100, the program knows, from the previous calculations, that atoms with label 4 and 560 are indeed two of the ten nearest neighbors to 1100.
I had in mind something like:
file=readmatrix('File.txt');
counter=0;
for i=1:length(file(:,1))
indices=file(file(:,4)~=file(i,4));
atoms=file(indices,1:3);
for j=1:length(indices)
distances(j)=norm(atoms(j,2:3)-file(i,2:3)); % Angstroms
end
sorted_distances=unique(distances); % Angstroms
interacting_distances=sorted_distances(1:10); % Angstroms
interacting_atoms=find(distances<=max(interacting_distances));
for j=1:length(interacting_atoms(:,1))
counter=counter+1;
interaction_list(counter,1)=file(i,1); % i-th atom id
interaction_list(counter,2)=interacting_atoms(j,1); % j-th atom id
interaction_list(counter,3)=somevalue;
end
clear indices distances sorted_distances interacting_distances interacting_atoms
end
Obviously, this approach does not seem to be very fast/efficient. Also, this won't meet one of the requirements I explained above: that if one finds that i-th atom interacts with j-th atom, one knows that j-th atom interacts with i-th atom (double counting). In that case, I would need to save in the interaction_list(counter,3)=-somevalue. If, for example, during the looping atom with atom id 30500 it was already identified as one of the ten nearest interlayer neighbors of ten different atoms, it would be desirable that, when i=30500+1, MatLab does proceed to look for the ten nearest interlayer neighbors, since they have been already identified before!
Any idea?

Answers (1)

John D'Errico
John D'Errico on 31 Jul 2024
Edited: John D'Errico on 31 Jul 2024
No. It is NOT true that in a set, if atom x is one of the k closest atoms to atom y, that the converse is also true, i.e., that atom y is one of the k closest atoms to atom x. You seem to want to make that assumption, and it is patently false. We need only consider a trivial counter-exmple. This works in 1 dimension as a counter-example, but it would also apply to higher dimensions.
x = [0,10:15];
In the vector x, what are the three nearest neighbors to 0? Clearly that are the elements [10,11,12]. But 0 is NOT one of the three nearest neighbors to any of the numbers 10, 11, or 12.
Nearest neighbor-ness does not commute.
Anyway, the simple solution would seem to be to use tools like knnsearch, or, better yet, KDtreeSearcher. Then you will be able to create a kd-tree, and then search it repeatedly.
help KDTreeSearcher
KDTreeSearcher Neighbor search object using a kd-tree. A KDTreeSearcher object performs KNN (K-nearest-neighbor) search or radius search using a kd-tree. You can create a KDTreeSearcher object based on X using either of the following syntaxes: CREATENS function: NS = CREATENS(X,'NSMethod','kdtree') KDTreeSearcher constructor: NS = KDTreeSearcher(X) Rows of X correspond to observations and columns correspond to variables. When one of the above methods creates a KDTreeSearcher object, it creates and saves a kd-tree based on X. The KNNSEARCH method and RANGESEARCH use the kd-tree to find requested points in X to the query points. KDTreeSearcher properties: X - Data used to create the KDTreeSearcher object Distance - The distance metric DistParameter - The additional parameter for the distance metric BucketSize - Bucket size of each leaf node in the kd-tree KDTreeSearcher methods: KDTreeSearcher - Construct a KDTreeSearcher object knnsearch - Find nearest neighbors for query points rangesearch - Find points within a given radius for query points. Example: % Create a KDTreeSearcher object for data X with the 'euclidean' % distance. X = randn(100,5); ns = createns(X,'nsmethod','kdtree'); % Find 5 nearest neighbors in X and the corresponding distance % values for each point in Y. Y = randn(25, 5); [idx, dist] = knnsearch(ns,Y,'k',5); % Find the points in X whose distance are not greater than 1.3 to % the points in Y. idx = rangesearch(ns,Y,1.3); See also CREATENS, ExhaustiveSearcher, STATS/KNNSEARCH, STATS/RANGESEARCH, hnswSearcher. Documentation for KDTreeSearcher doc KDTreeSearcher
  2 Comments
Richard Wood
Richard Wood on 1 Aug 2024
Hi @John D'Errico, thank you for your answer, you are totally right, my bad. The method that you propose seems very interesting, do you have in mind an efficient way of restrict the search for neighbors between atoms of different layers (neglecting the search for intralayer ones, that it is, with the same height along the z-th spatial direction)?
John D'Errico
John D'Errico on 1 Aug 2024
If you want to search within a plane, then build a series of kd-trees, one for each plane. You may need to reduce your data then to 2-dimensions, but at that point, the problem is essenitally already 2-dimensional.

Sign in to comment.

Categories

Find more on Particle & Nuclear Physics in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!