View source for Molecular dynamics

{{Short description|Computer simulations to discover and understand chemical properties}}
{{Use American English|date = February 2019}}
{{cs1 config|name-list-style=vanc|display-authors=6}}
[[File:Cudeposition.gif|thumb|upright=1.35|Example of a molecular dynamics simulation in a simple system: deposition of one [[copper]] (Cu) [[atom]] on a cold crystal of copper ([[Miller index]] (001) [[Surface science|surface]]). Each circle represents the position of one atom. The kinetic energy of the atom approaching from the top is redistributed among the other atoms, so instead of bouncing off it remains attached due to attractive forces between the atoms.]]
[[File:MD water.gif|thumb|upright=1.35|Molecular dynamics simulations are often used to study biophysical systems. Depicted here is a 100 ps simulation of water.]]
[[File:Molecular dynamics algorithm.png|thumb|upright=1.8|A simplified description of the standard molecular dynamics simulation algorithm, when a predictor-corrector-type integrator is used. The forces may come either from classical [[interatomic potential]]s (described mathematically as <math> F= - \nabla V(\vec r) </math>) or quantum mechanical (described mathematically as <math> F=F( \Psi(\vec r) ) </math>) methods. Large differences exist between different integrators; some do not have exactly the same highest-order terms as indicated in the flow chart, many also use higher-order time derivatives, and some use both the current and prior time step in variable-time step schemes.]]

