Efficient management of interacting atoms based on their mutual distance

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?

Answer 1

[John D'Errico]

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

[Roderick]

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]

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.