First look into (particle) simulations#

Learning goals:

To learn the basic terms and terminology of particle simulations
To have a basic understanding of molecular topologies
To understand the difference between particle and continuum models
To understand the concept of force field

Keywords: Particle simulations, molecular dynamics, Monte Carlo, continuum simulations, molecular topology, force field, interactions

It is finally the time to discuss some of the basic concepts of (particle) simulations. As this is the first description, it does not aim to be complete. Rather, the aim is to introduce some of basic concepts and terminology. This information will be needed and elaborated later, and this knowledge will help in setting up simulations and in data analysis.

Particle vs continuum simulations#

Although particle simulations seems almost too obvious a term, it needs some attention. The most common basic object in a particle simulation is an atom or, rather, a generalized atom. In this course we are mostly dealing with atoms in the sense that they don’t have electrons. The overall properties of atoms at the classical level are included in the force fields as we will discuss.

Force field: The parameters and functions that describe interactions in simulations.

These electron-less atoms do, however, have many other properties allowing for comparisons between the simulation results and experiments as long as one is in the same atomistic scales (for example energies, diffusion coefficients and such). The important concept here is that by particle simulations we refer to a system which is composed of discrete entities.

Atom: From the Greek word atomos, which means "uncuttable".

These discrete entities may also be coarse-grained atoms (sometimes also called superatoms). The concept of a coarse-grained atom can be thought of as follows: Assume that there is some method that allows us to average over some chosen physical properties of a group of atoms, for example a group of atoms in a polymer, and combine them to a single larger entity. That new entity is then a coarse-grained atom. At first sight such a procedure may sound like heresy or nonsense, but there are some well-justified and rigorous methods for such procedures. Such methods will not be discussed here. For coarse-grained systems and methods for coarse-graining, see for example references[5, 11, 13, 14, 16]. An example of such an approach is the MARTINI force field[9] in which four ‘atomistic’ atoms are a coarse-grained into one larger unit; Martini is the most popular coarse-grained approach at the moment.

If in some cataclysm all scientiﬁc knowledge were to be destroyed and only one sentence passed onto the next generation of creatures,what statement would contain the most information in the fewest words? I believe it is the atomic hypothesis that all things are made of atoms–little particles that move around in perpetual motion,attracting each other when they area little distance apart,but repelling upon being squeezed into one another. In that one sentence,you will see there is an enormous amount of information about the world, of just a little imagination and thinking is applied.

—Richard Feynman in Feynman Lectures in Physics

Simulations of even larger non-atomistic entities such as planets are in the category of particle simulations since they are discrete entities. In such simulations, the internal structure of planets or celestial bodies is not included in the simulated properties. This is similar in spirit to not taking electrons into account in atomic-level molecular dynamics simulations. Other examples that may not be so obvious at first sight are entities such as lattice defects in solid state systems. Defects can be treated as discrete entities as they have properties such as size, and they interact over distances. Here, we use ‘atomistic’ atoms and the term atomistic is used for describing systems at the atomic scale of Ångströms and nanometers.

In contrast to particle simulations, in continuum modeling one often uses a fixed underlying lattice or grid in which one solves finite difference equations. In such simulations the material (say, a liquid, gas or a solid) is described by a continuous field instead of an atom. This means that, unlike in an atomistic simulation, the material no longer has internal structure but instead is described by a field or several fields. The field(s) is (are) placed on a discrete grid and the time evolution(s) of the field(s) is (are) solved on a computer; the corresponding equations of motion are solved on a lattice or using finite element methods (FEM). Such approaches are commonly used in the fields of materials research and reaction-diffusion systems. As a conceptual point, it is worth noticing that electronic structure calculations should not be categorized as particle simulations. In such calculations the nuclei of the atoms are fixed, and the electronic structure is calculated by a minimization process: Minimization is not an actual simulation in the sense of the above.

Simulation methods for particle systems#

The main simulation methods for particle systems are the Monte Carlo method and Molecular Dynamics (MD).

Monte Carlo#

In this course the focus is on the molecular dynamics method, but for completeness, let’s briefly discuss the Monte Carlo (MC) method since it is one the most widely used simulation techniques (if not the most widely used) with applications ranging from social sciences, to optimizing public transport networks to particle simulations.

