Precise calculation of the local pressure tensor in Cartesian and spherical coordinates in LAMMPS

doi:10.1016/j.cpc.2014.11.017

Computer Physics Communications

Volume 190, May 2015, Pages 120-128

https://doi.org/10.1016/j.cpc.2014.11.017 Get rights and content

Abstract

An accurate and efficient algorithm for calculating the 3D pressure field has been developed and implemented in the open-source molecular dynamics package, LAMMPS. Additionally, an algorithm to compute the pressure profile along the radial direction in spherical coordinates has also been implemented. The latter is particularly useful for systems showing a spherical symmetry such as micelles and vesicles. These methods yield precise pressure fields based on the Irving–Kirkwood contour integration and are particularly useful for biomolecular force fields. The present methods are applied to several systems including a buckled membrane and a vesicle.

Program summary

Program title: Compute stress spatial

Catalogue identifier: AEVH_v1_0

Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEVH_v1_0.html

Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland

Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html

No. of lines in distributed program, including test data, etc.: 7799

No. of bytes in distributed program, including test data, etc.: 149739

Distribution format: tar.gz

Programming language: C++.

Computer: All.

Operating system: All.

Supplementary material: The input and expected output for the examples given in the manuscript can be downloaded here.

Classification: 7.7.

Nature of problem:

A precise calculation of the pressure (stress) field requires the implementation of the contour integration [2] for each pair in the cluster potentials.

Solution method:

We have implemented the method for calculating the local-volume average of the pressure field expressed by a contour integration with a delta function for the given molecular configuration obtained by molecular dynamics simulation [3].

Restrictions:

Because the definition of pressure field is based on the contour integration, the long-range interactions represented as Ewald summation are not included. In addition, the potential represented as a combination of pair and angle potential, such as Tersoff and Stillinger-Weber potential, and the interaction as a geometric constraint (shake, rigid body, etc.) are not available for current version.

Running time:

The examples provided take between 3 and 6 min each to run.

References:

[1] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1–19.

[2] P. Schofield and J.R. Henderson, Proc. R. Soc. Lond. A. 379 (1982) 231.

[3] T. Nakamura, W. Shinoda, and T. Ikeshoji, J. Chem. Phys. 135 (2011) 094106.

Introduction

Using molecular simulations, mechanical properties of materials and molecular assemblies have been widely investigated within a framework of elastic theory. In particular, by analyzing the local pressure (stress) distribution, several mechanical properties of heterogeneous, nano-structured systems can be well characterized. For example, local elastic compliances are represented by local fluctuation of the pressure tensor [1], an interfacial tension at the phase boundary is represented by the moment of pressure along the surface normal [2], a stress intensity factor of a crack is represented by the J-integral [3], a saddle-splay modulus and a spontaneous curvature for the fluid membrane are represented as a moment of pressure along the surface normal [4], and a lipid protein interaction is also characterized by a pressure profile [5]. Therefore the calculation of the local pressure distribution from molecular dynamic (MD) simulations is a highly promising technique for various systems in material science.

In MD simulations, the pressure tensor distribution is defined through the line integral between two points as a consequence of two- or $m (> 2)$ -body interactions between atoms [6]. Therefore, an algorithm to calculate the pressure tensor distribution is more complicated than that for calculating the mass-density or energy-density which can be obtained simply by plotting a histogram over the system.

As a simple approximation, one might define the pressure tensor distribution as a histogram by assigning a virial to each atom with a weight of $1 / m$ for the $m$ -cluster potential. This approach is referred to as “IK1” method [7], [8]. Even though the “IK1” method is intuitively reasonable and is simple to implement, it unfortunately yields an artificial atomic-scale wavy profile around the phase boundary [7], [8]. Therefore, in order to obtain a reasonable pressure tensor distribution at atomic-scale resolution, an accurate calculation method based on contour integration is required.

