Microcanonical Molecular Dynamics Simulations Research Articles

Unification of Ewald and shifted force methods to calculate Coulomb interactions in molecular simulations.

Three new Ewald series are derived using a new strategy that does not start with a proposed charge spreading function. Of these, the Ewald series produced using shifted potential interactions for the Ewald real space series converges relatively slowly, while the corresponding expression using a shifted force (SF) interaction does not converge. A comparison is made between several approximations of the Ewald method and the SF route to include Coulomb interactions in molecular dynamics (MD) computer simulations. MD simulations of a model bulk molten salt and water were carried out. The recently derived α' variant of Ewald, by K. D. Hammonds and D. M. Heyes [J. Chem. Phys. 157, 074108 (2022)], has been developed analytically and found to be more accurate and computationally efficient than SF in part due to the smaller real space truncation distance that can be used. In addition, with α', the number of reciprocal lattice vectors required is reduced considerably compared with the standard Ewald implementations to give the same accuracy. The invention of the α' method shifts the computational balance back toward using an Ewald construction. The SF method shows greater errors in the Coulomb pressure and time dependent fluctuation properties compared to α'. It does not conserve the shadow Hamiltonian in a microcanonical MD simulation, whereas the α' method does, which facilitates long time stability and insignificant drift of properties over time. The speed of the Ewald computer code is improved by using a new lookup table method.

Mechanisms of glycine formation in cold interstellar media: a theoretical study.

The possibility of the formation of glycine (Gly) from fundamental gas molecules in cold interstellar media was studied using quantum chemical methods, transition state theory and microcanonical molecular dynamics simulations with surface hopping dynamics (NVE-MDSH). This theoretical study emphasized five photochemical pathways in the lowest singlet-excited (S 1) state, thermochemical processes after non-radiative S 1→S 0 relaxations, and photo-to-thermal energy conversion in the NVE ensemble. The optimized reaction pathways suggested that to generate a reactive singlet dihydroxy carbene (HOCOH) intermediate, photochemical pathways involving the H2O…CO van der Waals and H2O-OC hydrogen bond precursors (Ch (1)_Step (1)) possess considerably lower energy barriers than the S 0 state pathways. The Gibbs free energy barriers (∆G ǂ ) calculated after the non-radiative S 1 →S 0 relaxations indicated higher spontaneous temperatures (T s) for the formation of the HOCOH intermediate (Ch (1)_Step (1)) than for Gly formation (Ch (1)_Step (2) and Ch (4)). Although the termolecular reaction in Ch (4) possesses a low energy barrier, and is thermodynamically favourable, the high exothermic S 1 →S 0 relaxation energy leads to the separation of the weakly associated H2O…CH2NH…CO complex into single molecules. The NVE-MDSH results also confirmed that the molecular processes after the S 1 →S 0 relaxations are thermally selective, and because the non-radiative S 1 →S 0 relaxation temperatures are exceedingly higher than T s, the formation of Gly on consecutive reaction pathways is non-synergistic with low yields and several side products. Based on the theoretical results, photo-to-thermal control strategies to promote desirable photochemical products are proposed. They could be used as guidelines for future theoretical and experimental research on photochemical reactions.

Photoluminescence mechanisms of BF2-formazanate dye sensitizers: a theoretical study.

