Molecular dynamics | EPFL Graph Search

Molecular dynamics (MD) is a computer simulation method for analyzing the physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic "evolution" of the system. In the most common version, the trajectories of atoms and molecules are determined by numerically solving Newton's equations of motion for a system of interacting particles, where forces between the particles and their potential energies are often calculated using interatomic potentials or molecular mechanical force fields. The method is applied mostly in chemical physics, materials science, and biophysics. Because molecular systems typically consist of a vast number of particles, it is impossible to determine the properties of such complex systems analytically; MD simulation circumvents this problem by using numerical methods. However, long MD simulations are mathematically ill-conditioned, generating cumulative errors in numerical int

Molecular Dynamics Simulations of Helium Atoms Clustering in bcc Iron

Ning Gao

Helium atoms are known to have a significant impact on materials used in fission and fusion reactors. In particular, the presence of helium atoms can change the mechanical properties and degrade the lifetime of reactors. In order to develop the helium-resistance materials, the underlying interactions between helium atoms and matrix atoms have to be understood. In the past, atomistic simulations have contributed a lot to the basic energy states of helium atoms in different defects of materials. Especially, insights into the formation energies of helium atoms at simple defects have been calculated by ab initio and molecular dynamics. The substitional helium has lower formation energy than the tetrahedral and octahedral interstitial helium atom. Additionally, the binding energies of helium with small HenVm cluster have also been found to be all positive. However, all these simulations have been performed without showing the kinetic formation process of such helium defects. The interactions of defect with larger number of helium atoms are also not explored in detail. In the first part of this thesis, a new relaxation-fitting method has been developed for cross potential of Fe-He system based on more energy data of helium-defect clusters from ab initio calculations and by including the migration energy of single free helium atom in the fitting process. Besides, the recently fitted magnetic potentials of Fe have been used to describe the Fe-Fe interactions for such cross potential. The new fitted Fe-He potentials have predicted well not only the properties of simple helium defect, but also the properties of helium clusters and diffusion process of helium atoms. Based on the new Fe-He potentials, the interactions between helium atoms and Fe with the effect of single vacancy have been simulated utilizing the classical molecular dynamics (MD) method in the single bcc Fe lattice at room temperatures by helium atoms up to 18. The binding energies of helium atom with HenV cluster are computed. By increasing number of helium atoms, the binding energy decreases first, then slowly increases up to n = 4. It is then almost constant between n = 5 to 15 but increases strongly at n = 16 to form a peak. This last increase is related to the athermal self-interstitial atom (SIA) emission in the form of dumbbell which relaxes the restrained system to He16V2 cluster. The second binding energy peak occurs at n = 18 by the formation of the second SIA to form He18V3 cluster. Close inspection of the local configuration shows the formed SIA combines with the helium-defect cluster. Calculation of the pressure of the helium-defect cluster shows local peak normal stress and shear stress values up to 9 GPa and 4 GPa, respectively. The local configurations of helium-defect cluster suggest that with increasing helium content, some symmetrical structures can be formed. By increasing the number of helium atoms inserted in single bcc Fe lattice, the clustering process of the formed SIAs with presence of helium cluster has been simulated with MD methods at room temperatures with the new fitted Fe-He potential. The first SIA is formed after 6 He atoms being inserted without vacancy. The simulations indicate that the SIA has mainly two configurations as either a dumbbell or a crowdion, which can be exchanged in the interface between helium bubble and matrix atoms. The SIA cluster can change direction and absorb more SIAs to grow. Small dislocation loops have been observed to form with Burgers vector b = 1/2 a0 on a {111} habit plane. The kinetic process of formation and growth of dislocation loop explored by MD simulations includes: (1) the formation of a self-interstitial atom by the high pressure of the helium cluster/bubble; (2) the diffusion of these self-interstitial atoms on the interface between helium bubble and matrix; (3) the clustering of the self-interstitial atoms in the form of a dislocation loop and its diffusion away from bubble when the number of SIAs in loop reaches a critical number. Moreover, the effect of grain boundary (GB) on the formation of helium clusters has been explored with MD simulations in symmetrical GBs and normal nano-grain boundary samples. Helium atoms can be trapped by GBs at free spaces. The local atomic excess free volume (LAEFV) is found to the one of the main parameters not only related with the energy state of GB but also with the nucleation and growth of He atoms in the GB.

