7 The pressure in molecular dynamics can be computed from the kinetic energy and
13 Whether or not pressure coupling is used within a simulation, the pressure
14 value for the simulation box will oscillate significantly. Instantaneous
15 pressure is meaningless, and not well-defined. Over a picosecond time scale it
16 usually will not be a good indicator of the true pressure. This variation is
17 entirely normal due to the fact that pressure is a macroscopic property and can
18 only be measured properly as time average, while it is being measured and/or
19 adjusted with pressure coupling on the microscopic scale. How much it varies
20 and the speed at which it does depends on the number of atoms in the system,
21 the type of pressure coupling used and the value of the coupling constants.
22 Fluctuations of the order of hundreds of bar are typical. For a box of 216
23 waters, fluctuations of 500-600 bar are standard. Since the fluctuations go
24 down with the square root of the number of particles, a system of 21600 water
25 molecules (100 times larger) will still have pressure fluctuations of 50-60 bar.
29 Periodic boundary conditions
30 ----------------------------
32 Periodic boundary conditions (PBC) are used in molecular dynamics simulations
33 to avoid problems with boundary effects caused by finite size, and make the
34 system more like an infinite one, at the cost of possible periodicity effects.
36 Beginners visualizing a trajectory sometimes think they are observing a problem
39 * the molecule(s) does not stay in the centre of the box, or
40 * it appears that (parts of) the molecule(s) diffuse out of the box, or
41 * holes are created, or
42 * broken molecules appear, or
43 * their unit cell was a rhombic dodecahedron or cubic octahedron but it looks
44 like a slanted cube after the simulation, or
45 * crazy bonds all across the simulation cell appear.
47 This is not a problem or error that is occuring, it is what you should expect.
49 The existence of PBC means that any atom that leaves a simulation box by, say,
50 the right-hand face, then enters the simulation box by the left-hand face. In
51 the example of a large protein, if you look at the face of the simulation box
52 that is opposite to the one from which the protein is protruding, then a hole
53 in the solvent will be visible. The reason that the molecule(s) move from where
54 they were initially located within the box is (for the vast majority of
55 simulations) they are free to diffuse around. And so they do. They are not held
56 in a magic location of the box. The box is not centered around anything while
57 performing the simulation. Molecules are not made whole as a matter of course.
58 Moreover, any periodic cell shape can be expressed as a parallelepiped (a.k.a.
59 triclinic cell), and |Gromacs| does so internally regardless of the initial
62 These visual issues can be fixed after the conclusion of the simulation by
63 judicious use of the optional inputs to :ref:`gmx trjconv` to process the
64 trajectory files. Similarly, analyses such as RMSD of atomic positions can be
65 flawed when a reference structure is compared with a structure that needs
66 adjusting for periodicity effects, and the solution with :ref:`gmx trjconv`
67 follows the same lines. Some complex cases needing more than one operation will
68 require more than one invocation of :ref:`gmx trjconv` in order to work.
70 For further information, see `Reference Manual`_.
75 Fixing periodicity effects with :ref:`gmx trjconv` to suit visualization or
76 analysis can be tricky. Multiple invocations can be necessary. You may need to
77 create custom index groups (e.g. to keep your ligand with your protein)
78 Following the steps below in order (omitting those not required) should help
79 get a pleasant result. You will need to consult ``gmx trjconv -h`` to find out
80 the details for each step. That's deliberate -- there is no magic "do what I
81 want" recipe. You have to decide what you want, first. :-)
83 #. First make your molecules whole if you want them whole.
84 #. Cluster your molecules/particles if you want them clustered.
85 #. If you want jumps removed, extract the first frame from the trajectory to
86 use as the reference, and then use ``-pbc nojump`` with that first
88 #. Center your system using some criterion. Doing so shifts the system, so
89 don't use ``-pbc nojump`` after this step.
90 #. Perhaps put everything in some box with the other ``-pbc`` or ``-ur``
92 #. Fit the resulting trajectory to some (other) reference structure (if
93 desired), and don't use any PBC related option afterwards.
95 With point three, the issue is that :ref:`gmx trjconv` removes the jumps from
96 the first frame using the reference structure provided with -s. If the reference
97 structure (run input file) is not clustered/whole, using ``-pbc nojump``
98 will undo steps 1 and 2.
105 Thermostats are designed to help a simulation sample from the correct ensemble
106 (i.e. NVT or NPT) by modulating the temperature of the system in some fashion.
107 First, we need to establish what we mean by temperature. In simulations, the
108 "instantaneous (kinetic) temperature" is usually computed from the kinetic
109 energy of the system using the equipartition theorem. In other words, the
110 temperature is computed from the system's total kinetic energy.
112 So, what's the goal of a thermostat? Actually, it turns out the goal is not to
113 keep the temperature constant, as that would mean fixing the total kinetic
114 energy, which would be silly and not the aim of NVT or NPT. Rather, it's to
115 ensure that the average temperature of a system be correct.
117 To see why this is the case, imagine a glass of water sitting in a room.
118 Suppose you can look very closely at a few molecules in some small region of
119 the glass, and measure their kinetic energies. You would not expect the kinetic
120 energy of this small number of particles to remain precisely constant; rather,
121 you'd expect fluctuations in the kinetic energy due to the small number of
122 particles. As you average over larger and larger numbers of particles, the
123 fluctuations in the average get smaller and smaller, so finally by the time you
124 look at the whole glass, you say it has "constant temperature".
126 Molecular dynamics simulations are often fairly small compared to a glass of
127 water, so we have bigger fluctuations. So it's really more appropriate here to
128 think of the role of a thermostat as ensuring that we have
130 (a) the correct average temperature, and
131 (b) the fluctuations of the correct size.
133 See `Reference Manual`_ for details on how temperature coupling is applied and
134 the types currently available.
136 .. _gmx-thermostats-do:
141 Some hints on practices that generally are a good idea:
143 * Preferably, use a thermostat that samples the correct distribution of
144 temperatures (for examples, see the corresponding manual section), in addition
145 to giving you the correct average temperature.
146 * At least: use a thermostat that gives you the correct average temperature,
147 and apply it to components of your system for which they are justified (see
148 the first bullet in `What not to do`_). In some cases, using
149 ``tc-grps = System`` may lead to the "hot solvent/cold solute" problem
150 described in the 3rd reference in `Further reading`_.
152 .. _gmx-thermostats-dont:
157 Some hints on practices that generally not a good idea to use:
159 * Do not use separate thermostats for every component of your system. Some
160 molecular dynamics thermostats only work well in the thermodynamic limit. A
161 group must be of sufficient size to justify its own thermostat. If you use one
162 thermostat for, say, a small molecule, another for protein, and another for
163 water, you are likely introducing errors and artifacts that are hard to
164 predict. In particular, do not couple ions in aqueous solvent in a separate
165 group from that solvent. For a protein simulation, using ``tc-grps = Protein
166 Non-Protein`` is usually best.
167 * Do not use thermostats that work well only in the limit of a large number of
168 degrees of freedom for systems with few degrees of freedom. For example, do
169 not use Nosé-Hoover or Berendsen thermostats for types of free energy
170 calculations where you will have a component of the system with very few
171 degrees of freedom in an end state (i.e. a noninteracting small molecule).
176 #. `Cheng, A. & Merz, K. M. Application of the nosé- hoover chain algorithm to
177 the study of protein dynamics. *J. Phys. Chem.* **100** (5), 1927–1937
178 (1996). <http://pubs.acs.org/doi/abs/10.1021/jp951968y>`__
179 #. `Mor, A., Ziv, G. & Levy, Y. Simulations of proteins with inhomogeneous
180 degrees of freedom: the effect of thermostats. *J. Comput. Chem.* **29**
181 (12), 1992–1998 (2008). <http://dx.doi.org/10.1002/jcc.20951>`__
182 #. `Lingenheil, M., Denschlag, R., Reichold, R. & Tavan, P. The
183 "hot-solvent/cold-solute" problem revisited. *J. Chem. Theory Comput.* **4**
184 (8), 1293–1306 (2008). <http://pubs.acs.org/doi/abs/10.1021/ct8000365>`__
189 In principle, a molecular dynamics simulation should conserve the total energy,
190 the total momentum and (in a non-periodic system) the total angular momentum. A
191 number of algorithmic and numerical issues make that this is not always the
194 * Cut-off treatment and/or long-range electrostatics treatment (see Van Der
195 Spoel, D. & van Maaren, P. J. The origin of layer structure artifacts in
196 simulations of liquid water. *J. Chem. Theor. Comp.* **2**, 1–11 (2006).)
197 * Treatment of pair lists,
198 * Constraint algorithms (see e.g. Hess, B. P-LINCS: A parallel linear constraint
199 solver for molecular simulation. *J. Chem. Theor. Comp.* **4**, 116–122
201 * The integration timestep,
202 * Temperature coupling and pressure coupling,
203 * Round-off error (in particular in single precision), for example subtracting
204 large numbers (Lippert, R. A. et al. A common, avoidable source of error in
205 molecular dynamics integrators. *J. Chem. Phys.* **126**, 046101 (2007).)
206 * The choice of the integration algorithm (in |Gromacs| this is normally
208 * Removal of center of mass motion: when doing this in more than one group the
209 conservation of energy will be violated.
214 Various |Gromacs| utilities can compute average structures. Presumably the idea
215 for this comes from something like an ensemble-average NMR structure. In some
216 cases, it makes sense to calculate an average structure (as a step on the way
217 to calculating root-mean-squared fluctuations (RMSF), for example, one needs
218 the average position of all of the atoms).
220 However, it's important to remember that an average structure isn't necessarily
221 meaningful. By way of analogy, suppose I alternate holding a ball in my left
222 hand, then in my right hand. What's the average position of the ball? Halfway
223 in between -- even though I always have it either in my left hand or my right
224 hand. Similarly, for structures, averages will tend to be meaningless anytime
225 there are separate metastable conformational states. This can happen on a
226 sidechain level, or for some regions of backbone, or even whole helices or
227 components of the secondary structure.
229 Thus, if you derive an average structure from a molecular dynamics simulation,
230 and find artifacts like unphysical bond lengths, weird structures, etc., this
231 doesn't necessarily mean something is wrong. It just shows the above: an
232 average structure from a simulation is not necessarily a physically meaningful
240 *Blowing up* is a highly technical term used to describe a common sort of
241 simulation failure. In brief, it describes a failure typically due to an
242 unacceptably large force that ends up resulting in a failure of the integrator.
244 To give a bit more background, it's important to remember that molecular
245 dynamics numerically integrates Newton's equations of motion by taking small,
246 discrete timesteps, and using these timesteps to determine new velocities and
247 positions from velocities, positions, and forces at the previous timestep. If
248 forces become too large at one timestep, this can result in extremely large
249 changes in velocity/position when going to the next timestep. Typically, this
250 will result in a cascade of errors: one atom experiences a very large force one
251 timestep, and thus goes shooting across the system in an uncontrolled way in
252 the next timestep, overshooting its preferred location or landing on top of
253 another atom or something similar. This then results in even larger forces the
254 next timestep, more uncontrolled motions, and so on. Ultimately, this will
255 cause the simulation package to crash in some way, since it can't cope with
256 such situations. In simulations with constraints, the first symptom of this
257 will usually be some LINCS or SHAKE warning or error -- not because the
258 constraints are the source of the problem, but just because they're the first
259 thing to crash. Similarly, in simulations with domain decomposition, you may
260 see messages about particles being more than a cell length out of the domain
261 decomposition cell of their charge group, which are symptomatic of your
262 underlying problem, and not the domain decomposition algorithm itself. Likewise
263 for warnings about tabulated or 1-4 interactions being outside the distance
264 supported by the table. This can happen on one computer system while another
265 resulted in a stable simulation because of the impossibility of numerical
266 reproducibility of these calculations on different computer systems.
268 Possible causes include:
270 * you didn't minimize well enough,
271 * you have a bad starting structure, perhaps with steric clashes,
272 * you are using too large a timestep (particularly given your choice of
274 * you are doing particle insertion in free energy calculations without using
276 * you are using inappropriate pressure coupling (e.g. when you are not in
277 equilibrium, Berendsen can be best while relaxing the volume, but you will
278 need to switch to a more accurate pressure-coupling algorithm later),
279 * you are using inappropriate temperature coupling, perhaps on inappropriate
281 * your position restraints are to coordinates too different from those present
283 * you have a single water molecule somewhere within the system that is
284 isolated from the other water molecules, or
285 * you are experiencing a bug in :ref:`gmx mdrun`.
287 Because blowing up is due, typically, to forces that are too large for a
288 particular timestep size, there are a couple of basic solutions:
290 * make sure the forces don't get that large, or
291 * use a smaller timestep.
293 Better system preparation is a way to make sure that forces don't get large, if
294 the problems are occurring near the beginning of a simulation.
296 .. _system-diagnosis:
298 Diagnosing an unstable system
299 -----------------------------
301 Troubleshooting a system that is blowing up can be challenging, especially for
302 an inexperienced user. Here are a few general tips that one may find useful
303 when addressing such a scenario:
305 #. If the crash is happening relatively early (within a few steps), set
306 ``nstxout`` (or ``nstxout-compressed``) to 1, capturing all possible frames.
307 Watch the resulting trajectory to see which atoms/residues/molecules become
309 #. Simplify the problem to try to establish a cause:
311 * If you have a new box of solvent, try minimizing and simulating a single
312 molecule to see if the instability is due to some inherent problem with
313 the molecule's topology or if instead there are clashes in your starting
315 * If you have a protein-ligand system, try simulating the protein alone in
316 the desired solvent. If it is stable, simulate the ligand in vacuo to see
317 if its topology gives stable configurations, energies, etc.
318 * Remove the use of fancy algorithms, particularly if you haven't
319 equilibrated thoroughly first
321 #. Monitor various components of the system's energy using :ref:`gmx energy`.
322 If an intramolecular term is spiking, that may indicate improper bonded
323 parameters, for example.
324 #. Make sure you haven't been ignoring error messages (missing atoms when
325 running :ref:`gmx pdb2gmx`, mismatching names when running :ref:`gmx grompp`,
326 etc.) or using work-arounds (like using ``gmx grompp -maxwarn`` when you
327 shouldn't be) to make sure your topology is intact and being interpreted
329 #. Make sure you are using appropriate settings in your :ref:`mdp` file for the
330 force field you have chosen and the type of system you have. Particularly
331 important settings are treatment of cutoffs, proper neighbor searching
332 interval (``nstlist``), and temperature coupling. Improper settings can lead
333 to a breakdown in the model physics, even if the starting configuration of
334 the system is reasonable.
336 When using no explict solvent, starting your equilibration with a smaller time
337 step than your production run can help energy equipartition more stably.
339 There are several common situations in which instability frequently arises,
340 usually in the introduction of new species (ligands or other molecules) into
341 the system. To determine the source of the problem, simplify the system (e.g.
342 the case of a protein-ligand complex) in the following way.
344 #. Does the protein (in water) minimize adequately by itself? This is a test of
345 the integrity of the coordinates and system preparation. If this fails,
346 something probably went wrong when running :ref:`gmx pdb2gmx` (see below), or
347 maybe :ref:`gmx genion` placed an ion very close to the protein (it is
349 #. Does the ligand minimize in vacuo? This is a test of the topology. If it
350 does not, check your parameterization of the ligand and any implementation of
351 new parameters in force field files.
352 #. (If previous item is successful) Does the ligand minimize in water, and/or
353 does a short simulation of the ligand in water succeed?
355 Other sources of possible problems are in the biomolecule topology itself.
357 #. Did you use ``-missing`` when running :ref:`gmx pdb2gmx`? If so, don't.
358 Reconstruct missing coordinates rather than ignoring them.
359 #. Did you override long/short bond warnings by changing the lengths? If so,
360 don't. You probably have missing atoms or some terrible input geometry.
362 .. _reference manual: gmx-manual-parent-dir_
369 .. TODO refer to reference manual definitions
376 Force fields are sets of potential functions and parametrized interactions that can be used to study
377 physical systems. A general introduction to their history, function and use is beyond the scope of this
378 guide, and the user is asked to consult either the relevant literature or
379 try to start at the relevant `Wikipedia page`_.
381 .. _Wikipedia page: https://en.wikipedia.org/wiki/Force_field_(chemistry)