Single PBC distance calculation

[richardjgowers]

Originally started in #1188 , there's currently not a way to calculate a single pairwise distance, with PBC taken into account. (For many distances there's lib.distances.distance_array)

We need:

• MDAnalysis.lib.distances.distance(a, b, [box]) - which calculates distance between a and b
• box can be either an orthogonal or triclinic box
• Rather than place this in the above function, make a dedicated MDAnalysis.lib.distances.minimum_image(vector, box) which takes a vector and an orthogonal or triclinic box, and returns the vector with periodic boundary conditions applied
• then make MDAnalysis.lib.mdamath.angle also be able to use the minimum_image function

For now, this can just use a Python (numpy) implementation

[Shtkddud123]

Hi. I've got a semi-working version of this at the moment. It's nowhere near pretty, but it'll hopefully work.

Just one question - my knowledge of MD might be lacking here - when we are using a cubic box, we can just use the format::

if (periodic_x) then
dx = x(j) - x(i)
if (dx > x_size * 0.5) dx = dx - x_size
if (dx <= -x_size * 0.5) dx = dx + x_size
endif

But is there any reason why this wouldn't work on the triclinic/orthrhombic box?

Thanks.

[richardjgowers]

Hi Sang!

So your version works some of the time, but not always! If the vector is 2.2 x_size, you need to take away 2 box lengths.

Triclinic boxes are a bit tricky, but they can be expressed as a triangular matrix, so for the first dimension you have to alter all three dimensions to correct the first dimension, then alter 2 dimensions for the second etc. There's a link to the C code which we use for this below

https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/lib/include/calc_distances.h#L33

Don't worry about pretty code yet, we can fix that in review

[Shtkddud123]

Hi Richard,

Thanks for your prompt reply.

So, just another query - when you are asking for a single molecule distance, rather than an array, I've identified what atom a (reference) and atom b (target) are using the 'bynum' (index in the gromacs/namd/etc .. file). I figured that a person concerned with this would rather want to identify the specific atom rather than molecule name, but I was just passing this by you and ask whether this would be acceptable? Or in the input of the function, would you need me to extensively give an option for which molecule/residue etc?

Thanks

Sang

[richardjgowers]

@Shtkddud123 the signature should be something like this:

def distance(a, b):
    return ???

a1 = u.atoms[0]
a2 = u.atoms[4]

separation = distance(a1.position, a2.position, u.dimensions)

So you input 2 coordinates and a box, and it returns a distance

[Shtkddud123]

Dear Richard,

Ok, that seems to be pretty much what I had in mind. Ok, I'll get back to you.

Thanks
Sang

[Shtkddud123]

Hi, so, I've gotten a premilimary version, which seems to have some issues getting the right values out, which I presume is due to a bug somewhere in the code.

The format goes: distance(atom a, atom b, box):

• You just need to type in the index of the atom from a .gro file for atoms a and b, and depending on the universe dimensions, the code will decide what type of box it is.

• I've explicitly addressed the triclinic scenario, because it is the most complicated one, but I've put a fairly generic correction (or minimum image convention) for the orthogonal version as well. This hasn't been completely tested yet but as mentioned before, I just wanted to put it out here.

• The code loops over the x,y z of atom a and b and appliies the MIC.

In any case, any feedback and suggestions for corrections would be much appreciated.

