Investigating various many-body force fields for their ability to predict reduction in elastic modulus due to vacancies using molecular dynamics simulations

ABSTRACT

Molecular dynamics simulations are more frequently being utilised to predict macroscale mechanical properties as a result of atomistic defects. However, the interatomic force field can significantly affect the resulting mechanical properties. While several studies exist which demonstrate the ability of various force fields to predict mechanical properties, the investigation into which is most accurate for the investigation of vacancies is limited. To obtain meaningful predictions of mechanical properties, a clear understanding of force field parameterisation is required. As such, the current study evaluates various many-body force fields to demonstrate the reduction in mechanical properties of iron and iron–chromium due to the presence of vacancies while undergoing room temperature atomistic uniaxial tension. Reduction was normalised in each case with the zero-vacancy elastic modulus, removing the need to predict an accurate nominal elastic modulus. Comparisons were made to experimental data and an empirical model from literature. It was demonstrated that accurate fitting to vacancy formation and migration energy allowed for accurate predictions. In addition, bond-order based force fields showed enhanced predictions regardless of fitting procedure. Overall, these findings highlight the need to understand capabilities and limitations of available force fields, as well as the need for enhanced parameterisation of force fields.     

Bridging microscopic and mesoscopic simulations of lipid bilayers

A lipid bilayer is modeled using a mesoscopic model designed to bridge atomistic bilayer simulations with macro-scale continuum-level simulation. Key material properties obtained from detailed atomistic-level simulations are used to parameterize the meso-scale model. The fundamental length and time scale of the meso-scale simulation are at least an order of magnitude beyond that used at the atomistic level. Dissipative particle dynamics cast in a new membrane formulation provides the simulation methodology. A meso-scale representation of a dimyristoylphosphatidylcholine membrane is examined in the high and low surface tension regimes. At high surface tensions, the calculated modulus is found to be slightly less than the atomistically determined value. This result agrees with the theoretical prediction that high-strain thermal undulations still persist, which have the effect of reducing the value of the atomistically determined modulus. Zero surface tension simulations indicate the presence of strong thermal undulatory modes, whereas the undulation spectrum and the calculated bending modulus are in excellent agreement with theoretical predictions and experiment.     

Molecular dynamics analysis of the transient temperature increase at void locations in shocked materials: RDX and Cu

Molecular dynamics (MD) simulations of high velocity impact (1–6 km/s) of RDX crystal with a nanometer-sized void, has been carried out to understand the mechanism of increase in temperature at void locations under shock loading. Similar simulations are then carried out on single-crystal copper for better interpretation of the results. A reactive potential that can simulate chemical reactions (ReaxFF) has been used for RDX, whereas an EAM potential has been used for Cu. Increased temperature at the void locations are observed under shock loading. The atomic motion, temperature, average potential energy per atom (PE), and average kinetic energy per atom (KE) in and around the voids are closely monitored in order to understand the reason for temperature increase. We compare our results with existing proposed mechanisms and show that some of the proposed mechanisms are not necessary for increased temperature at a void location. It is shown that the directed particle velocity is efficiently is converted into randomized velocity due to the presence of voids thereby increasing the local temperature transiently. In this initial stage (few picoseconds) of the shock, chemical reactions of energetic materials do not play a part in the temperature rise.     

Unphysical nucleation of diamond in the extended cutoff Tersoff potential

In simulations of carbon materials it is common practice to view the coefficients of the cutoff function as free parameters which can be optimised according to the system of interest. This work examines a common modification to the widely used Tersoff potential in which the coefficient of the upper cutoff is increased to improve the properties of amorphous carbon. Using molecular dynamics simulations, we show that this so-called extended cutoff Tersoff model leads to nucleation of diamond nanocrystals during annealing of amorphous carbon. By varying the density of the system, and examining the radial distribution function in conjunction with the modified cutoff function, we demonstrate that this behaviour is unphysical, and does not represent a new pathway for synthesising diamond. Viewed from a broader perspective, this observation provides a cautionary tale against altering the parameters of empirical potentials without fully considering the wider implications.     

Coarse-Grained Molecular Models of Water: A Review

