Empirical Molecular Dynamics Simulations Research Articles

Alkali Metal Ion Recognition by 18-Crown-6 in Aqueous Solutions: Evidence from Local Structures

The underlying recognition mechanisms of alkali metal ions by crown ethers in aqueous solutions are yet to be fully understood at the molecular level. We report direct experimental and theoretical evidence for the structure and recognition sequence of alkali metal ions (Li+, Na+, K+, Rb+, and Cs+) by 18-crown-6 in aqueous solutions by wide-angle X-ray scattering combined with an empirical potential structure refinement modeling and ab initio molecular dynamics simulation. Li+, Na+, and K+ are located in the negative potential cavity of 18-crown-6, with Li+ and Na+ deviating from the centroid of 18-crown-6 by 0.95 and 0.35 Å, respectively. Rb+ and Cs+ lie outside the 18-crown-6 ring and deviate from the centroid of 18-crown-6 by 0.05 and 1.35 Å, respectively. The formation of the 18-crown-6/alkali metal ion complexes is dominated by electrostatic attraction between the cations and the oxygen atoms (Oc) of 18-crown-6. Li+, Na+, K+, and Rb+ form the H2O···18-crown-6/cation···H2O "sandwich" hydrates, while water molecules only hydrate with Cs+ of the 18-crown-6/Cs+ complex on the same side of Cs+. Based on the local structure, the recognition sequence of 18-crown-6 for alkali metal ions in an aqueous solution follows K+ > Rb+ >Na+ >Li+, which is completely different from that (Li+ > Na+ > K+ > Rb+ > Cs+) in the gas phase, confirming that the solvation medium seriously affects the cation recognition of crown ethers. This work provides atomic insights into understanding the host-guest recognition and solvation behavior of crown ether/cation complexes.

Computed Pre-reactive Complex Association Lifetimes Explain Trends in Experimental Reaction Rates for Peroxy Radical Recombinations

The lifetimes of pre-reactive complexes, although implicitly part of the equations used to model many gas-phase bimolecular reactions, have seldom been included in quantitative calculations of rate coefficients. Here, we demonstrate the application of empirical molecular dynamics simulations of collisions between peroxy radicals to model association lifetimes. With the exception of the methyl peroxy–acetyl peroxy system, measurements of the lifetimes based on a phenomenological model are shown to correlate well with available experimental data for recombination reactions of peroxy radicals in cases where the rate-limiting transition state lies below the reactants in energy. Further, we predict reaction rates for larger α-pinene-derived peroxy radicals, and we interpret our results in tandem with available experimental data on these systems, which are of great relevance to improve our understanding of atmospheric aerosol formation.

Insight into Atomic-Scale Adhesion at the C–Cu Interface During the Initial Stage of Nanoindentation

Adhesion is a common phenomenon in nanomachining which affects processing accuracy and repeatability. As material removal approaches the atomic or close-to-atomic scale, quantum mechanics becomes the dominant principle behind the atomic-level interaction. However, atomic-scale effects cannot be properly described by empirical potential function-based molecular dynamics simulations. This study uses a first-principles method to reveal the atomic-scale adhesion between a diamond tip and a copper slab during initial-stage nanoindentation. Using a simplified tip and slab model, adhesion energy, electronic distribution, and density of states are analyzed based on quantum chemistry calculation. Results show that atomic adhesion is primarily due to the covalent bonding interaction between C and Cu atoms, which can induce structural changes to the diamond tip and copper slab. The effects of tip position and angles on adhesion are further studied through a series of simulations. The results show that adhesion between the tip and slab is sensitive to the lattice structure and a variant in angstroms is enough to cause different adhesion and structural changes. The actual determinants of adhesion can only be the atomic and electronic structures at the tip–slab interface. Bond rotation and breakage are observed during simulation and their effects on adhesion are further discussed. To conclude, the first-principles method is important for the analysis of an atomic-scale interaction system, even if only as an aid to describing adhesion at atomic and electronic scales.