Algorithms to calculate accurate local pressures have been proposed for several different coordinate systems. For a planar surface, the pressure tensor profile along the surface normal is typically calculated [9], [10] based on the Irving–Kirkwood contour, though a similar algorithm has been implemented in the NAMD packages [11], based on the Harashima contour [12]. The choice of the integration contour does affect the calculated pressure profile [13]. Even though the Harashima contour is often used for convenience, it is known to be problematic at least for the case of spherical coordinate [14]. Recently, for an analysis of curved surfaces, an algorithm to calculate the pressure tensor profile along the radial direction of the spherical coordinate was developed based on the Irvin–Kirkwood contour [15]. Moreover, a method to calculate the pressure tensor over a 3-dimensional grid with small mesh size in Cartesian coordinate has also been developed [16].

In this work, we have implemented methods for calculating the local pressure tensor under three geometrical conditions into LAMMPS [17]. Our implementation is based on the algorithm proposed in Ref. [15] for a radial pressure profile in spherical coordinates and the algorithm proposed in Ref. [10] for the pressure profile along a certain direction in Cartesian coordinates. In addition, we propose a new algorithm to calculate the 3-dimensional local pressure distribution in Cartesian coordinates. This new algorithm, in principle, yields the local pressure at higher precision than the previous algorithm [16]. Our implementation allows us to calculate the pressure tensor distribution either in the course of a MD simulation or as post-processing using a MD trajectory.

This paper is organized as follows: In Section 2, we give a definition of the local pressure tensor $P_{α β} (x)$ from the equations of motion of $N_{f}$ particles. We calculate the local pressure tensor averaged over a sliced or grid microvolume, ${[P]}_{α β} (V)$ (Section 2.1), of which size determines the resolution of the calculated pressure profile. We also explain how to calculate and average the local pressure tensor for each geometry of the microvolume element, $V$ . Here we consider three different geometries: orthogonal volume elements divided into a 3-dimensional grid pattern in Section 2.2, planar bins sliced along a one-dimensional axis in Section 2.3, and spherical bins sliced along the radial axis from a chosen center in Section 2.4. In Section 3, it will be shown how we have implemented the algorithms for practical use into the LAMMPS MD code. After we demonstrate several examples of the actual computations of the pressure profile in Section 4, we give a summary of this paper in Section 6.

Section snippets

Microscopic definition of the local pressure tensor and its average over a local (micro)volume

Consider a $N_{f}$ -particles Hamilton system. The position and momentum of the $i$ th particle are represented by three dimensional vectors, $r_{i}$ and $p_{i} (1 \leq i \leq N_{f})$ , respectively. The equations of motion for the $i$ th particle are written as ${\dot{p}}_{i} = - \nabla_{i} Φ ({r_{i}}_{i = 1}^{N_{f}})$ and ${\dot{r}}_{i} = p_{i} / m_{i}$ . Here, $m_{i}$ is the mass of the $i$ th particle, $Φ ({r_{i}}_{i = 1}^{N_{f}})$ is the total potential energy, and $\nabla_{i}$ is a gradient operator with respect to the $i$ th particle position $r_{i}$ . When no one-body potential is included in $Φ$ , the local pressure tensor $P$ at

Installation instructions

To compute pressure profile by LAMMPS, we have prepared a patch file on the website [19] for the stable version “lammps-1Feb14”. To install the patch, download the patch file “Add_Stress_Spatial.patch” from [19], then copy it to the directory where the src directory exists. In the directory, usually the “lammps-1Feb14” directory, type a command of the kind:

Next, move to the src directory and activate the files that were modified as follows:

Then, upon building LAMMPS in the usual way, the

Applications

We demonstrate here several examples of the actual calculation of the pressure tensor distribution in different geometries.

Benchmark

