Computing the Local Stress Tensor in Classical MD Simulations

Custom GROMACS version to compute the local stress tensor in 3D from a molecular dynamics simulation. This program computes the local stress based on the Hardy stress definition (see Ch.8 in Tadmor & Miller (2011) Modeling Materials: Continuum, Atomistic and Multiscale Techniques). The atomic virial stress (A. P. Thompson et al., J. Chem. Phys. 131, 154107 (2009)) definition was also recently added. Methodological details of the implementation can be found in: 
Importance of Force Decomposition for Local Stress Calculations in Biomembrane Molecular Simulations by J. M. Vanegas, A. Torres-Sanchez, and M. Arroyo. J. Chem. Theor. Comput.  10, 691 (2014). 
Examining the Mechanical Equilibrium of Microscopic Stresses in Molecular Simulations by  A. Torres-Sánchez, J. M. Vanegas, and M. Arroyo. Physical Review Letters, 114, 258102 (2015). 
We have patched the GROMACS source code (v4.5.5) at various locations in the routines that calculate particles forces and velocities. Our patches rely heavily on the code structure and functions implemented in a previous local pressure code (obtained from and This previous code implements the methodology outlined in 
3D Pressure Field in Lipid Membranes and Membrane-Protein Complexes by S. Ollila et al. Phys. Rev. Lett. 102, 078101 (2009).  
Highlights of GROMACS-LS:  
- Decomposition of multi-body potential forces using:
    1. Covariant central force decomposition (cCFD) (Torres-Sanchez A. Vanegas, J. M. and Arroyo, M. (Phys. Rev. Lett., 2015)).
    2. Non-covariant central force decomposition (nCFD) (N. C. Admal and E. B. Tadmor; J. Elast. 100, 63-143, 2010).
    3. Goetz-Lipowsky decomposition (GLD) (R. Goetz and R. J. Lipowsky; J. Chem. Phys. 108, 7397-7409, 1998.). 
    4. Method of planes (Heinz, A.  Paul W. and Binder K. Phys. Rev. E. 72 066704 (2005)).
- Ability to output the total or individual contributions to the local stress such as those from vdw, electrostatics, angles, and others.
- Finite discretization over a rectangular grid using tri-linear weight functions, which result in smoother stress fields and also make the discretization exact regardless of the grid size.
- Consistent treatment of forces arising from bond constraints (LINCS, SETTLE and SHAKE algorithms), which produces constant profiles for the normal component (perpendicular to the membrane plane) of the stress. This is necessary from an equilibrium stand-point and was an unresolved problem for stress profiles obtained from atomistic membrane simulations.
- Virial stress per atom as defined in A. P. Thompson et al., J. Chem. Phys. 131, 154107 (2009). 
If you publish results using GROMACS-LS, we kindly ask you to cite the papers by Vanegas et al., Torres-Sanchez et al., and Ollila et al., which we recommend you to read before continuing with the rest of the manual. In these papers we show the main differences between the different definitions of the stress. We show that the virial stress per atom and the stress from the method of planes do not satisfy balance of linear momentum, while the stress from cCFD, nCFD and GLD do. We also show that the GLD and the stress from the decomposition on geometric centers lead to non-symmetric stresses, which therefore do not satisfy balance of angular momentum. The two flavors of the CFD do satisfy both balance of linear and angular momentum by construction. However, for potentials beyond 4-body, such as CMAP, the nCFD leads to unphysical stresses. We therefore suggest that the cCFD definition is preferred in any case.
The types of interaction potentials that are included in stress calculations and have been tested are the following:
- Bonds
- Angles
- Dihedrals (proper, improper, and Ryckaert-Bellemans)
- Constraints (SETTLE and LINCS, SHAKE is included but has not been tested)
- Van der Waals
- Coulomb (plain cut-off and reaction-field) 
As such, if you use this code to analyze a system that may use other types of potentials, we cannot guarantee that it will be correct. Please thoroughly check the output of the program and compare to known quantities such as the total pressure. For such comparisons you should first obtain the total virial pressure in double precision. To do this, reanalyze the trajectory (i.e. mdrun_d -rerun) using a double precision version of GROMACS and recompute the pressure with g_energy_d. 
Contact information:
Juan M. Vanegas:
Alejandro Torres-Sánchez:
Binary Data gromacs-4.5.5-ls-5.0.tar.gz17.7 MB
PDF icon Local_stress.pdf3.86 MB