[gmx-developers] Difference between the distance used for force and virial calculations in umbrella pull code
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Fri Aug 8 07:29:43 CEST 2008
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.
More information about the gromacs.org_gmx-developers
mailing list