[gmx-developers] Difference between the distance used for force and virial calculations in umbrella pull code

[chris.neale at utoronto.ca]

Hello,

Regarding the umbrella pull code, I would appreciate if somebody could  
explain to me the difference between the distance used to calculate  
the force and the distance used to calculate the virial.

Specifically, I am currently attempting to modify gromacs 3.3.1 to  
allow pulling to a specified distance from the reference group. I  
realize that this is already in the cvs, but I'd rather avoid  
production runs with the cvs until it becomes gromacs 4.

I am stuck on the virial for the pull code. I have attempted to follow  
the logic in the cvs pull.c, but I am wondering if there has been some  
change in the way the virial is calculated for the pull code between  
gmx 3.3.3 and the current cvs.

 From what I can understand,

In 3.3.3, pull.c has do_umbrella() that determines a variable called  
dr in one way to get the force:
     d_pbc_dx(box, pull->grp[i].x_ref, pull->grp[i].x_unc, dr);
as a vector between the current and desired positions for the pulled group.

And then recalculates dr for use in the virial:
       d_pbc_dx(box, pull->ref.x_unc, pull->grp[i].x_unc, dr);
as a vector between the current position of the reference group and  
the current position of the pulled group.

In the cvs, however, pull.c has do_pull_pot() that calls  
get_pullgrp_distance(), which determines dev (a more advanced version  
of dr, I think) as:
   pull_d_pbc_dx(pull->npbcdim, box, pgrp->x, pref->x, dr);

I am unsure what is used there for the virial. Is it the variable V  
that is returned by the function? or is it still via dr by way of the  
bVir boolean?:

     if (bVir) {
       /* Add the pull contribution to the virial */
       for(j=0; j<DIM; j++)
         for(m=0;m<DIM;m++)
           vir[j][m] += 0.5*pgrp->f[j]*dr[m];
     }

#########

Here is what I have so far for a special version that only does  
umbrella pull code by a specified distance, in case it is useful:

static void do_umbrella(t_pull *pull, rvec *x,rvec *f, tensor vir, matrix box,
                         t_mdatoms *md, int start, int homenr)
{
   int i,ii=0,j,m,g;
   dvec cm;    /* center of mass displacement of reference */
   dvec dr;    /* distance from com to potential center */
   dvec drp;
   real desiredDist,actualDist;

   /* loop over the groups that are being umbrella sampled */
   for(i=0;i<pull->ngrp;i++) {

     /* Fix distance stuff
        pull->ref.x_unc[0] has position of reference group
        pull->x_ref is the coordinate that we would like to be at
        pull->spring has the corrected position of the pulled group
     */

     if(pull->bCyl){
       fprintf(stderr,"error: Dynamic Case not allowed as it is not  
coded for this source code modification\n");
       exit(1);
     }

     /* Find vector from reference group to pulled group */
     d_pbc_dx(box, pull->grp[i].x_unc, pull->ref.x_unc, drp);
     /* Normalize this vector such that it has the same distance as
        pull->grp[i].UmbPos[XX] (hack). While using this hack, one must set
        UmbPos[xx] to the desired distance and yy,zz to zero
        in the input .ppa file.
     */
     actualDist=sqrt(drp[XX]*drp[XX]+drp[YY]*drp[YY]+drp[ZZ]*drp[ZZ]);
     desiredDist=pull->grp[i].UmbPos[XX];
     drp[XX]*=desiredDist/actualDist;
     drp[YY]*=desiredDist/actualDist;
     drp[ZZ]*=desiredDist/actualDist;

     /* Set desired position in x_ref[i] */
     pull->grp[i].x_ref[XX]=pull->ref.x_unc[XX]+drp[XX];
     pull->grp[i].x_ref[YY]=pull->ref.x_unc[YY]+drp[YY];
     pull->grp[i].x_ref[ZZ]=pull->ref.x_unc[ZZ]+drp[ZZ];

     /* Find vector between current and desired positions */
     d_pbc_dx(box, pull->grp[i].x_ref, pull->grp[i].x_unc, dr);

     /* Select the components we want */
     for(m=DIM-1;m>=0;m--) {
       dr[m] *= pull->dims[m];
     }
     copy_dvec(dr,pull->grp[i].spring);

     /* f = um_width*dx */
     dsvmul(pull->grp[i].UmbCons,dr,pull->grp[i].f);

     /* VIRIAL NOT CHANGED YET */
     d_pbc_dx(box, pull->ref.x_unc, pull->grp[i].x_unc, dr);

     for(m=0; m<DIM; m++)
       dr[m] *= pull->dims[m];
     for(j=0; j<DIM; j++)
       for(m=0;m<DIM;m++)
         vir[j][m] += 0.5*pull->grp[i].f[j]*dr[m];
   }
   apply_forces(pull, md, f, start, homenr);
}

#######

where I am sure that the virial is not calculated correctly after I  
have made my changes.

Thanks,
Chris.