Coarse-grained (CG) models have proven to be very effective tools in the study of phenomena or systems that involve large time- and length-scales. By decreasing the degrees of freedom in the system and using softer interactions than seen in atomistic models, larger timesteps can be used and much longer simulation times can be studied. CG simulations are widely used to study systems of biological importance that are beyond the reach of atomistic simulation, necessitating a computationally efficient and accurate CG model for water. In this review, we discuss the methods used for developing CG water models and the relative advantages and disadvantages of the resulting models. In general, CG water models differ with regards to how many waters each CG group or bead represents, whether analytical or tabular potentials have been used to describe the interactions, and how the model incorporates electrostatic interactions. Finally, how the models are parameterized depends on their application, so, while some are fitted to experimental properties such as surface tension and density, others are fitted to radial distribution functions extracted from atomistic simulations.     

Stability Curve Prediction of Homologous Proteins Using Temperature-Dependent Statistical Potentials

The unraveling and control of protein stability at different temperatures is a fundamental problem in biophysics that is substantially far from being quantitatively and accurately solved, as it requires a precise knowledge of the temperature dependence of amino acid interactions. In this paper we attempt to gain insight into the thermal stability of proteins by designing a tool to predict the full stability curve as a function of the temperature for a set of 45 proteins belonging to 11 homologous families, given their sequence and structure, as well as the melting temperature () and the change in heat capacity () of proteins belonging to the same family. Stability curves constitute a fundamental instrument to analyze in detail the thermal stability and its relation to the thermodynamic stability, and to estimate the enthalpic and entropic contributions to the folding free energy. In summary, our approach for predicting the protein stability curves relies on temperature-dependent statistical potentials derived from three datasets of protein structures with targeted thermal stability properties. Using these potentials, the folding free energies () at three different temperatures were computed for each protein. The Gibbs-Helmholtz equation was then used to predict the protein''s stability curve as the curve that best fits these three points. The results are quite encouraging: the standard deviations between the experimental and predicted ''s, ''s and folding free energies at room temperature () are equal to 13 , 1.3 ) and 4.1 , respectively, in cross-validation. The main sources of error and some further improvements and perspectives are briefly discussed.     

Mechanical properties of a slow muscle in the cockroach

The mechanical properties of the metacoxal muscle, 177d, in the cockroach, Periplaneta americana, was investigated. The muscle exhibited a mean resting tension of 2.6 ± 1.3g SD. Neurally evoked tension summed with the resting tension and the relaxation phase of the evoked tension varied from less than 1 s to several minutes. This residual tension varied not only in duration but also in amplitude. Stimulation of inhibitory axons increased the rate of relaxation and thereby abolished the residual tension. However, inhibitory stimulation never reduced the resting tension. Stimulation of the main leg nerve at several times the threshold of the inhibitory axons could evoke residual tension. Recording of synaptic potentials from the two histochemically different fiber types (dorsal and ventral groups) revealed large hyperpolarizations in the ventral fibers and decreased duration and amplitude of excitatory potentials in the dorsal fibers. These results suggest that there are a variety of ways in which tension can be evoked, maintained, and controlled in these muscles.     

Tensile and compressive properties of healthy and osteoarthritic human articular cartilage

Osteoarthritis (OA) is a disease affecting articular cartilage and the underlying bone, resulting from many biological and mechanical interacting factors which change the extracellular matrix (ECM) and cells and lead to increasing levels of cartilage degeneration, like softening, fibrillation, ulceration and cartilage loss. The early diagnosis of the disease is fundamental to prevent pain, further tissue degeneration and reduce hospital costs. Although morphological modifications can be detected by modern non-invasive diagnostic techniques, they may not be evident in the early stages of OA. The mechanical properties of articular cartilage are related to its composition and structure and are sensitive to even small changes in the ECM that could occur in early OA. The aim of the present study was to compare the mechanical properties of healthy and OA cartilage using a combined experimental-numerical approach. Experimental assessments consisted of step wise confined and unconfined compression and tension stress relaxation tests on disks (for compression) or strips (for tension) of cartilage obtained from human femoral heads discarded from the operating room after total hip replacement. The numerical model was based on the biphasic theory and included the tension-compression non-linearity. Considering OA samples vs normal samples, the static compressive modulus was 55-68% lower, the permeability was 60-80% higher, the dynamic compressive modulus was 59-64% lower, the static tension modulus was 72-83% lower. The model successfully simulated the experimental tests performed on healthy and OA cartilage and was used in combination with the experimental tests to evaluate the role of different ECM components in the mechanical response of normal and OA cartilage.     

GPU-accelerated molecular dynamics simulation of solid covalent crystals