Photodynamic therapy (PDT) is an alternative, minimally invasive treatment for human diseases such as cancer. PDT uses a photosensitizer to transfer photon energy directly to cellular 3O2 to generate 1O2 (Type II), the toxicity of which leads to cancer cell death. In this work, the photoluminescence mechanisms of a BF2-formazanate dye sensitizer (BF2-FORM) and its iodinated derivative (BF2-FORM-D) were studied using complementary theoretical approaches; the photoluminescence pathways in the S1 and T1 states were studied using density functional theory (DFT) and time-dependent (TD)-DFT methods, the kinetic and thermodynamic properties of the pathways using the transition state theory (TST), and the time evolution and dynamics of key processes using non-adiabatic microcanonical molecular dynamics simulations with surface-hopping dynamics (NVE-MDSH). Evaluation of the potential energy surfaces (PESs) in terms of the rotations of the phenyl rings suggested a pathway for the S1 → S0 transition for the perpendicular structure, whereas two pathways were anticipated for the T1 → S0 transition, namely, [T1 → S0]1 occurring immediately after the S1/T1 intersystem crossing (ISC) and [T1 → S0]2 occurring after the S1/T1 ISC and T1 equilibrium structure relaxation, with the T1 → S0 energy gap being comparable to the energy required for 3O2 → 1O2. The PESs also showed that because of the heavy-atom effect, BF2-FORM-D possessed a significantly smaller S1/T1 energy gap than BF2-FORM. The TST results revealed that at room temperature, BF2-FORM-D was thermodynamically more favorable than the parent molecule. Analysis of the NVE-MDSH results suggested that the librational motions of the phenyl rings play an important role in the internal conversion (IC) and ISC, and the S1/T1 ISC and T1 → S0 transitions could be enhanced by varying the irradiation wavelength and controlling the temperature. These findings can be used as guidelines to improve and/or design photosensitizers for PDT.

Long-Time Dynamics of Selected Molecular-Motor Components Using a Physics-Based Coarse-Grained Approach.

Molecular motors are essential for the movement and transportation of macromolecules in living organisms. Among them, rotatory motors are particularly efficient. In this study, we investigated the long-term dynamics of the designed left-handed alpha/alpha toroid (PDB: 4YY2), the RBM2 flagellum protein ring from Salmonella (PDB: 6SD5), and the V-type Na+-ATPase rotor in Enterococcus hirae (PDB: 2BL2) using microcanonical and canonical molecular dynamics simulations with the coarse-grained UNRES force field, including a lipid-membrane model, on a millisecond laboratory time scale. Our results demonstrate that rotational motion can occur with zero total angular momentum in the microcanonical regime and that thermal motions can be converted into net rotation in the canonical regime, as previously observed in simulations of smaller cyclic molecules. For 6SD5 and 2BL2, net rotation (with a ratcheting pattern) occurring only about the pivot of the respective system was observed in canonical simulations. The extent and direction of the rotation depended on the initial conditions. This result suggests that rotatory molecular motors can convert thermal oscillations into net rotational motion. The energy from ATP hydrolysis is required probably to set the direction and extent of rotation. Our findings highlight the importance of molecular-motor structures in facilitating movement and transportation within living organisms.

Dynamics and mechanisms of nonradiative photodissociation of glycine: A theoretical study

Photochemistry of small biomolecules is fundamental to understand the photochemical reactions of larger molecules, such as proteins. In this work, the dynamics and mechanisms of photodissociations of glycine conformers were studied using the DFT, TD-DFT and CASPT2 methods with the aug-cc-pVDZ basis set and microcanonical molecular dynamics simulations with surface hopping dynamics (NVE-MDSH). This study emphasized unimolecular dissociations and isomerization, nonradiative S1 → S0 relaxations that generate molecular products and the interplay between photoenergies and thermal energies. Four nonplanar conformers were suggested by the TD-DFT/B3LYP/aug-cc-pVDZ and CASPT2(8,14)/aug-cc-pVDZ methods to be stable in the S1 state. The potential energy curves obtained from the TD-DFT/B3LYP/aug-cc-pVDZ and CASPT2(8,14)/aug-cc-pVDZ methods and SA-CPMCSCF/6-31G(d) conical intersection optimizations suggested S0/S1 intersections for distorted glycine structures, Cα–N and NH dissociations and NH → OC isomerization. The NVE-MDSH results showed that although no dissociated molecular product was observed in the S1 state, isomerization-mediated dissociation and exothermic S1 → S0 relaxation temperature could lead to eight molecular product channels, in which the vertically excited Wigner-sampled structures generated eleven reactive and stable molecular products in different temperature ranges. These thermally selective reactions favored the lowest energy conformer (structure Ip) as the precursor, from which reactive and stable molecular products were formed above the threshold exothermic relaxation temperature. These theoretical results showed in detail for the first time how unimolecular photodissociations and nonradiative S1 → S0 relaxations create molecules in cold interstellar media without heat transfer from the surroundings. The results in this work could be used in future theoretical and experimental studies on the same or similar systems.

