Large-scale Molecular Dynamics Research Articles

Shock-induced plasticity and phase transformation in single crystal magnesium: an interatomic potential and non-equilibrium molecular dynamics simulations

An effective and reliable Finnis–Sinclair (FS) type potential is developed for large-scale molecular dynamics (MD) simulations of plasticity and phase transition of magnesium (Mg) single crystals under high-pressure shock loading. The shock-wave profiles exhibit a split elastic–inelastic wave in the [0001]HCP shock orientation and a three-wave structure in the [10-10]HCP and [-12-10]HCP directions, namely, an elastic precursor, a followed plastic front, and a phase-transition front. The shock Hugoniot of the particle velocity (U p) vs the shock velocity (U s) of Mg single crystals in three shock directions under low shock strength reveals apparent anisotropy, which vanishes with increasing shock strength. For the [0001]HCP shock direction, the amorphization caused by strong atomic strain plays an important role in the phase transition and allows for the phase transition from an isotropic stressed state to the product phase. The reorientation in the shock directions [10-10]HCP and [-12-10]HCP, as the primary plasticity deformation, leads to the compressed hexagonal close-packed (HCP) phase and reduces the phase-transition threshold pressure. The phase-transition pathway in the shock direction [0001]HCP includes a preferential contraction strain along the [0001]HCP direction, a tension along [-12-10]HCP direction, an effective contraction and shear along the [10-10]HCP direction. For the [10-10]HCP and [-12-10]HCP shock directions, the phase-transition pathway consists of two steps: a reorientation and the subsequent transition from the reorientation hexagonal close-packed phase (RHCP) to the body-centered cubic (BCC). The orientation relationships between HCP and BCC are (0001)HCP ⟨-12-10⟩HCP // {110}BCC ⟨001⟩BCC. Due to different slipping directions during the phase transition, three variants of the product phase are observed in the shocked samples, accompanied by three kinds of typical coherent twin-grain boundaries between the variants. The results indicate that the highly concentrated shear stress leads to the crystal lattice instability in the elastic precursor, and the plasticity or the phase transition relaxed the shear stress.

Self-Consistent-Field Method for Correlated Many-Electron Systems with an Entropic Cumulant Energy

A self-consistent field method is presented within density matrix functional theory. The computational cost for a correlated many-electron calculation is reduced to that of the self-consistent-field Hartree-Fock method, while the accuracy still reaches that of sophisticated configuration interaction based methods. In this method, the two-electron cumulant energy is measured with an information entropy associated with the Fermi-Dirac distribution of the occupation numbers. An eigenvalue equation for the orbitals is obtained, with the eigenvalues (orbital energies) connected to the occupation numbers through the Fermi-Dirac distribution. The occupation numbers for the strongly occupied orbitals are very close to the natural orbital occupation numbers from wave function methods. It covers in a single scheme the nondynamical correlation in weak or breaking bonds as well as the dynamical correlation at all distances. The method is well suited to large-scale potential energy surface calculation and molecular dynamics simulation.

Crystallographic-orientation-dependence plasticity of niobium under shock compressions

Plastic deformation mechanism of metals under shock compressions is one of longstanding interests in compression science and materials related field. In this work, shock responses of Nb under the shock compressions along [001], [110] and [111] directions are investigated using large-scale nonequilibrium molecular dynamics (NEMD) simulations. A new reliable embedded atom method (EAM) potential is specially developed for purposes of studying the deformation twinning under high pressures. Our results indicate that the shock wave front exhibits split two-wave structures in all three shock direction. In contrast to traditional understandings, the overdriven pressure of [110] is smaller than that of [001], which is attributed to their different twinning threshold and growth speed. The deformation twins are considerably more pronounced for the shock along [001] and [110] directions than that along [111] direction, which take place along {112}〈111〉 systems. Twins nucleate at the shock front and rapidly grow accompanied by dislocation nucleation and multiplication between the twins. Different twinning mechanisms for the shock along [001] and [110] directions are identified. For the shock along [001], the deformation twin is formed through successive movements of atoms on the alternative (11¯2) plane along the [11¯0] and [11¯3¯] direction. For the shock along [110], the atoms on each (112¯) plane directly move along the [111] direction layer by layer. Using a lattice model combined with the transition state theory, we find that uniaxial compressions along specific directions play the key role for the twinning mechanism.