Molecular dynamics and density functional theory study on the potassium distribution and lattice thermal conductivity of KxRhO2

We present a theoretical study on the K ions' diffusion and lattice thermal conductivities in the layered KxRhO2 by the empirical molecular dynamics simulations and first-principles calculations. The calculated lattice thermal conductivities of the non-stoichiometric KxRhO2 are much smaller than that of the stoichiometric one, which implies that the disordered distribution of K ions has significantly suppressed the thermal conductivity. This study could inspire more theoretical works about the concept of ‘phonon liquid’ in other similar thermoelectric materials.

An internal energy-dependent model for the Grüneisen parameter of silicate liquids

We investigated the accuracy of the Mie-Grüneisen approximation, which treats the Grüneisen parameter (γ) as a one-parameter function of volume, for use in describing the thermal equation of state of a silicate liquid. For this study, we focused on a single composition: the diopside-anorthite eutectic, an Fe-free basalt analog that has been extensively studied by shock wave experiments. We tuned an empirical force-field to a small set of ab initio N-V-T molecular dynamics simulations to ensure that it reproduces pressure, heat capacity, and γ at high and low pressures. We then used empirical force-field molecular dynamics simulations in a larger system and for longer run times to ensure accurate extraction of γ at numerous N-V-E state points. To first order, the results show the expected volume-dependence for silicate liquids, with γ increasing as volume decreases. However, there are also significant and systematic variations of γ with internal energy (E) at constant volume. We propose a simple model form that captures the volume and E dependence of γ with only one more free parameter than a typical Mie-Grüneisen formulation. We demonstrate the utility of this new model for well-constrained fitting to sparse shock wave experiment data below 200 GPa, obtaining a marked improvement in the ability to simultaneously fit the pre-heated liquid Hugoniot and points substantially offset from this Hugoniot.

(Keynote) An Automatic Rapid Screening Tool for Fast Ionic Conductors

Identifying new materials that combine high ionic conductivity with structural and electrochemical stability so far remains a slow trial and error search process. To rationally accelerate materials design and exploit the opportunities in the materials genome a dependable rapid screening of materials is required that can pre-select structures that merit higher level computational as well as experimental characterization. Here we report on the progress of our “softBV” bond-valence site energy-based automated pathway analysis utilizing our new Bond Valence Pathway Analyzer (BVPA) – with fast bond valence site energy calculations to quickly obtain suggested candidate materials for fast ionic conduction. BVPA provides rapid and simplified visualization, in order to bridge the gap between experimentalists and simulation [1,2]. Calculation of ion transport pathways can be done extremely quickly on the order of seconds or minutes on desktop PCs providing a speedup factor of 3 to 5 orders of magnitude compared to DFT-based NEB methods. Combined with a graphical user interface our software suite (that can be downloaded from [2] and is free for academic use) should enable experimentalists to quickly identify candidate solid electrolyte materials. We also aim to integrate the pre-screening into an automated workflow for subsequent DFT characterization [3].Results will be benchmarked against both experimental and DFT NEB migration barriers. Besides the migration barriers the approach now also comprises an AI-based dopant predictor utilizing bond-valence-based crystal chemical descriptors to assist experimentalists in exploring favorable substitutional doping strategies. We will also compare the predictability of absolute room temperature conductivities from static energy landscape analysis, bond-valence based empirical MD simulations and ab initio molecular dynamics (AIMD) simulations. While for small fast-ion conductor structures at sufficiently high temperatures AIMD appears to be the gold standard, the less reliable but computationally empirical approaches have an advantage in modelling complex disordered interfaces at low temperatures over longer periods. This eliminates the hazards involved in extrapolations down to room temperature properties for the frequent cases of order-disorder phase transitions at intermediate temperatures.As an example we will discuss lithium and sodium compounds containing multiple anions, in particular the combination of thiophosphate and halide anions or various MS4 polyanions. Based on computational screening using our bond valence site approach and DFT studies several thiophosphate halides along the A3PS4-LiX (Cl, Br, I; A = Li, Na) tie line [4] and the Ax(MS4)y(M’S4)z phase space [5] have been explored and their properties discussed based on BVSE pathway models and molecular dynamics simulations in combination with experimental (X-ray and neutron) diffraction, solid state NMR and electrochemical characterisation. MD simulations e.g. show that Li5(PS4)Cl2 is found to undergo an order-disorder phase transition and thus should, contrasting to earlier predictions, not be a fast Li+ ion conductor. The newly predicted thermodynamically stable cubic solid electrolyte Li15(PS4)4Cl3 was successfully prepared and characterized. Though its conductivity does not reach a superionic level, it demonstrates that the computational approach can successfully predict a completely new classes of solid electrolytes and can predict its optimization by doping. The simplicity of the approach also facilitates the study of homogeneity ranges as exemplified for the solid solution systems Li4-xPS4Ix (0<x<0.67) and Na9+x(MS4)3-x(SnS4)x (M = P, Sb; x ≈ 2).

