Search

Biased and accelerated molecular simulations (BAMS) are widely used tools to observe relevant molecular phenomena occurring on time scales inaccessible to standard molecular dynamics, but evaluation of the physical time scales involved in the processes is not directly possible from them. For this reason, the problem of recovering dynamics from such kinds of simulations is the object of very active research due to the relevant theoretical and practical implications of dynamics on the properties of both natural and synthetic molecular systems. In a recent paper [A. Rapallo et al., J. Comput. Chem. 42, 586–599 (2021)], it has been shown how the coupling of BAMS (which destroys the dynamics but allows to calculate average properties) with Extended Diffusion Theory (EDT) (which requires input appropriate equilibrium averages calculated over the BAMS trajectories) allows to effectively use the Smoluchowski equation to calculate the orientational time correlation function of the head–tail unit vector defined over a peptide in water solution. Orientational relaxation of this vector is the result of the coupling of internal molecular motions with overall molecular rotation, and it was very well described by correlation functions expressed in terms of weighted sums of suitable time-exponentially decaying functions, in agreement with a Brownian diffusive regime. However, situations occur where exponentially decaying functions are no longer appropriate to capture the actual dynamical behavior, which exhibits persistent long time correlations, compatible with the so called subdiffusive regimes. In this paper, a generalization of EDT will be given, exploiting a fractional Smoluchowski equation (FEDT) to capture the non-exponential character observed in the relaxation of intramolecular distances and molecular radius of gyration, whose dynamics depend on internal molecular motions only. The calculation methods, proper to EDT, are adapted to implement the generalization of the theory, and the resulting algorithm confirms FEDT as a tool of practical value in recovering dynamics from BAMS, to be used in general situations, involving both regular and anomalous diffusion regimes. [ABSTRACT FROM AUTHOR]

The first paper of this series [J. Chem. Phys. 158, 034103 (2023)] demonstrated that excess entropy scaling holds for both fine-grained and corresponding coarse-grained (CG) systems. Despite its universality, a more exact determination of the scaling relationship was not possible due to the semi-empirical nature. In this second paper, an analytical excess entropy scaling relation is derived for bottom-up CG systems. At the single-site CG resolution, effective hard sphere systems are constructed that yield near-identical dynamical properties as the target CG systems by taking advantage of how hard sphere dynamics and excess entropy can be analytically expressed in terms of the liquid packing fraction. Inspired by classical equilibrium perturbation theories and recent advances in constructing hard sphere models for predicting activated dynamics of supercooled liquids, we propose a new approach for understanding the diffusion of molecular liquids in the normal regime using hard sphere reference fluids. The proposed "fluctuation matching" is designed to have the same amplitude of long wavelength density fluctuations (dimensionless compressibility) as the CG system. Utilizing the Enskog theory to derive an expression for hard sphere diffusion coefficients, a bridge between the CG dynamics and excess entropy is then established. The CG diffusion coefficient can be roughly estimated using various equations of the state, and an accurate prediction of accelerated CG dynamics at different temperatures is also possible in advance of running any CG simulation. By introducing another layer of coarsening, these findings provide a more rigorous method to assess excess entropy scaling and understand the accelerated CG dynamics of molecular fluids. [ABSTRACT FROM AUTHOR]

This paper series aims to establish a complete correspondence between fine-grained (FG) and coarse-grained (CG) dynamics by way of excess entropy scaling (introduced in Paper I). While Paper II successfully captured translational motions in CG systems using a hard sphere mapping, the absence of rotational motions in single-site CG models introduces differences between FG and CG dynamics. In this third paper, our objective is to faithfully recover atomistic diffusion coefficients from CG dynamics by incorporating rotational dynamics. By extracting FG rotational diffusion, we unravel, for the first time reported to our knowledge, a universality in excess entropy scaling between the rotational and translational diffusion. Once the missing rotational dynamics are integrated into the CG translational dynamics, an effective translation-rotation coupling becomes essential. We propose two different approaches for estimating this coupling parameter: the rough hard sphere theory with acentric factor (temperature-independent) or the rough Lennard-Jones model with CG attractions (temperature-dependent). Altogether, we demonstrate that FG diffusion coefficients can be recovered from CG diffusion coefficients by (1) incorporating "entropy-free" rotational diffusion with translation-rotation coupling and (2) recapturing the missing entropy. Our findings shed light on the fundamental relationship between FG and CG dynamics in molecular fluids. [ABSTRACT FROM AUTHOR]

Protein folding in confined media has attracted wide attention over the past 15 years due to its importance to both in vivo and in vitro applications. It is generally believed that protein stability increases by decreasing the size of the confining medium, if the medium’s walls are repulsive, and that the maximum folding temperature in confinement is in a pore whose size D0 is only slightly larger than the smallest dimension of a protein’s folded state. Until recently, the stability of proteins in pores with a size very close to that of the folded state has not received the attention it deserves. In a previous paper [L. Javidpour and M. Sahimi, J. Chem. Phys. 135, 125101 (2011)], we showed that, contrary to the current theoretical predictions, the maximum folding temperature occurs in larger pores for smaller α-helices. Moreover, in very tight pores, the free energy surface becomes rough, giving rise to a new barrier for protein folding close to the unfolded state. In contrast to unbounded domains, in small nanopores proteins with an α-helical native state that contain the β structures are entropically stabilized implying that folding rates decrease notably and that the free energy surface becomes rougher. In view of the potential significance of such results to interpretation of many sets of experimental data that could not be explained by the current theories, particularly the reported anomalously low rates of folding and the importance of entropic effects on proteins’ misfolded states in highly confined environments, we address the following question in the present paper: To what extent the geometry of a confined medium affects the stability and folding rates of proteins? Using millisecond-long molecular dynamics simulations, we study the problem in three types of confining media, namely, cylindrical and slit pores and spherical cavities. Most importantly, we find that the prediction of the previous theories that the dependence of the maximum folding temperature Tf on the size D of a confined medium occurs in larger media for larger proteins is correct only in spherical geometry, whereas the opposite is true in the two other geometries that we study. Also studied is the effect of the strength of the interaction between the confined media’s walls and the proteins. If the walls are only weakly or moderately attractive, a complex behavior emerges that depends on the size of the confining medium. [ABSTRACT FROM AUTHOR]

Stable crystalline structures of confined water can be different from bulk ice. In Paper I [T. Yagasaki et al., J. Chem. Phys. 151, 064702 (2019)] of this study, it was shown, using molecular dynamics (MD) simulations, that a zeolite-like ice structure forms in nanobrushes consisting of (6,6) carbon nanotubes (CNTs) when the CNTs are located in a triangle arrangement. The melting temperature of the zeolite-like ice structure is much higher than the melting temperature of ice Ih when the distance between the surfaces of CNTs is ∼0.94 nm, which is the best spacing for the bilayer structure of water. In this paper, we perform MD simulations of nanobrushes of CNTs that are different from (6,6) CNTs in radius. Several new porous ice structures form spontaneously in the MD simulations. A stable porous ice forms when the radius of its cavities matches the radius of the CNTs well. All cylindrical porous ice structures found in this study can be decomposed into a small number of structural blocks. We provide a new protocol to classify cylindrical porous ice crystals on the basis of this decomposition. [ABSTRACT FROM AUTHOR]

