Full Molecular Dynamics Simulations Research Articles

Markov state models and molecular alchemy

In recent years, Markov state models (MSMs) have attracted a considerable amount of attention with regard to modelling conformation changes and associated function of biomolecular systems. They have been used successfully, e.g. for peptides including time-resolved spectroscopic experiments, protein function and protein folding , DNA and RNA, and ligand-receptor interaction in drug design and more complicated multivalent scenarios. In this article, a novel reweighting scheme is introduced that allows to construct an MSM for certain molecular system out of an MSM for a similar system. This permits studying how molecular properties on long timescales differ between similar molecular systems without performing full molecular dynamics simulations for each system under consideration. The performance of the reweighting scheme is illustrated for simple test cases, including one where the main wells of the respective energy landscapes are located differently and an alchemical transformation of butane to pentane where the dimension of the state space is changed.

Sub-THz Phonon drag on dislocations by coarse-grained atomistic simulations

In this paper, dislocation–phonon interactions in fcc crystals are simulated using the concurrent atomistic-continuum (CAC) method. With significantly less computational cost than that by full molecular dynamics (MD) simulations, dislocation migration, phonon transport, and their interactions are explicitly modeled. Sub-THz Phonon drag coefficients on moving dislocations are predicted using CAC and compared with MD and experimental results. During the interactions between phonons and dislocations, the local temperature is found to rise due to energy dissipation, and dislocations are observed to exhibit “breathing modes” in which the separation between leading and trailing partials varies in time in a quasi-periodic fashion. Simulation results show that the phonon drag coefficient on dislocation migration increases with phonon wave packet magnitudes or sizes but is insensitive to the incident angles.

Investigation of multiscale fluid flow characteristics based on a hybrid atomistic–continuum method

A computer program based on a molecular dynamics–continuum hybrid method has been developed in which the Navier–Stokes equations are solved in the continuum region and the molecular dynamics in the atomistic region. The coupling between the atomistic and continuum is constructed through constrained dynamics within an overlap region where both molecular and continuum equations are solved simultaneously. The simulation geometries are solved in three dimensions and an overlap region is introduced in two directions to improve the choice of using the molecular region in smaller areas. The proposed method is used to simulate steady and start-up Couette flow showing quantitative agreement with results from analytical solutions and full molecular dynamics simulations. The prepared algorithm and the computer code are capable of modeling fluid flows in micro and nano-scale geometries.

Adaptive multiscale method for two-dimensional nanoscale adhesive contacts

There are two separate traditional approaches to model contact problems: continuum and atomistic theory. Continuum theory is successfully used in many domains, but when the scale of the model comes to nanometer, continuum approximation meets challenges. Atomistic theory can catch the detailed behaviors of an individual atom by using molecular dynamics (MD) or quantum mechanics, although accurately, it is usually time-consuming. A multiscale method coupled MD and finite element (FE) is presented. To mesh the FE region automatically, an adaptive method based on the strain energy gradient is introduced to the multiscale method to constitute an adaptive multiscale method. Utilizing the proposed method, adhesive contacts between a rigid cylinder and an elastic substrate are studied, and the results are compared with full MD simulations. The process of FE meshes refinement shows that adaptive multiscale method can make FE mesh generation more flexible. Comparison of the displacements of boundary atoms in the overlap region with the results from full MD simulations indicates that adaptive multiscale method can transfer displacements effectively. Displacements of atoms and FE nodes on the center line of the multiscale model agree well with that of atoms in full MD simulations, which shows the continuity in the overlap region. Furthermore, the Von Mises stress contours and contact force distributions in the contact region are almost same as full MD simulations. The method presented combines multiscale method and adaptive technique, and can provide a more effective way to multiscale method and to the investigation on nanoscale contact problems.

A hybrid molecular-continuum simulation method for incompressible flows in micro/nanofluidic networks

