Large-scale Molecular Dynamics Research Articles

Proton Collective Quantum Tunneling Induces Anomalous Thermal Conductivity of Ice under Pressure.

Proton tunneling is believed to be nonlocal in ice, but its range has been shown to be limited to only a few molecules. Here, we measured the thermal conductivity of ice under pressure up to 50GPa and found it increases with pressure until 20GPa but decreases at higher pressures. We attribute this nonmonotonic thermal conductivity to the collective tunneling of protons at high pressures, supported by large-scale quantum molecular dynamics simulations. The collective tunneling loops span several picoseconds in time and are as large as nanometers in space, which match the phonon periods and wavelengths, leading to strong phonon scattering at high pressures. Our results show direct evidence of global quantum motion existing in high-pressure ice and provide a new perspective to understanding the coupling between phonon propagation and atomic tunneling.

Orientation-dependent multi-spall performance of monocrystalline NiTi alloys under shock compression

The primary objective of this work was committed to large-scale non-equilibrium molecular dynamics modeling for the investigation on multiple spall performance, microstructural evolutions, and melting state of NiTi alloys along three typical crystal orientations. We discovered that there was a single-wave structure present in the [100] orientation, whereas elastic-plastic two-wave structures were found in [110] and [111], exhibiting a significant anisotropy in elastic-plastic wave velocity. Compared to [110] and [111] orientations, multi-spallation was more likely to be generated in the [100] direction. During accumulative dynamic damage, there was also a strong sensitivity of void evolution to orientation. For [100], voids tended to be densely distributed in the material and exhibited a strip-shaped manner. In contrast, a more scattered void distribution was observed for the other two orientations, showing a flaky pattern. Interestingly, dislocations in the spall region were almost annihilated before approaching void nucleation. According to void statistics, the peak number of voids, void nucleation rate, and spall region size were identified to be in a good uniformity along [110], [100], and [111]. For all the cases, the samples experienced partial melting in the released state and were followed by a recovery of structural order degree. However, there was an insignificant anisotropy in the melting state for each crystal orientation during multi-spallation.

Rational Design of Dual-Functionalized Gd@C82 Nanoparticles to Relieve Neuronal Cytotoxicity in Alzheimer's Disease via Inhibition of Aβ Aggregation.

The accumulation of amyloid-β (Aβ) peptides is a major hallmark of Alzheimer's disease (AD) and plays a crucial role in its pathogenesis. Particularly, the structured oligomeric species rich in β-sheet formations were implicated in neuronal organelle damage. Addressing this formidable challenge requires identifying candidates capable of inhibiting peptide aggregation or disaggregating preformed oligomers for effective antiaggregation-based AD therapy. Here, we present a dual-functional nanoinhibitor meticulously designed to target the aggregation driving force and amyloid fibril spatial structure. Leveraging the exceptional structural stability and facile tailoring capability of endohedral metallofullerene Gd@C82, we introduce desired hydrogen-binding sites and charged groups, which are abundant on its surface for specific designs. Impressively, these designs endow the resultant functionalized-Gd@C82 nanoparticles (f-Gd@C82 NPs) with high capability of redirecting peptide self-assembly toward disordered, off-pathway species, obstructing the early growth of protofibrils, and disaggregating the preformed well-ordered protofibrils or even mature Aβ fibrils. This results in considerable alleviation of Aβ peptide-induced neuronal cytotoxicity, rescuing neuronal death and synaptic loss in primary neuron models. Notably, these modifications significantly improved the dispersibility of f-Gd@C82 NPs, thus substantially enhancing its bioavailability. Moreover, f-Gd@C82 NPs demonstrate excellent cytocompatibility with various cell lines and possess the ability to penetrate the blood-brain barrier in mice. Large-scale molecular dynamics simulations illuminate the inhibition and disaggregation mechanisms. Our design successfully overcomes the limitations of other nanocandidates, which often overly rely on hydrophobic interactions or photothermal conversion properties, and offers a viable direction for developing anti-AD agents through the inhibition and even reversal of Aβ aggregation.

Micro-jetting: Areal density calculation from a triangular groove

