Aidan P. Thompson

General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions

Three distinct forms are derived for the force virial contribution to the pressure and stress tensor of a collection of atoms interacting under periodic boundary conditions. All three forms are written in terms of forces acting on atoms, and so are valid for arbitrary many-body interatomic potentials. All three forms are mathematically equivalent. In the special case of atoms interacting with pair potentials, they reduce to previously published forms. (i) The atom-cell form is similar to the standard expression for the virial for a finite nonperiodic system, but with an explicit correction for interactions with periodic images. (ii) The atom form is particularly suited to implementation in modern molecular dynamics simulation codes using spatial decomposition parallel algorithms. (iii) The group form of the virial allows the contributions to the virial to be assigned to individual atoms.

A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks.

The time evolution of species concentrations in biochemical reaction networks is often modeled using the stochastic simulation algorithm (SSA) [Gillespie, J. Phys. Chem. 81, 2340 (1977)]. The computational cost of the original SSA scaled linearly with the number of reactions in the network. Gibson and Bruck developed a logarithmic scaling version of the SSA which uses a priority queue or binary tree for more efficient reaction selection [Gibson and Bruck, J. Phys. Chem. A 104, 1876 (2000)]. More generally, this problem is one of dynamic discrete random variate generation which finds many uses in kinetic Monte Carlo and discrete event simulation. We present here a constant-time algorithm, whose cost is independent of the number of reactions, enabled by a slightly more complex underlying data structure. While applicable to kinetic Monte Carlo simulations in general, we describe the algorithm in the context of biochemical simulations and demonstrate its competitive performance on small- and medium-size networks, as well as its superior constant-time performance on very large networks, which are becoming necessary to represent the increasing complexity of biochemical data for pathways that mediate cell function.

Physical mechanism of anisotropic sensitivity in pentaerythritol tetranitrate from compressive-shear reaction dynamics simulations

We propose computational protocol (compressive shear reactive dynamics) utilizing the ReaxFF reactive force field to study chemical initiation under combined shear and compressive load. We apply it to predict the anisotropic initiation sensitivity observed experimentally for shocked pentaerythritol tetranitrate single crystals. For crystal directions known to be sensitive we find large stress overshoots and fast temperature increase that result in early bond-breaking processes whereas insensitive directions exhibit small stress overshoot, lower temperature increase, and little bond dissociation. These simulations confirm the model of steric hindrance to shear and capture the thermochemical processes dominating the phenomena of shear-induced chemical initiation.

Reactive Molecular Dynamics Simulations of Shock Through a Single Crystal of Pentaerythritol Tetranitrate

Large-scale molecular dynamics simulations and the reactive force field ReaxFF were used to study shock-induced initiation in crystalline pentaerythritol tetranitrate (PETN). In the calculations, a PETN single crystal was impacted against a wall, driving a shockwave back through the crystal in the [100] direction. Two impact speeds (4 and 3 km/s) were used to compare strong and moderate shock behavior. The primary difference between the two shock strengths is the time required to exhibit the same qualitative behaviors with the lower impact speed lagging behind the faster impact speed. For both systems, the shock velocity exhibits an initial deceleration due to onset of endothermic reactions followed by acceleration due to the onset of exothermic reactions. At long times, the shock velocity reaches a steady value. After the initial deceleration period, peaks are observed in the profiles of the density and axial stress with the strongly shocked system having sharp peaks while the weakly shocked system developed broad peaks due to the slower shock velocity acceleration. The dominant initiation reactions in both systems lead to the formation of NO(2) with lesser quantities of NO(3) and formaldehyde also produced.

Direct molecular simulation of gradient-driven diffusion

Recent work in the active area of grand canonical molecular dynamics methods is first briefly reviewed followed by an overview of the dual control volume grand canonical molecular dynamics (DCV-GCMD) method, designed to enable the dynamic simulation of a system with a steady-state chemical potential gradient. A short review of the methods and systems used to prototype the DCV-GCMD method and its parallel implementation follows. Finally a new, novel implementation of the DCV-GCMD method that enables the establishment of a steady-state chemical potential gradient in a multicomponent system without having to insert or delete one of the components is presented and discussed.

Nonequilibrium molecular dynamics simulation of electro-osmotic flow in a charged nanopore

