Programmer's Guide

main()

(mdrun.c) The main() function prepares a communication record cr, parses command line arguments, and opens new log files. Then it passes these information to mdrunner().

• cr = init_par() : initialize the communication record cr.
• parse_common_args(): parse program parameters.
• init_multisystem(cr) : modify the communication record cr for multiple-system simulations (e.g., replica exchange).
• prepare other parameters, and call mdrunner().

mdrunner()

(md.c)

The mdrunner() function first prepares an input record, a topology record, a force record, as well as runtime information, such as local topology and system configuration, then calls an integrator such as do_md() to perform simulation.

• init_parallel() or init_single(): MPI master loads the .tpr file, extracts the gmx_mtop_t and t_inputrec structures and scatters this data to all the clients
• integrator[ inputrec->eI ].func(): function pointer depending on inputrec->eI is called. (see const gmx_intp_t integrator[eiNR] in md.c).

This calls one of ({ {do_md}, {do_steep}, {do_cg}, {do_md}, {do_md}, {do_nm}, {do_lbfgs}, {do_tpi}, {do_tpi}, {do_md} }) depending on the type of integrator defined in the mdp file. From here on we will assume that do_md was selected!

do_md()

(md.c)

• dd_partition_system(): distribute atoms among clients (if MPI is enabled), creates an atom index
• atoms2md(): translates the atom type information from gmx_moltype_t into the runtime structures (t_mdatoms)
• update_mdatoms(): interpolates mass of perturbed atoms to the current lambda value, charges are done in PME code, currently I don't know where lj params are interpolated

while(step<nsteps)
{
 dd_partition_system(): redistribute atoms every neighbor search steps, calls atoms2md() with index
 do_force(): calculate forces (bonded and nonbonded) acting on the (home) atoms (no position update here)
 write_traj(): write current coordinates to trajectory file (gathers coordinates before hands if MPI-enabled)
 update(): tcoupl, constraints, second part of leap frog, update X coords
 step++;
}

do_force()

(sim_util.c)

• dd_move_x(): synchronize X coords with neighbor clients (if MPI-enabled)
• do_ns(): neighbor searching
• do_force_lowlevel(): calc bonded and non bonded forces acting on the (home) atoms
• dd_move_f(): synchronize forces with neighbor clients (if MPI-enabled)
• calc_viral(): calc viral tensor for home atoms based on local atom coordinate and force
• does periodic boundary conditions if domain decomposition is not enabled

do_force_lowlevel()

(force.c)

• do_walls(): calc forces and viral for walls
• do_nonbonded(): calc forces and viral from nonbonded interactions
• calc_bonds(): calc bonds
• gmx_pme_do(): particle mesh ewald (or similar function for different estatic type)

do_nonbonded()

(nonbonded.c)

• uses nb_kernel_list[nblists->nlist_sr[i].il_code] (function pointer to kernels, normally optimized for architecture) to calc forces
• nblists are constructed during neighbor search
• gromacs can be forced to use non assembly routines to do these calculations (easier to debug) by defining "NOASSEMBLYLOOPS=1" as environment variable (e.g. put "export NOASSEMBLYLOOPS=1" without the quotes in your /etc/profile)

atoms2md()

(mdatoms.c)

• copy the attributes of the (home) atoms from gmx_moltype_t to the runtime data structure t_mdatoms
• if an index is supplied (by user or domain decomposition) then the index contains the global atom index. So if an atom has offset i in the global topology (global index) then with index the atom is at the position where index[j]==i (see index, local index).
• here only the attributes in t_atom are copied, other interactions like lennard jones from gmx_ffparams_t are done in mk_nbfp

mk_nbfp()

(force.c)

• copys the lennard jones or the BHAM parameters from gmx_ffparams_t into a symmetric matrix of size atnr²
• interactions between type a and b can be found by matrix[b*atnr+a], code uses C6(),C12() and other macros
• these interactions are used at runtime in the [assembly] kernels
• the interactions are stored in the t_forcerec struct defined in mdrunner()

dd_partition_system()

(domdec.c)