Solvation of Na? in the Sodide Solution, LiNa?10MeNH2.

Alkalides, the alkali metals in their −1 oxidation state, represent some of the largest and most polarizable atomic species in condensed phases. This study determines the solvation environment around the sodide anion, Na–, in a system of co-solvated Li+. We present isotopically varied total neutron scattering experiments alongside empirical potential structure refinement and ab initio molecular dynamics simulations for the alkali–alkalide system, LiNa·10MeNH2. Both local coordination modes and the intermediate range liquid structure are determined, which demonstrate that distinct structural correlations between cation and anion in the liquid phase extend beyond 8.6 Å. Indeed, the local solvation around Na– is surprisingly well defined with strong solvent orientational order, in contrast to the classical description of alkalide anions not interacting with their environment. The ion-paired Li(MeNH2)4+·Na– species appears to be the dominant alkali–alkalide environment in these liquids, whereby Li+ and Na– share a MeNH2 molecule through the amine group in their primary solvation spheres.

On the hydration of DOPE in solution.

The atomic-scale hydration structure around the 1,2-dioleoyl-sn-glycero-3-phosphoethanolamine (DOPE) headgroup in a chloroform/water solution has been investigated using neutron diffraction enhanced by isotopic substitution and NMR, coupled with empirical potential structure refinement and molecular dynamics simulations. The results obtained show the preferential binding sites for water molecules on the DOPE headgroups, with the most predominant interactions being with the ammonium and phosphate groups. Interestingly, the level of hydration, as well as the association of DOPE molecules, varies according to the simulation method used. The results here suggest the presence of a tight water network around these lipid headgroups that could affect the permeability of the membrane for lipid-mediated diffusion.

Determining influential descriptors for polymer chain conformation based on empirical force-fields and molecular dynamics simulations

Many properties of bulk polymers are closely related to the molecular chain conformation. It is ideal if guidance exists to link molecular features to polymer chain conformation to significantly narrow down the design space. Here, we analyze a series of molecular descriptors derived from empirical force-fields and relate them to single chain conformation characterized by radius of gyration calculated from molecular dynamics simulations. Using data analyses and machine learning techniques, we identify that the weakest dihedral angle, characterized by the energy constants, along the chain backbone is the most influential descriptor determining the single chain radius of gyration.

Ab initio and classical atomistic modelling of structure and defects in crystalline orthorhombic polyethylene: Twin boundaries, slip interfaces, and nature of barriers

