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

[chris.neale at utoronto.ca]

My appologies, I wasn't doing mass weighting to determine the initial offset.
It works perfectly with a loop like this in g_com.c:

     clear_rvec(coms);
     coms_count=0.0;
     for(i=0; i<nidx; i++){
       coms[XX]+=fr.x[index[i]][XX]*top.atoms.atom[index[i]].m;
       coms[YY]+=fr.x[index[i]][YY]*top.atoms.atom[index[i]].m;
       coms[ZZ]+=fr.x[index[i]][ZZ]*top.atoms.atom[index[i]].m;
       coms_count+=top.atoms.atom[index[i]].m;
     }
     numfr++;

Chris.


^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.
^