The MC method does not solve any equations of motion, instead it is a stochastic optimization method using an energy function: No time integration is involved, but one uses random numbers for optimizing selected properties. In the case molecular systems, one typically optimizes (minimizes) the potential energy. This also means that, in its basic form, MC cannot measure dynamical properties (there are variants like the kinetic Monte Carlo or KMC, but even that doesn’t have time integration). Instead, MC uses ensemble averaging to compute the system’s properties.

The basic idea of MC simulation is that it uses pseudo-random numbers to overcome local energy barriers to reach the minimum energy state. Random numbers are used to give the system a small probability to move in the direction that is energetically unfavorable. That is, in a crude sense the system ‘gambles’ to find the global minimum state; the name Monte Carlo comes from the Monte Carlo Casino in the city of Monaco on the French Riviera or Côte d’Azur. In terms of an atomistic system, only the positions of the atoms are used and there are no velocities. Since the particles have mutual interactions described by a potential energy, minimization of potential energy yields the equilibrium configuration of the system. The part that sometimes causes confusion is that in MC, term Monte Carlo step is used. The term indicates the number of so-called Monte Carlo sweeps, and it is sometimes falsely equated to the time steps in MD simulations. In MD, there is a real time step (typically in terms of femtoseconds) but the MC step does not correspond to any physical time.

energy-minimum — Fig. 23 The concept behind the Monte Carlo method. Consider that a system’s state is described by the green sphere and the that the energy landscape is given by the dashed line. Deterministic optimization finds the closest minimum as indicated by the arrow in the left hand side figure. Once the minimum has been reached, the system stays there. Monte Carlo is a stochastic optimization method. It uses random numbers to enable the system to move uphill, that is, in the energetically unfavorable direction, as indicated by the arrow in the right hand side figure. Once the system crosses the saddle point, it can find the global minimum.#

Technically, a Monte Carlo simulation is the evaluation of a statistical mechanical configuration integral. We will discuss it later, but MC generates representative samples of the system. The key ingredients are the application of a probability distribution and the use of random numbers. In the case of physical systems, the equilibrium state has a unique distribution: the Boltzmann distribution.

Although it may not be clear from the above, one of the great advantages of the MC method is its flexibility. In terms of atoms or particles, it can be used to simulate classical and quantum mechanical systems (quantum Monte Carlo). It can also be applied to systems of very different nature such as optimizing train scheduling for railway networks and to numerically solve complex mathematical integrals that have no analytical solution.

Classical Molecular Dynamics#

Now that we have an idea of the Monte Carlo, let’s move on to our main interest, the MD method. We can define classical Molecular Dynamics the following way:

Definition of Molecular Dynamics (MD)

Molecular Dynamics is the evolution of Newton’s equations motion in time. In short, the problem is simply to solve numerically F=ma.

This definition has several important points:

We have a dynamical system that evolves in time. That is, it provides time-resolved information. This is in contrast to the Monte Carlo method.
The evolution of the system in the phase space is determined by solving Newton’s equations of motion for all the particles that comprise the system.
In the absence of external forces and fields, the system must evolve toward equilibrium, that is, the state of minimum free energy. Statistical mechanics tells us that equilibrium is described by a unique distribution, the Boltzmann distribution.

The first article using classical MD simulations titled “Phase Transition for a Hard Sphere System” was published by Alder and Wainwright in 1957[1]. They simulated systems of 32 hard particles and also performed shorter simulations with 108 particles. As already discussed, the first MC simulations were performed a few years before in 1952 by Metropolis et al.[10]. Other landmarks include the paper of Gibson et al., “Dynamics of Radiation Damage” which is the first time when continuous potentials were used in an MD simulation in the same spirit as we do today[6]. Another interesting point is the number of time steps taken in their simulation. They had 500 atoms and reported a total 219 steps which took over 3.5 hours. Keep this number in mind for the time when we set up our own simulations. Another landmark in MD simulations is the emergence of force fields and the first QM/MM simulations which can be traced to the articles of Warshel, Karplus and Lifson in the late 1960’s and the early 1970s[8, 18].

What is the output of a particle simulation?#

The basic result of a particle simulation is a trajectory, that is, the positions of atoms written at a requested output interval (the output interval is typically given as an input parameter for the simulation). In the case of an MD simulation, it is also possible to output the velocities and forces together with the positions. An MC simulation, however, does not have velocities or forces. Instead, it has energies. They can be written out at the requested interval. One should also notice that an MD simulation has time and hence the output interval is given in terms of time (or steps that can be easily converted to time), in a Monte Carlo simulation the MC step has no relation to real time.