Nonequilibrium molecular dynamics simulations were performed for Poiseuille and electro-osmotic flow in a charged cylindrical nanopore. The goal was to examine any deviations from continuum flow behavior and to compare and contrast the Poiseuille and electro-osmotic flow situations. The fluid was composed of cationic counterions and nonpolar monatomic solvent molecules. The cylindrical surface of the pore wall was represented by a stochastic scattering boundary condition. The lack of any surface roughness and the computational efficiency of the fluid model enabled the velocity profile near the wall to be measured at very high spatial resolution. The simulation results indicate that both Poiseuille and electro-osmotic flow conform to continuum transport theories except in the first monolayer of fluid at the pore wall. The apparent viscosity in this region was highly nonuniform and exhibited singularities. Despite this, the viscosity profiles obtained from Poiseuille and electro-osmotic flow were in good mu...

Effect of pressure, membrane thickness, and placement of control volumes on the flux of methane through thin silicalite membranes: A dual control volume grand canonical molecular dynamics study

The flux of methane through the straight channels of thin silicalite membranes is studied via dual control volume grand canonical molecular dynamics. The adsorption layers on the surfaces of the thin membranes are found to provide a significant resistance to the flux of methane. This strong surface effect for thin membranes requires that the control volumes (where insertions and deletions are performed) must be placed far enough away from the membrane surface that they do not overlap with the surface adsorption layer. The permeance (flux/pressure drop) of methane through the surface layer is shown to be insensitive to both the average pressure and the pressure drop. In contrast, the permeance through the interior of the membrane increases with decreasing average pressure. These results are explained using a model which treats the transport through the surface barrier as driven by the pressure gradient and transport through the zeolite as driven by the chemical potential gradient. A new force field named D...

Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials

We present a new interatomic potential for solids and liquids called Spectral Neighbor Analysis Potential (SNAP). The SNAP potential has a very general form and uses machine-learning techniques to reproduce the energies, forces, and stress tensors of a large set of small configurations of atoms, which are obtained using high-accuracy quantum electronic structure (QM) calculations. The local environment of each atom is characterized by a set of bispectrum components of the local neighbor density projected onto a basis of hyperspherical harmonics in four dimensions. The bispectrum components are the same bond-orientational order parameters employed by the GAP potential 1]. The SNAP potential, unlike GAP, assumes a linear relationship between atom energy and bispectrum components. The linear SNAP coefficients are determined using weighted least-squares linear regression against the full QM training set. This allows the SNAP potential to be fit in a robust, automated manner to large QM data sets using many bispectrum components. The calculation of the bispectrum components and the SNAP potential are implemented in the LAMMPS parallel molecular dynamics code. We demonstrate that a previously unnoticed symmetry property can be exploited to reduce the computational cost of the force calculations by more than one order of magnitude. We present results for a SNAP potential for tantalum, showing that it accurately reproduces a range of commonly calculated properties of both the crystalline solid and the liquid phases. In addition, unlike simpler existing potentials, SNAP correctly predicts the energy barrier for screw dislocation migration in BCC tantalum.

DIRECT MOLECULAR SIMULATION OF GRADIENT-DRIVEN DIFFUSION OF LARGE MOLECULES USING CONSTANT PRESSURE

Dual control volume grand canonical molecular dynamics (DCV-GCMD) is a boundary-driven nonequilibrium molecular-dynamics technique for simulating gradient-driven diffusion in multicomponent systems. Two control volumes are established at opposite ends of the simulation box. Constant temperature and chemical potential of diffusing species are imposed in the control volumes (i.e., constant-μ1⋯μn−1μnVT). This results in stable chemical potential gradients and steady-state diffusion fluxes in the region between the control volumes. We present results and detailed analysis for a new constant-pressure variant of the DCV-GCMD method in which one of the diffusing species for which a steady-state diffusion flux exists does not have to be inserted or deleted. Constant temperature, pressure, and chemical potential of all diffusing species except one are imposed in the control volumes (i.e., constant-μ1⋯μn−1NnPT). The constant-pressure method can be applied to situations in which insertion and deletion of large molec...

Topology of cyclo-octane energy landscape

Understanding energy landscapes is a major challenge in chemistry and biology. Although a wide variety of methods have been invented and applied to this problem, very little is understood about the actual mathematical structures underlying such landscapes. Perhaps the most general assumption is the idea that energy landscapes are low-dimensional manifolds embedded in high-dimensional Euclidean space. While this is a very mild assumption, we have discovered an example of an energy landscape which is nonmanifold, demonstrating previously unknown mathematical complexity. The example occurs in the energy landscape of cyclo-octane, which was found to have the structure of a reducible algebraic variety, composed of the union of a sphere and a Klein bottle, intersecting in two rings.