[gmx-developers] umbrella pull code modified to apply a force of constant magnitude

[chris.neale at utoronto.ca]

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);
}