Graphics processing unit (GPU) is becoming a powerful computational tool in science and engineering. In this paper, different from previous molecular dynamics (MD) simulation with pair potentials and many-body potentials, two MD simulation algorithms implemented on a single GPU are presented to describe a special category of many-body potentials – bond order potentials used frequently in solid covalent materials, such as the Tersoff potentials for silicon crystals. The simulation results reveal that the performance of GPU implementations is apparently superior to their CPU counterpart. Furthermore, the proposed algorithms are generalised, transferable and scalable, and can be extended to the simulations with general many-body interactions such as Stillinger–Weber potential and so on.     

Structural and dynamical properties predicted by reactive force fields simulations for four common pure fluids at liquid and gaseous non-reactive conditions

Four common pure fluids were chosen to elucidate the reliability of reactive force fields in estimating bulk properties of selected molecular systems: CH₄, H₂O, CO₂ and H₂. The pure fluids are not expected to undergo chemical reactions at the conditions chosen for these simulations. The ‘combustion’ ReaxFF was chosen as reactive force field. In the case of water, we also considered the ‘aqueous’ ReaxFF model. The results were compared to data obtained implementing popular classic force fields. In the gas phase, it was found that simulations conducted using the ‘combustion’ ReaxFF formalism yield structural properties in reasonable good agreement with classic simulations for CO₂ and H₂, but not for CH₄ and H₂O. In the liquid phase, ‘combustion’ ReaxFF simulations reproduce reasonably well the structure obtained from classic simulations for CH₄, degrade for CO₂ and H₂, and are rather poor for H₂O. In the gas phase, the simulation results are compared to experimental second virial coefficient data. The ‘combustion’ ReaxFF simulations yield second virial coefficients that are not sufficiently negative for CH₄ and CO₂, and slightly too negative for H₂. The ‘combustion’ ReaxFF parameterisation induces too strong an effective attraction between water molecules, while the ‘aqueous’ ReaxFF yields a second virial coefficient that is in reasonable agreement with experiments. The ‘combustion’ ReaxFF parameterisation yields acceptable self-diffusion coefficients for gas-phase properties of CH₄, CO₂ and H₂. In the liquid phase, the results are good for CO₂, while the self-diffusion coefficient predicted for liquid CH₄ is slower, and that predicted for liquid H₂ is about nine times faster than those expected based on classic simulations. The ‘aqueous’ ReaxFF parameterisation yields good results for both the structure and the diffusion of both liquid and vapour water.     

Molecular Dynamics Study of Zinc-Blende GaN,AIN and InN

Abstract

We present a result of the molecular dynamics calculations with used a three-body empirical Tersoff potential. The parameters of the Tersoff potential are determined for nitride compound semiconductors such as GaN, AlN and InN. The structural and thermodynamic properties of GaN, AlN and InN in zinc-blende structure are presented. We report the equilibrium lattice constants, the bulk moduli, the cubic clastic constants, thermal expansion coefficient and specific heat. Good agreement is obtained with recent experimental and theoretical results for all constants.     

An ultrathin carbon nanoribbon study as a component of nanoelectromechanical devices

We investigated a carbon nanoribbon (CNR) using atomistic simulations based on Tersoff–Brenner potential function. The CNR was obtained from a compressed (5,5) carbon nanotube (CNT). The obtained CNR had a cross-sectional view as a binocular telescope structure composed of both sp² and sp³ bonds. One carbon atom per ten carbon atoms had sp³ bond. For the optimized structures, the residual forces on the CNR were 3-order higher than that on the CNR and the lattice constant of the CNR was higher 0.0624?Å than that of the CNT along the tube axis. The Young's modulus of the CNR was the same as that of the CNT whereas the critical strain of the CNR was significantly lower than that of the CNT because the residual stresses on the CNR was very higher than those on the CNT. The tensile force curve vs. the strain of the CNT was slightly higher than that of the CNR.     

Alzheimer's Aβ(1-40) Amyloid Fibrils Feature Size-Dependent Mechanical Properties