We present a method for calculating the mass ejected during the reflection of a shock wave on a triangular groove. This calculation is based on the combination of two models taken from the literature, BMPT-2, on the one hand, for the calculation of the velocity and density of the jet, and fragmentation zone propagation (FZP ) on the other hand, for the calculation of the ejected mass, certain parameters required for FZP being determined by BMPT-2. Compared with previous work, FZP has been extended to deal with the various stages of jet formation. The approach was first evaluated on tin using a large-scale molecular dynamics simulation. This first step validated the overall phenomenology and the associated theoretical tools, and enabled us to propose a procedure for adjusting FZP. Next, we used the BMPT-2/FZP combination to analyze ejected mass measurement experiments using Asay foil. The areal mass curves are well reproduced with few parameters, showing that there is no inconsistency between BMPT-2/FZP and the experiments. Finally, a more detailed analysis of the results obtained enables us to set the limits of the jet thickness at the moment of rupture, and to propose a simple analytical form of its profile compatible with the model used.

Molecular Dynamics Simulation of the Interaction within the Carboxymethyl Cellulose-Modified Montmorillonite Lamellae at Aggressive CuCl2 Concentrations.

To elucidate the degradation mechanism of the CMC-modified MMT composite at aggressive Cu2+ concentrations, large scale molecular dynamics simulation was conducted for CuCl2 concentrations ranging from 0 to 800 mM. Both CMC and MMT followed the Langmuir isotherm for Cu2+ adsorption, and the adsorption capacity of CMC (8.75 mmol/g) was much higher than that of MMT (0.83 mmol/g). Despite the CMC mass ratio being only 4.1%, it adsorbed up to 34.3% of the total adsorbed Cu2+. The Cu2+ attraction ability hierarchy of oxygen-containing functional groups in the CMC is as follows: carboxylic oxygens > alcoholic oxygens > carbinolic oxygens > bridging oxygens > glucose oxygens. Carboxyls were the most effective in chelating and complexing with Cu2+, and they could be intentionally added in artificially synthesized polymer-MMT composites for Cu2+ containment. Formation of the Cu2+ cation bridge between CMC and MMT at aggressive CuCl2 concentrations contributed to the transition of CMC density distribution from unimodality to bimodality and enhanced resistance of polymer elution. As the CuCl2 concentration increased, the stoichiometric ratio between the chelated Cu2+ and carboxylic oxygens increased from 1:2 to 1:1, suggesting the evolution of the Cu2+ chelation mechanism.

Unraveling the Hall-Petch to inverse Hall-Petch transition in nanocrystalline high entropy alloys under shock loading

The transition from Hall-Petch (HP) to inverse Hall-Petch (IHP) behaviors associated with grain size reduction has been recognized for over two decades. However, the underlying mechanisms for such transition in high entropy alloys (HEAs) under dynamic loading, in which abundant deformation mechanisms could be activated either sequentially or simultaneously, remain unclear. Here, we investigate the HP to IHP transition in nanocrystalline CoCrFeMnNi HEAs under shock loading by examining their deformation mechanisms and flow stresses using large-scale molecular dynamics (MD) simulations. It is found that this transition is strongly dependent on the shock pressure as a result of the complex interplay among multiple competing deformation mechanisms, including the hardening mechanisms such as dislocations interactions and grain boundary (GB) blocking, as well as the softening mechanisms like phase formation, amorphization, GB thickening, and grain rotation. Moreover, there exists a critical shock pressure, which corresponds to the largest critical grain size for the HP-IHP transition. Below the critical shock pressure, the critical grain size increases with pressure due to a stronger hardening effect in grain interior (GIs), while above the critical pressure, the critical grain size first decreases and then undergoes a pressure-insensitive plateau before further decrease due to softening effects in GIs. A theoretical model that includes different deformation mechanisms is proposed for the first time to capture the shock pressure-dependent HP-IHP transition. Our work provides valuable guidance for optimizing the grain size of nanocrystalline HEAs for applications involving dynamic loadings.

Enhancing surface strength of tungsten by gradient nano-grained structure