EPFL2011

Solute energy based REMD

Molecular dynamics (MD) simulations have increasingly contributed to the understanding of biomolecular processes, allowing for predictions of thermodynamic and structural properties. Unfortunately, the holy grail of protein structure prediction was soon found to be severely hampered by the very rugged free energy surface of proteins, with small relative free energies separating native, folded protein conformations from unfolded states. These multiple minima frequently trap present-day protein MD simulations permanently. In order to allow the simulation to escape minima and explore wider portions of conformational space, enhanced sampling techniques were developed. One of the most popular ones, replica exchange molecular dynamics (REMD), is based on multiple parallel MD simulations that are performed with replicas of a system at increasing temperatures T1, T2, etc. Periodic Monte Carlo exchange moves are attempted, aiming to allow conformations to exchange temperature ensembles with a probability that depends on their potential energy and temperature difference. Thus, conformations are simulated at all temperatures and escape local minima with the kinetic energy provided at higher temperatures, while Boltzmann distributions are generated at all temperatures. REMD has been successful in ab-initio folding of a variety of small peptides and proteins (up to 20-30 residues). However, with larger proteins, the overlap of potential energy distributions diminishes, since the potential energy and its fluctuation scale with fkBT, respectively with √fkBT, where f is the number of degrees of freedom of the system, kB Boltzmann's constant and T the temperature. Consequently, the related Monte Carlo exchange probability and number of exchanges in a simulation are also diminished. This is generally compensated by choosing smaller temperature intervals between replicas and thereby increasing the number of necessary replicas (as well as the computational cost of the simulation) to cover a given temperature range. An additional problem relates to explicit solvent simulations, in which solvent to solvent interactions account for the largest part of the total potential energy. Consequently, explicit solvent REMD simulations almost exclusively sample solvent degrees of freedom. These two limitations have lead to the development of REMD protocols for large explicit solvent systems that are based on exchange probabilities computed with subsystem (e.g. protein only) potential energy functions, allowing for a targeted sampling of protein degrees of freedom and a reduction of the computational effort. In this thesis, this approximation is tested by implementing its simplest variation that entirely neglects solvent-solvent interaction as a new REMD protocol termed REM Dpe (Chapter 2). Possible REMD limitations for large explicit solvent systems are tested (Chapter 3) with REM Dpe, which is further applied to perform a thorough and comparative investigation of prion (Chapter 4) and doppel (Chapter 5) protein misfolding. In Chapter 2, the practical validity of the REM Dpe approximation is assessed with simulations of the prion protein, a system that is too large (i.e. 103 residues) to allow for efficient simulations using traditional REMD over the necessary temperature range. A first validation consists in testing whether protein and total potential energy distributions are consistent with their analogs from straightforward reference MD simulations. Second, the overlap pattern of the total potential energy distributions are characterized at different temperatures and show that the exchanges in the REM Dpe simulations are performed according to a Boltzmann weight. Native structures are found to have the lowest protein and total potential energies, as compared to higher energies found for various unfolded structures. Although no obvious bias is detected in the three validations, the conformational landscapes of the REM Dpe simulation at low temperatures progressively shift to non-native regions of the free energy surface. In Chapter 3, this phenomenon is quantified, and its origin identified in insufficient low temperature residence times required for refolding native-like structures. REMD is based on the assumption that systems have to be decorrelated between exchange attempts. Increasing inter-exchange times accordingly would allow for decorrelation and sufficient low temperature residence times but is practically impossible to achieve for large protein simulations, highlighting a major limitation and possible source of bias for present-day REMD simulations. We have chosen the prion as a test case because of its link to transmissible spongiform encephalopathies. Diseases of this category are believed to be caused by a rare prion protein (PrP) misfold leading from the cellular, monomeric, soluble, α-helical PrPC isoform to a pathogenic, aggregated, insoluble, β-rich PrPSc isoform of unknown structure. Gaining experimental knowledge of the PrPSc structure has remained elusive, and aroused interest in predictions supplied by computer simulations. REMD provides a powerful tool allowing to explore a diversity of misfolds and select stable ones that accumulate at lower temperatures. In Chapter 4, we describe a PrP REM Dpe simulation in which rare new β-strands are formed and arrange into a multitude of different β-sheets, reproducing the α-helix → β-sheet conversion observed with circular dichroism spectra. The α-helical and β-sheet propensities along the sequence can thus be computed. We develop and apply the β contact map clustering (bcmc) protocol to identify the most frequent β-sheet pattern defining β-rich folds. 10 new β-rich folds are found and compared to recent experimental data characterizing PrPSc, providing atomistically detailed models for putative monomeric precursors of PrPSc or β-oligomeric conformations. In Chapter 5, an analogous simulation is performed with doppel, a structural homolog of prion (with an identical three α-helix, two β-strand fold) originating from the same gene family, but characterized by a different sequence (only 25% sequence homology), expression pattern and physiological function. Unrelated to amyloid neurodegenerative diseases, doppel supplies the perfect test system to investigate the misfolding of a non-amyloidogenic protein. Prion and doppel misfolding are compared in their monomeric form in the quest to identify prion-specific features that might reveal the mechanism of conversion to PrPSc. In agreement with experiments, we find a lower thermal stability for doppel. Surprisingly, we also observe β-rich forms for doppel. However, the β-rich folds of the two proteins are very different. Moreover, a major difference is found in the free energy barriers leading from the native structure to such conformations as well as to non-native conformations in general: These barriers are low for prion and can already be crossed at 300K, while for doppel they are at least 3 times higher. This difference suggests an intrinsic misfolding and β-enrichment propensity for the monomeric form of prion as compared to doppel.