Amyloid fibrils are highly ordered protein aggregates that are associated with several pathological processes, including prion propagation and Alzheimer''s disease. A key issue in amyloid science is the need to understand the mechanical properties of amyloid fibrils and fibers to quantify biomechanical interactions with surrounding tissues, and to identify mechanobiological mechanisms associated with changes of material properties as amyloid fibrils grow from nanoscale to microscale structures. Here we report a series of computational studies in which atomistic simulation, elastic network modeling, and finite element simulation are utilized to elucidate the mechanical properties of Alzheimer''s Aβ(1-40) amyloid fibrils as a function of the length of the protein filament for both twofold and threefold symmetric amyloid fibrils. We calculate the elastic constants associated with torsional, bending, and tensile deformation as a function of the size of the amyloid fibril, covering fibril lengths ranging from nanometers to micrometers. The resulting Young''s moduli are found to be consistent with available experimental measurements obtained from long amyloid fibrils, and predicted to be in the range of 20–31 GPa. Our results show that Aβ(1-40) amyloid fibrils feature a remarkable structural stability and mechanical rigidity for fibrils longer than ≈100 nm. However, local instabilities that emerge at the ends of short fibrils (on the order of tens of nanometers) reduce their stability and contribute to their disassociation under extreme mechanical or chemical conditions, suggesting that longer amyloid fibrils are more stable. Moreover, we find that amyloids with lengths shorter than the periodicity of their helical pitch, typically between 90 and 130 nm, feature significant size effects of their bending stiffness due the anisotropy in the fibril''s cross section. At even smaller lengths (⪅50 nm), shear effects dominate lateral deformation of amyloid fibrils, suggesting that simple Euler-Bernoulli beam models fail to describe the mechanics of amyloid fibrils appropriately. Our studies reveal the importance of size effects in elucidating the mechanical properties of amyloid fibrils. This issue is of great importance for comparing experimental and simulation results, and gaining a general understanding of the biological mechanisms underlying the growth of ectopic amyloid materials.     

A simple embedded atom potential for Pd-H alloys

Metal Hydride systems are an important research topic in material science because of their many practical, industrial, and scientific applications. Therefore, the development of reliable and efficient interatomic potentials for metal hydrides systems, to be utilised in molecular simulations, can be of great value in accelerating the research in this field. In this paper, a fully analytical interatomic Embedded Atom Potential (EAM) for the Pd-H system has been developed, that can be easily extended to ternary Palladium-based hydride systems. The new potential has fewer fitting parameters than previously developed EAM Pd-H potentials and can better predict the cohesive energy, lattice constant, bulk modulus, elastic constants, and the stable alloy crystal structures during molecular dynamics (MD) simulations for various hydrogen concentrations. The EAM potential also predicts the miscibility gap, the separation of a Pd-H system to dilute (α) and concentrated (β) phases.     

A Discrete Electromechanical Model for Human Cardiac Tissue: Effects of Stretch-Activated Currents and Stretch Conditions on Restitution Properties and Spiral Wave Dynamics

We introduce an electromechanical model for human cardiac tissue which couples a biophysical model of cardiac excitation (Tusscher, Noble, Noble, Panfilov, 2006) and tension development (adjusted Niederer, Hunter, Smith, 2006 model) with a discrete elastic mass-lattice model. The equations for the excitation processes are solved with a finite difference approach, and the equations of the mass-lattice model are solved using Verlet integration. This allows the coupled problem to be solved with high numerical resolution. Passive mechanical properties of the mass-lattice model are described by a generalized Hooke''s law for finite deformations (Seth material). Active mechanical contraction is initiated by changes of the intracellular calcium concentration, which is a variable of the electrical model. Mechanical deformation feeds back on the electrophysiology via stretch-activated ion channels whose conductivity is controlled by the local stretch of the medium. We apply the model to study how stretch-activated currents affect the action potential shape, restitution properties, and dynamics of spiral waves, under constant stretch, and dynamic stretch caused by active mechanical contraction. We find that stretch conditions substantially affect these properties via stretch-activated currents. In constantly stretched medium, we observe a substantial decrease in conduction velocity, and an increase of action potential duration; whereas, with dynamic stretch, action potential duration is increased only slightly, and the conduction velocity restitution curve becomes biphasic. Moreover, in constantly stretched medium, we find an increase of the core size and period of a spiral wave, but no change in rotation dynamics; in contrast, in the dynamically stretching medium, we observe spiral drift. Our results may be important to understand how altered stretch conditions affect the heart''s functioning.     

Effect of temperature on tether extraction, surface protrusion, and cortical tension of human neutrophils