Strongly coupled Yukawa trilayer liquid: Structure and dynamics.

The equilibrium structure and the dispersion relations of collective excitations in trilayer Yukawa systems in the strongly coupled liquid regime are examined. The equilibrium correlations reveal a variety of structures in the liquid phase, reminiscent of the corresponding structures in the solid phase. At small layer separation substitutional disorder becomes the governing feature. Theoretical dispersion relations are obtained by applying the quasilocalized charge approximation (QLCA) formalism, while numerical data are generated by microcanonical molecular dynamics simulations. The dispersions and polarizations of the collective excitations obtained through both of these methods are compared and discussed in detail. We find that the QLCA method is, in general, very satisfactory, but that there are phenomena not covered by the QLCA. In particular, by analyzing the dynamical longitudinal and transverse current fluctuation spectra we discover the existence of a structure not related to the collective mode spectra. This also provides insight into the long-standing problem of the gap frequency discrepancy, observed in strongly coupled layered systems in earlier studies.

Extension of the UNRES Coarse-Grained Force Field to Membrane Proteins in the Lipid Bilayer.

The physics-based UNRES coarse-grained force field for the simulations of protein structure and dynamics has been extended to treat membrane proteins. The lipid bilayer has been modeled by introducing a continuous nonpolar phase with the water-interface region of appropriate thickness. The potentials for average electrostatic and correlation interactions of the peptide groups have been rescaled to account for the reduction of the dielectric permittivity compared to the water phase and new potentials for protein side-chain-side-chain interactions inside and across the lipid phase have been introduced. The model was implemented in the UNRES package for coarse-grained simulations of proteins, and the package with the new functionality was tested for total energy conservation and thermostat behavior in microcanonical and canonical molecular dynamics simulations runs, respectively. The method was validated by running unrestricted ab initio blind-prediction tests of 10 short α-helical membrane proteins, all runs started from the extended structures. The modified UNRES force field was able to predict correctly the overall folds of the membrane proteins studied.

Chemical dynamics characteristics of Kapton polyimide damaged by electron beam irradiation

Polyimides, particularly poly-(pyromellitic dianhydride-4,40-oxydianiline) or PMDA-ODA (Kapton-H®), are ubiquitous in the aircraft and aerospace industries, flexible printed circuits, and many other high-performance applications. For many of these applications, the polymer's stability during and after radiation damage is a serious concern. The effect of electron irradiation on the breaking of inter-molecular bonds and creation of free radicals in Kapton polyimide were studied using reactive molecular dynamics simulation with the polarizable ReaxFF method and experimentally on pristine and electron irradiated polyimide by vibrational spectroscopy, electron paramagnetic resonance, nuclear magnetic resonance, wide angle x-ray diffraction, small angle x-ray scattering, and dynamic mechanical analysis. In our simulations, the electron beam is approximated by a dense column of negatively charged particles, which is inserted at random positions in the polyimide matrix. In order to observe inter-molecular bond breaking during the time frame of the molecular simulations, each electron beam is switched on for 2 fs, causing an average energy of 21 eV per electron trajectory to be transferred to the polyimide through Coulomb interactions during the simulation. Experimentally, polyimide is aged in vacuum with a 90 keV electron beam with a maximum penetration depth greater than the material's thickness.Throughout the simulations, modifications in the chemical structure of the polyimide during and after the electron beam irradiations are monitored. During electron beam irradiation, hydrogen atoms separate from different parts of the polyimide through C–H bond cleavage. After turning off the electron beams and running a microcanonical molecular dynamics simulation for 0.5 ns, several chemical and structural rearrangements are observed in the polyimide structure. The chemical modifications mostly include ring opening and other changes to the imide rings. The computations are compared to infrared spectroscopy, electron paramagnetic resonance, nuclear magnetic resonance, x-ray scattering, and mechanical measurements on aged and pristine polyimide to both guide and validate the conclusions from the computational study. Apart from these changes in chemistry we also evaluate temperature evolution and change in mechanical response of the polyimide slab caused by electron beam exposure. Taken together, the direct measurements of the electronic, chemical, and mechanical properties of the irradiated Kapton films validate the key results of the model.

