[gmx-users] pull code yields non-zero separations at t=0

[chris.neale at utoronto.ca]

Using gromacs 3.3.1: I am currently trying to using the pull code to do
non-equilibrium simulations that will remove a ligand from a binding pocket in
a directionally unbiased manner. As far as I can tell, some type of harmonic
distance restraint is not possible using center  of mass that has a minimum
biasing energy at a given distance (other than 0.0). Therefore I am extracting
my ligand in a non-equilibrium manner using a negative force constant and the
umbrella option of the pull code.

Since the first motions away from the initial position are likely to determine
the direction in which an extraction will be attempted (although less so with
an inverted harmonic restraint like this) I am doing many repeated runs and am
relying on the fact that the runs all start *exactly* at the center of the
inverse restraint while using gen_vel=yes and unconstrained_start=yes in order
to get an ensemble of simulations that, together, attempt extraction in an
unbiased direction.

However, when running the pull code, my separations at t=0 are non-zero.
For me it is fine, since each run starts from a different conformation and
the rounding should be at random in a directional sense. I mention it first to
put it on the list and second in case those updating the pull code for
version 4 are interested in looking into this.

In my implementation, I script the determination of the exact offset of the
ligand from a convenient measure of the COM of the cognate receptor. In X/Y/Z:

Offset for run 0 is .015937 -.271066 .670338 based on
Receptor: 2.370000 2.510879 3.596662 and
Ligand:   2.385937 2.239813 4.267000

And then the output .pdo from the run indicates:

# UMBRELLA      3.0
# Component selection: 1 1 0
# nSkip 1
# Ref. Group  
'r_280-291_r_311-318_r_324-333_r_339-351_r_359-370_r_380-391_r_394-401_r_410-419_&_C_O_CA_N'
# Nr. of pull groups 1
# Group 1 'r_422'  Umb. Pos. 0.015937 -0.271066  Umb. Cons.  
-10000.000000 -10000.000000
#####
0.000000        0.002477        -0.000157
...

For which I had expected the output at t=0 to be 0,0.

I realize that the eventual crash with a negative force constant of this
magnitude is entirely unrelated to any other issue. I only used it in
a first run to be sure that a negative force constant would even work
(it does).

I also realize that the problem here could be with my determination of the COM
to start this run. I used my own g_com program to do that, but it is so simple
that I think it must be correct. The relevant lines of that program look like
this:

   /* This is the main loop over frames */
   numfr=0;
   tot_bend=0;


   resnr=top.atoms.atom[index[0]].resnr;
   warned=0;
   for(i=1; i<nidx; i++){
     if (top.atoms.atom[index[i]].resnr!=resnr && !warned){
       printf("Warning: all atoms not part of same residue\n");
       warned=1;
       //exit(1);
     }
   }
   fb=xvgropen(bendiness_file,(char *)NULL,(char *)NULL,(char *)NULL);

   fprintf(fb,"Frame\tXofCOM\tYofCOM\tZofCOM\n");
   do {
     /* coordinates are available in the vector fr.x
      * you can find this and all other structures in
      * the types directory under the gromacs include dir.
      * Note how flags determines wheter to read x/v/f!
      */
     copy_mat(box,box_pbc);
     if (bPBC) {
       rm_pbc(&top.idef,top.atoms.nr,box,fr.x,fr.x);
       set_pbc(&pbc,box_pbc);
     }

     clear_rvec(coms);
     coms_count=0;
     for(i=0; i<nidx; i++){
       rvec_inc(coms,fr.x[index[i]]);
       ++coms_count;
     }
     numfr++;
     if (coms_count==0){
       printf("Error: coms_count==0\n");
       fclose(fb);
       exit(1);
     }
     coms[ZZ]/=(double)coms_count;
     coms[YY]/=(double)coms_count;
     coms[XX]/=(double)coms_count;
     fprintf(fb,"%d\t%f\t%f\t%f\n",numfr,coms[XX],coms[YY],coms[ZZ]);
   } while(read_next_frame(status,&fr));

Thanks,
Chris.

As a PS note to developers: a distance-from-group pull-code option would be
awesome to have in gromacs 4 since this type of restraint flexibility
(and additional PBCs -- e.g. P2_1) are the only things that CHARMM still
has over gromacs.