# Optimize molecule distance analyzing code

I have a very large dataset (31552000 lines) of xyz coordinates in the following format

1 2 3
4 5 6
7 8 9
. . .


I have to take a distance using the special method below.

Distance[{a_, b_, c_}, {d_, e_, f_}] :=
Sqrt[(If[Abs[a - d] >= (40/2), Abs[a - d] - 40, Abs[a - d]])^2 + (If[
Abs[b - e] >= (40/2), Abs[b - e] - 40, Abs[b - e]])^2 + (If[
Abs[c - f] >= (40/2), Abs[c - f] - 40, Abs[c - f]])^2]


Then I import the data.

data = Partition[
Partition[ReadList["input.txt", {Real, Real, Real}], 16], 38];


The formatting is kind of strange. Every 16 rows is one molecule, and every 38 molecules is one timestep. I take the distance between the 16th atom of each molecule and the 5th atom of each molecule.Then I select the distances that are less than 5.55 and determine the length of the resulting list. This is repeated for each of the 29,000 timesteps.

analysis =
Flatten[
Table[
Table[
Length[
Select[
Table[
Distance[data[[r, y, 16]], data[[r, x, 5]]],
{x, 1, 38}],
# <= 5.55 &]
],
{y, 1, 38}],
{r, 1, 29000}]
];


This last section is my most computationally intensive part. For 29000 timesteps and 38 molecules, it takes 40 minutes to process fully. It also takes too much memory (16+ gigs per kernel) to parallelize. Is there any other method that will improve the performance? I have tried using compile, but I realized that Table, the biggest bottleneck, is already complied to machine code.

Below is an example of a dataset that takes my computer 2 minutes to complete with the analysis code. It is scalable to larger timesteps by changing 4000 to larger numbers.

data = Partition[
Partition[Partition[Table[RandomReal[{0, 40}], (3*16*38*4000)], 3],
16], 38]

• An example dataset could help to develop potential solutions. May 31, 2021 at 19:07
• I added the example dataset at the bottom of my post. Thanks! May 31, 2021 at 19:27

If we reduce the dataset to work with (you are only interested in the 5th an 16th atom and this will help your RAM problems) and compile the whole algorithm into something which is still readable but more compact, then we land an around x27 speed improvement.

dataDistances = Compile[{{dat, _Real, 4}}, Module[
{m, count, res, a, d, b, e, c, f, dist},
m = Length[dat];
res = Table[0, {m*38}];
Do[
count = 0;
Do[
d = dat[[r, x, 1, 1]];
e = dat[[r, x, 1, 2]];
f = dat[[r, x, 1, 3]];
a = dat[[r, y, 2, 1]];
b = dat[[r, y, 2, 2]];
c = dat[[r, y, 2, 3]];
dist = (If[Abs[a - d] >= (40/2), Abs[a - d] - 40,
Abs[a - d]])^2 + (If[Abs[b - e] >= (40/2), Abs[b - e] - 40,
Abs[b - e]])^2 + (If[Abs[c - f] >= (40/2),
Abs[c - f] - 40, Abs[c - f]])^2;
If[dist <= 5.55^2, count++];
, {x, 1, 38}];
res[[38*(r - 1) + y]] = count;
, {r, 1, m}, {y, 1, 38}];
res
],CompilationTarget -> "C", Parallelization -> True];

rdata = data[[All, All, {5, 16}, All]];


I'm pretty sure, if one would attempt to calculate the algorithm through the flat dataset, or just an $$n\times 3$$ dataset, that we can get another factor two or higher in improvement.

Edit: Inserted CompilationTarget -> "C", Parallelization -> True

Edit 2: Changed from ConstantArray to Table. See comments

• What is the purpose of the constant array? This code certainly cuts the time faster than I thought possible. I am confused on how its counting and the res[[38*(r - 1) + y]] portion of the code. Jun 1, 2021 at 2:31
• ConstantArray is not compilable, can be replace with Table
– I.M.
Jun 1, 2021 at 6:28
• @I.M. i noticed the MainEvaluate but it was only one per invocation so I figured it'd be fine. I am seriously surprised to see that Table compiles when ConstantArray does not. Looks like an oversight to me. However, thanks for noticing, I changed it. Jun 1, 2021 at 10:11
• @bandannadata The constant array will store the resulting solution and will be populated while the algorithm runs. The counting is done in the inner two Do-loops. It iterates over your 38x38 combinations and counts a counter up whenever your distance condition holds true. The res[[38*(r-1)+y]] is a flattened array access. I could have made the result array two dimensional $m\times 38$ and store the results as res[[r,y]] and flatten it afterward. But since we are only interested in the flattened part, i wrote the algorithm to already store it like that. Jun 1, 2021 at 10:15
• It is very nice approach (+1), even that this code is not convertible to other problems. Jun 1, 2021 at 12:04

We can reduce computational time by 6 times with using compiled version of Distanceas follows

disc = Compile[{{r1, _Real, 1}, {r2, _Real,
1}}, \[Sqrt](If[Abs[r1[[1]] - r2[[1]]] >= 20,
Abs[r1[[1]] - r2[[1]]] - 40, Abs[r1[[1]] - r2[[1]]]]^2 +
If[Abs[r1[[2]] - r2[[2]]] >= 20, Abs[r1[[2]] - r2[[2]]] - 40,
Abs[r1[[2]] - r2[[2]]]]^2 +
If[Abs[r1[[3]] - r2[[3]]] >= 20, Abs[r1[[3]] - r2[[3]]] - 40,
Abs[r1[[3]] - r2[[3]]]]^2), CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True]


For data in the form

data = Partition[
Partition[Partition[Table[RandomReal[{0, 40}], (3*16*38*4000)], 3],
16], 38];


analysis =
Flatten[Table[
Table[Length[
Select[Table[
Distance[data[[r, y, 16]], data[[r, x, 5]]], {x, 1, 38}], # <=
5.55 &]], {y, 1, 38}], {r, 1, 1000}]]; // AbsoluteTiming

Out[]= {16.8741, Null}


and for compiled version

Flatten[Table[
Table[Table[
Length[Select[
Table[disc[data[[r, y, 16]], data[[r, x, 5]]], {x, 1,
38}], # <= 5.55 &]], {y, 1, 38}], {r, s, s + 100}], {s,
10}]]; // AbsoluteTiming

Out[]= {2.51888, Null}


Pay attention that we divide last table into 10 parts to make computation more effective. I think that computation time can be also reduced in the part of Select.

• How did you make the disc function? Im new to programing, so are there any resources or methods to convert conventional functions to compiled ones? May 31, 2021 at 23:11
• @bandannadata First you can look tutorial in the Help section. Also this forum is a good one resource, search Compile in the Search on Mathematica window mathematica.stackexchange.com/questions/tagged/compile Jun 1, 2021 at 11:39