Classical Molecular Dynamics Research Articles

A Molecular Dynamics Study of the Influence of Low-Dosage Methanol on Hydrate Formation in Seawater and Pure Water Metastable Solutions of Methane

The behavior of low concentrations of methanol (0.5 and 1.0 wt% of water) as a promoter for hydrate formation in seawater or pure water metastable solutions of methane was investigated using the classical molecular dynamics method at moderate temperature and pressure. The influence of methanol on the dynamics of the re-arrangement of the hydrogen bond network in seawater and pure water solutions of methane was studied by calculating order parameters of the tetrahedral environment and intermolecular torsion angles for water molecules, as well as by calculating the number of hydrogen bonds, hydrate, and hydrate-like cavities. It was found that hydrate nucleation can be considered a collective process in which the rate of hydrate growth is faster in systems with low concentrations of methanol, and confident hydrate growth begins earlier in a metastable solution without sea salt with a small amount of methanol than in systems without methanol.

TH-graphyne: a new porous bidimensional carbon allotrope.

Graphyne and two-dimensional porous carbon-based materials have garnered significant attention due to their interesting structural characteristics and essential properties for new technological applications. Within this scope, this work investigates the structural, thermal, electronic, optical, and mechanical properties of a novel two-dimensional allotrope that combines triangular (T) and hexagonal (H) rings, connected by acetylenic linkages (graphyne-like), thus named TH-graphyne (TH-GY). This study comprehensively characterizes the proposed system's behavior using density functional theory, ab initio molecular dynamics, and classical reactive molecular dynamics simulations. Our results confirm the structural stability of TH-GY. AIMD simulations demonstrate the material's thermal stability at elevated temperatures, while phonon dispersions indicate its dynamical stability. Electronic band structure calculations show that the system is metallic. The analysis of optical properties reveals intense activity in the visible and UV regions, with pronounced anisotropy. A machine learning interatomic potentials model was developed for TH-GY and used to determine the mechanical behavior of the system, which exhibits Young's modulus ranging from 263 to 356 GPa, highlighting its flexibility. Classical reactive MD simulations elucidate the fracture behavior of TH-GY, revealing distinct fracture patterns and mechanical anisotropy.

Effect of Protein-Polarized Ligand Charges on Relative Protein Ligand Binding Affinities.

A major challenge in computer-aided drug design is predicting relative binding energies of different molecules to a target protein using fast and accurate free-energy calculation methods. Free-energy calculations are primarily computed by utilizing classical molecular dynamics simulations based on all-atom force fields (FF) to model the interactions in the system. The present standard classical all-atom FFs contain fixed partial charges on the atoms, and hence electrostatic interactions are modeled between them. The parametrization process to determine these partial charges usually relies on quantum mechanics or semiempirical calculations of the molecule in the gas phase or homogeneous water surrounding. These present standard parametrization schemes of the partial charges neglect, therefore, polarization effects from the protein surrounding. The absence of protein polarization effects can lead to significant errors in free-energy calculations in proteins. We present a parametrization scheme for the partial charges of ligands, named protein-induced polarization (PIP) charges, which account for the electrostatic polarization due to the protein surrounding. The scheme involves single-point quantum mechanics/molecular mechanics calculations of the ligand charges in the protein/water surrounding. Using PIP ligand partial charges, we have calculated the relative binding free energies (RBFEs) of well-studied protein-ligand systems. We show here that RBFEs computed with PIP charges are either significantly improved or at least comparable to those computed with nonpolarized standard GAFF charges. Overall, we present a simple-to-use parametrization scheme to include protein polarization in any type of binding free-energy calculations. The parametrization scheme increases the accuracy in RBFE calculations, while it does not add significant computation time to standard parametrization procedures.

Molecular Investigation of the Self-Assembly Mechanism Underlying Polydopamine Coatings: The Synergistic Effect of Typical Building Blocks Acting on Interfacial Adhesion.