Bulk viscosity of molecular fluids.

The bulk viscosity of molecular models of gases and liquids is determined by molecular simulations as a combination of a dilute gas contribution, arising due to the relaxation of internal degrees of freedom, and a configurational contribution, due to the presence of intermolecular interactions. The dilute gas contribution is evaluated using experimental data for the relaxation times of vibrational and rotational degrees of freedom. The configurational part is calculated using Green-Kubo relations for the fluctuations of the pressure tensor obtained from equilibrium microcanonical molecular dynamics simulations. As a benchmark, the Lennard-Jones fluid is studied. Both atomistic and coarse-grained force fields for water, CO2, and n-decane are considered and tested for their accuracy, and where possible, compared to experimental data. The dilute gas contribution to the bulk viscosity is seen to be significant only in the cases when intramolecular relaxation times are in the μs range, and for low vibrational wave numbers (<1000 cm-1); This explains the abnormally high values of bulk viscosity reported for CO2. In all other cases studied, the dilute gas contribution is negligible and the configurational contribution dominates the overall behavior. In particular, the configurational term is responsible for the enhancement of the bulk viscosity near the critical point.

Extension of coarse-grained UNRES force field to treat carbon nanotubes

Carbon nanotubes (CNTs) have recently received considerable attention because of their possible applications in various branches of nanotechnology. For their cogent application, knowledge of their interactions with biological macromolecules, especially proteins, is essential and computer simulations are very useful for such studies. Classical all-atom force fields limit simulation time scale and size of the systems significantly. Therefore, in this work, we implemented CNTs into the coarse-grained UNited RESidue (UNRES) force field. A CNT is represented as a rigid infinite-length cylinder which interacts with a protein through the Kihara potential. Energy conservation in microcanonical coarse-grained molecular dynamics simulations and temperature conservation in canonical simulations with UNRES containing the CNT component have been verified. Subsequently, studies of three proteins, bovine serum albumin (BSA), soybean peroxidase (SBP), and α-chymotrypsin (CT), with and without CNTs, were performed to examine the influence of CNTs on the structure and dynamics of these proteins. It was found that nanotubes bind to these proteins and influence their structure. Our results show that the UNRES force field can be used for further studies of CNT-protein systems with 3–4 order of magnitude larger timescale than using regular all-atom force fields.Graphical abstractBovine serum albumin (BSA), soybean peroxidase (SBP), and α-chymotrypsin (CT), with and without CNTsᅟ

How fast fast-folding proteins fold in silico

In reported microcanonical molecular dynamics simulations, fast-folding proteins CLN025 and Trp-cage autonomously folded to experimentally determined native conformations. However, the folding times of these proteins derived from the simulations were more than 4–10 times longer than their experimental values. This article reports autonomous folding of CLN025 and Trp-cage in isobaric–isothermal molecular dynamics simulations with agreements within factors of 0.69–1.75 between simulated and experimental folding times at different temperatures. These results show that CLN025 and Trp-cage can now autonomously fold in silico as fast as in experiments, and suggest that the accuracy of folding simulations for fast-folding proteins begins to overlap with the accuracy of folding experiments. This opens new prospects of developing computer algorithms that can predict both ensembles of conformations and their interconversion rates for a protein from its sequence for artificial intelligence on how and when a protein acts as a receiver, switch, and relay to facilitate various subcellular-to-tissue communications. Then the genetic information that encodes proteins can be better read in the context of intricate biological functions.