We present a hybrid molecular-continuum simulation method for modelling nano- and micro-flows in network-type systems. In these types of problem, a full molecular dynamics (MD) description of the macroscopic flow behaviour would be computationally intractable, or at least too expensive to be practical for engineering design purposes. Systems that exhibit multiscale traits, such as this, can instead be solved using a hybrid approach that distinguishes the problem into macroscopic and microscopic dynamics, modelled by their respective solvers. The technique presented in this study is an extension and addition to a hybrid method developed by Borg et al. (J Comput Phys 233:400–413, 2013) for high-aspect-ratio channel geometries, known as the internal-flow multiscale method (IMM). Computational savings are obtained by replacing long channels in the network, which are highly scale-separated, by much smaller, but representative, MD simulations, without a substantial loss of accuracy. On the other hand, junction components do not exhibit this length-scale separation, and so must be simulated in their entirety using MD. The current technique combines all network elements (junctions and channels) together in a coupled simulation using continuum conservation laws. For the case of steady, isothermal, incompressible, low-speed flows, we use the conservation of mass and momentum flux equations to derive a set of molecular-continuum constraints. An algorithm is presented here that computes at each iteration the new constraints on the pressure differences to be applied over individual MD micro-elements (channels and junctions), successively moving closer to macroscopic mass and momentum conservation. We show that hybrid simulations of some example network cases converge quickly, in only a few iterations, and compare very well to the corresponding full MD results, which are taken as the most accurate solutions. Major computational savings can be afforded by the IMM-type approximation in the channel components, but for steady-state solutions, even greater savings are possible. This is because the micro-elements are coupled to a steady-state continuum conservation expression, which greatly speeds up the relaxation of individual micro-components to steady conditions as compared to that of a full MD simulation. Unsteady problems with high temporal scale separation can also be simulated, but general transient problems are beyond the capabilities of the current technique.

Aqueous Guanidinium–Carbonate Interactions by Molecular Dynamics and Neutron Scattering: Relevance to Ion–Protein Interactions

Guanidinium carbonate was used in this study as a simple proxy for the biologically relevant arginine-carbonate interactions in water. Molecular dynamics (MD) simulations of guanidinium carbonate were performed with nonpolarizible water using two implementations of the ion force fields. In the first, the ions had full charges, while in the second, the ions had reduced charges in order to effectively account for electronic polarization effects of water. The results from the simulations were then compared to data from previous neutron scattering experiments. It was found that there were significant discrepancies between the full charge force field MD simulations and the experimental results due to excessive ion pairing and clustering in the former. In contrast, reducing the ionic charges yields a more regular solution with a simulated structure, which fits well the experimental data.

A multiscale method for micro/nano flows of high aspect ratio

We develop a new multiscale scheme for simulating micro/nano flows of high aspect ratio in the flow direction, e.g. within long ducts, tubes, or channels, of varying section. The scheme couples conventional hydrodynamic conservation equations for mass and momentum-flux with molecular dynamics (MD) in a unified framework. The method is very much different from common ‘domain-decomposition’ hybrid methods, and is more related to micro-resolution methods, such as the Heterogeneous Multiscale Method. We optimise the use of the computationally-costly MD solvers by applying them only at a limited number of streamwise-distributed cross-sections of the macroscale geometry. The greater the streamwise scale of the geometry, the more significant is the computational speed-up when compared to a full MD simulation. We test our new multiscale method on the case of a converging/diverging nanochannel conveying a simple Lennard–Jones liquid. We validate the results from our simulations by comparing them to a full MD simulation of the same test case.

Coarse-grained atomistic modeling and simulation of inelastic material behavior

This paper presents a new methodology for coarse-grained atomistic simulation of inelastic material behavior including phase transformations in ceramics and dislocation mediated plasticity in metals. The methodology combines an atomistic formulation of balance equations and a modified finite element method. With significantly fewer degrees of freedom than those of a fully atomistic model and without additional constitutive rules but the interatomic force field, the new coarse-grained (CG) method is shown to be feasible in predicting the nonlinear constitutive responses of materials and also reproducing atomic-scale phenomena such as phase transformations (diamond → β-Sn) in silicon and dislocation nucleation and migration, formation of dislocation loops and stacking faults ribbons in single crystal nickel. Direct comparisons between CG and the corresponding full molecular dynamics (MD) simulations show that the present methodology is efficient and promising in modeling and simulation of inelastic material behavior without losing the essential atomistic features. The potential applications and the limitations of the CG method are also discussed.