Polydopamine (PDA) is well known as a mussel-inspired adhesive material composed of oligomeric heteropolymers. However, the conventional eumelanin-like structural assumption of PDA seems deficient in explaining its interfacial adhesion. To determine the decisive mechanism of PDA coating formation, experiments and simulations were performed in this study. 5,6-Dihydroxyindole (DHI), the signature building block of eumelanin, was introduced as the control group. Various typical building blocks in PDA were quantified by physicochemical characterizations, and the polar-group-dominated interfacial interaction was evaluated by classic molecular dynamics and metadynamics methods. Aminoethyl has been proven to be the key functional group inducing the adsorption of PDA on the hydroxylated silica substrates, while DHI shows limited adhesion to the substrate due to the absence of aminoethyl as the catechol-indole structure of DHI exhibits poor affinity to the silica surface. Pyrrole carboxylic acid, as an oxidative product detected from PDA/DHI, is unfavorable for its adhesion to silica substrates. Overall, the coating formation and self-aggregating precipitation of PDA are two competitive aminoethyl-consuming paths; thus, the in situ oxidative coupling of dopamine is indispensable for the PDA coating preparation. The collected PDA precipitates can no longer present satisfactory coating forming behavior, resulting from a shortage of aminoethyl moieties.

The structure of water-ammonia mixtures from classical and abinitio molecular dynamics.

The structure of aqueous ammonia solutions is investigated through classical molecular dynamics (MD) and abinitio molecular dynamics (AIMD) simulations. We have preliminarily compared three well-known classical force fields for liquid water (SPC, SPC/E, and TIP4P) in order to identify the most accurate one in reproducing AIMD results obtained at the Generalized Gradient Approximation (GGA) and meta-GGA levels of theory. Liquid ammonia has been simulated by implementing an optimized force field recently developed by Chettiyankandy et al. [Fluid Phase Equilib. 511, 112507 (2020)]. Analysis of the radial distribution functions for different ammonia concentrations reveals that the three water force fields provide comparable estimates of the mixture structure, with the SPC/E performing slightly better. Although a fairly good agreement between MD and AIMD is observed for conditions close to the equimolarity, at lower ammonia concentrations, important discrepancies arise, with classical force fields underestimating the number and strength of H-bonds between water molecules and between water and ammonia moieties. Here, we prove that these drawbacks are rooted in a poor sampling of the configurational space spanned by the hydrogen atoms lying in the H-bonds of H2O⋯H2O and, more critically, H2O⋯NH3 neighbors due to the lack of polarization and charge transfer terms. This way, non-polarizable classical force fields underestimate the proton affinity of the nitrogen atom of ammonia in aqueous solutions, which plays a key role under realistic dilute ammonia conditions. Our results witness the need for developing more suited polarizable models that are able to take into account these effects properly.

Carbon Dioxide Capture on Oxygen- and Nitrogen-Containing Carbon Quantum Dots.

To address global climate change challenges, an effective strategy involves capturing CO2 directly at its source using a sustainable, low-cost adsorbent. Carbon quantum dots (CQDs), derived from lignin, are employed to modify the internal surface of an activated carbon adsorbent, enabling selective adsorption based on electrostatic interactions. By manipulating charge distribution on CQDs through either doping (nitrogen) or functionalization (amine, carboxyl, or hydroxyl groups), the study confirms, through classical molecular dynamics simulations, the potential to adjust binding strength, adsorption capacity, and selectivity for CO2 over N2 and O2. For simulations with a single component gas, maximum selectivities of 3.6 and 6.7 are shown for CO2/N2 and CO2/O2, respectively, at 300 K. Simulations containing a wet flue gas indicate that the presence of water increases the CO2/N2 and CO2/O2 selectivities. The highest CO2/H2O selectivity obtained from a CQD/graphite system is 4.3. A comparison of graphite and lignin-based carbon composite (LBCC) substrates demonstrated that LBCC has enhanced adsorptive capacity. The roughness of the LBCC substrate prevents the diffusion of the CQD on the surface. This computational study takes another step toward identifying optimal CQD atomic architecture, dimensions, doping, and functionalization for a large-scale CQD/AC adsorbent solution for CO2 capture.

Mechanical performance of amorphous diamond-like carbon nanowires

