[gmx-developers] adding external forces to all atoms

[Julius Su]

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:

http://lists.gromacs.org/pipermail/gmx-developers/2011-May/005193.html
http://lists.gromacs.org/pipermail/gmx-developers/2011-April/005141.html

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!

Sincerely,

Julius

-----
(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
             snew(f_global,state_global->natoms);

(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
         MPI_Barrier(cr->mpi_comm_mygroup);

(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 
functions:

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

----