Molecular Modelling of G Protein-Coupled Receptors Through the Web.

With the recent crystallization of several G Protein-Coupled receptors (GPCRs), homology modelling and all atom molecular dynamics (MD) simulations have proven their usefulness for exploring the structure and function of this superfamily of membrane receptors. Subsequently, automated computational protocols have been implemented as web-based servers in the recent years to produce reliable models of GPCRs, providing partial or global solutions for the structural characterization and molecular simulation of GPCRs. These dedicated modelling services represent an attractive tool for the broader community of public researchers and pharmaceutical companies, in order to assist in the structure-based drug design of GPCRs. We here collect and analyze the existing web servers, among which a previously unreported service, GPCR-ModSim, offers for the first time full atom MD simulations in the pipeline for GPCR molecular modelling.

Multiscale simulations of nanoindentation and nanoscratch of single crystal copper

In this paper, a three-dimensional multiscale method coupling molecular dynamics (MD) and finite element (FE) analysis has been employed to study nanoindentation and nanoscratch processes of copper thin film. During the nanoindentation process, no discontinuities can be observed in the displacement gradient, which demonstrates a seamless coupling between MD and FE regions. The multiscale model of nanoindentation is further validated by the fact that the load–depth curves obtained from both multiscale and full MD simulations match each other reasonably well. The nanoindentation results using multiscale model are also compared with those using part MD model to elucidate the advantages of multiscale method. The multiscale simulation results of nanoindentation show that the indentation velocity has a great effect on the load–depth curve and hardness of copper. The continuity of displacement field during the nanoscratch also demonstrates a seamless information exchange between MD and FE regions. The deformation of substrate and forces in the nanoscratch process are also investigated using the multiscale method. The multiscale method can greatly expand the system size to be explored without significant increase in computational cost.

A Molecular Dynamics-Continuum Coupled Model for Heat Transfer in Composite Materials

The heat transfer problem in composite materials containing nanoscale interfaces is investigated. A molecular dynamics-continuum coupled model is developed to study heat transport from the macro- to the microscales. The model includes four major steps: (1) A reverse nonequilibrium molecular dynamics (RNEMD) method is used to calculate some physical parameters such as the thermal conductivities on the interface. (2) The homogenization method is applied to compute the homogenized thermal conductivities of composite materials. (3) The temperature field in the global structure of composite materials is computed with the multiscale asymptotic method for the macroscopic heat transfer equation. (4) A molecular dynamics-continuum coupled model is developed to reevaluate the temperature field of composite materials, in particular, the local temperature field near the interface. Numerical results in one-, two-, and three-dimensional structures of composite materials including the nanoscale interface are given. Good agreement is achieved between the numerical results of the proposed coupled algorithm and those of the full molecular dynamics simulation, demonstrating the accuracy of the present method and its potential applications in thermal engineering of composite materials.

Testing the semi-explicit assembly solvation model in the SAMPL3 community blind test

We report here a test of the Semi-Explicit Assembly (SEA) model in the solvation free energy category of the SAMPL3 blind prediction event (summer 2011). We tested how dependent the SEA results are on the chosen force field by performing calculations with both the General Amber and OPLS force fields. We compared our SEA results with full molecular dynamics simulations in explicit solvent. Of the 20 submissions, our SEA/OPLS results gave the second smallest RMS errors in free energies compared to experiments. SEA gives results that are very similar to those of its underlying force field and explicit solvent model. Hence, while the SEA water modeling approach is much faster than explicit solvent simulations, its predictions appear to be just as accurate.

Static substrate deposition: Toward longer time scale deposition simulations

