[gmx-developers] adding external forces to all atoms

Julius Su jsu at caltech.edu
Fri Oct 14 23:50:45 CEST 2011

Hi everyone,

I would like to present a solution I've arrived at to implement a 
feature -- the ability to add a set of external forces F(x1 ... xN) to 
all the atoms of a parallel GROMACS simulation -- which has been 
previously discussed:


This could useful for applying biasing potentials that span 
decomposition domains, such as ones acting on collective coordinates, or 
for applying novel experimental restraints which include large groups of 
atoms in the system. There is existing functionality in GROMACS to do 
this sort of thing, but the code below makes it straightforward to try 
out different potentials, without needing to perform a more extensive 
integration with the existing code base. So it may be useful for 
prototyping purposes ..

Note that in the approach below, all the work for evaluating F(x1 ... 
xN) is localized to one processor (the Master), so it should not a time 
consuming potential to evaluate -- else the parallel benefits of GROMACS 
are negated. Also the synchronization required to evaluate and 
distribute the forces may cause the parallel performance to suffer slightly.

In any case, since there appears to be some (scattered) interest, I'll 
present my code modifications. They may be of use to others in the 
community. Also if I've done anything egregiously (or even slightly) 
wrong, please let me know!



(1) In kernel/md.c, modify such that f_global is always initialized:

         //if (DDMASTER(cr->dd) && ir->nstfout) {
         if (1) {                                 // change so that 
f_global is always initialized

(2) In kernel/md.c, add lines in

        /*  ################## END TRAJECTORY OUTPUT ################ */

        // Collect all
        if (DOMAINDECOMP(cr))
           dd_collect_vec(cr->dd, state, f, f_global);

         if (MASTER(cr))
            // Fetch positions of all atoms
           for (i = 0; i < state_global->natoms; i++)
             current_x[i] = state_global->x[i][0];
             current_y[i] = state_global->x[i][1];
             current_z[i] = state_global->x[i][2];

           // Compute biasing potential (ADD NEW POTENTIAL HERE)
           ... Given current_x, current_y, current_z --> returns 
delta_fx, delta_fy, delta_fz.

           // Apply biasing potential to all atoms
           for (i = 0; i < state_global->natoms; i++)
              f_global[i][0] += delta_fx[i];
              f_global[i][1] += delta_fy[i];
              f_global[i][2] += delta_fz[i];

         // Distribute modified f_global back to individual processors
         if (DOMAINDECOMP(cr))
           dd_distribute_vec(cr->dd, dd_charge_groups_global(cr->dd), 
f_global, f);

         // Make sure all processes are synchronized before continuing

(3) In mdlib/domdec.c, modify dd_distribute_vec to have external scope:

//static void dd_distribute_vec(gmx_domdec_t *dd,t_block *cgs,rvec 
*v,rvec *lv)
void dd_distribute_vec(gmx_domdec_t *dd,t_block *cgs,rvec *v,rvec *lv) 
// allow external functions access

(4) In include/domdec.h, expose dd_distribute_vec for use by other 

void dd_distribute_vec(gmx_domdec_t *dd,t_block *cgs,rvec *v,rvec *lv);


More information about the gromacs.org_gmx-developers mailing list