A gradient nano-grained (GNG) structure demonstrates satisfactory surface strength. However, the underlying mechanism responsible for its strengthening lacks sufficient research. To explain how gradient nano-grained structures improve surface strength in detail, large-scale parallel molecular dynamics simulations are utilized in this study to investigate the mechanical deformation behavior of BCC tungsten with varying grain sizes during spherical nanoindentation. The findings suggest that a well-designed gradient structure can promote rational plasticity and an appropriate distribution of internal atomic stress. The critical point of maximum stress and hardness is observed when the initial grain size is 4.5 nm, with an average grain size of 7.1 nm. The interaction between grain boundary slip and migration in small grains, along with the enhanced activity of grain boundary dislocations in large grains, collectively contributes to the enhancement of the strength and hardness of the GNG structure. Compared with a homogeneous nano-grained structure, the gradient nano-grained structure exhibits a more rational distribution of dislocations and stress relaxation effects to enhance strength. The present work utilizes the molecular dynamics nanoindentation method to study GNG materials, providing a methodology for investigating the surface strengthening effects of GNG structures at the atomic scale and effectively revealing potential mechanisms for resisting surface deformation in GNG structures.

Strain Engineering of Ion-Coordinated Nanochannels in Nanocellulose.

Expanding the interlayer spacing plays a significant role in improving the conductivity of a cellulose-based conductor. However, it remains a challenge to regulate the cellulose nanochannel expanded by ion coordination. Herein, starting from multiscale mechanics, we proposed a strain engineering method to regulate the interlayer spacing of the cellulose nanochannels. First-principles calculations were conducted to select the most suitable ions for coordination. Large-scale molecular dynamics simulations were performed to reveal the mechanism of interlayer spacing expansion by the ion cross-linking. Combining the shear-lag model, we established the relationship between interfacial cross-link density and interlayer spacing of an ion-coordinated cellulose nanochannel. Consequently, fast ion transport and current regulation were realized via the strain engineering of nanochannels, which provides a promising strategy for the current regulation of a cellulose-based conductor.

MolDStruct: Modeling the dynamics and structure of matter exposed to ultrafast x-ray lasers with hybrid collisional-radiative/molecular dynamics.

We describe a method to compute photon-matter interaction and atomic dynamics with x-ray lasers using a hybrid code based on classical molecular dynamics and collisional-radiative calculations. The forces between the atoms are dynamically determined based on changes to their electronic occupations and the formation of a free electron cloud created from the irradiation of photons in the x-ray spectrum. The rapid transition from neutral solid matter to dense plasma phase allows the use of screened potentials, reducing the number of non-bonded interactions. In combination with parallelization through domain decomposition, the hybrid code handles large-scale molecular dynamics and ionization. This method is applicable for large enough samples (solids, liquids, proteins, viruses, atomic clusters, and crystals) that, when exposed to an x-ray laser pulse, turn into a plasma in the first few femtoseconds of the interaction. We present four examples demonstrating the applicability of the method. We investigate the non-thermal heating and scattering of bulk water and damage-induced dynamics of a protein crystal using an x-ray pump-probe scheme. In both cases, we compare to the experimental data. For single particle imaging, we simulate the ultrafast dynamics of a methane cluster exposed to a femtosecond x-ray laser. In the context of coherent diffractive imaging, we study the fragmentation as given by an x-ray pump-probe setup to understand the evolution of radiation damage in the time range of hundreds of femtoseconds.

Low-index mesoscopic surface reconstructions of Au surfaces using Bayesian force fields

Metal surfaces have long been known to reconstruct, significantly influencing their structural and catalytic properties. Many key mechanistic aspects of these subtle transformations remain poorly understood due to limitations of previous simulation approaches. Using active learning of Bayesian machine-learned force fields trained from ab initio calculations, we enable large-scale molecular dynamics simulations to describe the thermodynamics and time evolution of the low-index mesoscopic surface reconstructions of Au (e.g., the Au(111)-‘Herringbone,’ Au(110)-(1 × 2)-‘Missing-Row,’ and Au(100)-‘Quasi-Hexagonal’ reconstructions). This capability yields direct atomistic understanding of the dynamic emergence of these surface states from their initial facets, providing previously inaccessible information such as nucleation kinetics and a complete mechanistic interpretation of reconstruction under the effects of strain and local deviations from the original stoichiometry. We successfully reproduce previous experimental observations of reconstructions on pristine surfaces and provide quantitative predictions of the emergence of spinodal decomposition and localized reconstruction in response to strain at non-ideal stoichiometries. A unified mechanistic explanation is presented of the kinetic and thermodynamic factors driving surface reconstruction. Furthermore, we study surface reconstructions on Au nanoparticles, where characteristic (111) and (100) reconstructions spontaneously appear on a variety of high-symmetry particle morphologies.