We report on the development of a new method for simulating the deposition of atoms onto surfaces, which allows for longer time scales and/or larger system sizes. The method involves simulating only the deposited material on a purely static substrate and corrects for the missing thermal fluctuations and energy exchange effects that surface atoms experience as a result of a static substrate. This method may be viewed as an intermediate approximation between the discrete kinetic Monte Carlo and the full molecular dynamics (MD) methods, because it retains the full trajectory dynamics for the deposited atoms. The method can achieve speed-up times of approximately two orders of magnitude for submonolayer depositions on large simulation substrates. We first apply the method to an ideal prototypical system involving the deposition of Pt material on Pt substrates with two-dimensional Ag monolayer substrate patterns, where we show that it produces nanostructures and small cluster diffusion behavior almost identical to that observed during full MD simulations. As an ideal second application of the method, we use it to deposit Pt and Au onto a reconstructed Au(111) surface representative of a herringbone reconstruction and comment on the findings. Finally, we comment on the limitations and possible future improvements.

A hybrid model for multiscale simulations of nanoindentation

A numerical model for multiscale simulations of nanoindentation has been established by coupling molecular dynamics (MD) simulations with finite element analysis. The multiscale simulations are performed for two-dimensional nanoindentation of a cylindrical indenter against a copper (1 1 1) substrate. Results show that the load–depth curves obtained from both multiscale and full MD simulations match each other reasonably well and no visible discontinuity of deformation is observed in the handshake region. In the beginning, the normal load increases gradually, corresponding to elastic substrate deformation, followed by sudden drops in load–depth curves, which are related to the nucleation and movement of dislocations. The simulations also reveal that indenter radius and indentation velocity significantly affect the nanoindentation process. By the use of multiscale method, the system size to be explored can be greatly expanded without increasing computational cost.

Roughness effect on flow and thermal boundaries in microchannel/nanochannel flow using molecular dynamics‐continuum hybrid simulation

SUMMARYOur molecular dynamics‐continuum hybrid simulation method is further validated in terms of flexibility in domain decomposition. The roughness effects on the flow and thermal boundaries in liquid channel flow is studied. The results indicate that the molecules in the wall‐neighboring area can be firmly confined in the concaves due to geometric structure and strong liquid–solid interaction and cause locking boundary in the velocity profile and linear gradient in the temperature profile. The locking boundary can further lead to negative slip length, which varies in power law with channel height. The linear temperature gradient, as well as nearly constant temperature jump, can lead to obviously increasing Kapitza length versus channel height. Analysis on flow friction shows that the confinement on the molecules will equivalently narrow the channel, where larger pressure gradient is needed. Therefore, the roughness should be strictly restricted within a range shrinking correspondingly with channel height if the flow condition is supposed to be maintained unchanged. Finally, our hybrid simulation is compared with full molecular dynamics simulation in terms of computational efficiency. Great advantages of the hybrid simulation, such as exclusively flexibility and combination characteristics, demonstrate its potential values and promising applications in the field of microfluidics/nanofluidics. Copyright © 2011 John Wiley & Sons, Ltd.

Rigidity analysis of HIV-1 protease

We present a rigidity analysis on a large number of X-ray crystal structures of the enzyme HIV-1 protease using the 'pebble game' algorithm of the software FIRST. We find that although the rigidity profile remains similar across a comprehensive set of high resolution structures, the profile changes significantly in the presence of an inhibitor. Our study shows that the action of the inhibitors is to restrict the flexibility of the β-hairpin flaps which allow access to the active site. The results are discussed in the context of full molecular dynamics simulations as well as data from NMR experiments.

Does ℏ Play a Role in Multidimensional Spectroscopy? Reduced Hierarchy Equations of Motion Approach to Molecular Vibrations