We study the stability of twin boundaries and slip in crystalline orthorhombic polyethylene by means of density functional theory (DFT), using a nonempirical, truly nonlocal density function, and by means of classical molecular dynamics (MD). The results show that, in accordance with experimental observations, there is a clear preference to chain slip over transverse slip for all considered slip planes. The activation energy for pure chain slip lies in the range 10–20 mJ/m2 while that for transverse slip corresponds to 40–280 mJ/m2. For the (11¯0)-slip plane the energy landscape is non-convex with multiple potential energy minima, indicating the presence of stable stacking faults. This suggests that dissociation of perfect dislocations into partials may occur. For the two low-energy twin boundaries considered in this work, {110} and {310}, we find that the former is more stable than the latter, with ground state energies corresponding to 8.9 and 28 mJ/m2, respectively. We have also evaluated how well the empirical MD simulations with the all-atom optimized potential for liquid MD simulations (OPLS-AA) and the coarse-grained united atom (UA) potential concur with the DFT results. It is found that an all-atom potential is necessary to partially capture the γ-surface energy landscapes obtained from the DFT calculations. The OPLS-AA predicts chain slip activation energies comparable with DFT data, while the transverse slip energy thresholds are low in comparison, which is attributed to weak close ranged monomer repulsion. Finally, we find that the H-H interaction dominates the slip activation. While not explicitly represented in the UA potential, its key role is revealed by correlating the DFT energy landscape with changes in the electron distributions and by MD simulations in which components of the OPLS-AA intermolecular potential are selectively silenced.

Diffusion and aggregation process of oxygen embedded around an amorphous/crystal interface of Si(001) studied by molecular dynamics simulation

I performed empirical molecular dynamics (MD) simulations to understand the peculiar migration behavior of oxygen embedded in an amorphous Si (a-Si) layer near the crystal/amorphous (c/a) Si interface and investigated the time evolution of the atomic configuration at high temperatures from 1200 to 1500 K. The previously proposed sweeping effect, which is demonstrated in terms of the oxygen migration and precipitation in silicon taking place along the moving c/a interface, was definitely confirmed in this MD simulation. [Hoshino et al., J. Phys. D: Appl. Phys. 49, 315106 (2016)] In the present study, I reproducibly found the theoretical evidence of the novel sweeping and aggregation phenomenon of oxygen occurring in the recrystallization process of a-Si. The temperature-dependence revealed that the relationship between the displacement velocity of the oxygen and the c/a interface plays an important role in interpreting the behavior. The oxide precipitations in the recrystallized Si as well as the sweeping effect were well reproduced in the simulation in which the systems containing several oxygen atoms were assumed. These facts significantly well explain and support my interpretation in the previous papers reported on the synthesis mechanism of the ultrathin silicon-on-insulator/buried oxide structure prepared by low-energy implantation followed by relatively low temperature annealing, compared to the ordinary separation by the implanted oxygen process.

Highly Dynamic Anion-Quadrupole Networks in Proteins.

The dynamics of anion-quadrupole (or anion-π) interactions formed between negatively charged (Asp/Glu) and aromatic (Phe) side chains are for the first time computationally characterized in RmlC (Protein Data Bank entry 1EP0 ), a homodimeric epimerase. Empirical force field-based molecular dynamics simulations predict anion-quadrupole pairs and triplets (anion-anion-π and anion-π-π) are formed by the protein during the simulated trajectory, which suggests that the anion-quadrupole interactions may provide a significant contribution to the overall stability of the protein, with an average of -1.6 kcal/mol per pair. Some anion-π interactions are predicted to form during the trajectory, extending the number of anion-quadrupole interactions beyond those predicted from crystal structure analysis. At the same time, some anion-π pairs observed in the crystal structure exhibit marginal stability. Overall, most anion-π interactions alternate between an "on" state, with significantly stabilizing energies, and an "off" state, with marginal or null stabilizing energies. The way proteins possibly compensate for transient loss of anion-quadrupole interactions is characterized in the RmlC aspartate 84-phenylalanine 112 anion-quadrupole pair observed in the crystal structure. A double-mutant cycle analysis of the thermal stability suggests a possible loss of anion-π interactions compensated by variations of hydration of the residues and formation of compensating electrostatic interactions. These results suggest that near-planar anion-quadrupole pairs can exist, sometimes transiently, which may play a role in maintaining the structural stability and function of the protein, in an otherwise very dynamic interplay of a nonbonded interaction network as well as solvent effects.

Investigation on gallium ions impacting monolayer graphene