Neutrophil rolling on endothelial cells, the initial stage of its migrational journey to a site of inflammation, is facilitated by tether extraction and surface protrusion. Both phenomena have been studied extensively at room temperature, which is considerably lower than human body temperature. It is known that temperature greatly affects cellular mechanical properties such as viscosity. Therefore, we carried out tether extraction, surface protrusion, and cortical tension experiments at 37 degrees C with the micropipette aspiration technique. The experimental temperature was elevated using a custom-designed microscope chamber for the micropipette aspiration technique. To evaluate the constant temperature assumption in our experiments, the temperature distribution in the whole chamber was computed with finite element simulation. Our simulation results showed that temperature variation around the location where our experiments were performed was less than 0.2 degrees C. For tether extraction at 37 degrees C, the threshold force required to pull a tether (40 pN) was not statistically different from the value at room temperature (51 pN), whereas the effective viscosity (0.75 pN.s/microm) decreased significantly from the value at room temperature (1.5 pN.s/microm). Surface protrusion, which was modeled as a linear deformation, had a slightly smaller spring constant at 37 degrees C (40 pN/microm) than it did at room temperature (56 pN/microm). However, the cortical tension at 37 degrees C (5.7+/-2.2 pN/microm) was substantially smaller than that at room temperature (23+/-8 pN/microm). These data clearly suggest that neutrophils roll differently at body temperature than they do at room temperature by having distinct mechanical responses to shear stress of blood flow.     

Study on the anisotropic response of condensed-phase RDX under repeated stress wave loading via ReaxFF molecular dynamics simulation

Anisotropic mechanical response and chemical reaction process of cyclotrimethylene trinitramine (RDX) along crystal orientations were studied with molecular dynamics simulations using ReaxFF potential under repeated stress wave loading. In the simulations, shocks were propagated along the [010], [001], [210], [100], [111], and [102] orientations of crystal RDX at initial particle velocity U_p in the range of 1～4 km/s. For shocks at U_p?≤?2 km/s, local stacking fault and molecular conformational change can only cause marginal temperature and pressure increase without molecular decomposition. As shocks increase to U_p?≥?2.5 km/s, rupture of N-NO₂ bond accompanied by partial HONO elimination dominates the main chemical reactions at the initial stage. The ordering of the follow-up consumption of NO₂ and ring-breaking rate is directly consistent with that of increasing rate in temperature and pressure. The (210) and (100) planes are more sensitive to shocks in temperature and pressure profiles than the (111) plane, which agrees well with experimental observations and theoretical results in the literature. Therefore, the repeated dynamic loading model in conjunction with MD simulation using ReaxFF potential for crystal RDX indicates that these methods can be applied to study the mechanical response and chemical reaction process of polymer bonded explosives that are commonly subjected to compressive and tensile stress waves observed in practice.     

Modelling of the mechanical and mass transport properties of auxetic molecular sieves: an idealised organic (polymeric honeycomb) host–guest system

Force field-based simulations have been employed to model the mechanical properties of a range of undeformed molecular polymeric honeycombs having conventional and re-entrant hexagon pores. The conventional and re-entrant hexagon honeycombs are predicted to display positive and negative in-plane Poisson's ratios, respectively, confirming previous simulations. The structure, and mechanical and mass transport properties of a layered re-entrant honeycomb ((2,8)-reflexyne) were studied in detail for a uniaxial load applied along the x ₂ direction. The mechanical properties are predicted to be stress- (strain-) dependent and the trends can be interpreted using analytical expressions from honeycomb theory. Transformation from negative to positive Poisson's ratio behaviour is predicted at an applied stress of σ₂ = 2 GPa. Simulations of the loading of C₆₀ and C₇₀ guest molecules into the deformed layered (2,8)-reflexyne host framework demonstrate the potential for tunable size selectivity within the host framework. The entrapment and release of guest molecules is attributed to changes in the size and shape of the pores in this host–guest system.     

TENSION AT THE HIGHLY STRETCHED SURFACE OF SEA-URCHIN EGGS*

Unfertilized eggs of the sea urchin, Paracentrotus lividus, were placed between two parallel plates and flattened by a definite force to 20% of their original diameter, with two-fold increase in their surface area. The resulting tension at their surface was calculated from the relation of force and deformation. In spite of this extensive stretching, the tension was found to be not more than 0.2 dyne/cm, while under conditions involving mild stretching (3%) the tension still amounted to 0.12 dyne/cm. These results do not support Mela's theory (7, 8), which predicts a transition of the mechanical properties of the egg surface from a ‘subelastic’ to ‘elastic’ state when the surface is stretched to beyond 34% of its initial area.