Molecular dynamic simulations of plasticity and phase transition in Mg polycrystalline under shock compression

We have investigated the shock-induced plasticity and phase transition in the hexagonal columnar nanocrystalline Mg by large-scale nonequilibrium molecular dynamics simulations. The preexisting grain boundaries (GBs) induce the nucleation of the {10–12} twins for the local stress relaxation. The twins grow up in grains leading to the orientation rotation. The phase transition from the hexagonal close-packed phase to the body-centered cubic phase begins when the migrating twin GBs meet in A- and C-type grains, and continues in the plastic deformation regions. The phase-transition pathway involves two steps: the reorientation and phase transformation.

Deep learning molecular dynamics simulation on microwave high-temperature dielectric function of silicon nitride

Silicon nitride (<i>β</i>-Si<sub>3</sub>N<sub>4</sub>) is a most promising thermal wave-transparent material. The accurate measurement of its high-temperature dielectric function is essential to solving the “black barrier” problem of hypersonic vehicles and accelerating the design of silicon nitride-based thermal wave-transparent materials. Direct experimental measurement at high temperature is a difficult job and the accuracy of classical molecular dynamics (CMD) simulations suffers the choice of empirical potential. In this work, we build a <i>β</i>-Si<sub>3</sub>N<sub>4</sub> model on a nanoscale, train the deep learning potential (DLP) by using first-principles data, and apply the deep potential molecular dynamics (DPMD) to simulate the polarization relaxation process. The predicted energy and force by DLP are excellently consistent with first-principles calculations, which proves the high accuracy of DLP. The RMSEs for <i>β</i>-Si<sub>3</sub>N<sub>4</sub> are quite low (0.00550 meV/atom for energy and 7.800 meV/Å for force). According to the Cole-Cole formula, the microwave dielectric function in the temperature range of 300–1000 K is calculated by using the deep learning molecular dynamics method. Compared with the empirical potential, the computational results of the DLP are consistent with the experimental results in the sense of order of magnitude. It is also found that the DPMD performs well in terms of computational speed. In addition, a mathematical model of the temperature dependence of the relaxation time is established to reveal the pattern of relaxation time varying with temperature. The high-temperature microwave dielectric function of silicon nitride is calculated by implementing large-scale and high-precision molecular dynamics simulations. It provides fundamental data for promoting the application of silicon nitride in high-temperature thermal transmission.

Structures and localized vibrational states of defects in graphite by tight-binding calculations.

The structural and vibrational properties of pristine graphite and point defects in graphite are studied by tight-binding (TB) calculations using a three-center TB potential model. We showed that the three-center TB potential without "ad hoc" van der Waals interaction corrections can accurately describe the inter-layer distance of graphite and the lowest-energy structures and stabilities of typical point defects in graphite. The results from our TB calculations are in good agreement with those from density-functional theory calculations with van der Waals interaction corrections. We also investigated the vibrational properties to gain better understanding on the localization of vibrational states induced by the point defects. Our calculation results show that although localized or quasi-localized vibrational modes can be found in all defected graphite, the localization induced by Frenkel pair, dual-vacancy, and dual-interstitial defects is much stronger. Atomic displacements associated with the localized vibrational modes induced by these three point defects are also analyzed.

SISSO-assisted prediction and design of mechanical properties of porous graphene with a uniform nanopore array.

Mechanical properties of porous graphene can be effectively tuned by tailoring the nanopore arrangement. Knowledge of the relationship between the porous structure and overall mechanical properties is thus essential for the wide potential applications, and the existing challenge is to efficiently predict and design the mechanical properties of porous graphene due to the diverse nanopore arrangements. In this work, we report on how the SISSO (Sure Independence Screening and Sparsifying Operator) algorithm can be applied to build a bridge between the mechanical properties of porous graphene and the uniform nanopore array. We first construct a database using the strength and work of fracture calculated by large-scale molecular dynamics simulations. Then the SISSO algorithm is adopted to train a predictive model and automatically derive the optimal fitting formulae which explicitly describe the nonlinear structure–property relationships. These expressions not only enable the direct and accurate prediction of targeted properties, but also serve as a convenient and portable tool for inverse design of the porous structure. Compared with other forecasting methods including several popular machine learning algorithms, the SISSO algorithm shows its advantages in both accuracy and convenience.