In this paper, the physical phenomena of gallium (Ga+) ion impacting monolayer graphene in the nanosculpting process are investigated experimentally, and the mechanisms are explained by using Monte Carlo (MC) and molecular dynamics (MD) simulations. Firstly, the MC method is employed to clarify the phenomena happened to the monolayer graphene target under Ga+ ion irradiation. It is found that substrate has strong influence on the damage mode of graphene. The mean sputtering yield of graphene under 30 keV Ga+ ion irradiation is 1.77 and the least ion dose to completely remove carbon atoms in graphene is 21.6 ion/nm2. Afterwards, the focused ion beam over 21.6 ion/nm2 is used for the irradiation on a monolayer graphene supported by SiO2 experimentally, resulting in the nanostructures, i.e., nanodot and nanowire array on the graphene. The performances of the nanostructures are characterized by atomic force microscopy and Raman spectrum. A plasma plume shielding model is put forward to explain the nanosculpting results of graphene under different irradiation parameters. In addition, two damage mechanisms are found existing in the fabrication process of the nanostructures by using empirical MD simulations. The results can help us open the possibilities for better control of nanocarbon devices.

Thermal transport across graphene/SiC interface: effects of atomic bond and crystallinity of substrate

The effect of interatomic interaction between graphene and 4H-SiC on their interfacial thermal transport is investigated by empirical molecular dynamics simulation. Two magnitudes of interfacial thermal conductance (ITC) improvement are observed for graphene/4H-SiC interface interacting through covalent bonds than through van der Waals interaction, which can be explained by the bond strength and the number of covalent bonds. Besides, it is found that the ITC of covalent graphene/C-terminated SiC is larger than that Si-terminated SiC, which is due to the stronger bond strength of C–C than that of C–Si. The effect of crystallinity of the substrate is studied, and the result shows that the ITC of graphene/a-SiC is higher than that of graphene/c-SiC. These results are crucial to the understanding of thermal transport across graphene interfaces, which are useful for thermal design in graphene-based transistors.

Dynamical processes of low-energy carbon ion collision with the graphene supported by diamond

The dynamical processes of a low-energy carbon ion collision with the graphene sheet supported by diamond at three impact positions are studied by using empirical potential molecular dynamics simulations. The energy transformation and the structural evolution have been studied. Five types of processes are observed: adsorption, hybridization, defects formation in diamond, atom emission and transmission. We find that the irradiation damage is closely related to the incident energy and impact position. In our simulations, as the projectile collides at a graphene atom, it transfers most of its energy to the primary knock-on atom, and defects are created in graphene. When the projectile moves perpendicular towards the center of a C-C bond in the graphene sheet, the energy transferred from the projectile to the atoms associated with the bond increases firstly and then decreases with the increasing incident energy, and the graphene sheet remains two-dimensional crystal structure after collision when the incident energy is larger than 360 eV. While the impact location is the center of a hexagonal ring on the graphene sheet, the energy transferred from the projectile to the atoms of the target ring is very small regardless of how large is the incident energy, and the graphene sheet is able to keep perfect crystal structure when the incident energy is larger than 34 eV.

Network analysis of proton transfer in liquid water.

Proton transfer in macromolecular systems is a fascinating yet elusive process. In the last ten years, molecular simulations have shown to be a useful tool to unveil the atomistic mechanism. Notwithstanding, the large number of degrees of freedom involved make the accurate description of the process very hard even for the case of proton diffusion in bulk water. Here, multi-state empirical valence bond molecular dynamics simulations in conjunction with complex network analysis are applied to study proton transfer in liquid water. Making use of a transition network formalism, this approach takes into account the time evolution of several coordinates simultaneously. Our results provide evidence for a strong dependence of proton transfer on the length of the hydrogen bond solvating the Zundel complex, with proton transfer enhancement as shorter bonds are formed at the acceptor site. We identify six major states (nodes) on the network leading from the "special pair" to a more symmetric Zundel complex required for transferring the proton. Moreover, the second solvation shell specifically rearranges to promote the transfer, reiterating the idea that solvation beyond the first shell of the Zundel complex plays a crucial role in the process.