Most simulation programs also output other quantities such as energies, pressures, the size of the simulation box and so on. However, the trajectory is the most fundamental outcome of a particle simulation and it is the basis of more detailed analyses.

Example output from an MD simulation

The actual output format varies depending on software, but in the most typical case the positions of all the atoms are stored. It is also possible to store velocities and forces in the case of MD simulations (but that is rarely done). The example below shows an output in the format of a Gromacs .gro file.

The output below has been extracted (5 water molecules out of 216 were selected) from the file spc216.gro that comes with Gromacs; the full file will be used when we set up simulations that include water. Notice also that most software, including Gromacs, can output trajectory data in multiple formats. The example below is the human-readable .gro format but during a simulation binary formats are used since they take much less space to store. File formats will be discussed in detail in later lectures.

The 1st line has comments. Arbitrary information can be included and this line is ignored by the analysis software. In this case is says: 5 water molecules, SPC model, temperature of 300K and box size of 1.86206.
The 2nd line: The total number of atoms. This information is required by the .gro format.
The lines starting with a number and SOL: 1SOL is the 1st solvent molecule. 2SOL: second solvent molecule and so on. Each solvent/water molecule consist of three atoms: OW, HW1, HW2.
- OW: Water oxygen
- HW1 and HW2: Water hydrogens 1 and 2.
The 3rd column: Numbering of individual atoms. Each atom must have a unique identifier.
Columns 4-6: x, y and z coordinates in nanometers.
The last line: the size of the simulation box in nanometers.

For example, VMD understand this format, and can read and visualize the data. Notice also that formatting is strict regarding the positioning of the columns. For the full specification of the .gro format see The gro format.

5H2O,SPC-MODEL,300K,BOX(M)=1.86206
  15
    1SOL     OW    1    .230    .628    .113
    1SOL    HW1    2    .137    .626    .150
    1SOL    HW2    3    .231    .589    .021
    2SOL     OW    4    .225    .275   -.866
    2SOL    HW1    5    .260    .258   -.774
    2SOL    HW2    6    .137    .230   -.878
    3SOL     OW    7    .019    .368    .647
    3SOL    HW1    8   -.063    .411    .686
    3SOL    HW2    9   -.009    .295    .584
    4SOL     OW   10    .569   -.587   -.697
    4SOL    HW1   11    .476   -.594   -.734
    4SOL    HW2   12    .580   -.498   -.653
    5SOL     OW   13   -.307   -.351    .703
    5SOL    HW1   14   -.364   -.367    .784
    5SOL    HW2   15   -.366   -.341    .623
     1.86206   1.86206   1.86206

The output, when visualized using VMD, looks like this:

5spc — Fig. 24 The visualization (done with VMD), shows 5 SPC (simple point-charge) waters [2] and the coordinate system.#

System sizes#

Let’s start with atoms. The size of an atom is in the scale of Ångströms (10\(^{-10}\) m). Thus, it is quite easy to see that describing anything in macroscopic length scales is simply impossible; just make a simple estimate of how many water molecules would it take to fill a endoscopic container of size 10 nm\(\times\)10 nm\(\times\)10 nm! For reference, if you look at the example .gro file from the dropdown above, the last line gives the dimensions of the simulation box in nanometers: the box length is 1.86206 nanometers. When using the more or less correct density of water, the box can fit a maximum of about 216 water molecules. Since each water has (using the SPC model) three atoms, the total number of atoms is 648.

One advantage of coarse-grained models such as Martini or continuum simulations using fields, is that it is possible to describe systems in much larger length scales since one is not limited to atoms. What one loses is the atomic specificity and its consequences. There are approaches that allow continuum methods to be more material specific but such techniques are outside the scope of this course. For those interested, such methods include phase-field[17] and phase-field crystal methods[15].

Summary#

This section introduced a number of concepts related particle modeling. In the next sections, we will make these concepts more quantitative and put them into practise when we set up and analyze simulations and simulation data.

Bibliography#