Amorphous carbon nanowires (NW) are a fundamental piece in the study of novel nanoarchitectures with unexpected mechanical properties. However, the failure mechanism of amorphous carbon NW is still missing. In this study, classical molecular dynamics was employed to conduct stress-strain tests to address the mechanical response of amorphous carbons NW with different radii. This research characterizes the mechanical properties of aC NW with different sp3 contents, including Youngs modulus and ultimate tensile stress. Our study reveals that plastic deformation is mediated by chemical modifications in carbon bonding, leading to transitions from sp2 to sp3 and vice versa. We observed that denser amorphous carbon materials undergo a transition from sp3 to sp2 bonding prior to failure, facilitated by the recombination of sp2 clusters. Additionally, plastic deformation in amorphous carbon NW is facilitated by the formation of shear transformation zones and nanovoids, with the deformation mechanism strongly dependent on the average coordination of carbon atoms. These insights provide valuable information for designing nanomaterials with tailored mechanical properties.

DMDA-PatA mediates RNA sequence-selective translation repression by anchoring eIF4A and DDX3 to GNG motifs

Small-molecule compounds that elicit mRNA-selective translation repression have attracted interest due to their potential for expansion of druggable space. However, only a limited number of examples have been reported to date. Here, we show that desmethyl desamino pateamine A (DMDA-PatA) represses translation in an mRNA-selective manner by clamping eIF4A, a DEAD-box RNA-binding protein, onto GNG motifs. By systematically comparing multiple eIF4A inhibitors by ribosome profiling, we found that DMDA-PatA has unique mRNA selectivity for translation repression. Unbiased Bind-n-Seq reveals that DMDA-PatA-targeted eIF4A exhibits a preference for GNG motifs in an ATP-independent manner. This unusual RNA binding sterically hinders scanning by 40S ribosomes. A combination of classical molecular dynamics simulations and quantum chemical calculations, and the subsequent development of an inactive DMDA-PatA derivative reveals that the positive charge of the tertiary amine on the trienyl arm induces G selectivity. Moreover, we identified that DDX3, another DEAD-box protein, is an alternative DMDA-PatA target with the same effects on eIF4A. Our results provide an example of the sequence-selective anchoring of RNA-binding proteins and the mRNA-selective inhibition of protein synthesis by small-molecule compounds.

Improving the reliability of classical molecular dynamics simulations in battery electrolyte design

Explorations into new electrolytes have highlighted the critical impact of solvation structure on battery performance. Classical molecular dynamics (CMD) using semi-empirical force fields has become an essential tool for simulating solvation structures. However, mainstream force fields often lack accuracy in describing strong ion-solvent interactions, causing disparities between CMD simulations and experimental observations. Although some empirical methods have been employed in some of the studies to address this issue, their effectiveness has been limited. Our CMD research, supported by quantum chemical calculations and experimental data, reveals that the solvation structure is influenced not only by the charge model but also by the polarization description. Previous empirical approaches that focused solely on adjusting ion-solvent interaction strengths overlooked the importance of polarization effects. Building on this insight, we propose integrating the Drude polarization model into mainstream force fields and verify its feasibility in carbonate, ether, and nitrile electrolytes. Our experimental results demonstrate that this approach significantly enhances the accuracy of CMD-simulated solvation structures. This work is expected to provide a more reliable CMD method for electrolyte design, shielding researchers from the pitfalls of erroneous simulation outcomes.

Binary phase separation in strongly coupled plasma

We investigated the two-dimensional binary phase separation process of plasma species using classical molecular dynamics in the strongly coupled regime. Both the plasma species interact via a pairwise screened Coulomb (Debye–Hückel) potential; however, the screening parameter κ is different for like- and unlike-species and is the cause for phase separation. We characterize the separation process by measuring the domain growth of equilibrium phases as a function of time—generally, the more significant the inhomogeneity in pairwise interaction, the faster the domain growth. Typically, the domain growth follows a power law in time with an exponent β characterizing the underlying coarsening mechanism. We demonstrate that the growth law exponent is β=1/2 for equal-number-density mixtures and 1/3 otherwise. Further, by comparing these with the corresponding growth laws in binary mixtures of viscous fluids, we show that the viscoelastic nature of plasma fluid modifies the coarsening dynamics, which in turn leads to the observed growth law exponents, notably in the unequal-number-density case.

Water transport through monolayer fullerene membrane