Elucidation of the conformational dynamics and assembly of Argonaute–RNA complexes by distinct yet coordinated actions of the supplementary microRNA

Argonaute (AGO) proteins, the core of RNA-induced silencing complex, are guided by microRNAs (miRNAs) to recognize target RNA for repression. The miRNA–target RNA recognition forms initially through pairing at the seed region while the additional supplementary pairing can enhance target recognition and compensate for seed mismatch. The extension of miRNA lengths can strengthen the target affinity when pairing both in the seed and supplementary regions. However, the mechanism underlying the effect of the supplementary pairing on the conformational dynamics and the assembly of AGO–RNA complex remains poorly understood. To address this, we performed large-scale molecular dynamics simulations of AGO–RNA complexes with different pairing patterns and miRNA lengths. The results reveal that the additional supplementary pairing can not only strengthen the interaction between miRNA and target RNA, but also induce the increased plasticity of the PAZ domain and enhance the domain connectivity among the PAZ, PIWI, N domains of the AGO protein. The strong community network between these domains tightens the mouth of the supplementary chamber of AGO protein, which prevents the escape of target RNA from the complex and shields it from solvent water attack. Importantly, the inner stronger matching pairs between the miRNA and target RNA can compensate for weaker mismatches at the edge of supplementary region. These findings provide guidance for the design of miRNA mimics and anti-miRNAs for both clinical and experimental use and open the way for further engineering of AGO proteins as a new tool in the field of gene regulation.

Probabilistic Occlusion Culling using Confidence Maps for High-Quality Rendering of Large Particle Data.

Achieving high rendering quality in the visualization of large particle data, for example from large-scale molecular dynamics simulations, requires a significant amount of sub-pixel super-sampling, due to very high numbers of particles per pixel. Although it is impossible to super-sample all particles of large-scale data at interactive rates, efficient occlusion culling can decouple the overall data size from a high effective sampling rate of visible particles. However, while the latter is essential for domain scientists to be able to see important data features, performing occlusion culling by sampling or sorting the data is usually slow or error-prone due to visibility estimates of insufficient quality. We present a novel probabilistic culling architecture for super-sampled high-quality rendering of large particle data. Occlusion is dynamically determined at the sub-pixel level, without explicit visibility sorting or data simplification. We introduce confidence maps to probabilistically estimate confidence in the visibility data gathered so far. This enables progressive, confidence-based culling, helping to avoid wrong visibility decisions. In this way, we determine particle visibility with high accuracy, although only a small part of the data set is sampled. This enables extensive super-sampling of (partially) visible particles for high rendering quality, at a fraction of the cost of sampling all particles. For real-time performance with millions of particles, we exploit novel features of recent GPU architectures to group particles into two hierarchy levels, combining fine-grained culling with high frame rates.

The mechanism of plasticity and phase transition in iron single crystals under cylindrically divergent shock loading

Cylindrically divergent shock loading significantly influences the plasticity and phase transition of iron. In this study, the mechanisms of plasticity and α→ε phase transition are investigated using large-scale nonequilibrium molecular dynamics simulation. We find that iron undergoes a local bcc–hcp–bcc–hcp (bcc: body-centered cubic; hcp: hexagonal close-packed) phase transition under cylindrically divergent shock loading with cylinder axis lying in [001] and [110] directions. A bcc lattice reorientation driven by resolved shear stress occurs during the hcp–bcc reverse transition and results in bcc twins, which transform into hcp twin structures with an unexpected c-axis, as compared with planar shock. Moreover, we analyze the coupling behavior between plasticity and phase transition and find that collective movement and rearrangement of atoms induced by plastic slip on the shuffle plane play a key role in the activation of phase transition, which significantly differs from that of planar shock.

Insights into lithium manganese oxide-water interfaces using machine learning potentials.

