Editing Molecular dynamics (section)

== 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.