'''Molecular dynamics''' ('''MD''') is a [[computer simulation]] method for analyzing the [[Motion (physics)|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 [[dynamics (mechanics)|dynamic]] "evolution" of the system. In the most common version, the [[trajectory|trajectories]] of atoms and molecules are determined by [[Numerical integration|numerically solving]] [[Newton's laws of motion|Newton's equations of motion]] for a system of interacting particles, where [[Force (physics)|forces]] between the particles and their [[potential energy|potential energies]] are often calculated using [[interatomic potential]]s or [[molecular mechanics|molecular mechanical]] [[Force field (chemistry)|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 analysis|numerical]] methods. However, long MD simulations are mathematically [[condition number|ill-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 the [[ergodic hypothesis]], the evolution of one molecular dynamics simulation may be used to determine the macroscopic [[thermodynamic]] properties of the system: the time averages of an ergodic system correspond to [[microcanonical ensemble]] averages. MD has also been termed "[[statistical mechanics]] by numbers" and "[[Laplace]]'s vision of [[Newtonian mechanics]]" of predicting the future by animating nature's forces<ref>{{cite book |doi=10.1007/978-1-4612-4066-2_13 |chapter=Pursuing Laplace's Vision on Modern Computers |title=Mathematical Approaches to Biomolecular Structure and Dynamics |series=The IMA Volumes in Mathematics and its Applications |year=1996 | vauthors = Schlick T |volume=82 |pages=219–247 |isbn=978-0-387-94838-6 }}</ref> and allowing insight into molecular motion on an atomic scale.

== History ==
MD was originally developed in the early 1950s, following earlier successes with [[Monte Carlo simulation]]s{{Em dash}}which themselves date back to the eighteenth century, in the [[Buffon's needle problem]] for example{{Em dash}}but was popularized for [[statistical mechanics]] at [[Los Alamos National Laboratory]] by [[Marshall Rosenbluth]] and [[Nicholas Metropolis]] in what is known today as the [[Metropolis–Hastings algorithm]]. Interest in the time evolution of [[N-body problem|N-body systems]] dates much earlier to the seventeenth century, beginning with [[Isaac Newton]], and continued into the following century largely with a focus on [[celestial mechanics]] and issues such as the [[Stability of the Solar System|stability 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, the [[Verlet integration]] algorithm, was used as early as 1791 by [[Jean 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 with [[Analog computer|analog 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:<ref>{{cite journal |title=The Bakerian Lecture, 1962 The structure of liquids |journal=Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences |date=January 1997 |volume=280 |issue=1382 |pages=299–322 |doi=10.1098/rspa.1964.0147 |bibcode=1964RSPSA.280..299B | vauthors = Bernal JD |s2cid=178710030 }}</ref><blockquote>... 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.</blockquote>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 of [[irreversibility]], [[Enrico Fermi]] proposed in 1953, and published in 1955,<ref name=fput>Fermi E., Pasta J., Ulam S., Los Alamos report LA-1940 (1955).</ref> the use of the early computer [[MANIAC I]], also at [[Los 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 the [[Fermi–Pasta–Ulam–Tsingou problem]]. The time evolution of the energy from the original work is shown in the figure to the right.

[[File:Time evolution of energy for FPUT N-body dynamics.jpg|thumb|upright=1.35|One of the earliest simulations of an N-body system was carried out on the MANIAC-I by Fermi and coworkers to understand the origins of irreversibility in nature. Shown here is the energy versus time for a 64-particle system.]]
In 1957, [[Berni Alder]] and Thomas Wainwright used an [[IBM 704]] computer to simulate perfectly [[Elastic collision|elastic collisions]] between [[hard spheres]].<ref name="a&w">{{cite journal |vauthors=Alder BJ, Wainwright T |date=August 1959 |title=Studies in Molecular Dynamics. I. General Method |journal=The Journal of Chemical Physics |volume=31 |issue=2 |pages=459–466 |bibcode=1959JChPh..31..459A |doi=10.1063/1.1730376}}</ref> In 1960, in perhaps the first realistic simulation of matter, J.B. Gibson ''et al''. simulated radiation damage of [[Native copper|solid copper]] by using a [[Born–Mayer equation|Born–Mayer]] type of repulsive interaction along with a [[Cohesion (chemistry)|cohesive]] surface force.<ref>{{cite journal |vauthors=Gibson JB, Goland AN, Milgram M, Vineyard G |year=1960 |title=Dynamics of Radiation Damage |journal=Phys. Rev. |volume=120 |issue=4 |pages=1229–1253 |doi=10.1103/PhysRev.120.1229 |bibcode=1960PhRv..120.1229G}}</ref> In 1964, [[Aneesur Rahman]] published simulations of liquid [[argon]] that used a [[Lennard-Jones potential]]; calculations of system properties, such as the coefficient of [[self-diffusion]], compared well with experimental data.<ref name="a.rahman">{{cite journal |author-link1=Aneesur Rahman |vauthors=Rahman A |date=19 October 1964 |title=Correlations in the Motion of Atoms in Liquid Argon |journal=Physical Review |volume=136 |issue=2A |pages=A405–A411 |bibcode=1964PhRv..136..405R |doi=10.1103/PhysRev.136.A405}}</ref> Today, the Lennard-Jones potential is still one of the most frequently used [[Interatomic potential|intermolecular potentials]].<ref>{{cite journal |vauthors=Stephan S, Thol M, Vrabec J, Hasse H |title=Thermophysical Properties of the Lennard-Jones Fluid: Database and Data Assessment |journal=Journal of Chemical Information and Modeling |volume=59 |issue=10 |pages=4248–4265 |date=October 2019 |pmid=31609113 |doi=10.1021/acs.jcim.9b00620 |s2cid=204545481 |url=https://depositonce.tu-berlin.de/handle/11303/10447}}</ref><ref>{{cite journal |vauthors=Wang X, Ramírez-Hinestrosa S, Dobnikar J, Frenkel D |title=The Lennard-Jones potential: when (not) to use it |journal=Physical Chemistry Chemical Physics |volume=22 |issue=19 |pages=10624–10633 |date=May 2020 |pmid=31681941 |doi=10.1039/C9CP05445F |arxiv=1910.05746 |s2cid=204512243 |bibcode=2020PCCP...2210624W}}</ref> It is used for describing simple substances (a.k.a. [[Lennard-Jones potential|Lennard-Jonesium]]<ref>{{Cite journal |vauthors=Mick J, Hailat E, Russo V, Rushaidat K, Schwiebert L, Potoff J |date=December 2013 |title=GPU-accelerated Gibbs ensemble Monte Carlo simulations of Lennard-Jonesium |journal=Computer Physics Communications |language=en |volume=184 |issue=12 |pages=2662–2669 |doi=10.1016/j.cpc.2013.06.020 |bibcode=2013CoPhC.184.2662M}}</ref><ref>{{Cite journal |vauthors=Chapela GA, Scriven LE, Davis HT |date=October 1989 |title=Molecular dynamics for discontinuous potential. IV. Lennard-Jonesium |url=http://aip.scitation.org/doi/10.1063/1.456811 |journal=The Journal of Chemical Physics |language=en |volume=91 |issue=7 |pages=4307–4313 |doi=10.1063/1.456811 |bibcode=1989JChPh..91.4307C |issn=0021-9606}}</ref><ref>{{cite journal | vauthors = Lenhard J, Stephan S, Hasse H | title = A child of prediction. On the History, Ontology, and Computation of the Lennard-Jonesium | journal = Studies in History and Philosophy of Science | volume = 103 | pages = 105–113 | date = February 2024 | pmid = 38128443 | doi = 10.1016/j.shpsa.2023.11.007 | s2cid = 266440296 }}</ref>) for conceptual and model studies and as a building block in many [[Force field (chemistry)|force fields]] of real substances.<ref>{{Cite journal |vauthors=Eggimann BL, Sunnarborg AJ, Stern HD, Bliss AP, Siepmann JI |date=2013-12-24 |title=An online parameter and property database for the TraPPE force field |journal=Molecular Simulation |volume=40 |issue=1–3 |pages=101–105 |doi=10.1080/08927022.2013.842994 |s2cid=95716947 |issn=0892-7022}}</ref><ref>{{Cite journal |vauthors=Stephan S, Horsch MT, Vrabec J, Hasse H |date=2019-07-03 |title=MolMod – an open access database of force fields for molecular simulations of fluids |journal=Molecular Simulation |language=en |volume=45 |issue=10 |pages=806–814 |doi=10.1080/08927022.2019.1601191 |s2cid=119199372 |issn=0892-7022|url=https://osf.io/m87cy/ |arxiv=1904.05206 }}</ref>

== Areas of application and limits ==

First used in [[theoretical physics]], the molecular dynamics method gained popularity in [[materials science]] soon afterward, and since the 1970s it has also been commonly used in [[biochemistry]] and [[biophysics]]. MD is frequently used to refine 3-dimensional structures of [[protein]]s and other [[macromolecule]]s based on experimental constraints from [[X-ray crystallography]] or [[Nuclear magnetic resonance spectroscopy|NMR spectroscopy]]. In physics, MD is used to examine the dynamics of atomic-level phenomena that cannot be observed directly, such as [[thin film]] growth and ion subplantation, and to examine the physical properties of [[nanotechnology|nanotechnological]] devices that have not or cannot yet be created. In biophysics and [[structural biology]], the method is frequently applied to study the motions of macromolecules such as proteins and [[nucleic acid]]s, which can be useful for interpreting the results of certain biophysical experiments and for modeling interactions with other molecules, as in [[ligand docking]]. In principle, MD can be used for [[De novo protein structure prediction|''ab initio'' prediction]] of [[protein structure]] by simulating [[protein folding|folding]] of the [[polypeptide chain]] from a [[random coil]]. MD can also be used to compute other thermodynamic properties such as drug solubilities and free energies of solvation <ref>{{Cite journal |last=Suarez |first=Adiran Garaizar |last2=Göller |first2=Andreas H. |last3=Beck |first3=Michael E. |last4=Gheta |first4=Sadra Kashef Ol |last5=Meier |first5=Katharina |date=2024-10-29 |title=Comparative assessment of physics-based in silico methods to calculate relative solubilities |url=https://link.springer.com/article/10.1007/s10822-024-00576-y |journal=Journal of Computer-Aided Molecular Design |language=en |volume=38 |issue=1 |pages=36 |doi=10.1007/s10822-024-00576-y |issn=1573-4951}}</ref>  including in polymers. <ref>{{Cite journal |last=Higginbotham |first=T. |last2=Meier |first2=K. |last3=Ramírez |first3=J. |last4=Garaizar |first4=A. |date=2025-02-03 |title=Predicting Drug-Polymer Compatibility in Amorphous Solid Dispersions by MD Simulation: On the Trap of Solvation Free Energies |url=https://pubs.acs.org/doi/10.1021/acs.molpharmaceut.4c00810 |journal=Molecular Pharmaceutics |volume=22 |issue=2 |pages=760–770 |doi=10.1021/acs.molpharmaceut.4c00810 |issn=1543-8384}}</ref>  

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 (biophysicist)|Michael Levitt]], who shared the [[Nobel 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 of [[molecular mechanics]], namely that [[energy minimization]] or molecular dynamics generally leads to a model that is less like the experimental structure".<ref name="Koehl">{{cite journal | vauthors = Koehl P, Levitt M | title = A brighter future for protein structure prediction | journal = Nature Structural Biology | volume = 6 | issue = 2 | pages = 108–111 | date = February 1999 | pmid = 10048917 | doi = 10.1038/5794 | author-link2 = Michael Levitt (biophysicist) | s2cid = 3162636 }}</ref> Improvements in computational resources permitting more and longer MD trajectories, combined with modern improvements in the quality of [[force field (chemistry)|force field]] parameters, have yielded some improvements in both structure prediction and [[homology model]] refinement, without reaching the point of practical utility in these areas; many identify force field parameters as a key area for further development.<ref>{{cite journal | vauthors = Raval A, Piana S, Eastwood MP, Dror RO, Shaw DE | title = Refinement of protein structure homology models via long, all-atom molecular dynamics simulations | journal = Proteins | volume = 80 | issue = 8 | pages = 2071–2079 | date = August 2012 | pmid = 22513870 | doi = 10.1002/prot.24098 | s2cid = 10613106 }}</ref><ref>{{cite journal | vauthors = Beauchamp KA, Lin YS, Das R, Pande VS | title = Are Protein Force Fields Getting Better? A Systematic Benchmark on 524 Diverse NMR Measurements | journal = Journal of Chemical Theory and Computation | volume = 8 | issue = 4 | pages = 1409–1414 | date = April 2012 | pmid = 22754404 | pmc = 3383641 | doi = 10.1021/ct2007814 }}</ref><ref>{{cite journal | vauthors = Piana S, Klepeis JL, Shaw DE | title = Assessing the accuracy of physical models used in protein-folding simulations: quantitative evidence from long molecular dynamics simulations | journal = Current Opinion in Structural Biology | volume = 24 | pages = 98–105 | date = February 2014 | pmid = 24463371 | doi = 10.1016/j.sbi.2013.12.006 | doi-access = free }}</ref>

MD simulation has been reported for [[pharmacophore]] development and [[drug design]].<ref name="pmid25751016">{{cite journal | vauthors = Choudhury C, Priyakumar UD, Sastry GN | title = Dynamics based pharmacophore models for screening potential inhibitors of mycobacterial cyclopropane synthase | journal = Journal of Chemical Information and Modeling | volume = 55 | issue = 4 | pages = 848–60 | date = April 2015 | pmid = 25751016 | doi = 10.1021/ci500737b }}</ref> For example, Pinto ''et al''. implemented MD simulations of [[Bcl-xL|Bcl-xL complexes]] to calculate average positions of critical [[Amino acid|amino acids]] involved in ligand binding.<ref name="pmid15143800">{{cite journal | vauthors = Pinto M, Perez JJ, Rubio-Martinez J | title = Molecular dynamics study of peptide segments of the BH3 domain of the proapoptotic proteins Bak, Bax, Bid and Hrk bound to the Bcl-xL and Bcl-2 proteins | journal = Journal of Computer-aided Molecular Design | volume = 18 | issue = 1 | pages = 13–22 | date = January 2004 | pmid = 15143800 | doi = 10.1023/b:jcam.0000022559.72848.1c | bibcode = 2004JCAMD..18...13P | s2cid = 11339000 }}</ref> Carlson ''et al''. implemented molecular dynamics simulations to identify compounds that complement a [[Receptor (biochemistry)|receptor]] 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. Spyrakis ''et al''. relied on a workflow of MD simulations, fingerprints for ligands and proteins (FLAP) and [[linear discriminant analysis]] (LDA) to identify the best ligand-protein conformations to act as pharmacophore templates based on retrospective [[Receiver operating characteristic|ROC]] analysis of the resulting pharmacophores. In an attempt to ameliorate structure-based drug discovery modeling, ''vis-à-vis'' the need for many modeled compounds, Hatmal ''et 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.<ref name="pmid27722817">{{cite journal | vauthors = Hatmal MM, Jaber S, Taha MO | title = Combining molecular dynamics simulation and ligand-receptor contacts analysis as a new approach for pharmacophore modeling: beta-secretase 1 and check point kinase 1 as case studies | journal = Journal of Computer-aided Molecular Design | volume = 30 | issue = 12 | pages = 1149–1163 | date = December 2016 | pmid = 27722817 | doi = 10.1007/s10822-016-9984-2  | bibcode = 2016JCAMD..30.1149H | s2cid = 11561853 }}</ref>

An important factor is intramolecular [[hydrogen bond]]s,<ref name="Myers">{{cite journal | vauthors = Myers JK, Pace CN | title = Hydrogen bonding stabilizes globular proteins | journal = Biophysical Journal | volume = 71 | issue = 4 | pages = 2033–2039 | date = October 1996 | pmid = 8889177 | pmc = 1233669 | doi = 10.1016/s0006-3495(96)79401-8 | bibcode = 1996BpJ....71.2033M }}</ref> which are not explicitly included in modern force fields, but described as [[Coulomb interaction|Coulomb interactions]] of atomic [[Point charge|point charges]].{{Citation needed|date=February 2024|reason=Needs a citation to show this is a general issue and not an issue with specific forcefields}} This is a crude approximation because hydrogen bonds have a partially [[quantum mechanical]] and [[Quantum chemistry|chemical]] nature. Furthermore, electrostatic interactions are usually calculated using the [[Vacuum permittivity|dielectric constant of a vacuum]], even though the surrounding [[aqueous solution]] has a much higher dielectric constant. Thus, using the [[macroscopic]] dielectric constant at short interatomic distances is questionable. Finally, [[Van der Waals force|van der Waals interactions]] in MD are usually described by [[Lennard-Jones potential]]s<ref>{{Cite journal | vauthors = Lenhard J, Stephan S, Hasse H |date=June 2024 |title=On the History of the Lennard-Jones Potential |url=https://onlinelibrary.wiley.com/doi/10.1002/andp.202400115 |journal=Annalen der Physik |language=en |volume=536 |issue=6 |doi=10.1002/andp.202400115 |issn=0003-3804}}</ref><ref>{{Cite journal | vauthors = Fischer J, Wendland M |date=October 2023 |title=On the history of key empirical intermolecular potentials |journal=Fluid Phase Equilibria |language=en |volume=573 |pages=113876 |doi=10.1016/j.fluid.2023.113876|bibcode=2023FlPEq.57313876F |doi-access=free }}</ref> based on the [[Fritz London]] theory that is only applicable in a vacuum.{{Citation needed|date=February 2024|reason=needs a citation to show this true in the general case and not just specific forcefields.}} However, all types of van der Waals forces are ultimately of electrostatic origin and therefore depend on [[Permittivity|dielectric properties of the environment]].<ref name="Israelachvili">{{cite book | vauthors = Israelachvili J | author-link1 = Jacob Israelachvili | date = 1992 | title = Intermolecular and surface forces. | publisher = Academic Press | location = San Diego }}</ref> The direct measurement of attraction forces between different materials (as [[Hamaker constant]]) shows that "the interaction between [[Hydrocarbon|hydrocarbons]] across water is about 10% of that across vacuum".<ref name="Israelachvili" /> The environment-dependence of van der Waals forces is neglected in standard simulations, but can be included by developing polarizable force fields.

== Design constraints ==

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 the [[Kinetics (physics)|kinetics]] 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<ref name="ReferenceB">{{cite journal | vauthors = Cruz FJ, de Pablo JJ, Mota JP | title = Endohedral confinement of a DNA dodecamer onto pristine carbon nanotubes and the stability of the canonical B form | journal = The Journal of Chemical Physics | volume = 140 | issue = 22 | pages = 225103 | date = June 2014 | pmid = 24929415 | doi = 10.1063/1.4881422 | arxiv = 1605.01317 | s2cid = 15149133 | bibcode = 2014JChPh.140v5103C }}</ref><ref name="ReferenceC">{{cite journal | vauthors = Cruz FJ, Mota JP |title=Conformational Thermodynamics of DNA Strands in Hydrophilic Nanopores |journal=J. Phys. Chem. C |volume=120 |year=2016 |issue=36 |pages=20357–20367 |doi=10.1021/acs.jpcc.6b06234}}</ref> use data from simulations spanning nanoseconds (10<sup>−9</sup> s) to microseconds (10<sup>−6</sup> s). To obtain these simulations, several [[CPU time|CPU-days]] to CPU-years are needed. [[Parallel algorithm|Parallel algorithms]] allow the load to be distributed among [[Central processing unit|CPUs]]; an example is the spatial or force decomposition algorithm.<ref>{{cite web | vauthors = Plimpton S | url = http://www.sandia.gov/~sjplimp/md.html | title = Molecular Dynamics - Parallel Algorithms | work = sandia.gov }}</ref>

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. In [[big O notation]], common molecular dynamics simulations [[analysis of algorithms|scale]] by <math>O(n^2)</math> if all pair-wise [[electrostatic]] and [[van der Waals forces|van der Waals interactions]] must be accounted for explicitly. This computational cost can be reduced by employing [[electrostatics]] methods such as particle mesh [[Ewald summation]] ( <math>O(n \log(n))</math> ), particle-particle-particle mesh ([[P3M|P<sup>3</sup>M]]), or good spherical cutoff methods ( <math>O(n)</math> ). {{citation needed |date=December 2017}}

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 avoid [[Discretization error|discretization 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&nbsp;femtosecond (10<sup>−15</sup> s). This value may be extended by using algorithms such as the SHAKE [[constraint 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.<ref name="Streett">{{cite journal | vauthors = Streett WB, Tildesley DJ, Saville G |year=1978 |title=Multiple time-step methods in molecular dynamics |journal=Mol Phys |volume=35 |issue=3 |pages=639–648 |doi=10.1080/00268977800100471|bibcode= 1978MolPh..35..639S}}</ref><ref name="Tuckerman1991">{{cite journal | vauthors = Tuckerman ME, Berne BJ, Martyna GJ |year=1991 |title=Molecular dynamics algorithm for multiple time scales: systems with long range forces |journal=J Chem Phys |volume=94 |issue=10 |pages=6811–6815 |doi=10.1063/1.460259|bibcode= 1991JChPh..94.6811T }}</ref><ref name="Tuckerman1992">{{cite journal | vauthors = Tuckerman ME, Berne BJ, Martyna GJ |year=1992 |title=Reversible multiple time scale molecular dynamics |journal=J Chem Phys |volume=97 |issue=3 |pages=1990–2001 |doi=10.1063/1.463137|bibcode= 1992JChPh..97.1990T |s2cid=488073 }}</ref>

For simulating molecules in a [[solvent]], a choice should be made between an [[Water model|explicit]] and [[implicit solvent]]. Explicit solvent particles (such as the [[TIP3P]], SPC/E and [[Flexible SPC water model|SPC-f]] water models) must be calculated expensively by the force field, while implicit solvents use a [[Mean-field particle methods|mean-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 reproduce [[chemical kinetics]].

In all kinds of molecular dynamics simulations, the simulation box size must be large enough to avoid [[boundary condition]] artifacts. Boundary conditions are often treated by choosing fixed values at the edges (which may cause artifacts), or by employing [[periodic boundary conditions]] in which one side of the simulation loops back to the opposite side, mimicking a bulk phase (which may cause artifacts too).

[[File:Sampling in Monte Carlo and molecular dynamics.png|thumb|upright=1.5|Schematic representation of the sampling of the system's potential energy surface with molecular dynamics (in red) compared to Monte Carlo methods (in blue)]]

=== Microcanonical ensemble (NVE) ===

In the [[microcanonical ensemble]], the system is isolated from changes in [[Mole (unit)|moles]] (N), volume (V), and energy (E). It corresponds to an [[adiabatic 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 of ''N'' particles with coordinates <math>X</math> and velocities <math>V</math>, the following pair of first order differential equations may be written in [[Newton's notation for differentiation|Newton's notation]] as

:<math>F(X) = -\nabla U(X)=M\dot{V}(t)</math>
:<math>V(t) = \dot{X} (t).</math>

The [[Potential energy|potential energy function]] <math>U(X)</math> of the system is a function of the particle coordinates <math>X</math>. It is referred to simply as the ''potential'' in physics, or the ''[[Force field (chemistry)|force field]]'' in chemistry. The first equation comes from [[Newton's laws of motion]]; the force <math>F</math> acting on each particle in the system can be calculated as the negative gradient of <math>U(X)</math>.

For every time step, each particle's position <math>X</math> and velocity <math>V</math> may be integrated with a [[symplectic integrator]] method such as [[Verlet integration]]. The time evolution of <math>X</math> and <math>V</math> is called a trajectory. Given the initial positions (e.g., from theoretical knowledge) and velocities (e.g., randomized [[Gaussian]]), we can calculate all future (or past) positions and velocities.

One frequent source of confusion is the meaning of [[temperature]] 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 the ''instantaneous temperature'', which is found by equating the kinetic energy of the system to ''nk<sub>B</sub>T''/2, where ''n'' 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 100 [[Electronvolt|eV]]. In the real world, the 100 eV from the deposited atom would rapidly be transported through and shared among a large number of atoms (<math>10^{10}</math> 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.

=== Canonical ensemble (NVT) ===

In the [[canonical 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 a [[thermostat]].

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 the [[canonical ensemble]]. Popular methods to control temperature include velocity rescaling, the [[Nosé–Hoover thermostat]], Nosé–Hoover chains, the [[Berendsen thermostat]], the [[Andersen thermostat]] and [[Langevin dynamics]]. The Berendsen thermostat might introduce the [[flying 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.

=== Isothermal–isobaric (NPT) ensemble ===

In the [[isothermal–isobaric ensemble]], amount of substance (N), pressure (P) and temperature (T) are conserved. In addition to a thermostat, a [[barostat]] is needed. It corresponds most closely to laboratory conditions with a flask open to ambient temperature and pressure.

In the simulation of [[Biological membrane|biological membranes]], [[isotropic]] pressure control is not appropriate. For [[Lipid bilayer|lipid bilayers]], pressure control occurs under constant membrane area (NPAT) or constant surface tension "gamma" (NPγT).

=== Generalized ensembles ===

The [[replica 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<ref>{{cite journal | vauthors = Sugita Y, Okamoto Y |title=Replica-exchange molecular dynamics method for protein folding |journal=Chemical Physics Letters |date=November 1999 |volume=314 |issue=1–2 |pages=141–151 |doi=10.1016/S0009-2614(99)01123-9 |bibcode= 1999CPL...314..141S }}</ref> tries to overcome the multiple-minima problem by exchanging the temperature of non-interacting replicas of the system running at several temperatures.

== Potentials in MD simulations ==
{{Main|Interatomic potential|Force field (chemistry)|l2=Force field|Comparison of force field implementations}}
A molecular dynamics simulation requires the definition of a [[potential 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 a [[force field (chemistry)|force field]] and in materials physics as an [[interatomic potential]]. Potentials may be defined at many levels of physical accuracy; those most commonly used in chemistry are based on [[molecular mechanics]] and embody a [[classical mechanics]] treatment of particle-particle interactions that can reproduce structural and [[conformational change]]s but usually cannot reproduce [[chemical reaction]]s.

The reduction from a fully quantum description to a classical potential entails two main approximations. The first one is the [[Born–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 on [[quantum mechanics]] are used; some methods attempt to create hybrid [[QM/MM|classical/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 ===

Empirical potentials used in chemistry are frequently called [[Force field (chemistry)|force fields]], while those used in materials physics are called [[interatomic potential]]s.

Most [[Force field (chemistry)|force fields]] in chemistry are empirical and consist of a summation of bonded forces associated with [[chemical bond]]s, bond angles, and bond [[dihedral angle|dihedrals]], and non-bonded forces associated with [[van der Waals force]]s and [[electrostatic charge]].<ref>{{cite journal | vauthors = Rizzuti B | title = Molecular simulations of proteins: From simplified physical interactions to complex biological phenomena | journal =   Biochimica et Biophysica Acta (BBA) - Proteins and Proteomics| year = 2022 | volume = 1870 | issue = 3 | pages = 140757 | pmid = 35051666 | doi = 10.1016/j.bbapap.2022.140757 | s2cid = 263455009 }}</ref> Empirical potentials represent quantum-mechanical effects in a limited way through ad hoc functional approximations. These potentials contain free parameters such as [[electrostatic charge|atomic charge]], van der Waals parameters reflecting estimates of [[atomic radius]], and equilibrium [[bond length]], angle, and dihedral; these are obtained by fitting against detailed electronic calculations (quantum chemical simulations) or experimental physical properties such as [[Young's modulus|elastic constants]], lattice parameters and [[spectroscopy|spectroscopic]] 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 field (chemistry)|force fields]] employ numerical approximations such as shifted cutoff radii, [[reaction field method|reaction field]] algorithms, particle mesh [[Ewald summation]], or the newer particle–particle-particle–mesh ([[P3M]]).

Chemistry force fields commonly employ preset bonding arrangements (an exception being ''[[ab initio quantum chemistry methods|ab 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 the [[bond order potential|bond order formalism]] can describe several different coordinations of a system and bond breaking.<ref>{{cite journal | vauthors = Sinnott SB, Brenner DW |author-link1=Susan Sinnott|year= 2012 |title= Three decades of many-body potentials in materials research |journal= MRS Bulletin |volume= 37 |issue= 5 |pages= 469–473 |doi=10.1557/mrs.2012.88|doi-access= free |bibcode=2012MRSBu..37..469S }}</ref><ref>{{cite journal | vauthors = Albe K, Nordlund K, Averback RS |year= 2002 |title= Modeling metal-semiconductor interaction: Analytical bond-order potential for platinum-carbon |journal= Phys. Rev. B |volume= 65 |issue= 19 |page= 195124 |doi=10.1103/physrevb.65.195124|bibcode= 2002PhRvB..65s5124A}}</ref> Examples of such potentials include the [[Brenner potential]]<ref>{{cite journal | vauthors = Brenner DW | title = Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films | journal = Physical Review B | volume = 42 | issue = 15 | pages = 9458–9471 | date = November 1990 | pmid = 9995183 | doi = 10.1103/physrevb.42.9458 | url = https://apps.dtic.mil/sti/pdfs/ADA230023.pdf | url-status = live | bibcode = 1990PhRvB..42.9458B | archive-url = https://web.archive.org/web/20170922092328/http://www.dtic.mil/get-tr-doc/pdf?AD=ADA230023 | archive-date = September 22, 2017 }}</ref> for hydrocarbons and its
further developments for the C-Si-H<ref>{{cite journal|doi=10.1080/01418619608240734|title= Empirical potentials for C-Si-H systems with application to C<sub>60</sub> interactions with Si crystal surfaces|journal= Philosophical Magazine A|volume= 74|issue= 6|pages= 1439–1466|year= 1996| vauthors = Beardmore K, Smith R |bibcode= 1996PMagA..74.1439B}}</ref> and C-O-H<ref>{{cite journal|title=A reactive empirical bond order (rebo) potential for hydrocarbon oxygen interactions|journal=Journal of Physics: Condensed Matter|volume=16|issue=41|pages=7261–7275|doi=10.1088/0953-8984/16/41/008|year=2004| vauthors = Ni B, Lee KH, Sinnott SB |bibcode= 2004JPCM...16.7261N|s2cid=250760409 }}</ref> systems. The [[ReaxFF]] potential<ref>{{cite journal | vauthors = Van Duin AC, Dasgupta S, Lorant F, Goddard WA |title=ReaxFF: A Reactive Force Field for Hydrocarbons |journal=The Journal of Physical Chemistry A |date=October 2001 |volume=105 |issue=41 |pages=9396–9409 |doi=10.1021/jp004368u |bibcode= 2001JPCA..105.9396V |citeseerx=10.1.1.507.6992}}</ref> can be considered a fully reactive hybrid between bond order potentials and chemistry force fields.

=== Pair potentials versus many-body potentials ===

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 popular [[force field (chemistry)|force 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-bonded [[Lennard-Jones potential]] (also termed the 6–12 potential), used for calculating van der Waals forces.

:<math>
U(r) = 4\varepsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right]
</math>

Another example is the Born (ionic) model of the ionic lattice. The first term in the next equation is [[Coulomb'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.<ref>{{cite journal | vauthors = Cruz FJ, Lopes JN, Calado JC, Minas da Piedade ME | title = A molecular dynamics study of the thermodynamic properties of calcium apatites. 1. Hexagonal phases | journal = The Journal of Physical Chemistry B | volume = 109 | issue = 51 | pages = 24473–24479 | date = December 2005 | pmid = 16375450 | doi = 10.1021/jp054304p }}</ref><ref>{{cite journal | vauthors = Cruz FJ, Lopes JN, Calado JC |title=Molecular dynamics simulations of molten calcium hydroxyapatite |journal=Fluid Phase Equilibria |date=March 2006 |volume=241 |issue=1–2 |pages=51–58 |doi=10.1016/j.fluid.2005.12.021 |bibcode=2006FlPEq.241...51C }}</ref> When ''n<sub>l</sub>'' = 6, this potential is also called the [[Buckingham potential|Coulomb–Buckingham potential]].

:<math>U_{ij}(r_{ij}) = \frac {z_i z_j}{4 \pi \epsilon_0} \frac {1}{r_{ij}} + A_l \exp \frac {-r_{ij}}{p_l} + C_l r_{ij}^{-n_l} + \cdots
</math>

In [[Many-body problem|many-body potentials]], the potential energy includes the effects of three or more particles interacting with each other.<ref name="ReferenceA">{{cite journal| vauthors = Justo JF, Bazant MZ, Kaxiras E, Bulatov VV, Yip S |title=Interatomic potential for silicon defects and disordered phases|journal=Phys. Rev. B|date=1998|volume=58|issue=5|pages=2539–2550|doi=10.1103/PhysRevB.58.2539|arxiv= cond-mat/9712058 |bibcode= 1998PhRvB..58.2539J|s2cid=14585375}}</ref> 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, the [[Tersoff potential]],<ref>{{cite journal | vauthors = Tersoff J | title = Modeling solid-state chemistry: Interatomic potentials for multicomponent systems | journal = Physical Review B| volume = 39 | issue = 8 | pages = 5566–5568 | date = March 1989 | pmid = 9948964 | doi = 10.1103/physrevb.39.5566 | bibcode = 1989PhRvB..39.5566T }}</ref> which was originally used to simulate [[carbon]], [[silicon]], and [[germanium]], 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 the [[embedded atom model|embedded-atom method]] (EAM),<ref>{{cite journal | vauthors = Daw MS, Foiles SM, Baskes MI |author1-link=Murray S. Daw|author2-link=Stephen M. Foiles|author3-link=Michael Baskes |title=The embedded-atom method: a review of theory and applications |journal=Materials Science Reports |date=March 1993 |volume=9 |issue=7–8 |pages=251–310 |doi=10.1016/0920-2307(93)90001-U |url=https://zenodo.org/record/1258631 |doi-access=free }}</ref> the EDIP,<ref name="ReferenceA"/> and the Tight-Binding Second Moment Approximation (TBSMA) potentials,<ref>{{cite journal | vauthors = Cleri F, Rosato V | title = Tight-binding potentials for transition metals and alloys | journal = Physical Review B | volume = 48 | issue = 1 | pages = 22–33 | date = July 1993 | pmid = 10006745 | doi = 10.1103/physrevb.48.22 | bibcode = 1993PhRvB..48...22C }}</ref> 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 ===

[[Semi-empirical quantum chemistry methods|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, termed [[Tight binding (physics)|tight-binding]] potentials, which vary according to the atoms being modeled.

=== Polarizable potentials ===

{{Main|Force field (chemistry)|l1=Force field}}
Most classical force fields implicitly include the effect of [[polarizability]], 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 as [[Drude particle]]s 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.<ref name="Lamoureux3">{{cite journal | vauthors = Lamoureux G, Harder E, Vorobyov IV, Roux B, MacKerell AD |year=2006 |title=A polarizable model of water for molecular dynamics simulations of biomolecules |journal=Chem Phys Lett |volume=418 |issue=1 |pages=245–249 |doi=10.1016/j.cplett.2005.10.135|bibcode= 2006CPL...418..245L }}</ref><ref name="Sokhan2015">{{cite journal | vauthors = Sokhan VP, Jones AP, Cipcigan FS, Crain J, Martyna GJ | title = Signature properties of water: Their molecular electronic origins | journal = Proceedings of the National Academy of Sciences of the United States of America | volume = 112 | issue = 20 | pages = 6341–6346 | date = May 2015 | pmid = 25941394 | pmc = 4443379 | doi = 10.1073/pnas.1418982112 | doi-access = free | bibcode = 2015PNAS..112.6341S }}</ref><ref name="Cipcigan">{{cite journal | vauthors = Cipcigan FS, Sokhan VP, Jones AP, Crain J, Martyna GJ | title = Hydrogen bonding and molecular orientation at the liquid-vapour interface of water | journal = Physical Chemistry Chemical Physics | volume = 17 | issue = 14 | pages = 8660–8669 | date = April 2015 | pmid = 25715668 | doi = 10.1039/C4CP05506C | doi-access = free | bibcode = 2015PCCP...17.8660C | hdl = 20.500.11820/0bd0cd1a-94f1-4053-809c-9fb68bbec1c9 | hdl-access = free }}</ref> Some promising results have also been achieved for proteins.<ref>{{cite journal | vauthors = Mahmoudi M, Lynch I, Ejtehadi MR, Monopoli MP, Bombelli FB, Laurent S | title = Protein-nanoparticle interactions: opportunities and challenges | journal = Chemical Reviews | volume = 111 | issue = 9 | pages = 5610–5637 | date = September 2011 | pmid = 21688848 | doi = 10.1021/cr100440g }}</ref><ref name=Patel2004b>{{cite journal | vauthors = Patel S, Mackerell AD, Brooks CL | title = CHARMM fluctuating charge force field for proteins: II protein/solvent properties from molecular dynamics simulations using a nonadditive electrostatic model | journal = Journal of Computational Chemistry | volume = 25 | issue = 12 | pages = 1504–1514 | date = September 2004 | pmid = 15224394 | doi = 10.1002/jcc.20077 | s2cid = 16741310 | doi-access = free }}</ref> However, it is still uncertain how to best approximate polarizability in a simulation.{{Citation needed|date=April 2009}} 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.<ref>{{cite journal | vauthors = Najla Hosseini A, Lund M, Ejtehadi MR | title = Electronic polarization effects on membrane translocation of anti-cancer drugs | journal = Physical Chemistry Chemical Physics | volume = 24 | issue = 20 | pages = 12281–12292 | date = May 2022 | pmid = 35543365 | doi = 10.1039/D2CP00056C | bibcode = 2022PCCP...2412281N | s2cid = 248696332 }}</ref>

=== Potentials in ''ab initio'' methods ===
{{Main|Quantum chemistry|List of quantum chemistry and solid state physics software}}

In classical molecular dynamics, one potential energy surface (usually the ground state) is represented in the force field. This is a consequence of the [[Born–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 as [[density functional theory]]. This is named ''Ab 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 initio]]'' [[quantum mechanical]] and [[Quantum chemistry|chemical]] methods may be used to calculate the [[potential energy surface|potential energy]] of a system on the fly, as needed for conformations in a trajectory. This calculation is usually made in the close neighborhood of the [[reaction 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 using ''ab 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 like [[mixed quantum-classical dynamics]].

=== Hybrid QM/MM ===
{{Main|QM/MM}}
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).<ref>The methodology for such methods was introduced by Warshel and coworkers. In the recent years have been pioneered by several groups including: [[Arieh Warshel]] ([[University of Southern California]]), Weitao Yang ([[Duke University]]), Sharon Hammes-Schiffer ([[The Pennsylvania State University]]), Donald Truhlar and Jiali Gao ([[University of Minnesota]]) and Kenneth Merz ([[University of Florida]]).</ref>

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(n<sup>2</sup>), 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(n<sup>2</sup>). 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 simplest ''ab initio'' calculations typically scale O(n<sup>3</sup>) or worse (restricted [[Hartree–Fock]] calculations have been suggested to scale ~O(n<sup>2.7</sup>)). 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 liver [[alcohol dehydrogenase]]. In this case, [[quantum tunneling]] is important for the hydrogen, as it determines the reaction rate.<ref>{{cite journal | vauthors = Billeter SR, Webb SP, Agarwal PK, Iordanov T, Hammes-Schiffer S | title = Hydride transfer in liver alcohol dehydrogenase: quantum dynamics, kinetic isotope effects, and role of enzyme motion | journal = Journal of the American Chemical Society | volume = 123 | issue = 45 | pages = 11262–11272 | date = November 2001 | pmid = 11697969 | doi = 10.1021/ja011384b }}</ref>

=== Coarse-graining and reduced representations ===

At the other end of the detail scale are [[Coarse-grained modeling|coarse-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 called [[Coarse-grained modeling|coarse-grained models]].<ref name=":0">{{cite journal | vauthors = Kmiecik S, Gront D, Kolinski M, Wieteska L, Dawid AE, Kolinski A | title = Coarse-Grained Protein Models and Their Applications | journal = Chemical Reviews | volume = 116 | issue = 14 | pages = 7898–7936 | date = July 2016 | pmid = 27333362 | doi = 10.1021/acs.chemrev.6b00163 | doi-access = free }}</ref>

Examples for coarse graining (CG) methods are discontinuous molecular dynamics (CG-DMD)<ref>{{cite journal | vauthors = Voegler Smith A, Hall CK | title = alpha-helix formation: discontinuous molecular dynamics on an intermediate-resolution protein model | journal = Proteins | volume = 44 | issue = 3 | pages = 344–360 | date = August 2001 | pmid = 11455608 | doi = 10.1002/prot.1100 | s2cid = 21774752 }}</ref><ref>{{cite journal | vauthors = Ding F, Borreguero JM, Buldyrey SV, Stanley HE, Dokholyan NV | title = Mechanism for the alpha-helix to beta-hairpin transition | journal = Proteins | volume = 53 | issue = 2 | pages = 220–228 | date = November 2003 | pmid = 14517973 | doi = 10.1002/prot.10468 | s2cid = 17254380 }}</ref> and Go-models.<ref>{{cite journal | vauthors = Paci E, Vendruscolo M, Karplus M | title = Validity of Gō models: comparison with a solvent-shielded empirical energy decomposition | journal = Biophysical Journal | volume = 83 | issue = 6 | pages = 3032–3038 | date = December 2002 | pmid = 12496075 | pmc = 1302383 | doi = 10.1016/S0006-3495(02)75308-3 | bibcode = 2002BpJ....83.3032P }}</ref> 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.<ref>{{cite journal | vauthors = Chakrabarty A, Cagin T |title=Coarse grain modeling of polyimide copolymers |journal=Polymer |date=May 2010 |volume=51 |issue=12 |pages=2786–2794 |doi=10.1016/j.polymer.2010.03.060 }}</ref> 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 both [[enthalpy|enthalpic]] and [[entropy|entropic]] contributions to free energy in an implicit way.<ref>{{cite journal | vauthors = Foley TT, Shell MS, Noid WG | title = The impact of resolution upon entropy and information in coarse-grained models | journal = The Journal of Chemical Physics | volume = 143 | issue = 24 | pages = 243104 | date = December 2015 | pmid = 26723589 | doi = 10.1063/1.4929836 | bibcode = 2015JChPh.143x3104F }}</ref> 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:
* [[protein folding]] and [[protein structure prediction]] studies are often carried out using one, or a few, pseudo-atoms per amino acid;<ref name=":0" />
* [[liquid crystal]] phase transitions have been examined in confined geometries and/or during flow using the [[Gay-Berne potential]], which describes anisotropic species;
* [[Polymer]] glasses during deformation have been studied using simple harmonic or [[FENE]] springs to connect spheres described by the [[Lennard-Jones potential]];
* [[Supercoiling|DNA supercoiling]] has been investigated using 1–3 pseudo-atoms per basepair, and at even lower resolution;
* Packaging of [[DNA|double-helical DNA]] into [[bacteriophage]] has been investigated with models where one pseudo-atom represents one turn (about 10 basepairs) of the double helix;
* RNA structure in the [[ribosome]] and other large systems has been modeled with one pseudo-atom per nucleotide.

The simplest form of coarse-graining is the ''united atom'' (sometimes called ''extended 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 CH<sub>3</sub> methyl group explicitly (or all three atoms of CH<sub>2</sub> 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 the [[CHARMM]] 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 ===
Machine Learning Force Fields] (MLFFs) represent one approach to modeling interatomic interactions in molecular dynamics simulations.<ref name="Unke_2021">{{cite journal | vauthors = Unke OT, Chmiela S, Sauceda HE, Gastegger M, Poltavsky I, Schütt KT, Tkatchenko A, Müller KR | title = Machine Learning Force Fields | journal = Chemical Reviews | volume = 121 | issue = 16 | pages = 10142–10186 | date = August 2021 | pmid = 33705118 | pmc = 8391964 | doi = 10.1021/acs.chemrev.0c01111 }}</ref> MLFFs can achieve accuracy close to that of [[Ab initio quantum chemistry methods|ab 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, including [[Vienna Ab initio Simulation Package|VASP]]<ref name="Hafner_2008">{{cite journal | vauthors = Hafner J | title = Ab-initio simulations of materials using VASP: Density-functional theory and beyond | journal = Journal of Computational Chemistry | volume = 29 | issue = 13 | pages = 2044–78 | date = October 2008 | pmid = 18623101 | doi = 10.1002/jcc.21057 }}</ref> and open-source libraries like DeePMD-kit<ref name = "Wang_2018">{{cite journal | vauthors = Wang H, Zhang L, Han J, Weinan E |title=DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics |journal=Computer Physics Communications |date=July 2018 |volume=228 |pages=178–184 |doi=10.1016/j.cpc.2018.03.016|arxiv=1712.03641 }}</ref><ref name="Zeng_2023">{{cite journal | vauthors = Zeng J, Zhang D, Lu D, Mo P, Li Z, Chen Y, Rynik M, Huang L, Li Z, Shi S, Wang Y, Ye H, Tuo P, Yang J, Ding Y, Li Y, Tisi D, Zeng Q, Bao H, Xia Y, Huang J, Muraoka K, Wang Y, Chang J, Yuan F, Bore SL, Cai C, Lin Y, Wang B, Xu J, Zhu JX, Luo C, Zhang Y, Goodall RE, Liang W, Singh AK, Yao S, Zhang J, Wentzcovitch R, Han J, Liu J, Jia W, York DM, E W, Car R, Zhang L, Wang H | title = DeePMD-kit v2: A software package for deep potential models | journal = The Journal of Chemical Physics | volume = 159 | issue = 5 | pages = | date = August 2023 | pmid = 37526163 | pmc = 10445636 | doi = 10.1063/5.0155600 }}</ref> and [https://schnetpack.readthedocs.io/en/latest/ SchNetPack].<ref name="Schütt_2019">{{cite journal | vauthors = Schütt KT, Kessel P, Gastegger M, Nicoli KA, Tkatchenko A, Müller KR | title = SchNetPack: A Deep Learning Toolbox For Atomistic Systems | journal = Journal of Chemical Theory and Computation | volume = 15 | issue = 1 | pages = 448–455 | date = January 2019 | pmid = 30481453 | doi = 10.1021/acs.jctc.8b00908 | arxiv = 1809.01072 }}</ref><ref name="Schütt_2023">{{cite journal | vauthors = Schütt KT, Hessmann SS, Gebauer NW, Lederer J, Gastegger M | title = SchNetPack 2.0: A neural network toolbox for atomistic machine learning | journal = The Journal of Chemical Physics | volume = 158 | issue = 14 | pages = 144801 | date = April 2023 | pmid = | doi = 10.1063/5.0138367  | arxiv = 2212.05517 }}</ref>

== Incorporating solvent effects ==
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.<ref>{{Cite book|title=Molecular Modelling: Principles and Applications| vauthors = Leach A |date=30 January 2001|publisher=Prentice Hall|isbn=9780582382107|edition= 2nd|location=Harlow|language=en|id={{ASIN|0582382106|country=uk}} }}</ref> 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 a [[Potential of mean force|potential 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.<ref>{{Cite book | vauthors = Leach AR |url=https://www.worldcat.org/oclc/45008511 |title=Molecular modelling : principles and applications |date=2001 |publisher=Prentice Hall |isbn=0-582-38210-6 |edition=2nd |location=Harlow, England |pages=320 |oclc=45008511}}</ref>

== Long-range forces ==
A long range interaction is an interaction in which the spatial interaction falls off no faster than <math>r^{-d}</math> where <math>d</math> 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.<ref>{{Cite book|title=Computer Simulation of Liquids| vauthors = Allen MP, Tildesley DJ |date=2017-08-22|publisher=Oxford University Press|isbn=9780198803201|edition= 2nd|pages=216|language=en|id={{ASIN|0198803206|country=uk}} }}</ref>

== Steered molecular dynamics (SMD) ==

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.<ref name="Nienhaus">{{cite book |title=Protein-ligand interactions: methods and applications |url=https://archive.org/details/proteinligandint00gulr |url-access=registration | vauthors = Nienhaus GU |year=2005 |pages=[https://archive.org/details/proteinligandint00gulr/page/54 54–56] |publisher=Humana Press |isbn=978-1-61737-525-5}}</ref>

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 the [[potential of mean force]].<ref name="Leszczyński">{{cite book |title=Computational chemistry: reviews of current trends, Volume 9 | vauthors = Leszczyński J |year=2005 |pages=54–56 | publisher = World Scientific |isbn=978-981-256-742-0}}</ref> A popular method of computing PMF is through the weighted histogram analysis method (WHAM), which analyzes a series of umbrella sampling simulations.<ref>{{cite journal | vauthors = Kumar S, Rosenberg JM, Bouzida D, Swendsen RH, Kollman PA |title=The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method |journal=Journal of Computational Chemistry |date=October 1992 |volume=13 |issue=8 |pages=1011–1021 |doi=10.1002/jcc.540130812 |s2cid=8571486 }}</ref><ref>{{cite journal | vauthors = Bartels C |title=Analyzing biased Monte Carlo and molecular dynamics simulations |journal=Chemical Physics Letters |date=December 2000 |volume=331 |issue=5–6 |pages=446–454 |doi=10.1016/S0009-2614(00)01215-X |bibcode=2000CPL...331..446B }}</ref>

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,<ref>{{cite journal | vauthors = Lemkul JA, Bevan DR | title = Assessing the stability of Alzheimer's amyloid protofibrils using molecular dynamics | journal = The Journal of Physical Chemistry B | volume = 114 | issue = 4 | pages = 1652–1660 | date = February 2010 | pmid = 20055378 | doi = 10.1021/jp9110794 }}</ref> to study the protein ligand interaction in cyclin-dependent kinase 5<ref>{{cite journal | vauthors = Patel JS, Berteotti A, Ronsisvalle S, Rocchia W, Cavalli A | title = Steered molecular dynamics simulations for studying protein-ligand interaction in cyclin-dependent kinase 5 | journal = Journal of Chemical Information and Modeling | volume = 54 | issue = 2 | pages = 470–480 | date = February 2014 | pmid = 24437446 | doi = 10.1021/ci4003574 }}</ref> and even to show the effect of electric field on thrombin (protein) and aptamer (nucleotide) complex<ref>{{cite journal | vauthors = Gosai A, Ma X, Balasubramanian G, Shrotriya P | title = Electrical Stimulus Controlled Binding/Unbinding of Human Thrombin-Aptamer Complex | journal = Scientific Reports | volume = 6 | issue = 1 | pages = 37449 | date = November 2016 | pmid = 27874042 | pmc = 5118750 | doi = 10.1038/srep37449 | bibcode = 2016NatSR...637449G }}</ref> among many other interesting studies.

== Examples of applications ==
[[File:MD rotor 250K 1ns.gif|thumb|upright=1.15|Molecular dynamics simulation of a [[synthetic molecular motor]] composed of three molecules in a nanopore (outer diameter 6.7 nm) at 250 K<ref>{{cite journal | vauthors = Palma CA, Björk J, Rao F, Kühne D, Klappenberger F, Barth JV | title = Topological dynamics in supramolecular rotors | journal = Nano Letters | volume = 14 | issue = 8 | pages = 4461–4468 | date = August 2014 | pmid = 25078022 | doi = 10.1021/nl5014162 | bibcode = 2014NanoL..14.4461P }}</ref>]]{{Primary sources|section|date=January 2024}}
Molecular dynamics is used in many fields of science.
* First MD simulation of a simplified biological folding process was published in 1975. Its simulation published in Nature paved the way for the vast area of modern computational protein-folding.<ref>{{cite journal | vauthors = Levitt M, Warshel A | title = Computer simulation of protein folding | journal = Nature | volume = 253 | issue = 5494 | pages = 694–698 | date = February 1975 | pmid = 1167625 | doi = 10.1038/253694a0 | s2cid = 4211714 | bibcode = 1975Natur.253..694L }}</ref>
* First MD simulation of a biological process was published in 1976. Its simulation published in Nature paved the way for understanding protein motion as essential in function and not just accessory.<ref>{{cite journal | vauthors = Warshel A | title = Bicycle-pedal model for the first step in the vision process | journal = Nature | volume = 260 | issue = 5553 | pages = 679–683 | date = April 1976 | pmid = 1264239 | doi = 10.1038/260679a0 | s2cid = 4161081 | bibcode = 1976Natur.260..679W }}</ref>
* MD is the standard method to treat [[collision cascade]]s in the heat spike regime, i.e., the effects that energetic [[neutron]] and [[ion irradiation]] have on solids and solid surfaces.<ref name="Smith">{{cite book |editor=Smith, R. |title=Atomic & ion collisions in solids and at surfaces: theory, simulation and applications |publisher=Cambridge University Press |location=Cambridge, UK |year=1997}}{{page needed|date=April 2020}}</ref>

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):
* MD simulation of the full ''[[satellite tobacco mosaic virus]]'' (STMV) (2006, Size: 1 million atoms, Simulation time: 50 ns, program: [[NAMD]]) This virus is a small, icosahedral plant virus that worsens the symptoms of infection by Tobacco Mosaic Virus (TMV). Molecular dynamics simulations were used to probe the mechanisms of [[Virus#Structure|viral assembly]]. The entire STMV particle consists of 60 identical copies of one protein that make up the viral [[capsid]] (coating), and a 1063 nucleotide single stranded RNA [[genome]]. One key finding is that the capsid is very unstable when there is no RNA inside. The simulation would take one 2006 desktop computer around 35 years to complete. It was thus done in many processors in parallel with continuous communication between them.<ref>{{cite web |url=http://www.ks.uiuc.edu/Research/STMV/ |title=Molecular dynamics simulation of the Satellite Tobacco Mosaic Virus (STMV) |vauthors=Freddolino P, Arkhipov A, Larson SB, McPherson A, Schulten K |work=Theoretical and Computational Biophysics Group |publisher=University of Illinois at Urbana Champaign}}</ref>
* Folding simulations of the [[Villin]] Headpiece in all-atom detail (2006, Size: 20,000 atoms; Simulation time: 500 μs= 500,000 ns, Program: [[Folding@home]]) This simulation was run in 200,000 CPU's of participating personal computers around the world. These computers had the Folding@home program installed, a large-scale distributed computing effort coordinated by [[Vijay Pande]] at Stanford University. The kinetic properties of the Villin Headpiece protein were probed by using many independent, short trajectories run by CPU's without continuous real-time communication. One method employed was the Pfold value analysis, which measures the probability of folding before unfolding of a specific starting conformation. Pfold gives information about [[Phi value analysis|transition state]] structures and an ordering of conformations along the [[Protein folding|folding pathway]]. Each trajectory in a Pfold calculation can be relatively short, but many independent trajectories are needed.<ref>{{cite journal | vauthors = Jayachandran G, Vishal V, Pande VS | title = Using massively parallel simulation and Markovian models to study protein folding: examining the dynamics of the villin headpiece | journal = The Journal of Chemical Physics | volume = 124 | issue = 16 | pages = 164902 | date = April 2006 | pmid = 16674165 | doi = 10.1063/1.2186317 | bibcode = 2006JChPh.124p4902J | doi-access = free }}</ref>
* Long continuous-trajectory simulations have been performed on [[Anton (computer)|Anton]], a massively parallel supercomputer designed and built around custom [[application-specific integrated circuit]]s (ASICs) and interconnects by [[D. E. Shaw Research]]. The longest published result of a simulation performed using Anton is a 1.112-millisecond simulation of NTL9 at 355 K; a second, independent 1.073-millisecond simulation of this configuration was also performed (and many other simulations of over 250 μs continuous chemical time).<ref name="DESRES-Science2011">{{cite journal | vauthors = Lindorff-Larsen K, Piana S, Dror RO, Shaw DE | title = How fast-folding proteins fold | journal = Science | volume = 334 | issue = 6055 | pages = 517–520 | date = October 2011 | pmid = 22034434 | doi = 10.1126/science.1208351 | s2cid = 27988268 | citeseerx = 10.1.1.1013.9290 | bibcode = 2011Sci...334..517L }}</ref> In ''How Fast-Folding Proteins Fold'', researchers Kresten Lindorff-Larsen, Stefano Piana, Ron O. Dror, and [[David E. Shaw]] discuss "the results of atomic-level molecular dynamics simulations, over periods ranging between 100 μs and 1 ms, that reveal a set of common principles underlying the folding of 12 structurally diverse proteins." Examination of these diverse long trajectories, enabled by specialized, custom hardware, allow them to conclude that "In most cases, folding follows a single dominant route in which elements of the native structure appear in an order highly correlated with their propensity to form in the unfolded state."<ref name="DESRES-Science2011" /> In a separate study, Anton was used to conduct a 1.013-millisecond simulation of the native-state dynamics of bovine pancreatic trypsin inhibitor (BPTI) at 300 K.<ref>{{cite journal | vauthors = Shaw DE, Maragakis P, Lindorff-Larsen K, Piana S, Dror RO, Eastwood MP, Bank JA, Jumper JM, Salmon JK, Shan Y, Wriggers W | title = Atomic-level characterization of the structural dynamics of proteins | journal = Science | volume = 330 | issue = 6002 | pages = 341–346 | date = October 2010 | pmid = 20947758 | doi = 10.1126/science.1187409 | s2cid = 3495023 | bibcode = 2010Sci...330..341S }}</ref>

Another important application of MD method benefits from its ability of 3-dimensional characterization and analysis of microstructural evolution at atomic scale.
* MD simulations are used in characterization of grain size evolution, for example, when describing wear and friction of nanocrystalline Al and Al(Zr) materials.<ref>{{cite journal | vauthors = Shi Y, Szlufarska I |title=Wear-induced microstructural evolution of nanocrystalline aluminum and the role of zirconium dopants |journal=Acta Materialia |date=November 2020 |volume=200 |pages=432–441 |doi=10.1016/j.actamat.2020.09.005 |bibcode=2020AcMat.200..432S |s2cid=224954349 |doi-access=free }}</ref> Dislocations evolution and grain size evolution are analyzed during the friction process in this simulation. Since MD method provided the full information of the microstructure, the grain size evolution was calculated in 3D using the Polyhedral Template Matching,<ref>{{cite journal | vauthors = Larsen PM, Schmidt S, Schiøtz J |title=Robust structural identification via polyhedral template matching |journal=Modelling and Simulation in Materials Science and Engineering |date=1 June 2016 |volume=24 |issue=5 |pages=055007 |doi=10.1088/0965-0393/24/5/055007|arxiv=1603.05143 |bibcode=2016MSMSE..24e5007M |s2cid=53980652 }}</ref> Grain Segmentation,<ref>{{cite journal | vauthors = Hoffrogge PW, Barrales-Mora LA |title=Grain-resolved kinetics and rotation during grain growth of nanocrystalline Aluminium by molecular dynamics |journal=Computational Materials Science |date=February 2017 |volume=128 |pages=207–222 |doi=10.1016/j.commatsci.2016.11.027|arxiv=1608.07615 |s2cid=118371554 }}</ref> and Graph clustering<ref>{{cite arXiv | vauthors = Bonald T, Charpentier B, Galland A, Hollocou A |title=Hierarchical Graph Clustering using Node Pair Sampling |date=22 June 2018 |class=cs.SI |eprint=1806.01664}}</ref> methods. In such simulation, MD method provided an accurate measurement of grain size. Making use of these information, the actual grain structures were extracted, measured, and presented. Compared to the traditional method of using SEM with a single 2-dimensional slice of the material, MD provides a 3-dimensional and accurate way to characterize the microstructural evolution at atomic scale.

== Molecular dynamics algorithms ==
* [[Screened Coulomb potentials implicit solvent model]]<!--not sure which section this belongs in-->

=== Integrators ===
* [[Symplectic integrator]]
* [[Verlet integration|Verlet–Stoermer integration]]
* [[Runge–Kutta integration]]
* [[Beeman's algorithm]]
* [[Constraint algorithm]]s (for constrained systems)

=== Short-range interaction algorithms ===
* [[Cell lists]]
* [[Verlet list]]
* Bonded interactions

=== Long-range interaction algorithms ===
* [[Ewald summation]]
* Particle mesh [[Ewald summation]] (PME)
* Particle–particle-particle–mesh ([[P3M]])
* [[Shifted force method]]

=== Parallelization strategies ===
* [[Domain decomposition method]] (Distribution of system data for [[parallel computing]])

=== Ab-initio molecular dynamics ===
* [[Car–Parrinello molecular dynamics]]

== Specialized hardware for MD simulations ==
* [[Anton (computer)|Anton]] – A specialized, massively parallel supercomputer designed to execute MD simulations
* [[MDGRAPE]] – A special purpose system built for molecular dynamics simulations, especially protein structure prediction

== Graphics card as a hardware for MD simulations ==
{{Excerpt|Molecular modeling on GPUs|paragraphs=1-2}}

== See also ==
{{div col|colwidth=27em}}
* [[Molecular modeling]]
* [[Computational chemistry]]
* [[Force field (chemistry)]]
* [[Comparison of force field implementations]]
* [[Monte Carlo method]]
* [[Molecular design software]]
* [[Molecular mechanics]]
* [[Multiscale Green's function]]
* [[Car–Parrinello method]]
* [[Comparison of software for molecular mechanics modeling]]
* [[Quantum chemistry]]
* [[Discrete element method]]
* [[Comparison of nucleic acid simulation software]]
* [[Molecule editor]]
* [[Mixed quantum-classical dynamics]]
{{div col end}}

== References ==
{{Reflist|30em}}

=== General references ===
{{refbegin}}
* {{cite book | vauthors = Allen MP, Tildesley DJ | date = 1989 | title = Computer simulation of liquids | publisher = Oxford University Press | isbn = 0-19-855645-4}}
* {{cite book | vauthors = McCammon JA, Harvey SC | date = 1987 | title = Dynamics of Proteins and Nucleic Acids | publisher = Cambridge University Press | isbn = 0-521-30750-3}}
* {{cite book | vauthors = Rapaport DC | date = 1996 | title = The Art of Molecular Dynamics Simulation | isbn = 0-521-44561-2 }}
* {{cite book |title=Numerical Simulation in Molecular Dynamics |publisher=Springer |location=Berlin, Heidelberg |year=2007 |isbn=978-3-540-68094-9 | vauthors = Griebel M, Knapek S, Zumbusch G |author-link1=Michael Griebel}}
* {{cite book | vauthors = Frenkel D, Smit B |author-link= Daan Frenkel  |title= Understanding Molecular Simulation : from algorithms to applications |publisher= Academic Press |orig-year= 2001 |location= San Diego |isbn= 978-0-12-267351-1 |year= 2002}}
* {{cite book | vauthors = Haile JM | date = 2001 | title = Molecular Dynamics Simulation: Elementary Methods | publisher = Wiley | isbn = 0-471-18439-X }}
* {{cite book | vauthors = Sadus RJ | title = Molecular Simulation of Fluids: Theory, Algorithms and Object-Orientation | date = 2002 | publisher = Elsevier | isbn = 0-444-51082-6 }}
* {{cite book | vauthors = Becker OM, Mackerell Jr AD, Roux B, Watanabe M | date = 2001 | title = Computational Biochemistry and Biophysics | publisher = Marcel Dekker | isbn = 0-8247-0455-X}}
* {{cite book | vauthors = Leach A | date = 2001 | title = Molecular Modelling: Principles and Applications | edition = 2nd | publisher = Prentice Hall | isbn = 978-0-582-38210-7}}
* {{cite book | vauthors = Schlick T | author-link1 = Tamar Schlick | date = 2002 | title = Molecular Modeling and Simulation | publisher = Springer | isbn = 0-387-95404-X}}
* {{cite book | vauthors = Hoover WB | author-link1 = William G Hoover | date = 1991 | title = Computational Statistical Mechanics | publisher = Elsevier | isbn = 0-444-88192-1}}
* {{cite book | vauthors = Evans DJ, Morriss G | date = 2008 | title = Statistical Mechanics of Nonequilibrium Liquids | edition = Second | publisher = Cambridge University Press | isbn = 978-0-521-85791-8}}
{{refend}}

== External links ==
{{Commons category|Molecular dynamics simulation}}
* [http://gpugrid.net/ The GPUGRID.net Project] ([[GPUGRID.net]])
* [https://web.archive.org/web/20110606092842/http://researchweb.watson.ibm.com/bluegene/ The Blue Gene Project] ([[IBM]]) JawBreakers.org
* [http://www.sklogwiki.org/SklogWiki/index.php/Materials_modelling_and_computer_simulation_codes Materials modelling and computer simulation codes]
* [http://www.layruoru.com/dokuwiki/doku.php/ A few tips on molecular dynamics]
* [https://www.youtube.com/watch?v=x8Atqz5YvzQ Movie of MD simulation of water (YouTube)]
{{Branches of chemistry}}
[[Category:Molecular dynamics| ]]
[[Category:Computational chemistry]]
[[Category:Molecular modelling]]
[[Category:Simulation]]