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

[Berk Hess]

Hi,

Gromacs 4 has this option, you can compare your results with the current 
CVS version.

Berk.

chris.neale at utoronto.ca wrote:
> 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);
> }
>
>
> _______________________________________________
> 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
>