i don't know if this is faster but i store trajectory coordinates in a nested Python list and then i go through the list with a Python script, item by item (frame by frame), and i do my measurements, such as distances between points, and i store my results into Python dictionary.

I use PLUMED DRIVER to do this kind of analysis. Check out the manual page of Coordination Number Collective variable (CV) here: (https://plumed.github.io/doc-v2.4/user-doc/html/_c_o_o_r_d_i_n_a_t_i_o_n.html). If the distance between two atoms or a pair of atoms is within the specified distance (e.g. 4.5 A), this function outputs 1 and if the distance is >4.5 A, it returns 0. Ultimately, you will have a file with 0 or 1 value for each frame. You may then sum up the entire coordination number column to find out the number of times the contact was established in your trajectory.

MDAnalysis (https://www.mdanalysis.org/) has a number of functions to do this, the lib.distances one for example (https://www.mdanalysis.org/docs/documentation_pages/analysis/distances.html). As Martin Klvana said above, I'd then store residues below the cut off as a Python dictionary.

The easiest way to do so, is using gmx select in Gromacs. The tricky part is when you specify the option -select

Try this one:

1-make an index with the groups you want to compare (if you need only the protein create only one group)

2-use this line in gromacs

gmx select -f -s -n -of -ofpdb -selrpos atom -seltype res_com -select 'group 1 and within 0.4 of group 20'

-ofpdb is for generate one pdb with the occupancy factor (0 to 1). You can visualize it in VMD coloring by occupancy

-of is for generate occupancy.xvg when you can see the occupancy for each residue

-the select part includes a line in VMD syntaxis , 0.4 for cutoff, you can change it to 0.45 or any number

If you have any issue with this write me to [email protected]

Claudia Mejías Márquez

Thank you for your detailed explanation.

I would like to use gmx select for the same purpose but I am having problems with the selection part.

Let's say I have 2 groups: one for protein residues and the second for CHOL molecules.

I would like to see how many times a CHOL molecule is within a distance of 0.55 and 1.4nm of any residue of the protein, using their COM.

I tried something like:

gmx select -f traj.xtc -s md.tpr -n index.ndx -of occupancy.xvg -ofpdb occupancy.pdb -selrpos whole_res_com -seltype whole_res_com -select 'group 1 and within 0.55 and 1.4 of group 15' but it's not working. :/

Any ideas how to select a range of distance? :)

If you want to do this in gmx, look up mdmat

If this is a process you will need to repeat though, MDAnalysis coupled with matplotlib for imaging is going to be more efficient. This is one of the more trivial analysis codes to write. Here is an example code I've written, but if you are going to be doing analysis in the future it would help to write your own and use this as a reference (attached c_a.py)

^ you also need to import from scipy

from scipy.spatial import distance

Hi Sichaib,

You can use gmx distance that allow you calculate center of mass between two molecules.