Importance of Force Decomposition for Local Stress Calculations in Biomembrane Molecular Simulations

Abstract

Local stress fields are routinely computed from molecular dynamics trajectories to understand the structure and mechanical properties of lipid bilayers. These calculations can be systematically understood with the Irving-Kirkwood-Noll theory. In identifying the stress tensor, a crucial step is the decomposition of the forces on the particles into pairwise contributions.However, such a decomposition is not unique in general, leading to an ambiguity in the definition of the stress tensor, particularly for multibody potentials. Furthermore, a theoretical treatment of constraints in local stress calculations has been lacking. Here, we present a new implementation of local stress calculations that systematically treats constraints and considers a privileged decomposition, the central force decomposition, that leads to a symmetric stress tensor by construction. We focus on biomembranes, although the methodology presented here is widely applicable. Our results show that some unphysical behavior obtained with previous implementations, e.g. non-constant normal stress profiles along an isotropic bilayer in equilibrium, is a consequence of an improper treatment of constraints. Furthermore, other valid force decompositions produce significantly different stress profiles, particularly in the presence of * To whom correspondence should be addressed † Universitat Politècnica de Catalunya-BarcelonaTech ‡ Contributed equally to this work 1 dihedral potentials. Our methodology reveals the striking effect of unsaturations on the bilayer mechanics, missed by previous stress calculation implementations.

Pressure Profile Calculation with Mesh Ewald Methods

“…7 in Ref. 20 ). Sonne and collaborators 11 showed that the pressure profile calculated using the Irving-Kirkwood approach and a cutoff radius of 2 nm is qualitatively similar to the one calculated using the Harasima path with long-range electrostatic contributions, although with an indetermination of roughly 100 to 200 bar, which makes it only a qualitative test.…”

A direct two-dimensional pressure formulation in molecular dynamics

“…In our chosen example (argon liquid film in vapor) and also for lipid bilayer [7], the inhomogeneity is in two dimensions (say, and axes) and there is no density variation along the third dimension ( axis). The computational domain is shown in figure. 2a.…”

“…However, current literature on local pressure estimation is based on 3D [3] or 1D [4][5][6] pressure estimation. The 2D pressure distribution is obtained by averaging over the 3D pressure data, and is extremely computationally expensive [7] as it involves a 3D convolution. A generalized method for 3D stress calculations which included temporal averaging weight functions was derived by Yang [8].…”

“…A generalized method for 3D stress calculations which included temporal averaging weight functions was derived by Yang [8]. Recently, Vanegas [7] and Sanchez et al [9] applied the modified Hardy versions of IK stress to lipid bilayers, coiled coil protein and graphene sheet to determine continuum level properties from atomistic simulations. Further, there exist a few IrvingKirkwood versions [10][11][12] of 1D pressure calculations for 1D inhomogeneous system.…”

Central-force decomposition of spline-based modified embedded atom method potential

“…The possibility of applying Hardy stress to three-body potentials has also been discussed [56]. Only recently applications of Hardy stress to systems described by four-body potentials have been reported [63,64]. To our best knowledge, there are no examples in the literature of the use of Hardy stress in systems described by many-body potentials more involved than the original EAM potential.…”