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

[Berk Hess]

Hi,

Gromacs 3.3.3 and CVS do exactly the same thing.
In 3.3.3 things are bit more messy, since the variable x_unc is used,
which is somewhat confusing for umbrella pulling, but it is the same
variable that is used for constraint pulling.
In the CVS version I have separated the two functionalities better.
But in the code below you can see that dr (from pull_d_pbc_dx) is used 
for the virial.
This is the distance between the two groups.

Berk.

chris.neale at utoronto.ca wrote:
> 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.
>
>
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.
>
> This email was Anti Virus checked by Astaro Security Gateway. 
> http://www.astaro.com
>