6 Computational Chemistry and Molecular Modeling
7 ----------------------------------------------
9 |Gromacs| is an engine to perform molecular dynamics simulations and
10 energy minimization. These are two of the many techniques that belong to
11 the realm of computational chemistry and molecular modeling.
12 *Computational chemistry* is just a name to indicate the use of
13 computational techniques in chemistry, ranging from quantum mechanics of
14 molecules to dynamics of large complex molecular aggregates. *Molecular
15 modeling* indicates the general process of describing complex chemical
16 systems in terms of a realistic atomic model, with the goal being to
17 understand and predict macroscopic properties based on detailed
18 knowledge on an atomic scale. Often, molecular modeling is used to
19 design new materials, for which the accurate prediction of physical
20 properties of realistic systems is required.
22 Macroscopic physical properties can be distinguished by
24 #. *static equilibrium properties*, such as the binding constant of an
25 inhibitor to an enzyme, the average potential energy of a system, or the
26 radial distribution function of a liquid, and
28 #. *dynamic or non-equilibrium properties*, such as the viscosity of a liquid,
29 diffusion processes in membranes, the dynamics of phase changes,
30 reaction kinetics, or the dynamics of defects in crystals.
33 technique depends on the question asked and on the feasibility of the
34 method to yield reliable results at the present state of the art.
35 Ideally, the (relativistic) time-dependent Schrödinger equation
36 describes the properties of molecular systems with high accuracy, but
37 anything more complex than the equilibrium state of a few atoms cannot
38 be handled at this *ab initio* level. Thus, approximations are
39 necessary; the higher the complexity of a system and the longer the time
40 span of the processes of interest is, the more severe the required
41 approximations are. At a certain point (reached very much earlier than
42 one would wish), the *ab initio* approach must be augmented or replaced
43 by *empirical* parameterization of the model used. Where simulations
44 based on physical principles of atomic interactions still fail due to
45 the complexity of the system, molecular modeling is based entirely on a
46 similarity analysis of known structural and chemical data. The QSAR
47 methods (Quantitative Structure-Activity Relations) and many
48 homology-based protein structure predictions belong to the latter
51 Macroscopic properties are always ensemble averages over a
52 representative statistical ensemble (either equilibrium or
53 non-equilibrium) of molecular systems. For molecular modeling, this has
54 two important consequences:
56 - The knowledge of a single structure, even if it is the structure of
57 the global energy minimum, is not sufficient. It is necessary to
58 generate a representative ensemble at a given temperature, in order
59 to compute macroscopic properties. But this is not enough to compute
60 thermodynamic equilibrium properties that are based on free energies,
61 such as phase equilibria, binding constants, solubilities, relative
62 stability of molecular conformations, etc. The computation of free
63 energies and thermodynamic potentials requires special extensions of
64 molecular simulation techniques.
66 - While molecular simulations, in principle, provide atomic details of
67 the structures and motions, such details are often not relevant for
68 the macroscopic properties of interest. This opens the way to
69 simplify the description of interactions and average over irrelevant
70 details. The science of statistical mechanics provides the
71 theoretical framework for such simplifications. There is a hierarchy
72 of methods ranging from considering groups of atoms as one unit,
73 describing motion in a reduced number of collective coordinates,
74 averaging over solvent molecules with potentials of mean force
75 combined with stochastic dynamics :ref:`9 <refGunsteren90>`, to
76 *mesoscopic dynamics* describing densities rather than atoms and
77 fluxes as response to thermodynamic gradients rather than velocities
78 or accelerations as response to forces \ :ref:`10 <refFraaije93>`.
80 For the generation of a representative equilibrium ensemble two methods
83 #. *Monte Carlo simulations* and
85 #. *Molecular Dynamics simulations*.
87 For the generation of non-equilibrium
88 ensembles and for the analysis of dynamic events, only the second method
89 is appropriate. While Monte Carlo simulations are more simple than MD
90 (they do not require the computation of forces), they do not yield
91 significantly better statistics than MD in a given amount of computer
92 time. Therefore, MD is the more universal technique. If a starting
93 configuration is very far from equilibrium, the forces may be
94 excessively large and the MD simulation may fail. In those cases, a
95 robust *energy minimization* is required. Another reason to perform an
96 energy minimization is the removal of all kinetic energy from the
97 system: if several “snapshots” from dynamic simulations must be
98 compared, energy minimization reduces the thermal noise in the
99 structures and potential energies so that they can be compared better.
101 Molecular Dynamics Simulations
102 ------------------------------
104 MD simulations solve Newton’s equations of motion for a system of
105 :math:`N` interacting atoms:
107 .. math:: m_i \frac{\partial^2 \mathbf{r}_i}{\partial t^2} = \mathbf{F}_i, \;i=1 \ldots N.
108 :label: eqnnewtonslaws
110 The forces are the negative derivatives of a potential function
111 :math:`V(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N)`:
113 .. math:: \mathbf{F}_i = - \frac{\partial V}{\partial \mathbf{r}_i}
116 The equations are solved simultaneously in small time steps. The system
117 is followed for some time, taking care that the temperature and pressure
118 remain at the required values, and the coordinates are written to an
119 output file at regular intervals. The coordinates as a function of time
120 represent a *trajectory* of the system. After initial changes, the
121 system will usually reach an *equilibrium state*. By averaging over an
122 equilibrium trajectory, many macroscopic properties can be extracted
123 from the output file.
125 It is useful at this point to consider the limitations of MD
126 simulations. The user should be aware of those limitations and always
127 perform checks on known experimental properties to assess the accuracy
128 of the simulation. We list the approximations below.
130 **The simulations are classical**
132 - Using Newton’s equation of motion automatically implies the use of
133 *classical mechanics* to describe the motion of atoms. This is all
134 right for most atoms at normal temperatures, but there are
135 exceptions. Hydrogen atoms are quite light and the motion of
136 protons is sometimes of essential quantum mechanical character.
137 For example, a proton may *tunnel* through a potential barrier in
138 the course of a transfer over a hydrogen bond. Such processes
139 cannot be properly treated by classical dynamics! Helium liquid at
140 low temperature is another example where classical mechanics
141 breaks down. While helium may not deeply concern us, the high
142 frequency vibrations of covalent bonds should make us worry! The
143 statistical mechanics of a classical harmonic oscillator differs
144 appreciably from that of a real quantum oscillator when the
145 resonance frequency :math:`\nu` approximates or exceeds
146 :math:`k_BT/h`. Now at room temperature the wavenumber
147 :math:`\sigma = 1/\lambda = \nu/c` at which :math:`h
148 \nu = k_BT` is approximately 200 cm\ :math:`^{-1}`. Thus, all
149 frequencies higher than, say, 100 cm\ :math:`^{-1}` may misbehave
150 in classical simulations. This means that practically all bond and
151 bond-angle vibrations are suspect, and even hydrogen-bonded
152 motions as translational or librational H-bond vibrations are
153 beyond the classical limit (see :numref:`Table %s <tab-vibrations>`)
156 .. |H2CX| replace:: H\ :math:`_2`\ CX
157 .. |OHO1| replace:: O-H\ :math:`\cdots`\ O
158 .. |INCM| replace:: :math:`\mathrm{cm}~^{-1}`
163 Typical vibrational frequencies (wavenumbers) in molecules and hydrogen-bonded
164 liquids. Compare :math:`kT/h = 200~\mathrm{cm}^{-1}` at 300 K.
168 +---------------+-------------+------------+
169 | | type of | wavenumber |
170 | type of bond | vibration | |INCM| |
171 +===============+=============+============+
172 | C-H, O-H, N-H | stretch | 3000--3500 |
173 +---------------+-------------+------------+
174 | C=C, C=O | stretch | 1700--2000 |
175 +---------------+-------------+------------+
176 | HOH | bending | 1600 |
177 +---------------+-------------+------------+
178 | C-C | stretch | 1400--1600 |
179 +---------------+-------------+------------+
180 | |H2CX| | sciss, rock | 1000--1500 |
181 +---------------+-------------+------------+
182 | CCC | bending | 800--1000 |
183 +---------------+-------------+------------+
184 | |OHO1| | libration | 400--700 |
185 +---------------+-------------+------------+
186 | |OHO1| | stretch | 50--200 |
187 +---------------+-------------+------------+
191 - Well, apart from real quantum-dynamical simulations, we can do one
194 (a) If we perform MD simulations using harmonic oscillators for
195 bonds, we should make corrections to the total internal energy
196 :math:`U = E_{kin} + E_{pot}` and specific heat :math:`C_V` (and
197 to entropy :math:`S` and free energy :math:`A` or :math:`G` if
198 those are calculated). The corrections to the energy and specific
199 heat of a one-dimensional oscillator with frequency :math:`\nu`
200 are: \ :ref:`11 <refMcQuarrie76>`
202 .. math:: U^{QM} = U^{cl} +kT \left( {\frac{1}{2}}x - 1 + \frac{x}{e^x-1} \right)
205 .. math:: C_V^{QM} = C_V^{cl} + k \left( \frac{x^2e^x}{(e^x-1)^2} - 1 \right)
208 where :math:`x=h\nu /kT`. The classical oscillator absorbs too
209 much energy (:math:`kT`), while the high-frequency quantum
210 oscillator is in its ground state at the zero-point energy level
211 of :math:`\frac{1}{2} h\nu`.
213 (b) We can treat the bonds (and bond angles) as
214 *constraints* in the equations of
215 motion. The rationale behind this is that a quantum oscillator in
216 its ground state resembles a constrained bond more closely than a
217 classical oscillator. A good practical reason for this choice is
218 that the algorithm can use larger time steps when the highest
219 frequencies are removed. In practice the time step can be made
220 four times as large when bonds are constrained than when they are
221 oscillators \ :ref:`12 <refGunsteren77>`. |Gromacs| has this
222 option for the bonds and bond angles. The flexibility of the
223 latter is rather essential to allow for the realistic motion and
224 coverage of configurational space \ :ref:`13 <refGunsteren82>`.
226 **Electrons are in the ground state**
227 In MD we use a *conservative* force field that is a function of
228 the positions of atoms only. This means that the electronic
229 motions are not considered: the electrons are supposed to adjust
230 their dynamics instantly when the atomic positions change (the
232 approximation), and remain in their ground state. This is really
233 all right, almost always. But of course, electron transfer
234 processes and electronically excited states can not be treated.
235 Neither can chemical reactions be treated properly, but there are
236 other reasons to shy away from reactions for the time being.
238 **Force fields are approximate**
241 They are not really a part of the simulation method and their
242 parameters can be modified by the user as the need arises or
243 knowledge improves. But the form of the forces that can be used in
244 a particular program is subject to limitations. The force field
245 that is incorporated in |Gromacs| is described in Chapter 4. In the
246 present version the force field is pair-additive (apart from
247 long-range Coulomb forces), it cannot incorporate
248 polarizabilities, and it does not contain fine-tuning of bonded
249 interactions. This urges the inclusion of some limitations in this
250 list below. For the rest it is quite useful and fairly reliable
251 for biologically-relevant macromolecules in aqueous solution!
253 **The force field is pair-additive**
254 This means that all *non-bonded* forces result from the sum of
255 non-bonded pair interactions. Non pair-additive interactions, the
256 most important example of which is interaction through atomic
257 polarizability, are represented by *effective pair potentials*.
258 Only average non pair-additive contributions are incorporated.
259 This also means that the pair interactions are not pure, *i.e.*,
260 they are not valid for isolated pairs or for situations that
261 differ appreciably from the test systems on which the models were
262 parameterized. In fact, the effective pair potentials are not that
263 bad in practice. But the omission of polarizability also means
264 that electrons in atoms do not provide a dielectric constant as
265 they should. For example, real liquid alkanes have a dielectric
266 constant of slightly more than 2, which reduce the long-range
267 electrostatic interaction between (partial) charges. Thus, the
268 simulations will exaggerate the long-range Coulomb terms. Luckily,
269 the next item compensates this effect a bit.
271 **Long-range interactions are cut off**
272 In this version, |Gromacs| always uses a
274 radius for the Lennard-Jones
275 interactions and sometimes for the Coulomb interactions as well.
276 The “minimum-image convention” used by |Gromacs| requires that only
277 one image of each particle in the periodic boundary conditions is
278 considered for a pair interaction, so the cut-off radius cannot
279 exceed half the box size. That is still pretty big for large
280 systems, and trouble is only expected for systems containing
281 charged particles. But then truly bad things can happen, like
282 accumulation of charges at the cut-off boundary or very wrong
283 energies! For such systems, you should consider using one of the
284 implemented long-range electrostatic algorithms, such as
285 particle-mesh Ewald \ :ref:`14 <refDarden93>`,
286 :ref:`15 <refEssmann95>`.
288 **Boundary conditions are unnatural**
289 Since system size is small (even 10,000 particles is small), a
290 cluster of particles will have a lot of unwanted boundary with its
291 environment (vacuum). We must avoid this condition if we wish to
292 simulate a bulk system. As such, we use periodic boundary
293 conditions to avoid real phase boundaries. Since liquids are not
294 crystals, something unnatural remains. This item is mentioned last
295 because it is the least of the evils. For large systems, the
296 errors are small, but for small systems with a lot of internal
297 spatial correlation, the periodic boundaries may enhance internal
298 correlation. In that case, beware of, and test, the influence of
299 system size. This is especially important when using lattice sums
300 for long-range electrostatics, since these are known to sometimes
301 introduce extra ordering.
303 Energy Minimization and Search Methods
304 --------------------------------------
306 As mentioned in sec. :ref:`Compchem`, in many cases energy minimization
307 is required. |Gromacs| provides a number of methods for local energy
308 minimization, as detailed in sec. :ref:`EM`.
310 The potential energy function of a (macro)molecular system is a very
311 complex landscape (or *hypersurface*) in a large number of dimensions.
312 It has one deepest point, the *global minimum* and a very large number
313 of *local minima*, where all derivatives of the potential energy
314 function with respect to the coordinates are zero and all second
315 derivatives are non-negative. The matrix of second derivatives, which is
316 called the *Hessian matrix*, has non-negative eigenvalues; only the
317 collective coordinates that correspond to translation and rotation (for
318 an isolated molecule) have zero eigenvalues. In between the local minima
319 there are *saddle points*, where the Hessian matrix has only one
320 negative eigenvalue. These points are the mountain passes through which
321 the system can migrate from one local minimum to another.
323 Knowledge of all local minima, including the global one, and of all
324 saddle points would enable us to describe the relevant structures and
325 conformations and their free energies, as well as the dynamics of
326 structural transitions. Unfortunately, the dimensionality of the
327 configurational space and the number of local minima is so high that it
328 is impossible to sample the space at a sufficient number of points to
329 obtain a complete survey. In particular, no minimization method exists
330 that guarantees the determination of the global minimum in any practical
331 amount of time. Impractical methods exist, some much faster than
332 others \ :ref:`16 <refGeman84>`. However, given a starting configuration,
333 it is possible to find the *nearest local minimum*. “Nearest” in this
334 context does not always imply “nearest” in a geometrical sense (*i.e.*,
335 the least sum of square coordinate differences), but means the minimum
336 that can be reached by systematically moving down the steepest local
337 gradient. Finding this nearest local minimum is all that |Gromacs| can do
338 for you, sorry! If you want to find other minima and hope to discover
339 the global minimum in the process, the best advice is to experiment with
340 temperature-coupled MD: run your system at a high temperature for a
341 while and then quench it slowly down to the required temperature; do
342 this repeatedly! If something as a melting or glass transition
343 temperature exists, it is wise to stay for some time slightly below that
344 temperature and cool down slowly according to some clever scheme, a
345 process called *simulated annealing*. Since no physical truth is
346 required, you can use your imagination to speed up this process. One
347 trick that often works is to make hydrogen atoms heavier (mass 10 or
348 so): although that will slow down the otherwise very rapid motions of
349 hydrogen atoms, it will hardly influence the slower motions in the
350 system, while enabling you to increase the time step by a factor of 3 or
351 4. You can also modify the potential energy function during the search
352 procedure, *e.g.* by removing barriers (remove dihedral angle functions
353 or replace repulsive potentials by *soft-core*
354 potentials \ :ref:`17 <refNilges88>`), but always take care to restore the correct
355 functions slowly. The best search method that allows rather drastic
356 structural changes is to allow excursions into four-dimensional
357 space \ :ref:`18 <refSchaik93>`, but this requires some extra programming
358 beyond the standard capabilities of |Gromacs|.
360 Three possible energy minimization methods are:
362 - Those that require only function evaluations. Examples are the
363 simplex method and its variants. A step is made on the basis of the
364 results of previous evaluations. If derivative information is
365 available, such methods are inferior to those that use this
368 - Those that use derivative information. Since the partial derivatives
369 of the potential energy with respect to all coordinates are known in
370 MD programs (these are equal to minus the forces) this class of
371 methods is very suitable as modification of MD programs.
373 - Those that use second derivative information as well. These methods
374 are superior in their convergence properties near the minimum: a
375 quadratic potential function is minimized in one step! The problem is
376 that for :math:`N` particles a :math:`3N\times 3N` matrix must be
377 computed, stored, and inverted. Apart from the extra programming to
378 obtain second derivatives, for most systems of interest this is
379 beyond the available capacity. There are intermediate methods that
380 build up the Hessian matrix on the fly, but they also suffer from
381 excessive storage requirements. So |Gromacs| will shy away from this
384 The *steepest descent* method, available in |Gromacs|, is of the second
385 class. It simply takes a step in the direction of the negative gradient
386 (hence in the direction of the force), without any consideration of the
387 history built up in previous steps. The step size is adjusted such that
388 the search is fast, but the motion is always downhill. This is a simple
389 and sturdy, but somewhat stupid, method: its convergence can be quite
390 slow, especially in the vicinity of the local minimum! The
391 faster-converging *conjugate gradient method* (see *e.g.*
392 :ref:`19 <refZimmerman91>`) uses gradient information from previous steps. In general,
393 steepest descents will bring you close to the nearest local minimum very
394 quickly, while conjugate gradients brings you *very* close to the local
395 minimum, but performs worse far away from the minimum. |Gromacs| also
396 supports the L-BFGS minimizer, which is mostly comparable to *conjugate
397 gradient method*, but in some cases converges faster.