• only called if MPI and domain decomposition are enabled
• does periodic boundary conditions
• divides the system into spatial cells and distributes the atoms among them
• creates the cell specific index and calls atoms2md(): on each client
• only neighbor cells communicate with each other (MPI)
• domain decomp data structures are in "cr->dd"
• to find an atom from its see here

write_traj()

(stat.c)

• first gathers all coordinate info from clients to master if MPI is enabled
• then dumps current coords into trajectory file

communication

Most of the communication going on requires all nodes which are part of a communication to perform the same (collective) communication operations in the same order. For example, all clients enter the write_traj() function after the same number of steps even though they don't actually write trajectories (only the master does). Within this function the clients prepare the data (atom coordinates) the master expects (what the master expects has been communicated beforehand) and call mpi_gather to send this data. The master enters write_traj() and sets up buffers to receive the different coordinate sets from the clients. Then it calls mpi_gather to receive the client data. So there is no dedicated communication protocol when sending or receiving data that tells the receiver what kind of data it is actually receiving. This information is implicit in that all nodes (master and clients) call the same collective MPI functions in the same order and frequency.

The client nodes (which own one of the domain decomposition cells per node) can communicate only with the master and neighboring nodes. So if data from a node is required that is more than one cell away then this fails!

The inter client communication takes place mostly using non blocking MPI_Irecv,MPI_Isend or MPI_Sendrecv calls. The synchronization of the master with the clients happens using gather and scatter.

Only the initial setup communication in mdrunner() uses blocking mpi calls.

local index

The local atom indices are constructed by the master in dd_partition_system (the atoms are distributed based on spatial subdivision of the simulation box among the available nodes). Only the master retains a global state which is updated from the local state on each node before writing the current atom positions to the trajectory file using gather. Each node calculates forces and positions for its home atoms and communicates these results with the neighboring cells (dd_move_x() and dd_move_f() calls in do_force). Each node constructs its local state from the index supplied by the master in atoms2md. In order to find the local coordinate of an atom given its global index (so its global offset into X) the atom must be looked up in the global to local atom lookup table:

int cell_id=cr->dd->ga2la[global_index].cell;
int local_index=cr->dd->ga2la[global_index].a;
if(0==cell_id)
{
 rvec* atom_position=&local_state->x[local_index];
}
else
{
 // atom is NOT on this node
 // if atom is on a neighbor node then its position can be acquired using special communication
}

Each node has access to a copy of the topology information, so that the global atom index can be found for each local atom.

home

The home atoms on a node are those atoms that the node is supposed to handle. On a single core run that will be all atoms. With MPI and domain decomposition each node only handles a few atoms (or only Ewald summation or similar). Which atoms are home atoms to a specific node is determined by the master in dd_partition_system() and the actual run time atom data is composed in the t_mdatoms structure in atoms2md(). So if you want to add another force term acting on the simulation atoms you must iterate over the local atoms like this (do_walls() in wall.c part of mdlib might be a good place to see how a simple additional force is handled):

for(i=md->start;i<md->start+md->homenr;i++) //iterate over home atoms
{
 rvec* atom_pos=&state->x[i];//get atom position
 float mass=md->massT[i];//get atom mass
 //must handle perturbation, how???
 /*
 float q=md->chargeA[i]
 int type=md->typeA[i];
 */
 f[i]+=F(atom_pos,q,mass,...);//calc force dependant on position and other stuff
 //f is the force array defined in mdrunner()
 //F() is a user defined function calculating the force
}

The viral contribution for forces in the force array are handled automatically in calc_viral().

If anyhow possible then don't modify atom positions/forces in places where the original code doesn't do it. Quite likely you will break some assumptions about charge group center of mass or domain decomposition which will cause weird crashes or worse simple numerical hick ups which are a nightmare to find.

gather

All clients communicate their current local state to the master using MPI_Gather.

scatter

The master distributes data, like the current domain decomposition layout, to all clients using MPI_Scatter. The amount of data transmitted for each individual client has been communicated before hands (setup_dd_communication() in domdec.c)