To investigate the role of quantum effects in vibrational spectroscopies, we have carried out numerically exact calculations of linear and nonlinear response functions for an anharmonic potential system nonlinearly coupled to a harmonic oscillator bath. Although one cannot carry out the quantum calculations of the response functions with full molecular dynamics (MD) simulations for a realistic system which consists of many molecules, it is possible to grasp the essence of the quantum effects on the vibrational spectra by employing a model Hamiltonian that describes an intra- or intermolecular vibrational motion in a condensed phase. The present model fully includes vibrational relaxation, while the stochastic model often used to simulate infrared spectra does not. We have employed the reduced quantum hierarchy equations of motion approach in the Wigner space representation to deal with nonperturbative, non-Markovian, and nonsecular system-bath interactions. Taking the classical limit of the hierarchy equations of motion, we have obtained the classical equations of motion that describe the classical dynamics under the same physical conditions as in the quantum case. By comparing the classical and quantum mechanically calculated linear and multidimensional spectra, we found that the profiles of spectra for a fast modulation case were similar, but different for a slow modulation case. In both the classical and quantum cases, we identified the resonant oscillation peak in the spectra, but the quantum peak shifted to the red compared with the classical one if the potential is anharmonic. The prominent quantum effect is the 1-2 transition peak, which appears only in the quantum mechanically calculated spectra as a result of anharmonicity in the potential or nonlinearity of the system-bath coupling. While the contribution of the 1-2 transition is negligible in the fast modulation case, it becomes important in the slow modulation case as long as the amplitude of the frequency fluctuation is small. Thus, we observed a distinct difference between the classical and quantum mechanically calculated multidimensional spectra in the slow modulation case where spectral diffusion plays a role. This fact indicates that one may not reproduce the experimentally obtained multidimensional spectrum for high-frequency vibrational modes based on classical molecular dynamics simulations if the modulation that arises from surrounding molecules is weak and slow. A practical way to overcome the difference between the classical and quantum simulations was discussed.

Ab initio electron correlated studies on the intracluster reaction of NO+(H2O)n → H3O+(H2O)n−2 (HONO) (n = 4 and 5)

Hydrated nitrosonium ion clusters NO(+)(H(2)O)(n) (n = 4 and 5) were investigated by using MP2/aug-cc-pVTZ level of theory to clarify isomeric reaction pathways for formation of HONO and fully hydrated hydride ions. We found some new isomers and transition state structures in each hydration number, whose lowest activation energies of the intracluster reactions were found to be 4.1 and 3.4 kcal mol(-1) for n = 4 and n = 5, respectively. These thermodynamic properties and full quantum mechanical molecular dynamics simulation suggest that product isomers with HONO and fully hydrated hydride ions can be obtained at n = 4 and n = 5 in terms of excess hydration binding energies which can overcome these activation barriers.

The solvent shell structure of aqueous iodide: X-ray absorption spectroscopy and classical, hybrid QM/MM and full quantum molecular dynamics simulations

The L 3 X-ray absorption spectrum of aqueous iodide is reported, and its EXAFS is compared to theoretical spectra reconstructed from the radial distribution function of the iodide hydration obtained from classical, hybrid Quantum Mechanics Molecular Mechanics, (QM/MM) and full quantum (density functional theory, DFT) molecular dynamics simulations. Since EXAFS is mainly sensitive to short distances around the iodide ion, it is a direct probe of the local solvation structure. The comparison shows that QM/MM simulations deliver a satisfactory description of the EXAFS signal, while nonpolarizable classical simulations are somewhat less satisfactory and DFT-based simulations perform poorly. We also identify a weak anisotropy of the water solvation shell around iodide, which may be of importance in electron photoejection experiments.

Effective pair potentials using anab initiovariational approach

We used a variational approach adapted to a quantum molecular-dynamics code to determine the best reference potential for warm dense aluminum. This ab initio variational approach was based on the Gibbs-Bogolyubov inequality. We used many-body reference systems interacting through inverse-power-law potentials, among which the Coulomb potential was a particular case defining the classical one-component plasma model. By comparisons with full quantum molecular-dynamics simulations, we found that the Coulomb potential was not always the best reference potential. We calculated the self-diffusion coefficient and the shear viscosity and discussed the results obtained using the Chisolm-Wallace relation in the warm dense matter regime.