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

[Chris Neale]

Dear Berk,

thanks for the explanation. Can you please confirm that the virial 
should be based on the distance between the two groups when the force is 
not actually applied along the vector connecting these two groups? For 
instance, while using epullgDIST, when the pulled group is exactly at 
the desired distance from the reference group the biasing force will be 
zero, and I would have thought that the contribution to the virial 
should also be zero, whereas I am not convinced that it is in the cvs.

Thanks,
Chris.

Berk Hess wrote:
> 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
>>
>
> _______________________________________________
> 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.