Blocking in AMDAT does not use all possible pairs of frames in trajectory

[pkawak]

Basically, AMDAT does not use all frame pairs possible for the latter parts of the trajectory, i.e., the pairs that are larger than the block size.

To demonstrate this, I have added some print statements in src/mean_square_displacement.cpp:

void Mean_Square_Displacement::postprocess_list()
{
    for(int timeii=0;timeii<n_times;timeii++)
    {
          msd[timeii] /= float(weighting[timeii]);
          std::cout << timeii << " " << weighting[timeii] << "\n";
    }
}

I ran this for a trajectory with 3 blocks with 3 frames per block and I targeted 3 atoms per molecule with 20 molecules for an MSD computation. For reference, here is the <time_scheme> AMDAT line:

exponential 3 3 1.2 0 0 1.0

The results from the print statement above are (with two added columns:

0 0 180 3
1 1 180 3
2 2 180 3
3 3 120 2
4 6 60 1
5 9 60 1

The first column is the time delay index, the second column is the actual time delay from the MSD file, the third column is the weighting, and the fourth column is the weighting divided by the number of atoms (3*20).

The weighting is expected to be equal to:

Weighting = Number of Pairs * Number of Atoms

Therefore, the third column is the number of pairs involved in the measurement at that time delay.

For the first time delay (0), the timestep pairs are (0, 0), (3, 3), (6, 6).
1: (0, 1), (3, 4), (6, 7).
2: (0, 2), (3, 5), (6, 8)
3: (0, 3), (3, 6)
4: (0, 6)
5: (0, 9)

So, it seems to me that we are simply not using the last timestep in the computation, except for the highest time delay. With the settings above (exponential 3 3 1.2 0 0 1.0), AMDAT prints out: Reading a 10 timestep trajectory. This last timestep could technically be used in 3 and 4 ((6, 9) and (3, 9), respectively) but its not and this is the reason that we don't have an extra pair there.

[pkawak]

Adding these print statements (below) confirms the above:

void Mean_Square_Displacement::list_displacementkernel(int timegapii,int thisii, int nextii)
{
    currenttime=thisii;
    nexttime=nextii;
    currenttimegap=timegapii;

    std::cout << timegapii << " " << currenttime << " " << nexttime << "\n";

    #pragma omp atomic
    weighting[timegapii]+=trajectory_list->show_n_trajectories(thisii);
    (trajectory_list[0]).listloop(this,timegapii, thisii, nextii);
}

Output:

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

[pkawak]

For linear scaling with:

linear 10 13600

We get these unique pairs for each time delay:

0 | 0,0 1,1 2,2 3,3 4,4 5,5 6,6 7,7 8,8
1 | 0,1 1,2 2,3 3,4 4,5 5,6 6,7 7,8
2 | 0,2 1,3 2,4 3,5 4,6 5,7 6,8
3 | 0,3 1,4 2,5 3,6 4,7 5,8
4 | 0,4 1,5 2,6 3,7 4,8
5 | 0,5 1,6 2,7 3,8
6 | 0,6 1,7 2,8
7 | 0,7 1,8
8 | 0,8
9 | 0,9