[gmx-developers] umbrella pull code modified to apply a force of constant magnitude
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Mon Jun 23 03:05:00 CEST 2008
Hello,
I am attempting to use the umbrella option of the pull code except
that I want the force to be of constant magnitude along the usual
vector, and not related to the current displacement along that vector.
I have made some modifications to src/mdlib/pull.c of version 3.3.1
(listed below with comments to the side of new lines /* Added for
constant force */). It compiles fine and appears to work fine. I was
just hoping to get some expert opinion as to whether I appear to have
done this correctly. Of course the final responsibility is mine.
Thanks,
Chris.
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 */
double normalizesum; /* Added for constant force */
/* 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
*/
/* Set desired position in x_ref[i] */
if(pull->bCyl) /* Dynamic Case */
dvec_add(pull->dyna[i].x_unc, pull->grp[i].UmbPos, pull->grp[i].x_ref);
else /* Normal Case */
dvec_add(pull->ref.x_unc, pull->grp[i].UmbPos, pull->grp[i].x_ref);
/* 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 */
normalizesum=0.0; /* Added
for constant force */
for(m=DIM-1;m>=0;m--) {
dr[m] *= pull->dims[m];
normalizesum+=dr[m]*dr[m]; /* Added
for constant force */
}
normalizesum=sqrt(normalizesum); /* Added
for constant force */
copy_dvec(dr,pull->grp[i].spring);
if(normalizesum>0.001){ /*
Added for constant force */
for(m=DIM-1;m>=0;m--) { /* Added
for constant force */
dr[m]/=normalizesum; /* Added
for constant force */
} /* Added
for constant force */
} /* Added
for constant force */
/* f = um_width*dx */
dsvmul(pull->grp[i].UmbCons,dr,pull->grp[i].f);
if (pull->bCyl) {
d_pbc_dx(box, pull->dyna[i].x_unc, pull->grp[i].x_unc, dr);
} else {
d_pbc_dx(box, pull->ref.x_unc, pull->grp[i].x_unc, dr);
}
normalizesum=0.0; /* Added
for constant force */
for(m=0; m<DIM; m++){
dr[m] *= pull->dims[m];
normalizesum+=dr[m]*dr[m]; /* Added
for constant force */
}
normalizesum=sqrt(normalizesum); /* Added
for constant force */
if(normalizesum>=0.001){ /*
Added for constant force */
for(m=DIM-1;m>=0;m--) { /* Added
for constant force */
dr[m]/=normalizesum; /* Added
for constant force */
} /* Added
for constant force */
} /* Added
for constant force */
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);
}
More information about the gromacs.org_gmx-developers
mailing list