Microcanonical molecular dynamics simulation of the vitreous phase transition of poly(binaphthylphosphazene)

Properties of poly(binaphthoxyphosphazene)s were investigated employing molecular dynamics simulations in the original microcanonical ensemble. In particular we have applied our methodology to the case of isotactic poly-(2,2´-dioxy-1,1´-binaphthylphosphazene) (P-DBNP) in order to estimate the vitreous transition temperature Tg, the energy barrier Δμβ for β-relaxation of the BNP unit, and the thermodynamic entropy of the system. According to our results, the transition is consistent with the Adam-Gibbs model and the specific rotation of α and β-relaxations, which varies significantly with the molecular weight Mw. Our results suggest a rapid interconversion between the different proportions of the chiral R and S repeating units (70% and 30%, respectively) of the non-isotactic copolymer below 523K and a slow atropisomerization of the DBNP units above 523K which, in agreement with recent predictions on glassy polymer matrices, becomes much faster as the temperature approaches Tg~573K.

Extended Lagrangian formulation of charge-constrained tight-binding molecular dynamics.

The extended Lagrangian Born-Oppenheimer molecular dynamics formalism [Niklasson, Phys. Rev. Lett., 2008, 100, 123004] has been applied to a tight-binding model under the constraint of local charge neutrality to yield microcanonical trajectories with both precise, long-term energy conservation and a reduced number of self-consistent field optimizations at each time step. The extended Lagrangian molecular dynamics formalism restores time reversal symmetry in the propagation of the electronic degrees of freedom, and it enables the efficient and accurate self-consistent optimization of the chemical potential and atomwise potential energy shifts in the on-site elements of the tight-binding Hamiltonian that are required when enforcing local charge neutrality. These capabilities are illustrated with microcanonical molecular dynamics simulations of a small metallic cluster using an sd-valent tight-binding model for titanium. The effects of weak dissipation on the propagation of the auxiliary degrees of freedom for the chemical potential and on-site Hamiltonian matrix elements that is used to counteract the accumulation of numerical noise during trajectories was also investigated.

Thermal phonons and defects in semiconductors: The physical reason why defects reduce heat flow, and how to control it

It is generally accepted that heat-carrying phonons in materials scatter off each other (normal or Umklapp scattering) as well as off defects. This assumes static defects, implies quasi-instantaneous interactions and at least some momentum transfer. However, when defect dynamics are explicitly included, the nature of phonon-defect interactions becomes more subtle. Ab initio microcanonical molecular-dynamics simulations show that (1) spatially localized vibrational modes (SLMs), associated with all types of defects in semiconductors, can trap thermal phonons; (2) the vibrational lifetimes of excitations in SLMs are one to two orders of magnitude longer (dozens to hundreds of periods of oscillation) than those of bulk phonons of similar frequency; (3) it is phonon trapping by defects (in SLMs) rather than bulk phonon scattering, which reduces the flow of heat; and (4) the decay of trapped phonons and therefore heat flow can be predicted and controlled—at least to some extent—by the use of carefully selected interfaces and δ layers.

ANALYSIS OF THE MELTING PROCESS OF THE 14-PARTICLE MORSE CLUSTER THROUGH THE ATOMIC EQUIVALENCE INDEXES

In this paper, by using the standard microcanonical molecular-dynamics simulations, the melting processes of the 14-particle Morse clusters are investigated. The three steps melting process found in 2004 by one of the authors of this paper are analyzed further by the atomic equivalence indexes. The properties of the structure of the cluster in different melting processes are explored by the atomic equivalence indexes. It is found that in the melting process for most of the time the cluster keeps the basic structure frame similar to that in the solid-like state while there are very short time periods when the particles of the cluster exchange their locations.

Heat Transfer Calculations for Decomposition of Structure I Methane Hydrates by Molecular Dynamics Simulation

