Some implementation details#
In this chapter we will present some implementation details. This is far from complete, but we deemed it necessary to clarify some things that would otherwise be hard to understand.
Single Sum Virial in GROMACS#
The virial
where
Here it is shown how it is possible to extract the virial calculation from the inner loop 177.
Virial#
In a system with periodic boundary conditions, the periodicity must be taken into account for the virial:
where
or in shorthand:
In a triclinic system, there are 27 possible images of
Virial from non-bonded forces#
Here the derivation for the single sum virial in the non-bonded force routine is given. There are a couple of considerations that are special to GROMACS that we take into account:
When calculating short-range interactions, we apply the minimum image convention and only consider the closest image of each neighbor - and in particular we never allow interactions between a particle and any of its periodic images. For all the equations below, this means
.In general, either the
or particle might be shifted to a neighbor cell to get the closest interaction (shift ). However, with minimum image convention there can be at most 27 different shifts for particles in the central cell, and for typical (very short-ranged) biomolecular interactions there are typically only a few different shifts involved for each particle, not to mention that each interaction can only be present for one shift.For the GROMACS nonbonded interactions we use this to split the neighborlist of each
particle into multiple separate lists, where each list has a constant shift for the partlcle. We can represent this as a sum over shifts (for which we use index ), with the constraint that each particle interaction can only contribute to one of the terms in this sum, and the shift is no longer dependent on the particles. For any sum that does not contain complex dependence on , this means the sum trivially reduces to just the sum over and/or .To simplify some of the sums, we replace sums over
with double sums over all particles (remember, ) and divide by 2.
Starting from the above definition of the virial, we then get
In the second-last stage, we have used the property that each shift
vector itself does not depend on the coordinates of particle
which is the total force on
we must, in the implementation, double the term containing the shift
This separates the total virial
The intra-molecular shift (mol-shift)#
For the bonded forces and SHAKE it is possible to make a mol-shift list, in which the periodicity is stored. We simple have an array mshift in which for each atom an index in the shiftvec array is stored.
The algorithm to generate such a list can be derived from graph theory, considering each particle in a molecule as a bead in a graph, the bonds as edges.
Represent the bonds and atoms as bidirectional graph
Make all atoms white
Make one of the white atoms black (atom
) and put it in the central boxMake all of the neighbors of
that are currently white, grayPick one of the gray atoms (atom
), give it the correct periodicity with respect to any of its black neighbors and make it blackMake all of the neighbors of
that are currently white, grayIf any gray atom remains, go to [5]
If any white atom remains, go to [3]
Using this algorithm we can
optimize the bonded force calculation as well as SHAKE
calculate the virial from the bonded forces in the single sum method again
Find a representation of the bonds as a bidirectional graph.
Virial from Covalent Bonds#
Since the covalent bond force gives a contribution to the virial, we have:
The virial contribution from the bonds then is:
Virial from SHAKE#
An important contribution to the virial comes from shake. Satisfying the constraints a force G that is exerted on the particles “shaken.” If this force does not come out of the algorithm (as in standard SHAKE) it can be calculated afterward (when using leap-frog) by:
This does not help us in the general case. Only when no periodicity is needed (like in rigid water) this can be used, otherwise we must add the virial calculation in the inner loop of SHAKE.
When it is applicable the virial can be calculated in the single sum way:
where
Optimizations#
Here we describe some of the algorithmic optimizations used in GROMACS, apart from parallelism.
Inner Loops for Water#
GROMACS uses special inner loops to calculate non-bonded interactions for water molecules with other atoms, and yet another set of loops for interactions between pairs of water molecules. There highly optimized loops for two types of water models. For three site models similar to SPC 80, i.e.:
There are three atoms in the molecule.
The whole molecule is a single charge group.
The first atom has Lennard-Jones (sec. The Lennard-Jones interaction) and Coulomb (sec. Coulomb interaction) interactions.
Atoms two and three have only Coulomb interactions, and equal charges.
These loops also works for the SPC/E 178 and TIP3P 128 water models. And for four site water models similar to TIP4P 128:
There are four atoms in the molecule.
The whole molecule is a single charge group.
The first atom has only Lennard-Jones (sec. The Lennard-Jones interaction) interactions.
Atoms two and three have only Coulomb (sec. Coulomb interaction) interactions, and equal charges.
Atom four has only Coulomb interactions.
The benefit of these implementations is that there are more floating-point operations in a single loop, which implies that some compilers can schedule the code better. However, it turns out that even some of the most advanced compilers have problems with scheduling, implying that manual tweaking is necessary to get optimum performance. This may include common-sub-expression elimination, or moving code around.