This work follows a companion article, which will be referred to as Paper I [Campeggio et al., J. Chem. Phys. 158, 244104 (2023)] in which a quantum-stochastic Liouville equation for the description of the quantum–classical dynamics of a molecule in a dissipative bath has been formulated in curvilinear internal coordinates. In such an approach, the coordinates of the system are separated into three subsets: the quantum coordinates, the classical relevant nuclear degrees of freedom, and the classical irrelevant (bath) coordinates. The equation has been derived in natural internal coordinates, which are bond lengths, bond angles, and dihedral angles. The resulting equation needs to be parameterized. In particular, one needs to compute the potential energy surfaces, the friction tensor, and the rate constants for the nonradiative jumps among the quantum states. While standard methods exist for the calculation of energy and dissipative properties, an efficient evaluation of the transition rates needs to be developed. In this paper, an approximated treatment is introduced, which leads to a simple explicit formula with a single adjustable parameter. Such an approximated expression is compared with the exact calculation of transition rates obtained via molecular dynamics simulations. To make such a comparison possible, a simple sandbox system has been used, with two quantum states and a single internal coordinate (together with its conjugate momentum). Results show that the adjustable parameter, which is an effective decoherence time, can be parameterized from the effective relaxation times of the autocorrelation functions of the conjugated momenta of the relevant nuclear coordinates. [ABSTRACT FROM AUTHOR]

Ninety years ago, Wigner derived the leading order expansion term in ℏ2 for the tunneling rate through a symmetric barrier. His derivation included two contributions: one came from the parabolic barrier, but a second term involved the fourth-order derivative of the potential at the barrier top. He left us with a challenge, which is answered in this paper, to derive the same but for an asymmetric barrier. A crucial element of the derivation is obtaining the ℏ2 expansion term for the projection operator, which appears in the flux-side expression for the rate. It is also reassuring that an analytical calculation of semiclassical transition state theory (TST) reproduces the anharmonic corrections to the leading order of ℏ2. The efficacy of the resulting expression is demonstrated for an Eckart barrier, leading to the conclusion that especially when considering heavy atom tunneling, one should use the expansion derived in this paper, rather than the parabolic barrier approximation. The rate expression derived here reveals how the classical TST limit is approached as a function of ℏ and, thus, provides critical insights to understand the validity of popular approximate theories, such as the classical Wigner, centroid molecular dynamics, and ring polymer molecular dynamics methods. [ABSTRACT FROM AUTHOR]