Unraveling the atomistic and the electronic structure of solid-liquid interfaces is the key to the designof new materials for many important applications, from heterogeneous catalysis to battery technology. Density functional theory (DFT) calculations can, in principle, provide a reliable description of such interfaces, but the high computational costs severely restrict the accessible time and length scales. Here, we report machine learning-driven simulations of various interfaces between water and lithium manganese oxide (LixMn2O4), an important electrode material in lithium ion batteries and a catalyst for the oxygen evolution reaction. We employ a high-dimensional neural network potential to compute the energies and forces several orders of magnitude faster than DFT without loss in accuracy. In addition, a high-dimensional neural network for spin prediction is utilized to analyze the electronic structure of the manganese ions. Combining these methods, a series of interfaces is investigated by large-scale molecular dynamics. The simulations allow us to gain insights into a variety of properties, such as the dissociation of water molecules, proton transfer processes, and hydrogen bonds, as well as the geometric and electronic structure of the solid surfaces, including the manganese oxidation state distribution, Jahn-Teller distortions, and electron hopping.

Heat transport across graphene/hexagonal-BN tilted grain boundaries from phase-field crystal model and molecular dynamics simulations

We study the interfacial thermal conductance of grain boundaries (GBs) between monolayer graphene and hexagonal boron nitride (h-BN) sheets using a combined atomistic approach. First, realistic samples containing graphene/h-BN GBs with different tilt angles are generated using the phase-field crystal model developed recently [P. Hirvonen et al., Phys. Rev. B 100, 165412 (2019)] that captures slow diffusive relaxation inaccessible to molecular dynamics (MD) simulations. Then, large-scale MD simulations using the efficient GPUMD package are performed to assess heat transport and rectification properties across the GBs. We find that lattice mismatch between the graphene and h-BN sheets plays a less important role in determining the interfacial thermal conductance as compared to the tilt angle. In addition, we find no significant thermal rectification effects for these GBs.

A gradient Eshelby force on twinning partial dislocations and associated detwinning mechanism in gradient nanotwinned metals

It is well known that the driving force on dislocation glide is the Peach-Koehler force which is proportional to the resolved shear stress on the slip plane. Here, we report a type of configurational force, referred to as the gradient Eshelby force, that can drive the motion of twinning partial dislocations on twin boundaries (TBs) in the absence of any resolved shear stress and cause detwinning in gradient nanotwinned (GNT) metals, an emerging class of multiscale metallic materials with exceptional mechanical properties and novel deformation mechanisms. Specifically, we consider the Eshelby-force-driven motion of twinning partial dislocations and associated detwinning mechanism in GNT metals made of preferentially oriented TBs and columnar grain boundaries (GBs) with spatially varying twin spacing and grain size. Large-scale molecular dynamics simulations validate the proposed Eshelby-driven detwinning mechanism, where twinning partial dislocations are nucleated from GBs and glide on TBs near the region with the steepest local gradient, leading to extensive TB migration and twin annihilation. This study demonstrates the important role of Eshelby force in controlling twinning partial glide and TB migration in GNT metals, which may have broad implications on plastic deformation in gradient nanostructured metals.

Cluster Fragments in Amorphous Phosphorus and their Evolution under Pressure.

Amorphous phosphorus (a-P) has long attracted interest because of its complex atomic structure, and more recently as an anode material for batteries. However, accurately describing and understanding a-P at the atomistic level remains a challenge. Here, it is shown that large-scale molecular-dynamics simulations, enabled by a machine-learning (ML)-based interatomic potential for phosphorus, can give new insights into the atomic structure of a-P and how this structure changes under pressure. The structural model so obtained contains abundant five-membered rings, as well as more complex seven- and eight-atom clusters. Changes in the simulated first sharp diffraction peak during compression and decompression indicate a hysteresis in the recovery of medium-range order. An analysis of cluster fragments, large rings, and voids suggests that moderate pressure (up to about 5GPa) does not break the connectivity of clusters, but higher pressure does. The work provides a starting point for further computational studies of the structure and properties of a-P, and more generally it exemplifies how ML-driven modeling can accelerate the understanding of disordered functional materials.

Understanding the phase transformation mechanisms that affect the dynamic response of Fe-based microstructures at the atomic scales