Water transport through nanoporous materials is important in water treatment, desalination, and nanofiltration. Two-dimensional (2D) membranes such as porous graphene have been explored for high-permeance water transport. However, water transport through a new class of 2D membranes based on two-dimensional covalently linked fullerene monolayers has not been fully explored. Here we use classical molecular dynamics simulations to investigate both vapor and liquid water transport through a monolayer fullerene membrane. We find that a quasi-tetragonal phase fullerene membrane possesses the right pore size and geometry that allows fast water vapor transport (∼ 50 g m−2 day−1 Pa−1) and water liquid transport (∼ 2.0 g m−2 day−1 Pa−1). Furthermore, simulation of sea water transport through the fullerene membrane shows 100 % salt rejection. The much faster vapor transport rate is attributed to the funnel-shaped pore and the optimal size that allows free rotation of water molecules permeating through, while the slower liquid transport is due to the need to desolvate a water molecule to break its hydrogen-bond network across the hydrophobic pore. This work shows the great potential of using monolayer fullerene membranes as 2D membranes for fast and selective water transport.

Molecular dynamics simulation of the diffusion behaviour in Ti–Ag using diffusion couple method

The diffusion behaviour of the Ti–Ag system was investigated using classical molecular dynamics by simulating a diffusion couple at various temperatures. One bulk of the diffusion couple consisted of titanium atoms arranged in a hexagonal close-packed (hcp) structure, i.e., α-Ti, and the other consisted of silver atoms arranged in a face-centred cubic (fcc) structure. The modern second nearest-neighbour modified embedded-atom method (2NN-MEAM) interatomic potential was used. The mean square displacement (MSD) of Ti atoms diffusing into Ag and Ag atoms into α-Ti was recorded during the simulation. Diffusion coefficients were determined from the slopes of lines fit to linear regions of MSD dependence on time. The Arrhenius relation was used to determine the activation energy and the frequency of attempts. The activation energy of Ti was 0.61eV and frequency of attempts 1.90⋅10−8m2s−1. The determined activation energy of Ag and frequency of attempts were not considered reliable because the Ag diffusion coefficient did not grow exponentially with temperature. Ti atoms penetrated deeper and greater amounts into Ag than Ag into α-Ti and also had a greater diffusion coefficient at all temperatures. The results were compared with the available experimental data, and for these purposes, self-diffusion in α-Ti and Ag was further investigated. Finally, it was determined that Ag diffuses in α-Ti by hopping over interstitial positions and Ti in Ag via vacancies.

Computational study on the bifurcation mechanism in the H2CO− + CH3Cl →CH3CH2O + Cl−/H2CO + CH3 + Cl− reaction: The importance of intramolecular vibrational redistributions

We have investigated the bifurcation dynamics of the H2CO− + CH3Cl →CH3CH2O + Cl− (C-substitution channel)/H2CO + CH3 + Cl− (electron transfer channel) reaction using three different molecular dynamics (MD) methods: quasi-classical trajectory (QCT), classical MD, and ring-polymer MD simulations. All MD calculations were performed using the on-the-fly approach at the ab initio molecular orbital theory level. We found that the C-substitution channel is dominant in the classical and ring-polymer MD calculations, while the electron transfer channel is dominant in the QCT calculations. The different bifurcation mechanisms in the types of MD simulations were investigated using supervised machine-learning classification. This difference is attributed to the varying efficiency of the intramolecular vibrational redistribution process among the different MD methods. Thus, in the bifurcation mechanism of the H2CO− + CH3Cl reaction, the bifurcation fractions are primarily determined by the vibrational dynamics occurring in the later stages of the reaction.

Bridging classical nucleation theory and molecular dynamics simulation for homogeneous ice nucleation.

Water freezing, initiated by ice nucleation, occurs widely in nature, ranging from cellular to global phenomena. Ice nucleation has been experimentally proven to require the formation of a critical ice nucleus, consistent with classical nucleation theory (CNT). However, the accuracy of CNT quantitative predictions of critical cluster sizes and nucleation rates has never been verified experimentally. In this study, we circumvent this difficulty by using molecular dynamics (MD) simulation. The physical properties of water/ice for CNT predictions, including density, chemical potential difference, and diffusion coefficient, are independently obtained using MD simulation, whereas the calculation of interfacial free energy is based on thermodynamic assumptions of CNT, including capillarity approximation among others. The CNT predictions are compared to the MD evaluations of brute-force simulations and forward flux sampling methods. We find that the CNT and MD predicted critical cluster sizes are consistent, and the CNT predicted nucleation rates are higher than the MD predicted values within three orders of magnitude. We also find that the ice crystallized from supercooled water is stacking-disordered ice with a stacking of cubic and hexagonal ices in four representative types of stacking. The prediction discrepancies in nucleation rate mainly arise from the stacking-disordered ice structure, the asphericity of ice cluster, the uncertainty of ice-water interfacial free energy, and the kinetic attachment rate. Our study establishes a relation between CNT and MD to predict homogeneous ice nucleation.