Orientation-dependent deformation mechanisms of alpha-uranium single crystals under shock compression

Large-scale non-equilibrium molecular dynamics (NEMD) simulations were employed to investigate the dynamic deformations of alpha-uranium (α-U) single crystals subjected to varying shock strengths along low-index crystallographic orientations. The pronounced anisotropy of α-U gives rise to a complex microstructural evolution under shock loading. In-depth microstructural analysis of post-shock specimens reveals the identification of multiple dynamic deformation mechanisms. Notably, when the shock loading direction aligns with the a-axis, dynamic deformation of the α-U single crystals is primarily dominated by lattice instability, which attributes to a crystalline-to-amorphous transition serving as the dominant shear stress relaxation pathway. On the other hand, shock loading along the b-axis results in an abundance of deformation twins, with twinning planes identified as (130) and (13¯0). During the twinning event, the α-U matrix undergoes a transition to a metastable intermediate phase, subsequently decomposing into a composite structure comprising α-U twins and matrix. This unconventional twinning mechanism significantly deviates from classical theories. Furthermore, upon loading along the c-axis, twinning and a phase transition from α-U to body-centered tetragonal phase (bct-U) occur in α-U single crystal samples. Given that the pressure threshold of this phase transition predicted by ab initio calculations is as high as ∼270 GPa, the phase transition from α-U to bct-U might be implausible. An alternative interatomic potential of uranium with the higher pressure threshold was employed to reinvestigate the shock response of α-U single crystals along the c-axis. The phase transition of α-U to bct-U disappears, and twinning dominates the plastic deformation, with the twinning orientation conforming to the {112} twinning. The strong anisotropy of the α-U lattice triggers a wealth of orientation-dependent dynamic deformation mechanisms. The activation of the twinning system is evidently associated with the loading direction, constituting the potential cause for the discovery of multiple twinning variants during the deformation in polycrystalline uranium.

Exploring the Interactions between two Ligands, UCB-J and UCB-F, and Synaptic Vesicle Glycoprotein 2 Isoforms.

In silico modeling was applied to study the efficiency of two ligands, namely, UCB-J and UCB-F, to bind to isoforms of the synaptic vesicle glycoprotein 2 (SV2) that are involved in the regulation of synaptic function in the nerve terminals, with the ultimate goal to understand the selectivity of the interaction between UCB-J and UCB-F to different isoforms of SV2. Docking and large-scale molecular dynamics simulations were carried out to unravel various binding patterns, types of interactions, and binding free energies, covering hydrogen bonding and nonspecific hydrophobic interactions, water bridge, π-π, and cation-π interactions. The overall preference for bonding types of UCB-J and UCB-F with particular residues in the protein pockets can be disclosed in detail. A unique interaction fingerprint, namely, hydrogen bonding with additional cation-π interaction with the pyridine moiety of UCB-J, could be established as an explanation for its high selectivity over the SV2 isoform A (SV2A). Other molecular details, primarily referring to the presence of π-π interactions and hydrogen bonding, could also be analyzed as sources of selectivity of the UCB-F tracer for the three isoforms. The simulations provide atomic details to support future development of new selective tracers targeting synaptic vesicle glycoproteins and their associated diseases.

Super-Suppression of Long-Wavelength Phonons in Constricted Nanoporous Geometries.

In a typical semiconductor material, the majority of the heat is carried by long-wavelength, long-mean-free-path phonons. Nanostructuring strategies to reduce thermal conductivity, a promising direction in the field of thermoelectrics, place scattering centers of size and spatial separation comparable to the mean free paths of the dominant phonons to selectively scatter them. The resultant thermal conductivity is in most cases well predicted using Matthiessen's rule. In general, however, long-wavelength phonons are not as effectively scattered as the rest of the phonon spectrum. In this work, using large-scale molecular-dynamics simulations, non-equilibrium Green's function simulations, and Monte Carlo simulations, we show that specific nanoporous geometries that create narrow constrictions in the passage of phonons lead to anticorrelated heat currents in the phonon spectrum. This effect results in super-suppression of long-wavelength phonons due to heat trapping and reductions in the thermal conductivity to values well below those predicted by Matthiessen's rule.