As written in Introduction, the present method yields the precise 3D pressure tensor distribution, since it is free from the approximation for the weight calculation of the pressure tensor to each bin as used in the previous method by Ollila et al. [16]. In this section, we demonstrate that the present method is beneficial even for the computational efficiency if we require a pressure estimation within a small error e.g., less than 10 atm. We used the buckled membrane system (as used in

Summary

We have developed a program code for calculating the pressure field over a simulated molecular system, which is incorporated into the open source molecular dynamics package, LAMMPS. The implemented method yields the pressure (stress) profile along an axis in Cartesian coordinates, the pressure profile along the radial axis in spherical coordinate, and also the three-dimensional pressure field at an arbitrary chosen resolution. Using the program, in principle, heterogeneous structure of

Acknowledgments

We thank Dr. Axel Kohlmeyer for helpful discussion. This work is supported by JSPS KAKENHI Grant Number 23350014, the Next Generation Super Computing Project, TCCI/CMSI in the Strategic Programs for Innovative Research, MEXT, Japan. The simulations were conducted partly using computational resources of the HPCI systems provided by T2K and COMA at University of Tsukuba, TSUBAME2.0 and TSUBAME2.5 at Tokyo Institute of Technology, and Cray XE6 at Kyoto University through the HPCI System Research

References (31)

S. Plimpton
J. Comput. Phys.
(1995)
W. Shinoda et al.
Curr. Opin. Struct. Biol.
(2012)
J.F. Lutsuko
J. Appl. Phys.
(1988)
J.S. Rowlinson et al.
Molecular Theory of Capillarity
(1982)
H. Inoue et al.
Int. J. Fract.
(1994)
I. Szleifer et al.
J. Chem. Phys.
(1990)
S. Cantor
Biochemistry
(1997)
P. Schofield et al.
Proc. R. Soc. Lond. Ser. A
(1982)
B.D. Todd et al.
Phys. Rev. E
(1995)
T. Ikeshoji et al.
Mol. Simul.
(2003)

R. Goetz et al.

J. Chem. Phys.

(1998)

E. Lindahl et al.

J. Chem. Phys.

(2000)

J.C. Phillips et al.

J. Comput. Chem.

(2005)

A. Harashima

Adv. Chem. Phys.

(1958)

W. Shinoda et al.

Soft Matter

(2011)

Cited by (29)

Interfacial behaviors of the H<inf>2</inf>O+CO<inf>2</inf>+CH<inf>4</inf>+C<inf>10</inf>H<inf>22</inf> system in three phase equilibrium: A combined molecular dynamics simulation and density gradient theory investigation
2023, Journal of Molecular Liquids
Citation Excerpt :
However, there is an ongoing debate about the correct form of pressure tensor profile to be used in molecular simulation given the nonunique definition of the local pressure tensor [62–66]. More importantly, it is challenging to calculate pressure tensor profile in the systems with long-range interaction [66–68]. Hence, we calculated the IFT using the pressure tensor of the entire simulation box of the 2-phase systems split from the 3-phase system based on the Kirkwood and Buff approach[60] where the pressure profile is not required.
The interfacial behaviors in H₂O+gas+oil 3-phase systems are critical for CO₂ near-miscible/immiscible flooding and sequestration processes but require further investigation. In this article, molecular dynamics (MD) simulation and density gradient theory (DGT) with PC-SAFT equation of state were simultaneously conducted to study the interfacial behaviors in the H₂O+CO₂+C₁₀H₂₂ 3-phase system and the effect of impurity gas CH₄ on the interfacial properties at different temperatures (323–423 K) and pressures (up to around 16 MPa). We found reasonable agreement between the estimations from MD and DGT. When the H₂O+CO₂+C₁₀H₂₂ 3-phase system is in contact with CH₄, the interfacial tensions (IFTs) of all three interfaces increase. The increment of IFT is pronounced at low temperatures and high pressures. Remarkably, CH₄ molecules accumulate in all three interfaces. However, the positive surface excesses of CH₄ are smaller than those of CO₂, which may explain the increment of IFT. Moreover, the spreading coefficients S in the H₂O+CO₂+C₁₀H₂₂ 3-phase system are negative indicating the existence of 3-phase contact. The behaviors of S in the H₂O+CO₂+CH₄+C₁₀H₂₂ 3-phase system are similar to those in the system without CH₄. Nevertheless, the effect of temperature on S is little due to the changes of IFT caused by CH₄. These insights could help to enhance the understanding of the effects of impurities on CO₂ enhanced oil recovery methods under geological conditions.
Study of interfacial properties of water + methane + oil three-phase systems by a simple molecular simulation protocol
2022, Journal of Molecular Liquids
Citation Excerpt :
Nevertheless, there is an ongoing debate about the correct form of pressure tensor profile to be used in molecular simulation given the non-unique definition of the local pressure tensor [31–35]. In addition, it is challenging to calculate pressure tensor profile in the systems with long-range interaction and many-body potential [35–37]. Hence, it is desirable to develop an approach for IFT calculation in the three-phase systems without calculating the pressure tensor profile.
A straightforward molecular simulation protocol has been proposed to study the interfacial properties of the three-phase fluid systems. With this protocol, the interfctcial tension (IFT) can be calculated based on the pressure tensor of the entire simulation box of the two-phase systems split from the three-phase system without calculating the debatable pressure tensor profile. Water + methane + oil (decane, hexadecane, and toluene) three-phase systems was studied at different temperatures (323–423 K) and pressures (up to around 20 MPa). Reasonable agreement was found among the results obtained from molecular simulation, density gradient theory with the cubic-plus-association equation of state, and available experimental data in the literature. The IFT of the aqueous phase + vapor phase in the three-phase systems is smaller than the IFT in water + methane two-phase systems. Importantly, the reduction of IFT of aqueous phase + vapor phase in the three-phase systems containing decane or hexadecane is moderate while that in the three-phase systems containing toluene is significant, which can be explained by the stronger enrichment of toluene in the interfacial region in contrast to that of decane or hexadecane. Meanwhile, methane accumulates in the interfacial region of the aqueous phase + decane/hexadecane-rich phase in the three-phase systems, which causes the reduction of IFT with pressure while the opposite pressure effect was reported in the water + decane/hexadecane two-phase systems. The interfacial properties of the oil-rich phase + vapor phase in the three-phase systems are hardly affected by water due to the small amount of dissolved water. Furthermore, the calculated spreading coefficient of different types of oil in contact with water + methane under three-phase conditions follows this order: toluene $>$ decane $>$ hexadecane.
Phase equilibria and stability boundaries in a two–component Lennard-Jones mixture
2022, Fluid Phase Equilibria
The mechanical and diffusion stability conditions of binary Lennard-Jones mixtures with respect to infinitesimal changes in state parameters are considered. In stable and metastable states, by molecular dynamics method, the pressure is determined as a function of temperature, density and composition of a Lennard-Jones mixture simulating the methane–nitrogen system. The phase equilibrium parameters are calculated. Based on the data obtained, a unified mixture state equation for liquid and vapor was created. The diffusion and mechanical spinodals of the mixture were determined.
Atomistic characterization of the dispersed liquid droplet in immiscible Al–Pb alloy
2021, Journal of Materials Research and Technology
Citation Excerpt :
In the configurational contribution term of the Eq. (4), α and β are two directions in spherical coordinates, r, φ, θ. fij is the pairwise force acting on particle i and j, l denotes a position parameter on the Irving–Kirkwood integration contour Cij from particle i to particle j. Following Nakamura et al. [30,31], the configurational pressure is calculated by assigning the fractional virial to each fine-scale spherical shells, according to how much of the Irving–Kirkwood integration contour Cij (defined as the straight line between particle i and j) penetrates the spherical shells, see in Fig. 2 for details. All the interfacial profiles are reported in Gibbs' surface thermodynamics framework relative to an imaginary surface separating the droplet and bulk phases.
This study develops an atomistic simulation-based methodology for studying the dispersed droplet phase in immiscible alloy. The methodology is applied to the immiscible Al–Pb alloy liquid droplet system to study the temperature and curvature effects on the droplet radial distributions of thermodynamic properties, the mutual miscibilities within the spherical droplets, and the pressure differences across the curved liquid–liquid interfaces. The pressure tensor components profiles along radial direction are computed by the Irving–Kirkwood method coupling the spherical coordination and providing the inner droplet pressure data with sufficient statistical precision to validate the generalized Laplace equation in a liquid immiscible alloy system (both Al-rich convex droplet interface and Pb-rich concave droplet interface). Noticeable changes in mutual miscibilities in the droplet phases are seen, which correspond to the shifts of the solvus lines toward lower solubilities as temperature increases or the droplet size decreases. Compared with spherical droplets under the same temperature, cylindrical droplets with the same radii own the same Laplace pressure dependence of the mutual miscibilities but a weaker capability to tune the nano-alloy droplet phase diagram. Current findings contradict the miscibility changes in the two-complexion equilibria within the monoatomic layer at the Al(111)–Pb(liquid) interface [Acta Mater 2018;143:329]. The combined simulation-characterization methodology proposed in this study applies to broader types of immiscible alloy. The calculated data could provide guides to design well-dispersed droplet phases in the immiscible alloys and their manufacturing processes, facilitate the development of the thermodynamic theory/modeling of the nano-sized alloy phase diagram.
Simulation study of effects of grain boundary and helium bubble on lattice thermal resistance of tungsten
2021, Fusion Engineering and Design
Tungsten (W) as plasma facing material (PFM) has been used in ITER divertor. Many experiments indicate that due to helium (He) ion irradiation the thermal conductivities (TCs) of W can be significantly reduced. In this work, the lattice thermal transport property is studied based on the non-equilibrium molecular dynamics (NEMD) method. Firstly, the Kapitza resistance of grain boundaries (GBs) in W is obtained. It can be confirmed that the Kapitza resistance is positively correlated with GB energies, which reflect the lattice distortion degree in a sense. The result indicates that the distortion energy may play a key role in the thermal properties. Further simulation is performed to study the thermal resistance of GBs with nano He bubble. When the number of He atoms in bubble is less than a value, the thermal resistances are stable and change little. It implies that the effect of GBs is more significant than He bubble to block phonons diffusion. With the number of He atoms increasing, the property of heat conduction deteriorates further because the coupling of He bubble and GBs change the lattice distortion degree.
Nanovoids in uniaxially elongated polymer network filled with polydisperse nanoparticles via coarse-grained molecular dynamics simulation and two-dimensional scattering patterns
2019, Polymer
Citation Excerpt :
It is expected that these sharp tips related to the increased normal stress in the elongation direction when there were crosslinks. Detailed investigations with precious measurement of local stress profiles [105] are in progress. To demonstrate the possible correspondence between experimental results and CGMD simulations, Fig. 8 compares a sliced snapshot from the CGMD simulation and an in situ TEM image recorded by Jinnai and co-workers [59].
In order to investigate the relationship between nanovoids and two-dimensional scattering patterns (2DSPs) in elongated polymer nanocomposites (PNCs) filled with spherical nanoparticles (NPs), coarse-grained molecular dynamics simulations were performed with polydisperse NPs for Poisson's ratios of ν = 0.46 < 0.5. The polydispersity of NP diameter impacts on creation and early growth of nanovoids. Study of the growth of nanovoids under the attractive NP–polymer interactions showed that large NPs have a role of starting points of the growth. The 2DSP I(q_x, q_yz) is also calculated as the circular average of scattering intensity I(q), which is obtained by three-dimensional Fourier transformation of the density field weighted by contrasts of the NPs and polymers. It was confirmed that polydispersity of NPs causes disappearance of a ring pattern that was seen for uniform NPs. For the case with attractive NP–polymer interactions, the magnitudes of the central peaks of the 2DSPs was confirmed to become larger for larger elongation ratios. On the other hand, for the case with repulsive NP–polymer interactions, the central peak of the 2DSPs is much weaker than that for the case with attractive NP-polymer interaction, because of the occurrence of interfacial fractures and fine dispersal of nanovoids at the interface between NPs and polymers. In addition, creation and annihilation of a nanovoid during stretching process were observed for the cases with repulsive NP-polymer interactions. For the case with attractive NP–polymer interactions, a kink behavior on a plot of scattering invariant versus the elongation ratio was observed. From this behavior, it was confirmed that the detection of nanovoids through 2DSPs is effective for attractive NP–polymer interaction.

View all citing articles on Scopus

^☆: This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect (http://www.sciencedirect.com/science/journal/00104655).

View full text

Precise calculation of the local pressure tensor in Cartesian and spherical coordinates in LAMMPS☆

Abstract

Program summary

Introduction

Section snippets

Microscopic definition of the local pressure tensor and its average over a local (micro)volume

Installation instructions

Applications

Benchmark

Summary

Acknowledgments

J. Comput. Phys.

Curr. Opin. Struct. Biol.

J. Appl. Phys.

Molecular Theory of Capillarity

Int. J. Fract.

J. Chem. Phys.

Biochemistry

Proc. R. Soc. Lond. Ser. A

Phys. Rev. E

Mol. Simul.

J. Chem. Phys.

J. Chem. Phys.

J. Comput. Chem.

Adv. Chem. Phys.

Soft Matter