New Phonon Dynamics in a Diamond Crystal during Post-Annealing

A new mesoscopic response in a crystal during post-annealing (PA) after ion implantation has been discovered. PA is an indispensable procedure for activating the implanted dopant center by retrieving the lost crystallinity of the host crystal, which usually elapses over periods longer than a few milliseconds. By using an empirical molecular-dynamics simulation, we examined what actually happened during extended PA for up to one nanosecond in case of sub-keV N2 beam implantation into bulk diamond at room temperature with annealing at 1000 K. While observing the changes in the long-range-order (LRO) parameter, a notable phase transition from the crystalline to the amorphous state (CA transition) occurred synergistically at approximately a few picoseconds (ps). After the CA transition, PA was initiated. A few tens of ps later, a sequence of restoration spikes was observed in LRO parameter with a regular time interval of 70–80 ps. Each sharp peak showed the spontaneous restoration of “amorphous -> crystal -> amorphous” transitions. This process may trigger a “critical state” that occurs during the slow and steady recovery of the lost crystallinity in the early stages of PA. It was found that annealing not only enhances the vibration amplitude of the individual atoms in the bulk but also causes phonon-assisted restoration in a mesoscopic space.

Ab Initio Based conformational study of the crystalline α‐chitin

The equilibrium structure including the network of hydrogen bonds of an α-chitin crystal is determined combining density-functional theory (DFT), self-consistent DFT-based tight-binding (SCC-DFTB), and empirical forcefield molecular dynamics (MD) simulations. Based on the equilibrium geometry several possible crystal conformations (local energy minima) have been identified and related to hydrogen bond patterns. Our results provide new insight and allow to resolve the contradicting α-chitin structural models proposed by various experiments.

Transport properties and equation of state of 1-bar eutectic melt in the system CaAl2Si2O8-CaMgSi2O6 by molecular dynamics simulation

Empirical potential molecular dynamics (EPMD) simulations of 1-bar eutectic composition liquid in the system CaAl2Si2O8-CaMgSi2O6 have been conducted using the interatomic pair-potential of Matsui (1998). Simulations using ~10000 atoms over a wide range of conditions (ρ: 2200–5000 kg/m3; T : 1600–5500 K; P : 0–170 GPa) were used to derive an equation of state, determine self-diffusivities for all atoms, calculate melt viscosity, and investigate melt structures by coordination statistics. EOS results compare well to laboratory shock wave data up to ~25 GPa, diverging at higher pressure. Based on simulations of the end-member compositions of the join using the same potential, non-ideality in the volume of mixing at pressures below 10 GPa disappears at higher pressures. Ideal volume mixing at elevated pressure is consistent with inferences from laboratory shock wave studies of liquids in this system. The non-ideal volume of mixing at low pressure is directly correlated to structural differences between the end-member liquids and the mixing of cation-anion coordination polyhedra of differing volume. Self-diffusivities show reasonable agreement with laboratory values, with activation energies and activation volumes in the range 90–100 kJ/mol and 1–3 cm3/mol, respectively. Shear viscosities at 3500 K span from 1.8 × 10−3 Pa·s at low P to ~4.4 × 10−3 Pa·s at ~14 GPa.

Melt-Growth Dynamics in CdTe Crystals

We use a new, quantum-mechanics-based bond-order potential (BOP) to reveal melt growth dynamics and fine scale defect formation mechanisms in CdTe crystals. Previous molecular dynamics simulations of semiconductors have shown qualitatively incorrect behavior due to the lack of an interatomic potential capable of predicting both crystalline growth and property trends of many transitional structures encountered during the melt→crystal transformation. Here, we demonstrate successful molecular dynamics simulations of melt growth in CdTe using a BOP that significantly improves over other potentials on property trends of different phases. Our simulations result in a detailed understanding of defect formation during the melt growth process. Equally important, we show that the new BOP enables defect formation mechanisms to be studied at a scale level comparable to empirical molecular dynamics simulation methods with a fidelity level approaching quantum-mechanical methods.