Structural and Electronic Properties of Novel Azothiophene Dyes: A Multilevel Study Incorporating Explicit Solvation Effects.

Among azobenzene derivatives, azothiophenes represent a relatively recent family of compounds that exhibit similar characteristics as dyes and photoreactive systems. Their technological applications are extensive thanks to the additional design flexibility conferred by the heteroaromatic ring. In this study, we present a comprehensive investigation of the structural and electronic properties of novel dyes derived from 3-thiophenamine, utilizing a multilevel approach. We thoroughly examined the potential energy surfaces of the E and Z isomers for three molecules, each bearing different substituents on the phenyl ring at the para position relative to the diazo group. This exploration was conducted through quantum chemistry calculations at various levels of theory, employing a continuum solvent model. Subsequently, we incorporated an explicit solvent (a dimethyl sulfoxide-water mixture) to simulate the most stable isomers using classical molecular dynamics, delivering a clear picture of the local solvation structure and intermolecular interactions. Finally, a hybrid quantum mechanics/molecular mechanics (QM/MM) approach was employed to accurately describe the evolution of the solutes' properties within their environment, accounting for finite temperature effects.

Diffuse scattering from dynamically compressed single-crystal zirconium following the pressure-induced α→ω phase transition

The prototypical α→ω phase transition in zirconium is an ideal test bed for our understanding of polymorphism under extreme loading conditions. After half a century of study, a consensus had emerged that the transition is realized via one of two distinct displacive mechanisms, depending on the nature of the compression path. However, recent dynamic-compression experiments equipped with diffraction diagnostics performed in the past few years have revealed new transition mechanisms, demonstrating that our understanding of the underlying atomistic dynamics and transition kinetics is in fact far from complete. We present classical molecular dynamics simulations of the α→ω phase transition in single-crystal zirconium shock compressed along the [0001] axis using a machine-learning-class potential. The transition is predicted to proceed primarily via a modified version of the two-stage Usikov-Zilberstein mechanism, whereby the high-pressure ω phase heterogeneously nucleates at boundaries between grains of an intermediate β phase. We further observe the fomentation of atomistic disorder at the junctions between β grains, leading to the formation of highly defective interstitial material between the ω grains. We directly compare synthetic x-ray diffraction patterns generated from our simulations with those obtained using femtosecond diffraction in recent dynamic-compression experiments, and show that the simulations produce the same unique, anisotropic diffuse scattering signal unlike any previously seen from an elemental metal. Our simulations suggest that the diffuse signal arises from a combination of thermal diffuse scattering, nanoparticlelike scattering from residual kinetically stabilized α and β grains, and scattering from interstitial defective structures. Published by the American Physical Society 2024

Origin of Manipulating Selectivity in Prins Cyclization by [Ga4L6]12- Supramolecular Catalysis through Host-Guest Interactions.

The host effect of the supramolecular [Ga4L6]12- tetrahedral metallocage on Prins cyclization reaction of the substrate by encapsulated citronellal has been investigated by means of molecular dynamics and quantum mechanics. The encapsulation process of the substrate into the [Ga4L6]12- cavity was simulated via attach-pull-release (APR) methods. Thermodynamic calculations and classical molecular dynamics simulations assessed the substrate's microenvironment inside the cavity, guiding DFT-level modeling of the reaction. DFT calculations show diol product predominance in acidic solution but high enol selectivity inside [Ga4L6]12-, consistent with experimental findings. [Ga4L6]12- alters the selectivity of the Prins cyclization reaction by inhibiting diol formation. The activation strain model-based decomposition analysis (ASM-DA) of the barrier difference among distortion and interaction terms indicates that the more positive interaction between a host and guest in the diol transition state than enol determines the product selectivity, particularly the fewer C-H···O and O-H···O hydrogen-bonding interactions. These theoretical insights could contribute to a deeper understanding of the nature of supramolecular catalysis and to further develop new supramolecular catalysts.