Erratum: "Physics-based, neural network force fields for reactive molecular dynamics: Investigation of carbene formation from [EMIM+][OAc-]" [J. Chem. Phys. We have recomputed the PMFs with AIMD, utilizing a consistent collective variable definition of Eq. (1). This was done using the plumed/ASE interface with a CP2K calculator.[[2], [5]] Corrected PMFs that employ the collective variable of Eq. (1) for all calculations are shown in the corrected Fig. [Extracted from the article]

Vibrational spectroscopy is a useful technique for probing chemical environments. The development of models that can reproduce the spectra of nitriles and azides is valuable because these probes are uniquely suited for investigating complex systems. Empirical vibrational spectroscopic maps are commonly employed to obtain the instantaneous vibrational frequencies during molecular dynamics simulations but often fail to adequately describe the behavior of these probes, especially in its transferability to a diverse range of environments. In this paper, we demonstrate several reasons for the difficulty in constructing a general-purpose vibrational map for methyl thiocyanate (MeSCN), a model for cyanylated biological probes. In particular, we found that electrostatics alone are not a sufficient metric to categorize the environments of different solvents, and the dominant features in intermolecular interactions in the energy landscape vary from solvent to solvent. Consequently, common vibrational mapping schemes do not cover all essential interaction terms adequately, especially in the treatment of van der Waals interactions. Quantum vibrational perturbation (QVP) theory, along with a combined quantum mechanical and molecular mechanical potential for solute–solvent interactions, is an alternative and efficient modeling technique, which is compared in this paper, to yield spectroscopic results in good agreement with experimental FTIR. QVP has been used to analyze the computational data, revealing the shortcomings of the vibrational maps for MeSCN in different solvents. The results indicate that insights from QVP analysis can be used to enhance the transferability of vibrational maps in future studies. [ABSTRACT FROM AUTHOR]

Protein–calixarenes binding plays an increasingly central role in many applications, spanning from molecular recognition to drug delivery strategies and protein inhibition. These ligands obey a specific bio-supramolecular chemistry, which can be revealed by computational approaches, such as molecular dynamics simulations. In this paper, we rely on all-atom, explicit-solvent molecular dynamics simulations to capture the electrostatically driven association of a phosphonated calix-[4]-arene with cytochome-C, which critically relies on surface-exposed paired lysines. Beyond two binding sites identified in direct agreement with the x-ray structure, the association has a larger structural impact on the protein dynamics. Then, our simulations allow a direct comparison to analogous calixarenes, namely, sulfonato, similarly reported as "molecular glue." Our work can contribute to a robust in silico predictive tool to assess binding sites for any given protein of interest for crystallization, with the specificity of a macromolecular cage whose endo/exo orientation plays a role in the binding. [ABSTRACT FROM AUTHOR]

In order to prepare highly heat-resistant packaging insulation materials, in this paper, bismaleimide/epoxy resin (BMI/EP55) composites with different contents of BMI were prepared by melt blending BMI into amino tetrafunctional and phenolic epoxy resin (at a ratio of 5:5). The microstructures and thermal and electrical properties of the composites were tested. The electrostatic potential distribution, energy level distribution, and molecular orbitals of BMI were calculated using Gaussian. The results showed that the carbonyl group in BMI is highly electronegative, implying that the carbonyl group has a strong electron trapping ability. The thermal decomposition temperature of the composites gradually increased with the increase of BMI content, and the 20% BMI/EP55 composites had the highest heat-resistance index, along with a glass transition temperature (Tg) of >250 °C. At different test temperatures, with increase in the BMI content, the conductivity of epoxy resin composites showed a tendency to first decrease and then increase, the breakdown field strength showed a tendency to first increase and then decrease, and the dielectric constant was gradually decreased. Two trap centers were present simultaneously in the composites, where the shallow trap energy level is the deepest in 20% BMI/EP composites and the deep trap energy level is the deepest in 10% BMI/EP55 composites. Correspondingly, the 10% BMI/EP55 composite had a slower charge decay rate, while the 20% BMI/EP55 had a faster charge decay rate. In summary, the BMI/EP55 composites with high heat resistance and insulating properties were prepared in this study, which provided ideas for preparing high-temperature packaging insulating materials. [ABSTRACT FROM AUTHOR]

The large discrepancy among the nucleation kinetics extracted from experimental measurements and computer simulations and the prediction of the classical nucleation theory (CNT) has stimulated intense arguments about its origin in the past decades, which is crucially relevant to the validity of the CNT. In this paper, we investigate the atomistic mechanism of the nucleation in liquid Al in contact with amorphous substrates with atomic-level smooth/rough surfaces, using molecular dynamics (MD) simulations. This study reveals that the slightly distorted local fcc/hcp structures in amorphous substrates with smooth surfaces can promote heterogeneous nucleation through a structural templating mechanism, and on the other hand, homogeneous nucleation will occur at a larger undercooling through a fluctuation mechanism if the surface is rough. Thus, some impurities, previously thought to be impotent, could be activated in the homogeneous nucleation experiments. We further find that the initial growth of the nucleus on smooth surfaces of amorphous substrates is one order of magnitude faster than that in homogeneous nucleation. Both these factors could significantly contribute to the discrepancy in the nucleation kinetics. This study is also supported by a recent study of the synthesis of high-entropy alloy nanoparticles assisted with the liquid metal Ga [Cao et al., Nature 619, 73 (2023)]. In this study, we established that the boundary existed between homogeneous and heterogeneous nucleation, i.e., the structural templating is a general mechanism for heterogeneous nucleation, and in its absence, homogeneous nucleation will occur through the fluctuation mechanism. This study provides an in-depth understanding of the nucleation theory and experiments. [ABSTRACT FROM AUTHOR]

The present paper is the second part of a series of papers aimed at assessing the accuracy of alchemical computational approaches based on nonequilibrium techniques for solvation free energy of organic molecules in the context of molecular dynamics simulations. In Paper I [Procacci, J. Chem. Phys. 151, 144113 (2019)], we dealt with bidirectional estimates of solvation free energies using nonequilibrium approaches. Here, we assess accuracy and precision of unidirectional estimates with the focus on the Gaussian and Jarzynski estimators. We present a very simple methodology to increase the statistics in the work distribution, hence boosting the accuracy and precision of the Jarzynski unidirectional estimates at no extra cost, exploiting the observed decorrelation between the random variables represented by the Lennard-Jones solute-solvent recoupling or decoupling work and by the electrostatic work due to the charging/discharging of the solute in the solvent. [ABSTRACT FROM AUTHOR]

Introducing nanorods into a polymeric matrix can enhance the physical and mechanical properties of the resulting material. In this paper, we focus on understanding the dispersion and orientation patterns of nanorods in an unentangled polymer melt, particularly as a function of nanorod concentration, using molecular dynamics simulations. The system is comprised of flexible polymer chains and multi-thread nanorods that are equilibrated in the NPT ensemble. All interactions are purely repulsive except for those between polymers and rods. Results with attractive vs repulsive polymer–rod interactions are compared and contrasted. The concentration of rods has a direct impact on the phase behavior of the system. At lower concentrations, rods phase separate into nematic clusters, whereas at higher concentrations more isotropic and less structured rod configurations are observed. A detailed examination of the conformation of the polymer chains near the rod surface shows extension of the chains along the director of the rods (especially within clusters). The dispersion and orientation of the nanorods are a result of the competition between depletion entropic forces responsible for the formation of rod clusters, the enthalpic effects that improve mixing of rods and polymer, and entropic losses of polymers interpenetrating rod clusters. [ABSTRACT FROM AUTHOR]

Fewest-switches surface hopping (FSSH) has emerged as one of the leading methods for modeling the quantum dynamics of molecular systems. While its original formulation was limited to adiabatic populations, the growing interest in the application of FSSH to coherent phenomena prompts the question of how one should construct a complete density matrix based on FSSH trajectories. A straightforward solution is to define adiabatic coherences based on wavefunction coefficients. In this paper, we demonstrate that inconsistencies introduced in the density matrix through such treatment may lead to a violation of positivity. We furthermore show that a recently proposed coherent generalization of FSSH results in density matrices that satisfy positivity while yielding improved accuracy throughout much (but not all) of parameter space. [ABSTRACT FROM AUTHOR]

Using an approach based on the diffusion analog of the Cattaneo–Vernotte differential model, we find the exact analytical solution to the corresponding time-dependent linear hyperbolic initial boundary value problem, describing irreversible diffusion-controlled reactions under Smoluchowski's boundary condition on a spherical sink. By means of this solution, we extend exact analytical calculations for the time-dependent classical Smoluchowski rate coefficient to the case that includes the so-called inertial effects, occurring in the host media with finite relaxation times. We also present a brief survey of Smoluchowski's theory and its various subsequent refinements, including works devoted to the description of the short-time behavior of Brownian particles. In this paper, we managed to show that a known Rice's formula, commonly recognized earlier as an exact reaction rate coefficient for the case of hyperbolic diffusion, turned out to be only its approximation being a uniform upper bound of the exact value. Here, the obtained formula seems to be of great significance for bridging a known gap between an analytically estimated rate coefficient on the one hand and molecular dynamics simulations together with experimentally observed results for the short times regime on the other hand. A particular emphasis has been placed on the rigorous mathematical treatment and important properties of the relevant initial boundary value problems in parabolic and hyperbolic diffusion theories. [ABSTRACT FROM AUTHOR]

In this paper, we investigate prenucleation and heterogeneous nucleation in the liquid Pb/solid Al system as an example of systems with large lattice misfit using molecular dynamics simulation. Solid Pb and Al have a large positive lattice misfit (f) of 18.2% along the densely packed [110] direction. This study reveals that prenucleation occurs at 600 K (an undercooling of 15 K), and a 2-dimensional (2D) ordered structure forms at the interface with a coincidence site lattice (CSL) between the first Pb and first Al layers. The CSL accommodates the major part of the f, and only a small residual lattice misfit (fr) of 1.9% remains. The formation of the CSL transforms the original substrate into a considerably potent nucleant, where the first Pb layer becomes the new surface layer of the substrate. At an undercooling of about 22 K, nucleation proceeds by merging 2D ordered structure through structural templating: the second Pb layer is epitaxial to the CSL Pb layer, the third Pb layer largely accommodates the fr, and the fourth Pb layer is a nearly perfect crystalline plane. Further analysis indicates that the interface with the CSL has a lower interfacial energy than with a cube-to-cube orientation relationship. For the first time, we established that the CSL was an effective mechanism to accommodate the f for systems with a large positive misfits. Heterogeneous nucleation is governed not by a single mechanism (misfit dislocations in Turnbull's model), but instead by various mechanisms depending on f. This study sheds new light on the atomistic mechanism of heterogeneous nucleation. [ABSTRACT FROM AUTHOR]

This paper presents a state-averaged complete active space self-consistent field (SA-CASSCF) in the atomic multipole optimized energetics for biomolecular application (AMOEBA) polarizable water model, which enables rigorous simulation of non-adiabatic molecular dynamics with nonequilibrium solvation effects. The molecular orbital and configuration interaction coefficients of the solute wavefunction, and the induced dipoles on solvent atoms, are solved by minimizing the state averaged energy variationally. In particular, by formulating AMOEBA water models and the polarizable continuum model (PCM) in a unified way, the algorithms developed for computing SA-CASSCF/PCM energies, analytical gradients, and non-adiabatic couplings in our previous work can be generalized to SA-CASSCF/AMOEBA by properly substituting a specific list of variables. Implementation of this method will be discussed with the emphasis on how the calculations of different terms are partitioned between the quantum chemistry and molecular mechanics codes. We will present and discuss results that demonstrate the accuracy and performance of the implementation. Next, we will discuss results that compare three solvent models that work with SA-CASSCF, i.e., PCM, fixed-charge force fields, and the newly implemented AMOEBA. Finally, the new SA-CASSCF/AMOEBA method has been interfaced with the ab initio multiple spawning method to carry out non-adiabatic molecular dynamics simulations. This method is demonstrated by simulating the photodynamics of the model retinal protonated Schiff base molecule in water. [ABSTRACT FROM AUTHOR]

Hydration and, in particular, the coordination number of a metal ion is of paramount importance as it defines many of its (bio)physicochemical properties. It is not only essential for understanding its behavior in aqueous solutions but also determines the metal ion reference state and its binding energy to (bio)molecules. In this paper, for divalent metal cations Ca 2 + , Cd 2 + , Cu 2 + , Fe 2 + , Hg 2 + , Mg 2 + , Ni 2 + , Pb 2 + , and Zn 2 + , we compare two approaches for predicting hydration numbers: (1) a mixed explicit/continuum DFT-D3//COSMO-RS solvation model and (2) density functional theory based ab initio molecular dynamics. The former approach is employed to calculate the Gibbs free energy change for the sequential hydration reactions, starting from [M( H 2 O)2]2+ aqua complexes to [M( H 2 O) 9]2+, allowing explicit water molecules to bind in the first or second coordination sphere and determining the most stable [M( H 2 O) n]2+ structure. In the latter approach, the hydration number is obtained by integrating the ion–water radial distribution function. With a couple of exceptions, the metal ion hydration numbers predicted by the two approaches are in mutual agreement, as well as in agreement with the experimental data. [ABSTRACT FROM AUTHOR]

A Markov state model is a powerful tool that can be used to track the evolution of populations of configurations in an atomistic representation of a protein. For a coarse-grained linear chain model with discontinuous interactions, the transition rates among states that appear in the Markov model when the monomer dynamics is diffusive can be determined by computing the relative entropy of states and their mean first passage times, quantities that are unchanged by the specification of the energies of the relevant states. In this paper, we verify the folding dynamics described by a diffusive linear chain model of the crambin protein in three distinct solvent systems, each differing in complexity: a hard-sphere solvent, a solvent undergoing multi-particle collision dynamics, and an implicit solvent model. The predicted transition rates among configurations agree quantitatively with those observed in explicit molecular dynamics simulations for all three solvent models. These results suggest that the local monomer–monomer interactions provide sufficient friction for the monomer dynamics to be diffusive on timescales relevant to changes in conformation. Factors such as structural ordering and dynamic hydrodynamic effects appear to have minimal influence on transition rates within the studied solvent densities. [ABSTRACT FROM AUTHOR]

The network structure within polymers significantly influences their mechanical properties, including their strength, toughness, and fatigue resistance. All-atom molecular dynamics (AAMD) simulations offer a method to investigate the energy dissipation mechanism within polymers during deformation and fracture; Such an approach is, however, computationally inefficient when used to analyze polymers with complex network structures, such as the common chemically double-networked hydrogels. Alternatively, coarse-grained molecular dynamics (CGMD) models, which reduce the computational degrees of freedom by concentrating a set of adjacent atoms into a coarse-grained bead, can be employed. In CGMD simulations, a coarse-grained force field (CGFF) is a critical factor affecting the simulation accuracy. In this paper, we proposed a data-based method for predicting the CGFF parameters to improve the simulation efficiency of complex cross-linked network in polymers. Here, we utilized a typical chemically double-networked hydrogel as an example. An artificial neural network was selected, and it was trained with the tensile stress–strain data from the CGMD simulations using different CGFF parameters. The CGMD simulations using the predicted CGFF parameters show good agreement with the AAMD simulations and are almost fifty times faster. The data-inspired CGMD model presented here broadens the applicability of molecular dynamics simulations to cross-linked polymers and has the potential to provide insights that will aid the design of polymers with desirable mechanical properties. [ABSTRACT FROM AUTHOR]

The fewest switches surface hopping (FSSH) method proposed by Tully in 1990 [Tully, J. Chem. Phys. 93, 1061 (1990)]—along with its many later variations—forms the basis for most practical simulations of molecular dynamics with electronic transitions in realistic systems. Despite its popularity, a rigorous formal derivation of the algorithm has yet to be achieved. In this paper, we derive the energy-conserving momentum jumps employed by FSSH from the perspective of quantum trajectory surface hopping (QTSH) [Martens, J. Phys. Chem. A 123, 1110 (2019)]. In the limit of localized nonadiabatic transitions, simple mathematical and physical arguments allow the FSSH algorithm to be derived from first principles. For general processes, the quantum forces characterizing the QTSH method provide accurate results for nonadiabatic dynamics with rigorous energy conservation, at the ensemble level, within the consistency of the underlying stochastic surface hopping without resorting to the artificial momentum rescaling of FSSH. [ABSTRACT FROM AUTHOR]

Polymer crystallization is a long-standing interesting problem both in fundamental polymer physics and in polymer manufacturing. Fundamentally, the connectivity of the macromolecules provides a unique feature for the study of nucleation and growth of crystals in contrast to the crystallization of smaller molecules. In addition, understanding the crystallization in polymers is industrially important due to the necessity of its control to achieve mechanically durable plastic materials. Molecular dynamics simulations offer a suitable way of studying this phenomenon due to their capability to probe the small time and length scales that are characteristic of nucleation. In this paper, we use a long alkane chain model to study nucleation and the growth of polyethylene crystals both within bulk and nano-droplets whose diameters range from about 6 to 11 nm. It is found that the droplets approach being more spherical during the nucleation whereas they deviate from this shape during the growth regime. Strikingly, a mean first passage time analysis indicates that the nucleation rate per unit volume decreases as the droplet size is increased. Finally, visual inspection of the simulation snapshots suggests that nuclei majorly emerged from the surface of the droplets. [ABSTRACT FROM AUTHOR]

Polymer crystallization has long been a fascinating problem and is still attracting many researchers. Most of the previous simulations are concentrated on clarifying the universal aspects of polymer crystallization using model linear polymers such as polyethylene. We are recently focusing on a nearly untouched but very interesting problem of chiral selecting crystallization in helical polymers. We previously proposed a stepwise approach using two kinds of helical polymers, simple "bare" helical polymers made of backbone atoms only such as polyoxymethylene (POM) and "general" helical polymers containing complicated side groups such as isotactic polypropylene. We have already reported on the crystallization in oligomeric POM-like helix but have observed only weak chiral selectivity during crystallization. In the present paper, we investigate the crystallization of sufficiently long POM-like polymer both from the isotropic melt and from the highly stretched melt. We find in both cases that the polymer shows a clear chiral selecting crystallization. Specifically, the observation of a single crystal growing from the isotropic melt is very illuminating. It shows that the crystal thickness and the crystal chirality are closely correlated; thicker crystals show definite chirality while thinner ones are mostly mixtures of the R- and the L-handed stems. The single crystal is found to have a marked lenticular shape, where the thinner growth front, since being made of the mixture, shows no chiral selectivity. The final chiral crystal is found to be completed through helix reversal processes within thicker regions. [ABSTRACT FROM AUTHOR]

Molecular dynamics (MD) is the method of choice for understanding the structure, function, and interactions of molecules. However, MD simulations are limited by the strong metastability of many molecules, which traps them in a single conformation basin for an extended amount of time. Enhanced sampling techniques, such as metadynamics and replica exchange, have been developed to overcome this limitation and accelerate the exploration of complex free energy landscapes. In this paper, we propose Vendi Sampling, a replica-based algorithm for increasing the efficiency and efficacy of the exploration of molecular conformation spaces. In Vendi sampling, replicas are simulated in parallel and coupled via a global statistical measure, the Vendi Score, to enhance diversity. Vendi sampling allows for the recovery of unbiased sampling statistics and dramatically improves sampling efficiency. We demonstrate the effectiveness of Vendi sampling in improving molecular dynamics simulations by showing significant improvements in coverage and mixing between metastable states and convergence of free energy estimates for four common benchmarks, including Alanine Dipeptide and Chignolin. [ABSTRACT FROM AUTHOR]

Interpreting high-dimensional data from molecular dynamics simulations is a persistent challenge. In this paper, we show that for a small peptide, deca-alanine, metastable states can be identified through a neural net based on structural information alone. While processing molecular dynamics data, dimensionality reduction is a necessary step that projects high-dimensional data onto a low-dimensional representation that, ideally, captures the conformational changes in the underlying data. Conventional methods make use of the temporal information contained in trajectories generated through integrating the equations of motion, which forgoes more efficient sampling schemes. We demonstrate that EncoderMap, an autoencoder architecture with an additional distance metric, can find a suitable low-dimensional representation to identify long-lived molecular conformations using exclusively structural information. For deca-alanine, which exhibits several helix-forming pathways, we show that this approach allows us to combine simulations with different biasing forces and yields representations comparable in quality to other established methods. Our results contribute to computational strategies for the rapid automatic exploration of the configuration space of peptides and proteins. [ABSTRACT FROM AUTHOR]

Finding a low dimensional representation of data from long-timescale trajectories of biomolecular processes, such as protein folding or ligand–receptor binding, is of fundamental importance, and kinetic models, such as Markov modeling, have proven useful in describing the kinetics of these systems. Recently, an unsupervised machine learning technique called VAMPNet was introduced to learn the low dimensional representation and the linear dynamical model in an end-to-end manner. VAMPNet is based on the variational approach for Markov processes and relies on neural networks to learn the coarse-grained dynamics. In this paper, we combine VAMPNet and graph neural networks to generate an end-to-end framework to efficiently learn high-level dynamics and metastable states from the long-timescale molecular dynamics trajectories. This method bears the advantages of graph representation learning and uses graph message passing operations to generate an embedding for each datapoint, which is used in the VAMPNet to generate a coarse-grained dynamical model. This type of molecular representation results in a higher resolution and a more interpretable Markov model than the standard VAMPNet, enabling a more detailed kinetic study of the biomolecular processes. Our GraphVAMPNet approach is also enhanced with an attention mechanism to find the important residues for classification into different metastable states. [ABSTRACT FROM AUTHOR]

Hybrid methods that combine molecular dynamics methods capable of analyzing dynamics with Monte Carlo (MC) methods that can efficiently treat thermodynamically stable states are valuable for understanding complex chemical processes in which an equilibrium state is reached through many elementary processes. The hybrid MC (HMC) method is one such promising method; however, it often fails to sample configurations properly from the canonical multimodal distribution due to the rugged potential energy surfaces. In this paper, we extend the HMC method to overcome this difficulty. The new method, which is termed potential scaling HMC (PS-HMC), makes use of an artificially modulated trajectory to propose a new configuration. The trajectory is generated from Hamilton's equations, but the potential energy surface is scaled to be gradually flattened and then recovered to the original surface, which facilitates barrier-crossing processes. We apply the PS-HMC method to three kinds of molecular processes: the thermal motion of argon particles, butane isomerization, and an atom transfer chemical reaction. These applications demonstrate that the PS-HMC method is capable of correctly constructing the canonical ensemble with a multimodal distribution. The sampling efficiency and accepted trajectories are examined to clarify the features of the PS-HMC method. Despite the potential scaling, many reactive atom transfer trajectories (elementary processes) pass through the vicinity of the minimum energy path. Furthermore, we demonstrate that the method can properly imitate the relaxation process owing to the inherent configurational continuity. By comparing the PS-HMC method with other relevant methods, we can conclude that the new method is a unique approach for studying both the dynamic and thermodynamic aspects of chemical processes. [ABSTRACT FROM AUTHOR]

The remarkable success of x-ray free-electron lasers and their ability to image biological macromolecules while outrunning secondary radiation damage due to photoelectrons, by using femtosecond pulses, raise the question of whether this can be done using pulsed high-energy electron beams. In this paper, we use excited state molecular dynamics simulations, with tabulated potentials, for rare gas solids to investigate the effect of radiation damage due to inelastic scattering (by plasmons, excitons, and heat) on the pair distribution function. We use electron energy loss spectra to characterize the electronic excitations responsible for radiation damage. [ABSTRACT FROM AUTHOR]

The mechanism of vibrational relaxation and dissociation in the O2–O system at elevated temperatures is investigated by means of molecular dynamics. The most recent O3 potential energy surfaces (PESs), obtained from the first principles quantum mechanical calculations [Varga et al., J. Chem. Phys. 147, 154312 (2017)], are used to derive a complete set of state-specific rate coefficients of vibrational energy transfer and dissociation. Unlike most of the previous efforts that utilize only the lowest and supposedly most reactive 11A′ O3 PES [A. Varandas and A. Pais, Mol. Phys. 65, 843 (1988)], this paper demonstrates the necessity to account for a complete ensemble of all excited O3 PESs that correlate with O2(X) and O(3P) when high-temperature kinetics is of interest. At the same time, it is found that the Varandas 11A′ O3 PES adequately describes vibrational energy transfer and dissociating dynamics when compared to the most recent 11A′ O3 PES by Varga et al. [J. Chem Phys. 147, 154312 (2017)]. The differences between this new dataset and previous rate coefficients are quantified by the master equation model. [ABSTRACT FROM AUTHOR]

In the context of molecular dynamics simulations, alchemical approaches based on nonequilibrium techniques are recently emerging as a powerful method for the computation of solvation free energy of druglike compounds. Here, we present a rigorous and extensive analysis of the accuracy and precision of the method as a function of the parameters qualifying the nonequilibrium alchemical protocol (e.g., number and length of the nonequilibrium trajectories and decoupling or recoupling alchemical schedule) on a selection of drug-size organic compounds characterized by a nontrivial conformational free energy landscape. The study is organized in two contributions. The first paper includes a detailed description of method and of the conformational behavior of molecular systems. Results are focused on the accuracy and precision bidirectional estimates of solvation free energy, notably those based on the so-called Bennett acceptance ratio. In the second paper, unidirectional estimates for solvation free energy are analyzed in depth. [ABSTRACT FROM AUTHOR]

The hydration structure of cellulose is very important for understanding the hydrolysis of cellulose at the molecular level. In this paper, we report a joint experimental and theoretical study on x-ray absorption spectroscopy (XAS) of aqueous cellobiose, a disaccharide unit of cellulose. In the experimental part, high resolution measurements of the carbon K-edge XAS spectra were taken. In the theoretical part, ab initio molecular dynamics simulations and ensemble calculations of electronic excited states were performed to obtain the continuous XAS spectra. The XAS spectra were found to have three characteristic peaks at 289.3, 290.7, and 293.6 eV, each representing the absorption by carbon atoms of the alcohol group, the hemiacetal group, and both of these functional groups. It was found that the peak heights in the spectrum change considerably over the temperature range of 25–60 °C, which is a reflection of the number of hydrogen bonds between cellobiose and water. We suggest that this spectral change could be useful information for identifying the hydration of cellulose in various environments. [ABSTRACT FROM AUTHOR]

We present an investigation by molecular dynamics (MD)-simulations of the coherent dynamic structure factor, S(Q, t) (Q: momentum transfer), of liquid water at the mesoscale (0.1 Å−1 ≤ Q ≤ Qmax) [Qmax ≈ 2 Å−1: Q-value of the first maximum of the static structure factor, S(Q), of water]. The simulation cell—large enough to address the collective properties at the mesoscale—is validated by direct comparison with recent results on the dynamic structure factor in the frequency domain obtained by neutron spectroscopy with polarization analysis [Arbe et al., Phys. Rev. Res. 2, 022015 (2020)]. We have not only focused on the acoustic excitations but also on the relaxational contributions to S(Q, t). The analysis of the MD-simulation results—including the self- and distinct contributions to the diffusive part of S(Q, t)—nicely explains why the relaxation process hardly depends on Q in the low Q-range (Q ≤ 0.4 Å−1) and how it crosses over to a diffusion-driven process at Q ≈ Qmax. Our simulations also give support to the main assumptions of the model used to fit the experimental data in the above mentioned paper. The application of such a model to the simulation S(Q, t) data delivers (i) results for the relaxation component of S(Q, t) in agreement with those obtained from neutron experiments and (ii) longitudinal and transverse hydrodynamic-like components with similar features than those identified in previous simulations of the longitudinal and transverse current spectra directly. On the other hand, in general, our MD-simulations results of S(Q, t) qualitatively agree with the viscoelastic transition framework habitually used to describe inelastic x-ray scattering results. [ABSTRACT FROM AUTHOR]

The development of highly efficient methods for the calculation of electronic coupling matrix elements between the electron donor and acceptor is an important goal in theoretical organic semiconductor research. In Paper I [F. Gajdos, S. Valner, F. Hoffmann, J. Spencer, M. Breuer, A. Kubas, M. Dupuis, and J. Blumberger, J. Chem. Theory Comput. 10, 4653 (2014)], we introduced the analytic overlap method (AOM) for this purpose, which is an ultrafast electronic coupling estimator parameterized to and orders of magnitude faster than density functional theory (DFT) calculations at a reasonably small loss in accuracy. In this work, we reparameterize and extend the AOM to molecules containing nitrogen, oxygen, fluorine, and sulfur heteroatoms using 921 dimer configurations from the recently introduced HAB79 dataset. We find again a very good linear correlation between the frontier orbital overlap, calculated ultrafast in an optimized minimum Slater basis, and DFT reference electronic couplings. The new parameterization scheme is shown to be transferable to sulfur-containing polyaromatic hydrocarbons in experimentally resolved dimeric configurations. Our extension of the AOM enables high-throughput screening of very large databases of chemically diverse organic crystal structures and the application of computationally intense non-adiabatic molecular dynamics methods to charge transport in state-of-the-art organic semiconductors, e.g., non-fullerene acceptors. [ABSTRACT FROM AUTHOR]

Semiclassical (SC) vibrational spectroscopy is a technique capable of reproducing quantum effects (such as zero-point energies, quantum resonances, and anharmonic overtones) from classical dynamics runs even in the case of very large dimensional systems. In a previous study [Conte et al. J. Chem. Phys. 151, 214107 (2019)], a preliminary sampling based on adiabatic switching has been shown to be able to improve the precision and accuracy of semiclassical results for challenging model potentials and small molecular systems. In this paper, we investigate the possibility to extend the technique to larger (bio)molecular systems whose dynamics must be integrated by means of ab initio "on-the-fly" calculations. After some preliminary tests on small molecules, we obtain the vibrational frequencies of glycine improving on pre-existing SC calculations. Finally, the new approach is applied to 17-atom proline, an amino acid characterized by a strong intramolecular hydrogen bond. [ABSTRACT FROM AUTHOR]

In liquid water, there is a natural tendency to form aggregates that consist of water molecules linked by hydrogen bonds. Such spontaneously formed aggregates are surrounded by a "sea" of disordered water molecules, with both forms remaining in equilibrium. The process of creating water aggregates also takes place in the solvation water of proteins, but in this case, the interactions of water molecules with the protein surface shift the equilibrium of the process. In this paper, we analyze the structural properties of the solvation water in antifreeze proteins (AFPs). The results of molecular dynamics analysis with the use of various parameters related to the structure of solvation water on the protein surface are presented. We found that in the vicinity of the active region responsible for the binding of AFPs to ice, the equilibrium is clearly shifted toward the formation of "ice-like aggregates," and the solvation water has a more ordered ice-like structure. We have demonstrated that a reduction in the tendency to create "ice-like aggregates" results in a significant reduction in the antifreeze activity of the protein. We conclude that shifting the equilibrium in favor of the formation of "ice-like aggregates" in the solvation water in the active region is a prerequisite for the biological functionality of AFPs, at least for AFPs having a well-defined ice binding area. In addition, our results fully confirm the validity of the "anchored clathrate water" concept, formulated by Garnham et al. [Proc. Natl. Acad. Sci. U. S. A. 108, 7363 (2011)]. [ABSTRACT FROM AUTHOR]

In this paper, we study the permeation of polyatomic gas molecules through 2D graphene membranes. Using equilibrium molecular dynamics simulations, we investigate the permeation of pure gas compounds (CH4, CO2, O2, N2, and H2) through nanoporous graphene membranes with varying pore sizes and geometries. Our simulations consider the recrossing mechanism, often neglected in previous studies, which has a significant effect on permeation for intermediate pore size to molecular diameter ratios. We find that the permeation process can be decoupled into two steps: the crossing process of gas molecules through the pore plane and the escaping process from the pore region to a neighboring adsorption site, which prevents recrossing. To account for these mechanisms, we use a permeance model expressed as the product of the permeance for the crossing process and the probability of molecule escape. This phenomenological model is extended to account for small polyatomic gas molecules and to describe permeation regimes ranging from molecular sieving to effusion. The proposed model captures the temperature dependence and provides insights into the key parameters of the gas/membrane interaction controlling the permeance of the system. This work lays the foundation for predicting gas permeance and exploring membrane separation factors in 2D materials such as graphene. [ABSTRACT FROM AUTHOR]

This paper presents the molecular dynamics simulation of 1 mol kg−1 LiPF6 in a binary solvent of ethylene carbonate (EC) and dimethylcarbonate, which is a representative electrolyte solution for lithium-ion batteries. The simulation successfully reproduced the diffusion coefficient, ionic conductivity, and shear viscosity as functions of EC content at 300 K, which had been experimentally determined in our previous study. The Yukawa potential was adopted to model intercharge interactions to reduce computational costs, which consequently allowed us to precisely calculate the conductivity and viscosity by directly integrating time-correlation functions without explicitly modeling the molecular polarization. Breaking down microscopic current correlation functions into components revealed that, whereas the cation–anion attractive interaction dominantly impedes the conduction when the EC content is low, it is the cation–cation and anion–anion repulsive interactions that reduce the conductivity at a high EC content. An analysis of the pressure correlations revealed that all components positively contribute to the viscosity in the binary solvent without the electrolyte. On the other hand, negative terms are observed in five out of six cross correlations in the presence of the electrolyte, implying that these correlations negatively contribute to the shear stress and entropy production, both of which are net positive. [ABSTRACT FROM AUTHOR]

As biomolecules vibrate and rotate in the terahertz band, the biological effects of terahertz electromagnetic fields have drawn considerable attention from the physiological and medical communities. Ion channels are the basis of biological electrical signals, so studying the effect of terahertz electromagnetic fields on ion channels is significant. In this paper, the effect of a terahertz electromagnetic field with three different frequencies, 6, 15, and 25 THz, on the Kv1.2 potassium ion channel was investigated by molecular dynamics simulations. The results show that an electromagnetic field with a 15 THz frequency can significantly enhance the permeability of the Kv1.2 potassium ion channel, which is 1.7 times higher than without an applied electric field. By analyzing the behavior of water molecules, it is found that the electromagnetic field with the 15 THz frequency shortens the duration of frozen and relaxation processes when potassium ions pass through the channel, increases the proportion of the direct knock-on mode, and, thus, enhances the permeability of the Kv1.2 potassium ion channel. [ABSTRACT FROM AUTHOR]

The collective coupling of an ensemble of molecules to a light field is commonly described by the Tavis–Cummings model. This model includes numerous eigenstates that are optically decoupled from the optically bright polariton states. Accessing these dark states requires breaking the symmetry in the corresponding Hamiltonian. In this paper, we investigate the influence of non-unitary processes on the dark state dynamics in the molecular Tavis–Cummings model. The system is modeled with a Lindblad equation that includes pure dephasing, as it would be caused by weak interactions with an environment, and photon decay. Our simulations show that the rate of pure dephasing, as well as the number of two-level systems, has a significant influence on the dark state population. [ABSTRACT FROM AUTHOR]

The MACE architecture represents the state of the art in the field of machine learning force fields for a variety of in-domain, extrapolation, and low-data regime tasks. In this paper, we further evaluate MACE by fitting models for published benchmark datasets. We show that MACE generally outperforms alternatives for a wide range of systems, from amorphous carbon, universal materials modeling, and general small molecule organic chemistry to large molecules and liquid water. We demonstrate the capabilities of the model on tasks ranging from constrained geometry optimization to molecular dynamics simulations and find excellent performance across all tested domains. We show that MACE is very data efficient and can reproduce experimental molecular vibrational spectra when trained on as few as 50 randomly selected reference configurations. We further demonstrate that the strictly local atom-centered model is sufficient for such tasks even in the case of large molecules and weakly interacting molecular assemblies. [ABSTRACT FROM AUTHOR]

Inelastic collisions of OH with an inert liquid perfluoropolyether (PFPE) surface have been studied experimentally. A pulsed molecular beam of OH with a kinetic energy distribution peaking at 35 kJ mol−1 was directed at a continually refreshed PFPE surface. OH molecules were detected state-selectively with spatial and temporal resolution by pulsed, planar laser-induced fluorescence. The scattered speed distributions were confirmed to be strongly superthermal, regardless of the incidence angle (0° or 45°). Angular scattering distributions were measured for the first time; their reliability was confirmed through extensive Monte Carlo simulations of experimental averaging effects, described in Paper II [A. G. Knight et al., J. Chem. Phys. 158, 244705 (2023)]. The distributions depend markedly on the incidence angle and are correlated with scattered OH speed, consistent with predominantly impulsive scattering. For 45° incidence, the angular distributions are distinctly asymmetric to the specular side but peak at sub-specular angles. This, along with the breadth of the distributions, is incompatible with scattering from a surface that is flat on a molecular scale. New molecular dynamics simulations corroborate the roughness of the PFPE surface. A subtle but unexpected systematic dependence of the angular distribution on the OH rotational state was found, which may be dynamical in origin. The OH angular distributions are similar to those for kinematically similar Ne scattering from PFPE and hence not strongly perturbed by OH being a linear rotor. The results here are broadly compatible with prior predictions from independent quasiclassical trajectory simulations of OH scattering from a model-fluorinated self-assembled monolayer surface. [ABSTRACT FROM AUTHOR]

The osmotic pressure of dilute electrolyte solutions containing charged macro-ions as well as counterions can be computed directly from the particle distribution via the well-known cell model. Originally derived within the Poisson–Boltzmann mean-field approximation, the cell model considers a single macro-ion centered into a cell, together with counterions needed to neutralize the total cell charge, while it neglects the phenomena due to macro-ion correlations. While extensively applied in coarse-grained Monte Carlo (MC) simulations of continuum solvent systems, the cell model, in its original formulation, neglects the macro-ion shape anisotropy and details of the surface charge distribution. In this paper, by comparing one-body and two-body coarse-grained MC simulations, we first establish an upper limit for the assumption of neglecting correlations between macro-ions, and second, we validate the approximation of using a non-spherical macro-ion. Next, we extend the cell model to all-atom molecular dynamics simulations and show that protein concentration-dependent osmotic pressures can be obtained by confining counterions in a virtual, spherical subspace defining the protein number density. Finally, we show the possibility of using specific interaction parameters for the protein–ion and ion–ion interactions, enabling studies of protein concentration-dependent ion-specific effects using merely a single protein molecule. [ABSTRACT FROM AUTHOR]

Model Hamiltonians constructed from quantum chemistry calculations and molecular dynamics simulations are widely used for simulating nonadiabatic dynamics in the condensed phase. The most popular two-state spin-boson model could be built by mapping the all-atom anharmonic Hamiltonian onto a two-level system bilinearly coupled to a harmonic bath using the energy gap time correlation function. However, for more than two states, there lacks a general strategy to construct multi-state harmonic (MSH) models since the energy gaps between different pairs of electronic states are not entirely independent and need to be considered consistently. In this paper, we extend the previously proposed approach for building three-state harmonic models for photoinduced charge transfer to the arbitrary number of electronic states with a globally shared bath and the system–bath couplings are scaled differently according to the reorganization energies between each pair of states. We demonstrate the MSH model construction for an organic photovoltaic carotenoid–porphyrin–C60 molecular triad dissolved in explicit tetrahydrofuran solvent. Nonadiabatic dynamics was simulated using mixed quantum-classical techniques, including the linearized semiclassical and symmetrical quasiclassical dynamics with the mapping Hamiltonians, mean-field Ehrenfest, and mixed quantum-classical Liouville dynamics in two-state, three-state, and four-state harmonic models of the triad system. The MSH models are shown to provide a general and flexible framework for simulating nonadiabatic dynamics in complex systems. [ABSTRACT FROM AUTHOR]

In this paper, we present a fast and adaptive correlation guided enhanced sampling method (CORE-MD II). The CORE-MD II technique relies, in part, on partitioning of the entire pathway into short trajectories that we refer to as instances. The sampling within each instance is accelerated by adaptive path-dependent metadynamics simulations. The second part of this approach involves kinetic Monte Carlo (kMC) sampling between the different states that have been accessed during each instance. Through the combination of the partition of the total simulation into short non-equilibrium simulations and the kMC sampling, the CORE-MD II method is capable of sampling protein folding without any a priori definitions of reaction pathways and additional parameters. In the validation simulations, we applied the CORE-MD II on the dialanine peptide and the folding of two peptides: TrpCage and TrpZip2. In a comparison with long time equilibrium Molecular Dynamics (MD), 1 µs replica exchange MD (REMD), and CORE-MD I simulations, we find that the level of convergence of the CORE-MD II method is improved by a factor of 8.8, while the CORE-MD II method reaches acceleration factors of ∼120. In the CORE-MD II simulation of TrpZip2, we observe the formation of the native state in contrast to the REMD and the CORE-MD I simulations. The method is broadly applicable for MD simulations and is not restricted to simulations of protein folding or even biomolecules but also applicable to simulations of protein aggregation, protein signaling, or even materials science simulations. [ABSTRACT FROM AUTHOR]

An accurate system-bath model to investigate the quantum dynamics of hydrogen atoms chemisorbed on graphene is presented. The system comprises a hydrogen atom and the carbon atom from graphene that forms the covalent bond, and it is described by a previously developed 4D potential energy surface based on density functional theory ab initio data. The bath describes the rest of the carbon lattice and is obtained from an empirical force field through inversion of a classical equilibrium correlation function describing the hydrogen motion. By construction, model building easily accommodates improvements coming from the use of higher level electronic structure theory for the system. Further, it is well suited to a determination of the system-environment coupling by means of ab initio molecular dynamics. This paper details the system-bath modeling and shows its application to the quantum dynamics of vibrational relaxation of a chemisorbed hydrogen atom, which is here investigated at T = 0 K with the help of the multi-configuration time-dependent Hartree method. Paper II deals with the sticking dynamics. [ABSTRACT FROM AUTHOR]

Computationally inexpensive particle-based coarse-grained (CG) models are essential for use in molecular dynamics (MD) simulations of mesoscopically slow cooperative phenomena, such as plastic deformations in solids. Molecular crystals possessing complex symmetry present enormous practical challenges for particle-based coarse-graining at molecularly resolved scales, when each molecule is in a single-site representation, and beyond. Presently, there is no published pairwise non-bonded single-site CG potential that is able to predict the space group and structure of a molecular crystal. In this paper, we present a successful coarse-graining at a molecular level from first principles of an energetic crystal, hexahydro-1,3,5-trinitro-s-triazine (RDX) in the alpha phase, using the force-matching-based multiscale coarse-graining (MSCG/FM) approach. The new MSCG/FM model, which implements an optimal pair decomposition of the crystal Helmholtz free energy potential in molecular center-of-mass coordinates, was obtained by force-matching atomistic MD simulations of liquid, amorphous, and crystalline states and in a wide range of pressures (up to 20 GPa). The MSCG/FM potentials for different pressures underwent top-down optimization to fine-tune the mechanical and thermodynamic properties, followed by consolidation into a transferable density-dependent model referred to as RDX-TC-DD (RDX True-Crystal Density-Dependent). The RDX-TC-DD model predicts accurately the crystal structure of α-RDX at room conditions and reproduces the atomistic reference system under isothermal (300 K) hydrostatic compression up to 20 GPa, in particular, the Pbca symmetry of α-RDX in the elastic regime. The RDX-TC-DD model was then used to simulate the plastic response of uniaxially ([100]) compressed α-RDX resulting in nanoscale shear banding, a key mechanism for plastic deformation and defect-free detonation initiation proposed for many molecular crystalline explosives. Additionally, a comparative analysis of the effect of core-softening of the RDX-TC-DD potential and the degree of molecular rigidity in the all-atom treatment suggests a stress-induced short-range softening of the effective intermolecular interaction as a fundamental cause of plastic instability in α-RDX. The reported RDX-TC-DD model and overall workflow to develop it open up possibilities to perform high quality simulation studies of molecular energetic materials under thermal and mechanical stimuli, including extreme conditions. [ABSTRACT FROM AUTHOR]

The dynamics of a molecule in a magnetic field is significantly different from its zero-field counterpart. One important difference in the presence of a field is the Lorentz force acting on the nuclei, which can be decomposed as the sum of the bare nuclear Lorentz force and a screening force due to the electrons. This screening force is calculated from the Berry curvature and can change the dynamics qualitatively. It is therefore important to include the contributions from the Berry curvature in molecular dynamics simulations in a magnetic field. In this work, we present a scheme for calculating the Berry curvature numerically using a finite-difference technique, addressing challenges related to the arbitrary global phase of the wave function. The Berry curvature is calculated as a function of bond distance for H2 at the restricted and unrestricted Hartree–Fock levels of theory and for CH+ as a function of the magnetic field strength at the restricted Hartree–Fock level of theory. The calculations are carried out using basis sets of contracted Gaussian functions equipped with London phase factors (London orbitals) to ensure gauge-origin invariance. In this paper, we also interpret the Berry curvature in terms of atomic charges and discuss its convergence in basis sets with and without London phase factors. The calculation of the Berry curvature allows for its inclusion in ab initio molecular dynamics simulations in a magnetic field. [ABSTRACT FROM AUTHOR]

Two-dimensional electronic–vibrational spectroscopy (2DEVS) is an emerging spectroscopic technique which exploits two different frequency ranges for the excitation (visible) and detection (infrared) axes of a 2D spectrum. In contrast to degenerate 2D techniques, such as 2D electronic or 2D infrared spectroscopy, the spectral features of a 2DEV spectrum report cross correlations between fluctuating electronic and vibrational energy gaps rather than autocorrelations as in the degenerate spectroscopies. The center line slope of the spectral features reports on this cross correlation function directly and can reveal specific electronic–vibrational couplings and rapid changes in the electronic structure, for example. The involvement of the two types of transition moments, visible and infrared, makes 2DEVS very sensitive to electronic and vibronic mixing. 2DEV spectra also feature improved spectral resolution, making the method valuable for unraveling the highly congested spectra of molecular complexes. The unique features of 2DEVS are illustrated in this paper with specific examples and their origin described at an intuitive level with references to formal derivations provided. Although early in its development and far from fully explored, 2DEVS has already proven to be a valuable addition to the tool box of ultrafast nonlinear optical spectroscopy and is of promising potential in future efforts to explore the intricate connection between electronic and vibrational nuclear degrees of freedom in energy and charge transport applications. [ABSTRACT FROM AUTHOR]

Photoisomerization in the retinal leads to a channel opening in rhodopsins that triggers translocation or pumping of ions/protons. Crystal structures of rhodopsins contain several structurally conserved water molecules. It has been suggested that water plays an active role in facilitating the ion pumping/translocation process by acting as a lubricant in these systems. In this paper, we systematically investigate the localization, structure, dynamics, and energetics of the water molecules along the channel for the resting/dark state of KR2 rhodopsin. By employing several microseconds long atomistic molecular dynamics simulation of this trans-membrane protein system, we demonstrate the presence of five distinct water containing pockets/cavities separated by gateways controlled by protein side-chains. There exists a strong hydrogen bonded network involving these buried water molecules and functionally important key residues. We present evidence of significant structural and dynamical heterogeneity in the water molecules present in these cavities, with very rare exchange between them. The exchange time scale of such buried water with the bulk has an extremely wide range, from tens of nanoseconds to >1.5 µs. The translational and rotational dynamics of buried water are found to be strongly dependent on the protein cavity size and local interactions with a classic signature of trapped diffusion and rotational anisotropy. [ABSTRACT FROM AUTHOR]