def distance(atom1,atom2,box=None):

    name1 = "bynum" + " " + str(atom1) # concatenate string for first selection 
    name2 = "bynum" + " " + str(atom2) # concatenate string for the second selection 
    
    selection1 = u.select_atoms(name1) # select from universe
    selection2 = u.select_atoms(name2) # ditto          

    selection1_vec = selection1.positions[0] # As there is only one atom, we only need to extract the first row 
    selection2_vec = selection2.positions[0] 
    
    selection1_selection2_array =  [selection1_vec, selection2_vec]
    
     # --- For a triclinic box --- # 
    if box is not None:
        box_type = box_check(box) ## Check type of box and store it in box_type variable

        if box_type == 'tri_box':            
            for i in range(len(selection1_selection2_array)):
                if abs(selection1_selection2_array[i][2]) > u.dimensions[5]:
                   # Z axis minimum image - follows the LAMMPS convention
                    if selection1_selection2_array[i][2] < 0.0:
                           selection1_selection2_array[i][2] += u.dimensions[5]
                           selection1_selection2_array[i][1] += u.dimensions[4]
                           selection1_selection2_array[i][0] += u.dimensions[3]
                    else:
                           selection1_selection2_array[i][2] -= u.dimensions[5]
                           selection1_selection2_array[i][1] -= u.dimensions[4]
                           selection1_selection2_array[i][0] -= u.dimensions[3]
                    # Y axis minimum image 
                if abs(selection1_selection2_array[i][1]) > u.dimensions[4]:
                    if selection1_selection2_array[i][1] < 0.0:
                           selection1_selection2_array[i][1] += u.dimensions[4]
                           selection1_selection2_array[i][0] += u.dimensions[3]
                    else:                       
                           selection1_selection2_array[i][1] -= u.dimensions[4]
                           selection1_selection2_array[i][0] -= u.dimensions[3]
                     # X axis minimum image 
                if abs(selection1_selection2_array[i][0]) > u.dimensions[3]:
                    if selection1_selection2_array[i][0] < 0.0:
                           selection1_selection2_array[i][1] += u.dimensions[3]
                    else:                       
                           selection1_selection2_array[i][1] -= u.dimensions[3]
                   
        else:       # For non-trclinical boxes     
             for i in range(len(selection1_selection2_array)):
                    # X axis
                if abs(selection1_selection2_array[i][0]) < -u.dimensions[3]/2:
                           selection1_selection2_array[i][0] =  selection1_selection2_array[i][0] + u.dimensions[3]/2
                elif abs(selection1_selection2_array[i][0]) >= u.dimensions[3]/2:
                           selection1_selection2_array[i][0] =  selection1_selection2_array[i][0] - u.dimensions[3]/2
                         
                if abs(selection1_selection2_array[i][1]) < -u.dimensions[4]/2:
                           selection1_selection2_array[i][1] =  selection1_selection2_array[i][1] + u.dimensions[4]/2
                elif abs(selection1_selection2_array[i][1]) >= u.dimensions[4]/2:
                           selection1_selection2_array[i][1] =  selection1_selection2_array[i][1] - u.dimensions[4]/2
                                                     
                if abs(selection1_selection2_array[i][2]) < -u.dimensions[5]/2:
                           selection1_selection2_array[i][2] =  selection1_selection2_array[i][2] + u.dimensions[5]/2
                elif abs(selection1_selection2_array[i][2]) >= u.dimensions[5]/2:
                           selection1_selection2_array[i][2] =  selection1_selection2_array[i][2] - u.dimensions[5]/2
                       
    
                           separation = np.power((np.power((selection1_selection2_array[0][0] - selection1_selection2_array[1][0]),2) + np.power((selection1_selection2_array[1][0] - selection1_selection2_array[1][1]),2) + np.power((selection1_selection2_array[2][0] - selection1_selection2_array[2][1]),2)),2)

       return separation 

[mnmelo]

So, you should be presenting this code as a pull-request, on which it is much easier for everyone to comment. It doesn't matter if it's an initial non-working version. @orbeckst gave you a brief summary on how to accomplish this; we can guide you through it if you need help (but please keep that discussion on the mailing list, not this issue's conversation).

On your code: I can already see a problem, in that your choice of passing atom numbers will be somewhat undefined because the function doesn't know on which Universe to look for them (in other words, u is undefined). I also find it a somewhat inelegant approach. I'd either have the function take two Atom objects or two (X, Y, Z) position arrays (my favorite option since it parallels the behavior of distance_array).

Afterwards I'd seriously consider looking into using numpy for the fixing of PBC. It can make your code much faster, more compact, and would bring it in line with the rest of the MDAnalysis approach. I don't have time now for more details on how this should be done, but do give it a try.

[ayushsuhane]

Hi, I am not sure whether this is still a open issue or a closed one. I see one code in the issue discussion but didn's see any PR. I wish to take up this issue and resolve it. I would need some direction though. Can anyone help me in this regard?

[dpadula85]

Hi everyone,

is there a simple way to compute distances taking into account PBC in 1-2 dimensions only? As far as I understand, the current implementation only allows for either no PBC or PBC in all 3 dimensions.

Thanks,
D.

[IAlibay]

@dpadula85 may I ask you to raise this query in the mailing list please (rather than a closed issue)? Should you wish to raise this as a feature request, then a new issue can be raised a as feature request too (unless such a pre-existing feature request already exists).