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 :math:`\Xi` can be written in full tensor form as: .. math:: \Xi~=~-\frac{1}{2}~\sum_{i < j}^N~\mathbf{r}_{ij}\otimes\mathbf{F}_{ij} :label: eqnvirialfulltensorform where :math:`\otimes` denotes the *direct product* of two vectors. [1]_ When this is computed in the inner loop of an MD program 9 multiplications and 9 additions are needed. [2]_ Here it is shown how it is possible to extract the virial calculation from the inner loop \ :ref:`177 `. Virial ~~~~~~ In a system with periodic boundary conditions, the periodicity must be taken into account for the virial: .. math:: \Xi~=~-\frac{1}{2}~\sum_{i < j}^{N}~\mathbf{r}_{ij}^n\otimes\mathbf{F}_{ij} :label: eqnvirialperiodicity where :math:`\mathbf{r}_{ij}^n` denotes the distance vector of the *nearest image* of atom :math:`i` from atom :math:`j`. In this definition we add a *shift vector* :math:`\delta_i` to the position vector :math:`\mathbf{r}_i` of atom :math:`i`. The difference vector :math:`\mathbf{r}_{ij}^n` is thus equal to: .. math:: \mathbf{r}_{ij}^n~=~\mathbf{r}_i+\delta_i-\mathbf{r}_j :label: eqnvirialdiffvector or in shorthand: .. math:: \mathbf{r}_{ij}^n~=~\mathbf{r}_i^n-\mathbf{r}_j :label: eqnvirialdiffvecshort In a triclinic system, there are 27 possible images of :math:`i`; when a truncated octahedron is used, there are 15 possible images. 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 :math:`i \neq j`. - In general, either the :math:`i` or :math:`j` particle might be shifted to a neighbor cell to get the closest interaction (shift :math:`\delta_{ij}`). 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 :math:`i` particle into multiple separate lists, where each list has a constant shift :math:`\delta_i` for the :math:`i` partlcle. We can represent this as a sum over shifts (for which we use index :math:`s`), 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 :math:`j` particles. For any sum that does not contain complex dependence on :math:`s`, this means the sum trivially reduces to just the sum over :math:`i` and/or :math:`j`. - To simplify some of the sums, we replace sums over :math:`j`, *i.e.*: #. There are three atoms in the molecule. #. The whole molecule is a single charge group. #. The first atom has Lennard-Jones (sec. :ref:`lj`) and Coulomb (sec. :ref:`coul`) interactions. #. Atoms two and three have only Coulomb interactions, and equal charges. These loops also works for the SPC/E \ :ref:`178 ` and TIP3P \ :ref:`128 ` water models. And for four site water models similar to TIP4P \ :ref:`128 `: #. There are four atoms in the molecule. #. The whole molecule is a single charge group. #. The first atom has only Lennard-Jones (sec. :ref:`lj`) interactions. #. Atoms two and three have only Coulomb (sec. :ref:`coul`) 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. .. raw:: latex \clearpage .. [1] Note that some derivations, an alternative notation :math:`\xi_{\mathrm{alt}} = v_{\xi} = p_{\xi}/Q` is used. .. [2] The calculation of Lennard-Jones and Coulomb forces is about 50 floating point operations.