Development of All-Atom Empirical Potentials for Phloroglucinol (Phg) and Hexachlorocyclotriphosphazene (HCCP) Based on their Crystal Phase Structures.

Two existing generic force fields have been augmented with partial charges and tuned in order to give intercompatible all-atom empirical potentials that can satisfactorily represent the known crystal phase structures of the organic phloroglucinol (Phg) (C6H6O3) and inorganic hexachlorocyclotriphosphazene (HCCP) (N3P3Cl6) molecules at several temperatures. It has been proposed that HCCP-Phg network polymers could act as efficient H2 barrier layers in hydrogen storage tanks for cars. However, essential requirements for modeling such networks are adequate representations of both monomers in their pure dense phases using a common form of force field. Tests of their ability to maintain stable crystal structures have been made using classical molecular dynamics (MD) simulations on large 800-molecule supercells. The force fields have been optimized to match the densities calculated from the experimental unit cell dimensions at ambient conditions as well as the intermolecular potential energy, as estimated from experimental enthalpies of sublimation. For Phg, the crystal structure is stabilized by a network of hydrogen bonds and the Coulombic interactions contribute to over 55% of the total intermolecular potential energy. In contrast, the crystal structure of HCCP is intrinsically stabilized by the van der Waals terms. Both optimized force fields reproduce very well the orthorhombic symmetry of their respective crystals under constant-pressure NPT conditions. The model parameters tuned at ambient temperature also give reasonable agreement with crystallographic data at lower temperatures.

Temperature effects on the branching dynamics in the model ambimodal (6 + 4)/(4 + 2) intramolecular cycloaddition reaction.

C14H20 (tetradecapentaene) is a simple model system exhibiting post transition-state behavior, wherein both the (6 + 4) and (4 + 2) cycloaddition products are formed from one ambimocal transition state structure. We studied the bifurcation dynamics starting from the two ambimodal transition state structures, the chair-form and boat-form, using the quasi-classical trajectory, classical molecular dynamics, and ring-polymer molecular dynamics methods on the parameter-optimized semiempirical GFN2-xTB potential energy surface. It was found that the calculated branching fractions differ between the chair-form and boat-form due to the different nature in the IRC pathways. We also investigated the effects of temperature on bifurcation dynamics and found that, at higher temperatures, trajectories stay longer in the intermediate region of the potential energy surface.

Deciphering the Decomposition Mechanisms of Ether and Fluorinated Ether Electrolytes on Lithium Metal Surfaces: Insights from CMD and AIMD Simulations.

The performance of lithium metal batteries can be significantly enhanced by incorporating fluorinated ether-based electrolytes, yet the solid electrolyte interphase (SEI) formation mechanism on lithium metal surfaces remains elusive. This study employs classical and ab initio molecular dynamics simulations to investigate the decomposition mechanisms of lithium bis(fluoromethanesulfonyl)imide (LiFSI) in 1,2-diethoxyethane (DEE) and its fluorinated analogues, F5DEE and F2DEE, when in contact with lithium metal. Our findings indicate that F5DEE-based electrolytes favor the formation of a FSI-rich primary solvation shell around Li+, while F2DEE-based electrolytes yield a solvent-rich environment. The normalized number density at the Li/electrolyte/Li interface shows a depletion of FSI anions in the electrochemical double layer (EDL) structure near the Li anode upon charging, with the distance between the first main peak of the FSI anion and Li anode following the order F5DEE < DEE < F2DEE. Analysis of the electronic projected density of states and charge transfer dynamics unveils the reductive dissociation pathways of FSI anions and fluorinated DEE solvents on the lithium metal surface, taking into account the influence of the EDL structure. DEE is identified as the most reduction-stable solvent, leading to the selective dissociation of FSI anions and the formation of an entirely inorganic SEI. In contrast, F2DEE displays a pronounced reduction tendency, forming an organic-rich SEI due to the solvent-dominated lowest unoccupied molecular orbital at the interface. F5DEE, competing with FSI anions for reduction, results in the formation of an inorganic-rich hybrid SEI with the highest LiF content. The simulation results correlate well with experimental observations and underscore the pivotal role of various fluorinated functional groups in the formation of EDL and SEI near the lithium metal surface.