[gmx-users] Testing nonbonded list update time

[chris.neale at utoronto.ca]

To estimate the largest atomic displacement between updates of the nonbonded
list (e.g. in order to verify that the rvdw/rlist options are appropriate) :

(1) do a short run (~100ps) and save the trajectory every 10steps or however
often the nb list is updated.

(2) trjconv -pbc nojump

(3) Use the program below to calculate the largest single atomic displacement
that occurs between frame i and frame i+1 for all i.

I have included here only the modifications that should be made to
gromacs-3.3/exec/share/gromacs/template/template.c (rename it before making changes)

Thanks to Erik who directed me to the template.c program in share/template.

To compile:
(a) copy the Makefile... that exists there to a new name
(b) change all occurrences of "template" to the name of your modified executable.
(c) You may want to add ../../../bin/ before the name of the executable to put
it in the bin directory which is probably already sourced on your system.


Modifications to template.c
Add ~line 39
#include <math.h>

Add ~line 87
int first,k;
real dx,dy,dz,dss,dss_max,ds_max,time_max;
real **last;
int atom_max=0;

Add ~line 105
last=(real **)malloc((top.atoms.nr+1)*sizeof(real *));
if(last==NULL){
  printf("Unable to allcoate memory for last\n");
  exit(1);
}
for (k=0;k<=top.atoms.nr;k++){
  last[k]=(real *)malloc(3*sizeof(real));
  if(last[k]==NULL){
    printf("Unable to allcoate memory for last[]\n");
    exit(1);
  }
  last[k][0]=last[k][1]=last[k][2]=0.0;
}

Add ~line 120 (before the loop)
first=1;
dss_max=0.0;
time_max=0.0;

Comment out line 127
//printf("Coordinates at t=%8.3f : %8.5f %8.5f
%8.5f\n",fr.time,fr.x[n][XX],fr.x[n][YY],fr.x[n][ZZ]);

Replace ~line 127 with this block (inside the loop)
if(first){
  first=0;
}else{
  for(k=1;k<=top.atoms.nr;k++){
    dx=fr.x[k][XX]-last[k][0];
    dy=fr.x[k][YY]-last[k][1];
    dz=fr.x[k][ZZ]-last[k][2];
    dss=(dx*dx)+(dy*dy)+(dz*dz);
    if(dss>dss_max){
      dss_max=dss;
      time_max=fr.time;
      atom_max=k;
    }
  }
}
for(k=1;k<=top.atoms.nr;k++){
  last[k][0]=fr.x[k][XX];
  last[k][1]=fr.x[k][YY];
  last[k][2]=fr.x[k][ZZ];
}

Add ~line 129 (just after the loop)
ds_max=sqrt(dss_max);
printf("Atom=%d_%s MaxD=%8.3f
time=%8.3f\n",atom_max,*(top.atoms.atomname[atom_max]),dss_max,time_max);