In Bacterial Membranes Lipid II Changes the Stability of Pores Formed by the Antimicrobial Peptide Nisin.

Resistance to available antibiotics poses a growing challenge to modern medicine, as this often disallows infections to be controlled. This problem can only be alleviated by the development of new drugs. Nisin, a natural lantibiotic with broad antimicrobial activity, has shown promise as a potential candidate for combating antibiotic-resistant bacteria. However, nisin is poorly soluble and barely stable at physiological pH, which despite attempts to address these issues through mutant design has restricted its use as an antibacterial drug. Therefore, gaining a deeper understanding of the antimicrobial effectiveness, which relies in part on its ability to form pores, is crucial for finding innovative ways to manage infections caused by resistant bacteria. Using large-scale molecular dynamics simulations, we find that the bacterial membrane-specific lipid II increases the stability of pores formed by nisin and that the interplay of nisin and lipid II reduces the overall integrity of bacterial membranes by changing the local thickness and viscosity.

Evolutionary trends of respiratory syncytial viruses: Insights from large-scale surveillance and molecular dynamics of G glycoprotein

Human respiratory syncytial virus (RSV) is an underlying cause of lower respiratory illnesses in children, elderly and immunocompromised adults. RSV contains multiple structural and non-structural proteins with two major glycoproteins that control the initial phase of infection, fusion glycoprotein and the attachment (G) glycoprotein. G protein attaches to the ciliated cells of airways initiating the infection. The hypervariable G protein plays a vital role in evolution of RSV strains. We employed multiple bioinformatics tools on systematically accessed large-scale data to evaluate mutations, evolutionary history, and phylodynamics of RSV. Mutational analysis of central conserved region (CCR) on G protein-coding sequences between 163 and 189 positions revealed frequent mutations at site 178 in human RSV (hRSV) A while arginine to glutamine substitutions at site 180 positions in hRSV B, remained prevalent from 2009 to 2014. Phylogenetic analysis indicates multiple signature mutations within G protein responsible for diversification of clades. The USA and China have highest number of surveillance records, followed by Kenya. Markov Chain Monte Carlo Bayesian skyline plot revealed that RSV A evolved steadily from 1990 to 2000, and rapidly between 2003 and 2005. Evolution of RSV B continued from 2003 to 2022, with a high evolution stage from 2016 to 2020. Throughout evolution, cysteine residues maintained their strict conserved states while CCR has an entropy value of 0.0039(±0.0005). This study concludes the notion that RSV G glycoprotein is continuously evolving while the CCR region of G protein maintains its conserved state providing an opportunity for CCR-specific monoclonal antibodys (mAbs) and inhibitors as potential candidates for immunoprophylaxis.

Exploring Li-Ion Transport Properties of Li3TiCl6: A Machine Learning Molecular Dynamics Study

We performed large-scale molecular dynamics simulations based on a machine-learning force field (MLFF) to investigate the Li-ion transport mechanism in cation-disordered Li3TiCl6 cathode at six different temperatures, ranging from 25°C to 100°C. In this work, deep neural network method and data generated by ab − initio molecular dynamics (AIMD) simulations were deployed to build a high-fidelity MLFF. Radial distribution functions, Li-ion mean square displacements (MSD), diffusion coefficients, ionic conductivity, activation energy, and crystallographic direction-dependent migration barriers were calculated and compared with corresponding AIMD and experimental data to benchmark the accuracy of the MLFF. From MSD analysis, we captured both the self and distinct parts of Li-ion dynamics. The latter reveals that the Li-ions are involved in anti-correlation motion that was rarely reported for solid-state materials. Similarly, the self and distinct parts of Li-ion dynamics were used to determine Haven’s ratio to describe the Li-ion transport mechanism in Li3TiCl6. Obtained trajectory from molecular dynamics infers that the Li-ion transportation is mainly through interstitial hopping which was confirmed by intra- and inter-layer Li-ion displacement with respect to simulation time. Ionic conductivity (1.06 mS/cm) and activation energy (0.29eV) calculated by our simulation are highly comparable with that of experimental values. Overall, the combination of machine-learning methods and AIMD simulations explains the intricate electrochemical properties of the Li3TiCl6 cathode with remarkably reduced computational time. Thus, our work strongly suggests that the deep neural network-based MLFF could be a promising method for large-scale complex materials.

