Molecular dynamics (MD) is acomputer simulation method for analyzing thephysical movements ofatoms andmolecules. The atoms and molecules are allowed to interact for a fixed period of time, giving a view of thedynamic "evolution" of the system. In the most common version, thetrajectories of atoms and molecules are determined bynumerically solvingNewton's equations of motion for a system of interacting particles, whereforces between the particles and theirpotential energies are often calculated usinginteratomic potentials ormolecular mechanicalforce fields. The method is applied mostly inchemical physics,materials science, andbiophysics.
Because molecular systems typically consist of a vast number of particles, it is impossible to determine the properties of suchcomplex systems analytically; MD simulation circumvents this problem by usingnumerical methods. However, long MD simulations are mathematicallyill-conditioned, generating cumulative errors in numerical integration that can be minimized with proper selection of algorithms and parameters, but not eliminated.
For systems that obey theergodic hypothesis, the evolution of one molecular dynamics simulation may be used to determine the macroscopicthermodynamic properties of the system: the time averages of an ergodic system correspond tomicrocanonical ensemble averages. MD has also been termed "statistical mechanics by numbers" and "Laplace's vision ofNewtonian mechanics" of predicting the future by animating nature's forces[1] and allowing insight into molecular motion on an atomic scale.
MD was originally developed in the early 1950s, following earlier successes withMonte Carlo simulations—which themselves date back to the eighteenth century, in theBuffon's needle problem for example—but was popularized forstatistical mechanics atLos Alamos National Laboratory byMarshall Rosenbluth andNicholas Metropolis in what is known today as theMetropolis–Hastings algorithm. Interest in the time evolution ofN-body systems dates much earlier to the seventeenth century, beginning withIsaac Newton, and continued into the following century largely with a focus oncelestial mechanics and issues such as thestability of the Solar System. Many of the numerical methods used today were developed during this time period, which predates the use of computers; for example, the most common integration algorithm used today, theVerlet integration algorithm, was used as early as 1791 byJean Baptiste Joseph Delambre. Numerical calculations with these algorithms can be considered to be MD done "by hand".
As early as 1941, integration of the many-body equations of motion was carried out withanalog computers. Some undertook the labor-intensive work of modeling atomic motion by constructing physical models, e.g., using macroscopic spheres. The aim was to arrange them in such a way as to replicate the structure of a liquid and use this to examine its behavior.J.D. Bernal describes this process in 1962, writing:[2]
... I took a number of rubber balls and stuck them together with rods of a selection of different lengths ranging from 2.75 to 4 inches. I tried to do this in the first place as casually as possible, working in my own office, being interrupted every five minutes or so and not remembering what I had done before the interruption.
Following the discovery of microscopic particles and the development of computers, interest expanded beyond the proving ground of gravitational systems to the statistical properties of matter. In an attempt to understand the origin ofirreversibility,Enrico Fermi proposed in 1953, and published in 1955,[3] the use of the early computerMANIAC I, also atLos Alamos National Laboratory, to solve the time evolution of the equations of motion for a many-body system subject to several choices of force laws. Today, this seminal work is known as theFermi–Pasta–Ulam–Tsingou problem. The time evolution of the energy from the original work is shown in the figure to the right.
In 1957,Berni Alder and Thomas Wainwright used anIBM 704 computer to simulate perfectlyelastic collisions betweenhard spheres.[4] In 1960, in perhaps the first realistic simulation of matter, J.B. Gibsonet al. simulated radiation damage ofsolid copper by using aBorn–Mayer type of repulsive interaction along with acohesive surface force.[5] In 1964,Aneesur Rahman published simulations of liquidargon that used aLennard-Jones potential; calculations of system properties, such as the coefficient ofself-diffusion, compared well with experimental data.[6] Today, the Lennard-Jones potential is still one of the most frequently usedintermolecular potentials.[7][8] It is used for describing simple substances (a.k.a.Lennard-Jonesium[9][10][11]) for conceptual and model studies and as a building block in manyforce fields of real substances.[12][13]
First used intheoretical physics, the molecular dynamics method gained popularity inmaterials science soon afterward, and since the 1970s it has also been commonly used inbiochemistry andbiophysics. MD is frequently used to refine 3-dimensional structures ofproteins and othermacromolecules based on experimental constraints fromX-ray crystallography orNMR spectroscopy. In physics, MD is used to examine the dynamics of atomic-level phenomena that cannot be observed directly, such asthin film growth and ion subplantation, and to examine the physical properties ofnanotechnological devices that have not or cannot yet be created. In biophysics andstructural biology, the method is frequently applied to study the motions of macromolecules such as proteins andnucleic acids, which can be useful for interpreting the results of certain biophysical experiments and for modeling interactions with other molecules, as inligand docking. In principle, MD can be used forab initio prediction ofprotein structure by simulatingfolding of thepolypeptide chain from arandom coil.
The results of MD simulations can be tested through comparison to experiments that measure molecular dynamics, of which a popular method is NMR spectroscopy. MD-derived structure predictions can be tested through community-wide experiments in Critical Assessment of Protein Structure Prediction (CASP), although the method has historically had limited success in this area.Michael Levitt, who shared theNobel Prize partly for the application of MD to proteins, wrote in 1999 that CASP participants usually did not use the method due to "... a central embarrassment ofmolecular mechanics, namely thatenergy minimization or molecular dynamics generally leads to a model that is less like the experimental structure".[14] Improvements in computational resources permitting more and longer MD trajectories, combined with modern improvements in the quality offorce field parameters, have yielded some improvements in both structure prediction andhomology model refinement, without reaching the point of practical utility in these areas; many identify force field parameters as a key area for further development.[15][16][17]
MD simulation has been reported forpharmacophore development anddrug design.[18] For example, Pintoet al. implemented MD simulations ofBcl-xL complexes to calculate average positions of criticalamino acids involved in ligand binding.[19] Carlsonet al. implemented molecular dynamics simulations to identify compounds that complement areceptor while causing minimal disruption to the conformation and flexibility of the active site. Snapshots of the protein at constant time intervals during the simulation were overlaid to identify conserved binding regions (conserved in at least three out of eleven frames) for pharmacophore development. Spyrakiset al. relied on a workflow of MD simulations, fingerprints for ligands and proteins (FLAP) andlinear discriminant analysis (LDA) to identify the best ligand-protein conformations to act as pharmacophore templates based on retrospectiveROC analysis of the resulting pharmacophores. In an attempt to ameliorate structure-based drug discovery modeling,vis-à-vis the need for many modeled compounds, Hatmalet al. proposed a combination of MD simulation and ligand-receptor intermolecular contacts analysis to discern critical intermolecular contacts (binding interactions) from redundant ones in a single ligand–protein complex. Critical contacts can then be converted into pharmacophore models that can be used for virtual screening.[20]
An important factor is intramolecularhydrogen bonds,[21] which are not explicitly included in modern force fields, but described asCoulomb interactions of atomicpoint charges.[citation needed] This is a crude approximation because hydrogen bonds have a partiallyquantum mechanical andchemical nature. Furthermore, electrostatic interactions are usually calculated using thedielectric constant of a vacuum, even though the surroundingaqueous solution has a much higher dielectric constant. Thus, using themacroscopic dielectric constant at short interatomic distances is questionable. Finally,van der Waals interactions in MD are usually described byLennard-Jones potentials[22][23] based on theFritz London theory that is only applicable in a vacuum.[citation needed] However, all types of van der Waals forces are ultimately of electrostatic origin and therefore depend ondielectric properties of the environment.[24] The direct measurement of attraction forces between different materials (asHamaker constant) shows that "the interaction betweenhydrocarbons across water is about 10% of that across vacuum".[24] The environment-dependence of van der Waals forces is neglected in standard simulations, but can be included by developing polarizable force fields.
The design of a molecular dynamics simulation should account for the available computational power. Simulation size (n = number of particles), timestep, and total time duration must be selected so that the calculation can finish within a reasonable time period. However, the simulations should be long enough to be relevant to the time scales of the natural processes being studied. To make statistically valid conclusions from the simulations, the time span simulated should match thekinetics of the natural process. Otherwise, it is analogous to making conclusions about how a human walks when only looking at less than one footstep. Most scientific publications about the dynamics of proteins and DNA[25][26] use data from simulations spanning nanoseconds (10−9 s) to microseconds (10−6 s). To obtain these simulations, severalCPU-days to CPU-years are needed.Parallel algorithms allow the load to be distributed amongCPUs; an example is the spatial or force decomposition algorithm.[27]
During a classical MD simulation, the most CPU intensive task is the evaluation of the potential as a function of the particles' internal coordinates. Within that energy evaluation, the most expensive one is the non-bonded or non-covalent part. Inbig O notation, common molecular dynamics simulationsscale by if all pair-wiseelectrostatic andvan der Waals interactions must be accounted for explicitly. This computational cost can be reduced by employingelectrostatics methods such as particle meshEwald summation ( ), particle-particle-particle mesh (P3M), or good spherical cutoff methods ( ).[citation needed]
Another factor that impacts total CPU time needed by a simulation is the size of the integration timestep. This is the time length between evaluations of the potential. The timestep must be chosen small enough to avoiddiscretization errors (i.e., smaller than the period related to fastest vibrational frequency in the system). Typical timesteps for classical MD are on the order of 1 femtosecond (10−15 s). This value may be extended by using algorithms such as the SHAKEconstraint algorithm, which fix the vibrations of the fastest atoms (e.g., hydrogens) into place. Multiple time scale methods have also been developed, which allow extended times between updates of slower long-range forces.[28][29][30]
For simulating molecules in asolvent, a choice should be made between anexplicit andimplicit solvent. Explicit solvent particles (such as theTIP3P, SPC/E andSPC-f water models) must be calculated expensively by the force field, while implicit solvents use amean-field approach. Using an explicit solvent is computationally expensive, requiring inclusion of roughly ten times more particles in the simulation. But the granularity and viscosity of explicit solvent is essential to reproduce certain properties of the solute molecules. This is especially important to reproducechemical kinetics.
In all kinds of molecular dynamics simulations, the simulation box size must be large enough to avoidboundary condition artifacts. Boundary conditions are often treated by choosing fixed values at the edges (which may cause artifacts), or by employingperiodic boundary conditions in which one side of the simulation loops back to the opposite side, mimicking a bulk phase (which may cause artifacts too).
In themicrocanonical ensemble, the system is isolated from changes inmoles (N), volume (V), and energy (E). It corresponds to anadiabatic process with no heat exchange. A microcanonical molecular dynamics trajectory may be seen as an exchange of potential and kinetic energy, with total energy being conserved. For a system ofN particles with coordinates and velocities, the following pair of first order differential equations may be written inNewton's notation as
Thepotential energy function of the system is a function of the particle coordinates. It is referred to simply as thepotential in physics, or theforce field in chemistry. The first equation comes fromNewton's laws of motion; the force acting on each particle in the system can be calculated as the negative gradient of.
For every time step, each particle's position and velocity may be integrated with asymplectic integrator method such asVerlet integration. The time evolution of and is called a trajectory. Given the initial positions (e.g., from theoretical knowledge) and velocities (e.g., randomizedGaussian), we can calculate all future (or past) positions and velocities.
One frequent source of confusion is the meaning oftemperature in MD. Commonly we have experience with macroscopic temperatures, which involve a huge number of particles, but temperature is a statistical quantity. If there is a large enough number of atoms, statistical temperature can be estimated from theinstantaneous temperature, which is found by equating the kinetic energy of the system tonkBT/2, wheren is the number of degrees of freedom of the system.
A temperature-related phenomenon arises due to the small number of atoms that are used in MD simulations. For example, consider simulating the growth of a copper film starting with a substrate containing 500 atoms and a deposition energy of 100eV. In the real world, the 100 eV from the deposited atom would rapidly be transported through and shared among a large number of atoms ( or more) with no big change in temperature. When there are only 500 atoms, however, the substrate is almost immediately vaporized by the deposition. Something similar happens in biophysical simulations. The temperature of the system in NVE is naturally raised when macromolecules such as proteins undergo exothermic conformational changes and binding.
In thecanonical ensemble, amount of substance (N), volume (V) and temperature (T) are conserved. It is also sometimes called constant temperature molecular dynamics (CTMD). In NVT, the energy of endothermic and exothermic processes is exchanged with athermostat.
A variety of thermostat algorithms are available to add and remove energy from the boundaries of an MD simulation in a more or less realistic way, approximating thecanonical ensemble. Popular methods to control temperature include velocity rescaling, theNosé–Hoover thermostat, Nosé–Hoover chains, theBerendsen thermostat, theAndersen thermostat andLangevin dynamics. The Berendsen thermostat might introduce theflying ice cube effect, which leads to unphysical translations and rotations of the simulated system.
It is not trivial to obtain a canonical ensemble distribution of conformations and velocities using these algorithms. How this depends on system size, thermostat choice, thermostat parameters, time step and integrator is the subject of many articles in the field.
In theisothermal–isobaric ensemble, amount of substance (N), pressure (P) and temperature (T) are conserved. In addition to a thermostat, abarostat is needed. It corresponds most closely to laboratory conditions with a flask open to ambient temperature and pressure.
In the simulation ofbiological membranes,isotropic pressure control is not appropriate. Forlipid bilayers, pressure control occurs under constant membrane area (NPAT) or constant surface tension "gamma" (NPγT).
Thereplica exchange method is a generalized ensemble. It was originally created to deal with the slow dynamics of disordered spin systems. It is also called parallel tempering. The replica exchange MD (REMD) formulation[31] tries to overcome the multiple-minima problem by exchanging the temperature of non-interacting replicas of the system running at several temperatures.
A molecular dynamics simulation requires the definition of apotential function, or a description of the terms by which the particles in the simulation will interact. In chemistry and biology this is usually referred to as aforce field and in materials physics as aninteratomic potential. Potentials may be defined at many levels of physical accuracy; those most commonly used in chemistry are based onmolecular mechanics and embody aclassical mechanics treatment of particle-particle interactions that can reproduce structural andconformational changes but usually cannot reproducechemical reactions.
The reduction from a fully quantum description to a classical potential entails two main approximations. The first one is theBorn–Oppenheimer approximation, which states that the dynamics of electrons are so fast that they can be considered to react instantaneously to the motion of their nuclei. As a consequence, they may be treated separately. The second one treats the nuclei, which are much heavier than electrons, as point particles that follow classical Newtonian dynamics. In classical molecular dynamics, the effect of the electrons is approximated as one potential energy surface, usually representing the ground state.
When finer levels of detail are needed, potentials based onquantum mechanics are used; some methods attempt to create hybridclassical/quantum potentials where the bulk of the system is treated classically but a small region is treated as a quantum system, usually undergoing a chemical transformation.
Empirical potentials used in chemistry are frequently calledforce fields, while those used in materials physics are calledinteratomic potentials.
Mostforce fields in chemistry are empirical and consist of a summation of bonded forces associated withchemical bonds, bond angles, and bonddihedrals, and non-bonded forces associated withvan der Waals forces andelectrostatic charge.[32] Empirical potentials represent quantum-mechanical effects in a limited way through ad hoc functional approximations. These potentials contain free parameters such asatomic charge, van der Waals parameters reflecting estimates ofatomic radius, and equilibriumbond length, angle, and dihedral; these are obtained by fitting against detailed electronic calculations (quantum chemical simulations) or experimental physical properties such aselastic constants, lattice parameters andspectroscopic measurements.
Because of the non-local nature of non-bonded interactions, they involve at least weak interactions between all particles in the system. Its calculation is normally the bottleneck in the speed of MD simulations. To lower the computational cost,force fields employ numerical approximations such as shifted cutoff radii,reaction field algorithms, particle meshEwald summation, or the newer particle–particle-particle–mesh (P3M).
Chemistry force fields commonly employ preset bonding arrangements (an exception beingab initio dynamics), and thus are unable to model the process of chemical bond breaking and reactions explicitly. On the other hand, many of the potentials used in physics, such as those based on thebond order formalism can describe several different coordinations of a system and bond breaking.[33][34] Examples of such potentials include theBrenner potential[35] for hydrocarbons and itsfurther developments for the C-Si-H[36] and C-O-H[37] systems. TheReaxFF potential[38] can be considered a fully reactive hybrid between bond order potentials and chemistry force fields.
The potential functions representing the non-bonded energy are formulated as a sum over interactions between the particles of the system. The simplest choice, employed in many popularforce fields, is the "pair potential", in which the total potential energy can be calculated from the sum of energy contributions between pairs of atoms. Therefore, these force fields are also called "additive force fields". An example of such a pair potential is the non-bondedLennard-Jones potential (also termed the 6–12 potential), used for calculating van der Waals forces.
Another example is the Born (ionic) model of the ionic lattice. The first term in the next equation isCoulomb's law for a pair of ions, the second term is the short-range repulsion explained by Pauli's exclusion principle and the final term is the dispersion interaction term. Usually, a simulation only includes the dipolar term, although sometimes the quadrupolar term is also included.[39][40] Whennl = 6, this potential is also called theCoulomb–Buckingham potential.
Inmany-body potentials, the potential energy includes the effects of three or more particles interacting with each other.[41] In simulations with pairwise potentials, global interactions in the system also exist, but they occur only through pairwise terms. In many-body potentials, the potential energy cannot be found by a sum over pairs of atoms, as these interactions are calculated explicitly as a combination of higher-order terms. In the statistical view, the dependency between the variables cannot in general be expressed using only pairwise products of the degrees of freedom. For example, theTersoff potential,[42] which was originally used to simulatecarbon,silicon, andgermanium, and has since been used for a wide range of other materials, involves a sum over groups of three atoms, with the angles between the atoms being an important factor in the potential. Other examples are theembedded-atom method (EAM),[43] the EDIP,[41] and the Tight-Binding Second Moment Approximation (TBSMA) potentials,[44] where the electron density of states in the region of an atom is calculated from a sum of contributions from surrounding atoms, and the potential energy contribution is then a function of this sum.
Semi-empirical potentials make use of the matrix representation from quantum mechanics. However, the values of the matrix elements are found through empirical formulae that estimate the degree of overlap of specific atomic orbitals. The matrix is then diagonalized to determine the occupancy of the different atomic orbitals, and empirical formulae are used once again to determine the energy contributions of the orbitals.
There are a wide variety of semi-empirical potentials, termedtight-binding potentials, which vary according to the atoms being modeled.
Most classical force fields implicitly include the effect ofpolarizability, e.g., by scaling up the partial charges obtained from quantum chemical calculations. These partial charges are stationary with respect to the mass of the atom. But molecular dynamics simulations can explicitly model polarizability with the introduction of induced dipoles through different methods, such asDrude particles or fluctuating charges. This allows for a dynamic redistribution of charge between atoms which responds to the local chemical environment.
For many years, polarizable MD simulations have been touted as the next generation. For homogenous liquids such as water, increased accuracy has been achieved through the inclusion of polarizability.[45][46][47] Some promising results have also been achieved for proteins.[48][49] However, it is still uncertain how to best approximate polarizability in a simulation.[citation needed] The point becomes more important when a particle experiences different environments during its simulation trajectory, e.g. translocation of a drug through a cell membrane.[50]
In classical molecular dynamics, one potential energy surface (usually the ground state) is represented in the force field. This is a consequence of theBorn–Oppenheimer approximation. In excited states, chemical reactions or when a more accurate representation is needed, electronic behavior can be obtained from first principles using a quantum mechanical method, such asdensity functional theory. This is namedAb Initio Molecular Dynamics (AIMD). Due to the cost of treating the electronic degrees of freedom, the computational burden of these simulations is far higher than classical molecular dynamics. For this reason, AIMD is typically limited to smaller systems and shorter times.
Ab initioquantum mechanical andchemical methods may be used to calculate thepotential energy of a system on the fly, as needed for conformations in a trajectory. This calculation is usually made in the close neighborhood of thereaction coordinate. Although various approximations may be used, these are based on theoretical considerations, not on empirical fitting.Ab initio calculations produce a vast amount of information that is not available from empirical methods, such as density of electronic states or other electronic properties. A significant advantage of usingab initio methods is the ability to study reactions that involve breaking or formation of covalent bonds, which correspond to multiple electronic states. Moreover,ab initio methods also allow recovering effects beyond the Born–Oppenheimer approximation using approaches likemixed quantum-classical dynamics.
QM (quantum-mechanical) methods are very powerful. However, they are computationally expensive, while the MM (classical or molecular mechanics) methods are fast but suffer from several limits (require extensive parameterization; energy estimates obtained are not very accurate; cannot be used to simulate reactions where covalent bonds are broken/formed; and are limited in their abilities for providing accurate details regarding the chemical environment). A new class of method has emerged that combines the good points of QM (accuracy) and MM (speed) calculations. These methods are termed mixed or hybrid quantum-mechanical and molecular mechanics methods (hybrid QM/MM).[51]
The most important advantage of hybrid QM/MM method is the speed. The cost of doing classical molecular dynamics (MM) in the most straightforward case scales O(n2), where n is the number of atoms in the system. This is mainly due to electrostatic interactions term (every particle interacts with every other particle). However, use of cutoff radius, periodic pair-list updates and more recently the variations of the particle-mesh Ewald's (PME) method has reduced this to between O(n) to O(n2). In other words, if a system with twice as many atoms is simulated then it would take between two and four times as much computing power. On the other hand, the simplestab initio calculations typically scale O(n3) or worse (restrictedHartree–Fock calculations have been suggested to scale ~O(n2.7)). To overcome the limit, a small part of the system is treated quantum-mechanically (typically active-site of an enzyme) and the remaining system is treated classically.
In more sophisticated implementations, QM/MM methods exist to treat both light nuclei susceptible to quantum effects (such as hydrogens) and electronic states. This allows generating hydrogen wave-functions (similar to electronic wave-functions). This methodology has been useful in investigating phenomena such as hydrogen tunneling. One example where QM/MM methods have provided new discoveries is the calculation of hydride transfer in the enzyme liveralcohol dehydrogenase. In this case,quantum tunneling is important for the hydrogen, as it determines the reaction rate.[52]
At the other end of the detail scale arecoarse-grained and lattice models. Instead of explicitly representing every atom of the system, one uses "pseudo-atoms" to represent groups of atoms. MD simulations on very large systems may require such large computer resources that they cannot easily be studied by traditional all-atom methods. Similarly, simulations of processes on long timescales (beyond about 1 microsecond) are prohibitively expensive, because they require so many time steps. In these cases, one can sometimes tackle the problem by using reduced representations, which are also calledcoarse-grained models.[53]
Examples for coarse graining (CG) methods are discontinuous molecular dynamics (CG-DMD)[54][55] and Go-models.[56] Coarse-graining is done sometimes taking larger pseudo-atoms. Such united atom approximations have been used in MD simulations of biological membranes. Implementation of such approach on systems where electrical properties are of interest can be challenging owing to the difficulty of using a proper charge distribution on the pseudo-atoms.[57] The aliphatic tails of lipids are represented by a few pseudo-atoms by gathering 2 to 4 methylene groups into each pseudo-atom.
The parameterization of these very coarse-grained models must be done empirically, by matching the behavior of the model to appropriate experimental data or all-atom simulations. Ideally, these parameters should account for bothenthalpic andentropic contributions to free energy in an implicit way.[58] When coarse-graining is done at higher levels, the accuracy of the dynamic description may be less reliable. But very coarse-grained models have been used successfully to examine a wide range of questions in structural biology, liquid crystal organization, and polymer glasses.
Examples of applications of coarse-graining:
The simplest form of coarse-graining is theunited atom (sometimes calledextended atom) and was used in most early MD simulations of proteins, lipids, and nucleic acids. For example, instead of treating all four atoms of a CH3 methyl group explicitly (or all three atoms of CH2 methylene group), one represents the whole group with one pseudo-atom. It must, of course, be properly parameterized so that its van der Waals interactions with other groups have the proper distance-dependence. Similar considerations apply to the bonds, angles, and torsions in which the pseudo-atom participates. In this kind of united atom representation, one typically eliminates all explicit hydrogen atoms except those that have the capability to participate in hydrogen bonds (polar hydrogens). An example of this is theCHARMM 19 force-field.
The polar hydrogens are usually retained in the model, because proper treatment of hydrogen bonds requires a reasonably accurate description of the directionality and the electrostatic interactions between the donor and acceptor groups. A hydroxyl group, for example, can be both a hydrogen bond donor, and a hydrogen bond acceptor, and it would be impossible to treat this with one OH pseudo-atom. About half the atoms in a protein or nucleic acid are non-polar hydrogens, so the use of united atoms can provide a substantial savings in computer time.
Machine Learning Force Fields] (MLFFs) represent one approach to modeling interatomic interactions in molecular dynamics simulations.[59] MLFFs can achieve accuracy close to that ofab initio methods. Once trained, MLFFs are much faster than direct quantum mechanical calculations. MLFFs address the limitations of traditional force fields by learning complex potential energy surfaces directly from high-level quantum mechanical data. Several software packages now support MLFFs, includingVASP[60] and open-source libraries like DeePMD-kit[61][62] andSchNetPack.[63][64]
In many simulations of a solute-solvent system the main focus is on the behavior of the solute with little interest of the solvent behavior particularly in those solvent molecules residing in regions far from the solute molecule.[65] Solvents may influence the dynamic behavior of solutes via random collisions and by imposing a frictional drag on the motion of the solute through the solvent. The use of non-rectangular periodic boundary conditions, stochastic boundaries and solvent shells can all help reduce the number of solvent molecules required and enable a larger proportion of the computing time to be spent instead on simulating the solute. It is also possible to incorporate the effects of a solvent without needing any explicit solvent molecules present. One example of this approach is to use apotential mean force (PMF) which describes how the free energy changes as a particular coordinate is varied. The free energy change described by PMF contains the averaged effects of the solvent.
Without incorporating the effects of solvent simulations of macromolecules (such as proteins) may yield unrealistic behavior and even small molecules may adopt more compact conformations due to favourable van der Waals forces and electrostatic interactions which would be dampened in the presence of a solvent.[66]
A long range interaction is an interaction in which the spatial interaction falls off no faster than where is the dimensionality of the system. Examples include charge-charge interactions between ions and dipole-dipole interactions between molecules. Modelling these forces presents quite a challenge as they are significant over a distance which may be larger than half the box length with simulations of many thousands of particles. Though one solution would be to significantly increase the size of the box length, this brute force approach is less than ideal as the simulation would become computationally very expensive. Spherically truncating the potential is also out of the question as unrealistic behaviour may be observed when the distance is close to the cut off distance.[67]
Steered molecular dynamics (SMD) simulations, or force probe simulations, apply forces to a protein in order to manipulate its structure by pulling it along desired degrees of freedom. These experiments can be used to reveal structural changes in a protein at the atomic level. SMD is often used to simulate events such as mechanical unfolding or stretching.[68]
There are two typical protocols of SMD: one in which pulling velocity is held constant, and one in which applied force is constant. Typically, part of the studied system (e.g., an atom in a protein) is restrained by a harmonic potential. Forces are then applied to specific atoms at either a constant velocity or a constant force.Umbrella sampling is used to move the system along the desired reaction coordinate by varying, for example, the forces, distances, and angles manipulated in the simulation. Through umbrella sampling, all of the system's configurations—both high-energy and low-energy—are adequately sampled. Then, each configuration's change in free energy can be calculated as thepotential of mean force.[69] A popular method of computing PMF is through the weighted histogram analysis method (WHAM), which analyzes a series of umbrella sampling simulations.[70][71]
A lot of important applications of SMD are in the field of drug discovery and biomolecular sciences. For e.g. SMD was used to investigate the stability of Alzheimer's protofibrils,[72] to study the protein ligand interaction in cyclin-dependent kinase 5[73] and even to show the effect of electric field on thrombin (protein) and aptamer (nucleotide) complex[74] among many other interesting studies.
 | This sectionrelies excessively onreferences toprimary sources. Please improve this section by addingsecondary or tertiary sources. Find sources: "Molecular dynamics" – news ·newspapers ·books ·scholar ·JSTOR(January 2024) (Learn how and when to remove this message) |
Molecular dynamics is used in many fields of science.
The following biophysical examples illustrate notable efforts to produce simulations of a systems of very large size (a complete virus) or very long simulation times (up to 1.112 milliseconds):
Another important application of MD method benefits from its ability of 3-dimensional characterization and analysis of microstructural evolution at atomic scale.
Molecular modeling on GPU is the technique of using agraphics processing unit (GPU) for molecular simulations.[87]
In 2007,Nvidia introduced video cards that could be used not only to show graphics but also for scientific calculations. These cards include many arithmetic units (as of 2016[update], up to 3,584 in Tesla P100) working in parallel. Long before this event, the computational power of video cards was purely used to accelerate graphics calculations. The new features of these cards made it possible to develop parallel programs in a high-levelapplication programming interface (API) namedCUDA. This technology substantially simplified programming by enabling programs to be written inC/C++. More recently,OpenCL allowscross-platform GPU acceleration.
Branches ofchemistry | |
|---|---|
| Analytical | |
| Theoretical | |
| Physical | |
| Inorganic | |
| Organic | |
| Biological | |
| Interdisciplinarity | |
| See also | |