Microcanonical ensemble molecular dynamics simulations of structure I methane hydrate is presented in this work to study the endothermic decomposition process. The mechanism of decomposition of methane hydrate as a function of time was explained at the molecular level. The initial temperature and pressure of the simulation were chosen so as to depict the natural gas hydrate in conditions of oceanic sediments. A more realistic strategy was developed to perform the microcanonical ensemble simulation of solid–liquid interface of hydrate and amorphous water. Two water models, SPC/E and TIP4P, were used for the simulations, and the results of the simulations were compared. Heat transfer calculations were performed on the adiabatic system, and an attempt has been made to fit the MD simulation results to the heat balance equations derived from the heat transfer calculations. Estimates of the properties at the macroscopic scale, like the equilibrium temperature of methane hydrate and rate of supply of hot water f...

Fast method for quantum mechanical molecular dynamics

As the processing power available for scientific computing grows, first-principles Born-Oppenheimer molecular dynamics simulations are becoming increasingly popular for the study of a wide range of problems in materials science, chemistry, and biology. Nevertheless, the computational cost of Born-Oppenheimer molecular dynamics still remains prohibitively large for many potential applications. Here we show how to avoid a major computational bottleneck: the self-consistent-field optimization prior to force calculations. The optimization-free quantum mechanical molecular dynamics method gives trajectories that are almost indistinguishable from an ``exact'' microcanonical Born-Oppenheimer molecular dynamics simulation even when low-prefactor linear scaling sparse matrix algebra is used. Our findings show that the computational gap between classical and quantum mechanical molecular dynamics simulations can be significantly reduced.

Molecular simulation of the phase behavior of fluids and fluid mixtures using the synthetic method

A new molecular simulation procedure is reported for determining the phase behavior of fluids and fluid mixtures, which closely follows the experimental synthetic method. The simulation procedure can be implemented using Monte Calro or molecular dynamics in either the microcanonical or canonical statistical ensembles. Microcanonical molecular dynamics simulations are reported for the phase behavior of both the pure Lennard-Jones fluid and a Lennard-Jones mixture. The vapor pressures for the pure fluid are in good agreement with Monte Carlo Gibbs ensemble and Gibbs-Duhem calculations. The Lennard-Jones mixture is composed of equal size particles, with dissimilar energy parameters (ε(2)/ε(1) = 1/2, ε(12)/ε(1) = 1/√2). The binary Lennard-Jones mixture exhibits liquid-liquid equilibria at high pressures and the simulation procedure allows us to estimate the coordinates of the high-pressure branch of the critical curve.

Anomalous Behavior of Gear¡?s Predictor-corrector Algorithm in a Molecular Dynamics Simulation of Amorphous Silicon

An anomalous behavior of Gear’s predictor-corrector algorithm is observed during a microcanonical molecular dynamics simulation of amorphous silicon. While the amorphous silicon network generated by quenching was being annealed, an anomaly is started: both the total energy and the temperature of the system spontaneously and indefinitely increased, in violation of energy conservation. No such phenomenon was observed when the integrator was switched to the Verlet algorithm. To remove this anomaly with the Gear’s algorithm, it was necessary to reduce the size of the time step to 1/32 of the “usual value ( 1/100 of the fastest vibration period in the system)” at least. Our results show that Gear’s algorithm may become very unstable in a realistic simulation when the time step is not small enough, that is, orders of magnitude smaller than conventional value.

Nanoresonator Based on Relative Vibrations of the Walls of Carbon Nanotubes

A new type of ultrahigh frequency resonator based on the relative vibrations of carbon nanotubes walls is proposed. Microcanonical molecular dynamics simulations performed for a model resonator based on the (9,0)@(18,0) double-walled carbon nanotube with the movable outer wall give the value of frequency equal to 130 GHz and the quality factor up to 500. Possible applications of the resonator are discussed. The estimated mass sensitivity of the mass sensor based on the proposed resonator reaches the atomic mass limit at liquid helium temperature. The possibility of resonator based on relative vibrations of graphene layers is also considered.