Large-scale molecular dynamics (MD) simulations were carried out to investigate the shock-induced evolution of microstructure in Fe-based systems comprising single-crystal and layered Cu/Fe alloys with a distribution of interfaces. The shock compression of pure single-crystal Fe oriented along [110] above a threshold pressure results in a BCC (α) → HCP (ɛ) phase transformation behavior that generates a distribution of ε phase variants in the phase transformed region of the microstructure behind the shock front. The propagation of the release wave through a phase transformed ε phase causes a reverse ɛ → α phase transformation and renders a distribution of twins for the [110] oriented Fe that serve as void nucleation sites during spall failure. The simulations reveal that the α → ɛ → α transformation-induced twinning for shock loading along the [110] direction is due to a dominant ɛ phase variant formed during compression that rotates on the arrival of the release wave followed by a reverse phase transformation to twins in the α phase. The modifications in the evolution of the ɛ phase variants and twins in Fe behavior are also studied for Cu–Fe layered microstructures due to the shock wave interactions with the Cu/Fe interfaces using a newly constructed Cu–Fe alloy potential. The MD simulations suggest that interfaces affect the observed variants during shock compression and, hence, distributions of twins during shock release that affects the void nucleation stresses in the Fe phase of Cu/Fe microstructures.

Novel inhibitors to ADP ribose phosphatase of SARS-CoV-2 identified by structure-based high throughput virtual screening and molecular dynamics simulations

The outbreak of a new coronavirus (SARS-CoV-2) was first identified in Wuhan, People's Republic of China, in 2019, which has led to a severe, life-threatening form of pneumonia (COVID-19). Research scientists all around the world have been trying to find small molecule drugs to treat COVID-19. In the present study, a conserved macrodomain, ADP Ribose phosphatase (ADRP), of a critical non-structural protein (Nsp3) in all coronaviruses was probed using large-scale Molecular Dynamics (MD) simulations to identify novel inhibitors. In our virtual screening workflow, the recently-solved X-ray complex structure, 6W6Y, with a substrate-mimics was used to screen 17 million ZINC15 compounds using drug property filters and Glide docking scores. The top twenty output compounds each underwent 200 ns MD simulations (i.e. 20 × 200 ns) to validate their individual stability as potential inhibitors. Eight out of the twenty compounds showed stable binding modes in the MD simulations, as well as favorable drug properties from our predctions. Therefore, our computational data suggest that the resulting top eight out of twenty compounds could potentially be novel inhibitors to ADRP of SARS-CoV-2.

Overall mechanism of JP-10 pyrolysis unraveled by large-scale reactive molecular dynamics simulation

This paper reports the overall reaction mechanism of JP-10 pyrolysis obtained in large-scale reactive molecular dynamics simulations employing the force field of ReaxFF CHO-2008. The C<5 products obtained from the simulations are basically consistent with what detected using varied experimental techniques reported in literature. Particularly, the product evolution tendency of methane, ethane, ethylene, propylene, acetylene, allene,1−butene, propyne, 1,3-butadiene, cyclopentadiene with temperature is in good agreement with the results of single-pulse shock tube experiment performed in this work. The simulation results indicate that JP-10 pyrolysis reactions can be divided into three stages on basis of the evolution of total radicals, namely, the initial ring-openning through CC bond homolysis in stage I, the initiation and growth of chain radicals through β-scission reactions as well as some of the chain propagation via CH bond scission in stage II, the chain propagation through CH bond cleavage and the chain termination reactions in stage III. The C5H7• and C3H5• radicals were found playing dominant role in the reaction propagation. The obtained initial reaction sites of JP-10 pyrolysis are ranked as: B3>B4>B2>B5>B1>B6>B7 that will not be significantly affected by temperature. The obtained aromatic hydrocarbons in simulations and the increasing C/H ratio with temperature and time during the three pyrolysis stages exhibit the coking tendency during JP-10 pyrolysis.

Self-Assembly of NaOL-DDA Mixtures in Aqueous Solution: A Molecular Dynamics Simulation Study.