EPFL2008

A DNA Coarse-Grain Rigid Base Model and Parameter Estimation from Molecular Dynamics Simulations

Daiva Petkeviciute

Sequence dependent mechanics of DNA is believed to play a central role in the functioning of the cell through the expression of genetic information. Nucleosome positioning, gene regulation, DNA looping and packaging within the cell are only some of the processes that are believed to be at least partially governed by mechanical laws. Therefore it is important to understand how the sequence of DNA affects its mechanical properties. For exploring the mechanical properties of DNA, various discrete and continuum models have been, and continue to be, developed. A large family of these models, including the model considered in this work, assume that bases or base pairs of DNA are rigid bodies. The most standard are rigid base pair models, with parameters either obtained directly from experimental data or from Molecular Dynamics (MD) simulations. The drawback of current experimental data, such as crystal structures, is that only small ensembles of configurations are available for a small number of sequences. In contrast, MD simulations allow a much more detailed view of a larger number of DNA sequences. However, the drawback is that the results of these simulations depend on the choice of the simulation protocol and force field parameters. MD simulations also have sequence length limitations and are currently too intensive for (linear) molecules longer than a few tens of base pairs. The only way to simulate longer sequences is to construct a coarse-grain model. The goal of this work is to construct a small parameter set that can model a sequence- dependent equilibrium probability distribution for rigid base configurations of a DNA oligomer with any given sequence of any length. The model parameter sets previously available were for rigid base pair models ignoring all the couplings beyond nearest neighbour interactions. However it was shown in previous work, that this standard model of rigid base pair nearest neighbour interactions is inconsistent with a (then) large scale MD simulation of a single oligomer [36]. In contrast we here show that a rigid base nearest neighbour, dimer sequence dependent model is a quite good fit to many MD simulations of different duration and se- quence. In fact a hierarchy of rigid base models with different interaction range and length of sequence-dependence is discussed, and it is concluded that the nearest neighbour, dimer based model is a good compromise between accuracy and complexity of the model. A full parameter set for this model is estimated. An interesting feature is that despite the dimer dependence of the parameter set, due to the phenomenon of frustration, our model predicts non local changes in the oligomer shape as a function of local changes in the sequence, down to the level of a point mutation.

EPFL2012