[gmx-developers] Difference between the distance used for force and virial calculations in umbrella pull code
hessb at mpip-mainz.mpg.de
hessb at mpip-mainz.mpg.de
Tue Aug 12 08:47:33 CEST 2008
Hi,
Sorry, but I don't really understand your problem.
The formula (and the code) is very simple:
vir = 0.5 dr x f
If the force f is zero vir is zero.
The force ALWAYS works along vector dr connecting the two groups.
Berk.
> 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.
>
> _______________________________________________
> 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
>
>
More information about the gromacs.org_gmx-developers
mailing list