The self-assembly behaviors of sodium oleate (NaOL), dodecylamine (DDA), and their mixtures in aqueous solution were systematically investigated by large-scale molecular dynamics simulations, respectively. The interaction mechanisms between the surfactants, as well as the surfactants and solvent, were revealed via the radial distribution function (RDF), cluster size, solvent-accessible surface area (SASA), hydrogen bond, and non-bond interaction energy. Results showed that the molecules more easily formed aggregates in mixed systems compared to pure systems, indicating higher surface activity. The SASA values of DDA and NaOL decreased significantly after mixing, indicating a tighter aggregation of the mixed surfactants. The RDF results indicated that DDA and NaOL strongly interacted with each other, especially in the mixed system with a 1:1 molar ratio. Compared to van der Waals interactions, electrostatic interactions between the surfactant molecules were the main contributors to the improved aggregation in the mixed systems. Besides, hydrogen bonds were found between NaOL and DDA in the mixed systems. Therefore, the aggregates in the mixed systems were much more compact in comparison with pure systems, which contributed to the reduction of the repulsive force between same molecules. These findings indicated that the mixed NaOL/DDA surfactants had a great potential in application of mineral flotation.

Ionic screening in bulk and under confinement.

Recent experiments have shown that the repulsive force between atomically flat, like-charged surfaces confining room-temperature ionic liquids or concentrated electrolytes exhibits an anomalously large decay length. In our previous publication [J. Zeman, S. Kondrat, and C. Holm, Chem. Commun. 56, 15635 (2020)], we showed by means of extremely large-scale molecular dynamics simulations that this so-called underscreening effect might not be a feature of bulk electrolytes. Herein, we corroborate these findings by providing additional results with more detailed analyses and expand our investigations to ionic liquids under confinement. Unlike in bulk systems, where screening lengths are computed from the decay of interionic potentials of mean force, we extract such data in confined systems from cumulative charge distributions. At high concentrations, our simulations show increasing screening lengths with increasing electrolyte concentration, consistent with classical liquid state theories. However, our analyses demonstrate that-also for confined systems-there is no anomalously large screening length. As expected, the screening lengths determined for ionic liquids under confinement are in good quantitative agreement with the screening lengths of the same ionic systems in bulk. In addition, we show that some theoretical models used in the literature to relate the measured screening lengths to other observables are inapplicable to highly concentrated electrolytes.

Principles governing control of aggregation and dispersion of aqueous graphene oxide

Controlling the structure of graphene oxide (GO) phases and their smaller analogues, graphene (oxide) quantum dots (GOQDs), is vitally important for any of their widespread intended applications: highly ordered arrangements of nanoparticles for thin-film or membrane applications of GO, dispersed nanoparticles for composite materials and three-dimensional porous arrangements for hydrogels. In aqueous environments, it is not only the chemical composition of the GO flakes that determines their morphologies; external factors such as pH and the coexisting cations also influence the structures formed. By using accurate models of GO that capture the heterogeneity of surface oxidation and very large-scale coarse-grained molecular dynamics that can simulate the behaviour of GO at realistic sizes of GOQDs, the driving forces that lead to the various morphologies in aqueous solution are resolved. We find the morphologies are determined by a complex interplay between electrostatic, {pi }–{pi } and hydrogen bonding interactions. Assembled morphologies can be controlled by changing the degree of oxidation and the pH. In acidic aqueous solution, the GO flakes vary from fully aggregated over graphitic domains to partial aggregation via hydrogen bonding between hydroxylated domains, leading to the formation of planar extended flakes at high oxidation ratios and stacks at low oxidation ratios. At high pH, where the edge carboxylic acid groups are deprotonated, electrostatic repulsion leads to more dispersion, but a variety of aggregation behaviour is surprisingly still observed: over graphitic regions, via hydrogen bonding and “face-edge” interactions. Calcium ions cause additional aggregation, with a greater number of “face-face” and “edge-edge” aggregation mechanisms, leading to irregular aggregated structures. “Face-face” aggregation mechanisms are enhanced by the GO flakes possessing distinct domains of hydroxylated and graphitic regions, with {pi }–{pi } and hydrogen bonding interactions prevalent between these regions on aggregated flakes respectively. These findings furnish explanations for the aggregation characteristics of GO and GOQDs, and provide computational methods to design directed synthesis routes for self-assembled and associated applications.