Molecular Dynamics Simulations of Different Nanoparticles at Substrates.

We report the results of large-scale molecular dynamics simulations of adsorption nanoparticles on solid surfaces. The particles were modeled as stiff aggregates of spherical segments. Three types of particles were studied: rods, rectangles, and triangles built of the same number of segments. We show how the particle shape affects the adsorption, the structure of the surface layer, and the degree of the removal of particles from the solvent. The systems with different segment-segment and segment-surface interactions and different concentrations of particles were investigated. The ordered structures formed in adsorption monolayers were also analyzed. The results are consistent with experimental observations.

GENESIS CGDYN: large-scale coarse-grained MD simulation with dynamic load balancing for heterogeneous biomolecular systems.

Residue-level coarse-grained (CG) molecular dynamics (MD) simulation is widely used to investigate slow biological processes that involve multiple proteins, nucleic acids, and their complexes. Biomolecules in a large simulation system are distributed non-uniformly, limiting computational efficiency with conventional methods. Here, we develop a hierarchical domain decomposition scheme with dynamic load balancing for heterogeneous biomolecular systems to keep computational efficiency even after drastic changes in particle distribution. These schemes are applied to the dynamics of intrinsically disordered protein (IDP) droplets. During the fusion of two droplets, we find that the changes in droplet shape correlate with the mixing of IDP chains. Additionally, we simulate large systems with multiple IDP droplets, achieving simulation sizes comparable to those observed in microscopy. In our MD simulations, we directly observe Ostwald ripening, a phenomenon where small droplets dissolve and their molecules redeposit into larger droplets. These methods have been implemented in CGDYN of the GENESIS software, offering a tool for investigating mesoscopic biological processes using the residue-level CG models.

Turbulent spot formation in three-dimensional yukawa liquids using large-scale molecular dynamics simulation—effect of system size

We have performed classical ‘first principles’ 3D Molecular Dynamics (MD) simulation of plane-Couette flow (PCF) in a 3D Yukawa liquid. The main aim of the present investigation is to study the effect of stream-wise (flow-direction) and span-wise (transverse direction to the flow) width on the subcritical transition to turbulence in plane Couette flow (PCF), separately. In the past, taking very large-aspect ratio systems, subcritical transition to turbulence in PCF has been studied in great detail. However, the effect of stream-wise and span-wise width on the turbulent dynamics separately, have not been studied in detail. In the present work, we have investigated the effect of stream-wise and span-wise width or length and found that the stream-wise length enhances the large-scale energy, which gives rise to strong large-scale flow, and span-wise width enhances the small-scale energy, which gives rise to strong small-scale structures. In other words, topology of the turbulent spot is observed to change with the change in system size. The spectral separation between the modes governing the large and small-scale dynamics improves with the increase in system sizes. The connection between the stream-wise vortices and the stream-wise velocity streaks is investigated in details. We have found that the number stream-wise velocity streaks is one unit larger than the number of stream-wise vortices. A qualitative comparison between the results obtained from the MD simulation and the hydrodynamics experiment is presented in this article. The behavior of the system was observed to vary with the range of interactions among the dust grains.

Effects of shear strain on shock response in single crystal iron

With large-scale non-equilibrium molecular dynamics simulations and in situ x-ray diffraction analysis, we conducted a systematic investigation into the effects of pre-existing shear strain (γxy) on the shock response of single crystal iron. Our findings reveal significant effects of γxy on the deformation of the crystal structure during shock loading, leading to noticeable alterations in the propagation of shock waves. Specifically, during the elastic stage, the presence of γxy results in a reduction of shock strength, consequently diminishing the magnitude of elastic lattice strain (εe). In the plastic stage, γxy stimulates the α–ε phase transformation, and structure deformation undergoes a transition from the sequential activity of dislocation-to-transformation to the synchronous activity of dislocation and transformation. This transition inhibits the propagation of plastic waves and consequently broadens the elastic regime. Additionally, the introduction of γxy activates different slip systems, as it alters the corresponding resolved shear stress. Concurrently, the presence of γxy triggers the activation of different high-pressure phase variants. Our investigation sheds light on the fundamental physics of iron under shock compression and the influence of pre-existing shear strain on its behavior.