1: B J Alder and T E Wainwright. Phase transition for a hard sphere system. J. Chem. Phys., 27(5):1208, 1957. doi:10.1063/1.1743957.
2: H J C Berendsen, J P M Postma, W F van Gunsteren, and J Hermans. Interaction models for water in relation to protein hydration. In Bernard Pullman, editor, Intermolecular Forces: Proceedings of the Fourteenth Jerusalem Symposium on Quantum Chemistry and Biochemistry Held in Jerusalem, Israel, April 13–16, 1981, pages 331–342. Springer Netherlands, Dordrecht, 1981. doi:10.1007/978-94-015-7658-1\_21.
3: David M Ceperley and Stephen B Libby. Berni julian alder, theoretical physicist and inventor of molecular dynamics, 1925-2020. Proc. Natl. Acad. Sci. U. S. A., March 2021. doi:10.1073/pnas.2024252118.
4: Leticia F Cugliandolo. Coarsening phenomena. C. R. Phys., 16(3):257–266, April 2015. doi:10.1016/j.crhy.2015.02.005.
5: Dominik Fritz, Konstantin Koschke, Vagelis A Harmandaris, Nico F A van der Vegt, and Kurt Kremer. Multiscale modeling of soft matter: scaling of dynamics. Phys. Chem. Chem. Phys., 13(22):10412–10420, June 2011. doi:10.1039/c1cp20247b.
6: J B Gibson, A N Goland, M Milgram, and G H Vineyard. Dynamics of radiation damage. Phys. Rev., 120(4):1229–1253, November 1960. doi:10.1103/PhysRev.120.1229.
7: Teemu Leppänen, M Karttunen, Kimmo Kaski, Rafael A Barrio, and Limei Zhang. A new dimension to turing patterns. Physica D, 168-169:35–44, August 2002. doi:10.1016/S0167-2789(02)00493-1.
8: S. Lifson and A. Warshel. Consistent force field for calculations of conformations, vibrational spectra, and enthalpies of cycloalkane andn-alkane molecules. The Journal of Chemical Physics, 49(11):5116–5129, dec 1968. doi:10.1063/1.1670007.
9(1,2): Siewert J Marrink and D Peter Tieleman. Perspective on the Martini model. Chem. Soc. Rev., 42(16):6801–6822, August 2013. doi:10.1039/c3cs60093a.
10: Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller. Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6):1087, 1953. doi:10.1063/1.1699114.
11: Teemu Murtola, Alex Bunker, Ilpo Vattulainen, Markus Deserno, and Mikko Karttunen. Multiscale modeling of emergent materials: biological and soft matter. Phys. Chem. Chem. Phys., 11(12):1869–1892, March 2009. doi:10.1039/b818051b.
12: J E Pearson. Complex patterns in a simple system. Science, 261(5118):189–192, July 1993. doi:10.1126/science.261.5118.189.
13: Raffaello Potestio, Christine Peter, and Kurt Kremer. Computer simulations of soft matter: linking the scales. Entropy, 16(8):4199–4245, July 2014. doi:10.3390/e16084199.
14: Matej Praprotnik, Luigi Delle Site, and Kurt Kremer. Multiscale simulation of soft matter: from scale bridging to adaptive resolution. Annu. Rev. Phys. Chem., 59(1):545–571, 2008. doi:10.1146/annurev.physchem.59.032607.093707.
15: N Provatas, J A Dantzig, B Athreya, P Chan, P Stefanovic, N Goldenfeld, and K R Elder. Using the phase-field crystal method in the multi-scale modeling of microstructure evolution. JOM, 59(7):83–90, July 2007. doi:10.1007/s11837-007-0095-3.
16: Sereina Riniker, Jane R Allison, and Wilfred F van Gunsteren. On developing coarse-grained models for biomolecular simulation: a review. Phys. Chem. Chem. Phys., 14(36):12423–12430, September 2012. doi:10.1039/c2cp40934h.
17: Ingo Steinbach. Phase-field models in materials science. Modell. Simul. Mater. Sci. Eng., 17(7):073001, October 2009. doi:10.1088/0965-0393/17/7/073001.
18: A Warshel and M Karplus. Calculation of ground and excited state potential surfaces of conjugated molecules. i. formulation and parametrization. J. Am. Chem. Soc., 94(16):5612–5625, August 1972. doi:10.1021/ja00771a014.

Molecular Dynamics Simulations

First look into (particle) simulations

Contents

First look into (particle) simulations#

Particle vs continuum simulations#

Simulation methods for particle systems#

Monte Carlo#

Classical Molecular Dynamics#

What is the output of a particle simulation?#

System sizes#

Summary#

Bibliography#