Introduction to Molecular Dynamics Simulations: Applications in Hard and Soft Condensed Matter Physics

Molecular Dynamics is a two-volume compendium of the ever-growing applications of molecular dynamics simulations to solve a wider range of scientific and engineering challenges. The contents illustrate the rapid progress on molecular dynamics simulations in many fields of science and technology, such as nanotechnology, energy research, and biology, due to the advances of new dynamics theories and the extraordinary power of today's computers. This second book begins with an introduction of molecular dynamics simulations to macromolecules and then illustrates the computer experiments using molecular dynamics simulations in the studies of synthetic and biological macromolecules, plasmas, and nanomachines. Coverage of this book includes: Complex formation and dynamics of polymers Dynamics of lipid bilayers, peptides, DNA, RNA, and proteins Complex liquids and plasmas Dynamics of molecules on surfaces Nanofluidics and nanomachines

the computer experiment. Some traditionalists working in theory, who have not followed the modern developments of computer science and its applications in the realm of physics, biology, chemistry and many more scientific fields, still repudiate the term "experiment" in the context of computer simulations. However, this term is most certainly fully justified! In a computer experiment, a model is still provided by theory, but the calculations are carried out by the machine by following a series of instructions (the algorithm) usually coded in some high-level language and translated (compiled) into assembler commands which provide instructions how to manipulate the contents of processor registers. The results of computer simulations are just numbers, data which have to be interpreted by humans, either in the form of graphical output, as tables or as function plots. By using a machine to carry out the calculations necessary for solving a model, more complexity can be introduced and more realistic systems can be investigated.
Simulation is seen sometimes as theory, sometimes as experiment. On the one side, one is still dealing with models, not with "real systems" 2 On the other side, the procedure of verifying a model by computer simulation resembles an experiment very closely: One performs a run, then analyzes the results in pretty much the same way as an experimental physicist does. Simulations can come very close to experimental conditions which allows for interpreting and understanding the experiments at the microscopic level, but also for studying regions of systems which are not accessible in "real" experiments 3 ,tooexpensiveto perform, or too dangerous. In addition, computer simulations allow for performing thought experiments, which are impossible to do in reality, but whose outcome greatly increases our understanding of fundamental processes or phenomena. Imagination and creativity, just like in mathematics 4 , physics and other scientific areas, are very important qualities of a computer simulator! From a principal point of view, theory is traditionally based on the reductionist approach: one deals with complexity by reducing a system to simpler subsystems, continuing until the subsystems are simple enough to be represented with solvable models. From this perspective one can regard simulation as a convenient tool to verify and test theories and the models associated with them in situations which are too complex to be handled analytically 5 . Here, one implies that the model represents the level of the theory where the interest is focused.
However, it is important to notice that simulation can play a more important role than just being a tool to be used as an aid to reductionism because it can be considered as an alternative to it. Simulation increases considerably the threshold of complexity which separates solvable und unsolvable models. One can take advantage of this threshold shift and move up several levels in our description of physical systems. Thanks to the presence of simulation, we do not need to work with models as simple as those used in the past. This gives the researcher an additional freedom for exploration. As an example, one could mention the interatomic potentials which, in the past, were obtained by two-body potentials with simple Introduction to Molecular Dynamics Simulations: Applications in Hard and Soft Condensed Matter Physics analytical form, such as Morse or Lennard-Jones. Today, the most accurate potentials contain many-body terms and are determined numerically by reproducing as closely as possible the forces predicted by ab-initio methods. These new potentials could not exist without simulation, so simulation is not only a connecting link between theory and experiment, but it is also a powerful tool to make progress in new directions. Readers, interested in these more "philosophical" aspects of computational science will be able to find appropriate discussions in the first chapters of refs. (Haile, 1992;Steinhauser, 2008;2012). In the following, we focus on classical molecular dynamics (MD) simulations, i.e. a variant of MD which neglects any wave character of discrete atomic particles, describing them as classical in the Newtonian sense and not referring to any quantum mechanical considerations.

The objective of molecular dynamics simulations
Molecular dynamics computer experiments are done in an attempt of understanding the properties of assemblies of molecules in terms of their structure and the microscopic interactions between them. We provide a guess at the interactions between molecules, and obtain exact predictions of bulk properties. The predictions are "exact" in the sense that they can be made as accurate as we like, subject to the limitations imposed by our computer budget. At the same time, the hidden dynamic details behind bulk measurements can be revealed. An example is the link between the diffusion coefficient and the velocity autocorrelation function, with the latter being very hard to measure experimentally, but the former being very easy to measure. Ultimately one wants to make direct comparison with experimental measurements made on specific materials, in which case a good model of molecular interactions is essential. The aim of so-called ab-inito MD is to reduce the amount of guesswork and fitting in this process to a minimum. On the other hand, sometimes one is merely interested in phenomena of a rather generic nature, or one wants to discriminate between bad and good theories. In this case it is not necessary to have a perfectly realistic molecular model, but one that contains the essential physics may be quite suitable.

Molecular interactions
Classical MD simulation consists of the numerical, step-by-step solution of the classical Newtonian equations of motion, which for a simple atomic system may be written as where the vector symbol in the partial derivative is a physicsts' abbreviatory notation for the derivate of each individual component. To solve Eq. 1 one needs to calculate the forces F i acting on the atoms, which are usually derived from a potential energy function Φ( r N ),where r N =( r 1 , r 2 , ..., r N ) represents the complete set of 3N atomic coordinates.

Bonded interactions
Using the notion of intermolecular potentials acting between the particles of a system one cannot only model fluids made of simple spherically symmetric particles but also more complex molecules with internal degrees of freedom (due to their specific monomer connectivity). If one intends to incorporate all aspects of the chemical bond in complex molecules one has to treat the system with quantum chemical methods. Usually, one considers the internal degrees of freedom of polymers and biomacromolecules by using generic potentials that describe bond lengths l i ,b o n da n g l e sθ and torsion angles φ.W h e n neglecting the fast electronic degrees of freedom, often bond angles and bond lengths can be assumed to be constants. In this case, the potential includes lengths l 0 and angles θ 0 , φ 0 at equilibrium about which the molecules are allowed to oscillate, and restoring forces which ensure that the system attains these equilibrium values on average. Hence the bonded interactions Φ bonded for polymeric macromolecular systems with internal degrees of freedom can be treated by using some or all parts of the following potential term: Here, the summation indices sum up the number of bonds i at positions r i , the number of bond angles k between consecutive monomers along a macromolecular chain and the number of torsion angles m along the polymer chain. A typical value of κ = 5000 ensures that the fluctuations of bond angles are very small (below 1%). The terms l 0 , θ 0 and φ 0 are the equilibrium distance, bond angle and torsion angle, respectively.
In particular in polymer physics, very often a Finitely Extensible Non-linear Elastic (FENE) potential is used which -in contrast to a harmonic potential -restricts the maximum bond length of a polymer bond to a prefixed value R 0 : The FENE potential can be used instead of the first term on the right hand side of the bonded potential in Eq. 2. Figure 1 illustrates the different parameters which are used in the description of bonded interactions in Eq. 2. Further details on the use of potentials in macromolecular biology and polymer physics may be found in (Feller, 2000;Schlenkrich et al., 1996;Siu et al., 2008;Steinhauser, 2008).

Non-bonded interactions
Various physical properties are determined by different regions of the potential hypersurface of interacting particles. Thus, for a complete determination of potential curves, widespread experiments are necessary. For a N−body system the total energy Φ nb , i.e. the potential hypersurface of the non-bonded interactions can be written as (Allen & Tildesly, 1991) where φ 1 , φ 2 , φ 3 , ... are the interaction contributions due to external fields (e.g. the effect of container walls) and due to pair, triple and higher order interactions of particles. In classical MD one often simplifies the potential by the hypothesis that all interactions can be described by pairwise additive potentials. Despite this reduction of complexity, the efficiency of a MD algorithm taking into account only pair interactions of particles is rather low (of order O(N 2 )) and several optimization techniques are needed in order to improve the runtime behavior to O(N).
The simplest general Ansatz for a non-bonded potential for spherically symmetric systems, i.e. Φ( r)=Φ(r) with r = | r i − r j | is a potential of the following form: Parameters C 1 and C 2 are parameters of the attractive and repulsive interaction and the electrostatic energy Φ Coulomb (r) between the particles with position vectors r i and r j is given by: The constant k = 1 in the cgs-system of units and ǫ is the dielectric constant of the medium, for example ǫ air = 1forair ,ǫ prot = 4 for proteins or ǫ H 2 0 = 82 for water. The variables z i and z j denote the charge of individual monomers in the macromolecule and e is the electric charge of an electron.
The probably most commonly used form of the potential of two neutral atoms which are only bound by Van-der-Waals interactions, is the Lennard-Jones (LJ), or (a-b) potential which has the form (Haberland et al., 1995) where The most often used LJ-(6-12) potential for the interaction between two particles with a distance r = | r i − r j | then reads (cf. Eq. 5): Parameter ε determines the energy scale and σ 0 the length scale. In simulations one uses dimensionless reduced units which tend to avoid numerical errors when processing very small numbers, arising e.g. from physical constants such as the Boltzmann constant k B = 1.38 × 10 −23 J/K. In these reduced (simulation) units, one MD timestep is measured in units ofτ = (mσ 2 /ε) 1/2 ,w h e r em is the mass of a particle and ε and σ 0 are often simply set to σ 0 = ε = k B T = 1. Applied to real molecules, for example to Argon with m = 6.63 × 10 −23 kg, σ 0 ≈ 3.4 × 10 −10 m and ε/k B ≈ 120K one obtains a typical MD timestep ofτ ≈ 3.1 × 10 −13 s.
Using an exponential function instead of the repulsive r −12 term, one obtains the Buckingham potential (Buckingham, 1938): This potential however has the disadvantage of using a numerically very expensive exponential function and it is known to be unrealistic for many substances at small distances r where it has to be modified accordingly.
For reasons of efficiency, a classical MD potential should be short-ranged in order to keep the number of force calculations between interacting particles at a minimum. Therefore, instead of using the original form of the potential in Eq. 9, which approaches 0 at infinity, it is common to use a modified form, where the potential is simply cut off at its minimum value r = r min = 6 √ 2 and shifted to positive values by ε such that it is purely repulsive and smooth at r = r cut = 6 √ 2: Another extension of the potential in Eq. 9 is proposed in  where a smooth attractive part is introduced again, in order to allow for including different solvent qualities of the solvent surrounding the polymer: Φ cos (r)= 1 2 · cos(αr 2 + β)+γ ε.
This additional term adds an attractive part to the potential of Eq. 11 and at the same time -by appropriately choosing parameters α, β and γ -keeps the potential cutoff at r cut smooth. The parameters α, β and γ are determined analytically such that the potential tail of Φ cos has zero derivative at r = 2 1/6 and at r = r cut , while it is zero at r = r cut and has value γ at r = 2 1/6 , where γ is the depth of the attractive part. Further details can be found in . When setting r cut = 1.5 one sets γ = − and obtains α and β as solutions of the linear set of equations The total unbounded potential can then be written as: where λ is a new parameter of the potential which determines the depth of the attractive part. Instead of varying the solvent quality in the simulation by changing temperature T directly (and having to equilibrate the particle velocities accordingly), one can achieve a phase transition in polymer behavior by changing λ accordingly, cf. Fig. 2. Fig. 2. Graph of the total unbounded potential of Eq. 15 which allows for modeling the effects of different solvent qualities.
Using coarse-grained models in the context of lipids and proteins, where each amino acid of the protein is represented by two coarse-grained beads, it has become possible to simulate lipoprotein assemblies and protein-lipid complexes for several microseconds (Shih et al., 2006).
The assumption of a short ranged interaction is usually fulfilled very well for all (uncharged) polymeric fluids. However, as soon as charged systems are involved this assumption breaks down and the calculation of the Coulomb force requires special numerical treatment due to its infinite range.

Calculation of forces
The most crucial part of a MD simulation is the force calculation. At least 95% of a MD code is spent with the force calculation routine which uses a search algorithm that determines interacting particle pairs. Therefore this is the task of a MD program which has to be optimized first and foremost. We will review a few techniques that have become standard in MD simulations which enhance the speed of force calculations considerably and speed up the algorithm from O(N 2 ) run time to O(N) run time. Starting from the original LJ potential between two particles i and j with distance r = | r i − r j | of Eq. 7, one obtains the potential function for N interacting particles as the following double sum over all particles: Introduction to Molecular Dynamics Simulations: Applications in Hard and Soft Condensed Matter Physics www.intechopen.com The corresponding force F i exerted on particle i by particle j is given by the gradient with respect to r i as: where r ij =( r i − r j ) is the direction vector between particles i and j at positions r i and r j ,and r = | r i − r j |. Hence, in general, the force F i on particle i is the sum over all forces F ij := −∇ r i Φ between particle i and all other particles j: The least favorable method of looking for interacting pairs of particles and for calculating the double sum in Eqs. 16 and 17 is the "brute force" method that simply involves taking a double loop over all particles in the (usually) cubic simulation box, thus calculating 1 2 N(N − 1) interactions with a N 2 efficiency. This algorithm becomes extremely inefficient for systems of more than a few thousand particles, cf. Fig. 3(a).

The MD algorithm
The last decade has seen a rapid development in our understanding of numerical algorithms which have been summarized in a recent book (Steinhauser, 2008) that presents the current state of the field.
When introducing an N-dimensional position vector r N =( r 1 , r 2 , ..., r N ), the potential energy Φ( r N ) and the momenta p N =( p 1 , p 2 , ..., p N ), in terms of which the kinetic energy may be written as K( p N ), then the total energy H of a classical conservative system is given by H = Φ + K. The equations of motion determining all dynamics of the particles can be written as This is a system of coupled ordinary differential equations. Many methods exist to solve this set of equations numerically, among which the so-called velcity Verlet-algorithm is the one that is the most used. This algorithm integrates the equations of motion by performing the following four steps, where r i , v i , a i = T i /m i are the position, velocity and acceleration of the i-th particle, respectively: Derive a i (t + δt) from the interaction potential using r(t + δt), Calculate Further details about this standard algorithm can be found elsewhere (Steinhauser, 2008). It is exactly time reversible, symplectic, low order in time (hence permitting large timesteps), and it requires only one expensive force calculation per timestep.

Neighbor lists
In general, in molecular systems, the potential as well as the corresponding force decays very fast with the distance r between the particles. Thus, for reasons of efficiency, in molecular simulations one often uses the modified LJ potential of Eq. 11 which introduces a cutoff r cut for the potential. The idea here is to neglect all contributions in the sums in Eqs. 16 and 17 that are smaller than the threshold r cut which characterizes the range of the interaction. Thus, in this case the force F i on particle i is approximated by Contributions to the force on particle i that stem from particles j with r ≤ r cut are neglected.
This introduces a small error in the computation of the forces and the total energy of the system, but it reduces the overall computational effort from O(N 2 ) to O(N). For systems with short-ranged or rapidly decaying potentials, a very efficient algorithm for the search of potentially interacting particles, i.e. those particles that are within the cutoff distance r cut of a particle i, has been developed (Hockney, 1970). In MD this algorithm can be implemented most efficiently by geometrically dividing the volume of the (usually cubic) simulation box into small cubic cells whose sizes are slightly larger than the interaction range r cut of particles, cf. Fig. 3b. The particles are then sorted into these cells using the linked-cell algorithm (LCA). The LCA owes its name to the way in which the particle data are arranged in computer memory, namely as linked list for each cell. For the calculation of the interactions it is then sufficient to calculate the distances between particles in neighboring cells only, since cells which are further than one cell apart are by construction beyond the interaction range. Thus, the number of distance calculations is restricted to those particle pairs of neighboring cells only which means that the sums in Eq. 18 are now split into partial sums corresponding to the decomposition of the simulation domain into cells. For the force F i on particle i in cell number n one obtains a sum of the form where Ω(n) denotes cell n itself together with all cells that are direct neighbors of cell n.T h e linked-cell algorithm is a simple loop over all cells of the simulation box. For each cell there is a linked list which contains a root pointer that points to the first particle in the respective cell which then points to the next particle within this particular cell, until the last particle is reached which points to zero, indicating that all particles in this cell have been considered. Then the algorithm switches to the root pointer of the next cell and the procedure is repeated until all interacting cells have been considered, cf. Fig. 3.
Assuming the average particle density in the simulation box as ρ then the number of particles in each one of the subcells is ρ r 3 cut . The total number of subcells is N/ ρ r 3 cut and The linked-cell algorithm combined with neighbor lists which further reduces the search effort by using a list of potentially interacting neighbor particles which can be used for several timesteps before it has to be updated. In this 2D representation, the radius of the larger circle is r cut + r skin and the inner circle, which contains actually interacting particles, has radius r cut .
the total number of neighbor cells of each subcell is 26 in a cubic lattice in three dimensions (3D). Due to Newton's third law only half of the neighbors actually need to be considered. Hence, the order to which the linked-cell algorithm reduces the search effort is given by: For this method to function, the size of the simulation box has to be at least 3r cut ,c f . Fig. 3. For simulations of dense melts with many particles, this requirement is usually met. Consequently, by this method, the search-loop effort is reduced to O(N), but with a pre-factor that still can be very large, depending on the density of particles ρ and the interaction range r cut .

Boundary conditions
In a MD simulation only a very small number of particles can be considered. To avoid the (usually) undesired artificial effects of surface particles which are not surrounded by neighboring particles in all directions and thus are exerted to non-isotropic forces, one introduces periodic boundary conditions. Using this technique, one measures the "bulk" properties of the system, due to particles which are located far away from surfaces. As a rule, one uses a cubic simulation box were the particles are located. This cubic box is periodically repeated in all directions. If, during a simulation run, a particle leaves the central simulation box, then one of its image particles enters the central box from the opposite direction. Each of the image particles in the neighboring boxes moves in exactly the same way, cf. Fig. 4 for a two dimensional visualization.
The cubic box is used almost exclusively in simulations with periodic boundaries, mainly due to its simplicity, however also spherical boundary conditions have been investigated were the three-dimensional surface of the sphere induces a non-Euclidean metric. The use of periodic boundary conditions allows for the simulation of bulk properties of systems with a relatively small number of particles. Fig. 4. Two-dimensional schematic of periodic boundary conditions. The particle trajectories in the central simulation box are copied in every direction.

Minimum image convention
The question whether the measured properties with a small, periodically extended system are to be regarded as representative for the modeled system depends on the specific observable that is investigated and on the range of the intermolecular potential. For a LJ potential with cut-off as in Eq. 11 no particle can interact with one of its images and thus be exposed to the artificial periodic box structure which is imposed upon the system. For long range forces, also interactions of far away particles have to be included, thus for such systems the periodic box structure is superimposed although they are actually isotropic. Therefore, one only takes into account those contributions to the energy of each one of the particles which is contributed by a particles that lies within a cut-off radius that is at the most 1/2L B with boxlenth L B .T h i s procedure is called minimum image convention. Using the minimum image convention, each particle interacts with at the most (N − 1) particles. Particularly for ionic systems a cut-off has to be chosen such that the electro-neutrality is not violated. Otherwise, particles would start interacting with their periodic images which would render all calculations of forces and energies erroneous.

Complex formation of charged macromolecules
A large variety of synthetic and biological macromolecules are polyelectrolytes (Manning, 1969). The most well-known polyelectrolyte biopolymers, proteins, DNA and RNA, are responsible for functions in living systems which are incomparably more complex and diverse than the functions usually discussed for synthetic polymers present in the chemical industry. For example, polyacrylic acid is the main ingredient for diapers and dispersions of copolymers of acrylamide or methacrylamide and methacrylic acid are fundamental for cleaning water. In retrospect, during the past 30 years, despite the tremendous interest in polyelectrolytes, unlike neutral polymers (de Gennes, 1979;Flory, 1969), the general understanding of the behavior of electrically charged macromolecules is still rather poor. The contrast between our understanding of neutral and charged polymers results form the long range nature of the electrostatic interactions which introduce new length and time scales that render an analytical treatment beyond the Debye-Hückel approximation very complicated (Barrat & Joanny, 2007;Debye & Hückel, 1923). Here, the traditional separation of scales, which allows one to understand properties in terms of simple scaling arguments, does not work in many cases. Experimentally, often a direct test of theoretical concepts is not possible due to idealizing assumptions in the theory, but also because of a lack of detailed control over the experimental system, e.g. in terms of the molecular weight. Quite recently, there has been increased interest in hydrophobic polyelectrolytes which are water soluble, covalently bonded sequences of polar (ionizable) groups and hydrophobic groups which are not (Khoklov & Khalatur, 2005). Many solution properties are known to be due to a complex interplay between short-ranged hydrophobic attraction, long-range Coulomb effects, and the entropic degrees of freedom. Hence, such polymers can be considered as highly simplified models of biologically important molecules, e.g. proteins or lipid bilayers in cell membranes.
In this context, computer simulations are a very important tool for the detailed investigation of charged macromolecular systems. A comprehensive review of recent advances which have been achieved in the theoretical description and understanding of polyelectrolyte solutions can be found in (Holm et al., 2004).

Two oppositely charged macromolecules
The investigation of aggregates between oppositely charged macromolecules plays an important role in technical applications, particularly in biological systems. For example, DNA is associated with histone proteins to form the chromatin. Aggregates of DNA with cationic polymers or dendrimers are discussed in the context of their possible application as DNA vectors in gene therapies (Gössl et al., 2002;Yamasaki et al., 2001). Here, we present MD simulations of two flexible, oppositely charged polymer chains and illustrate the universal scaling properties of the resulting polyelectrolyte complexes that are formed when the chains collapse and build compact, cluster-like structures which are constrained to a small region in space (Steinauser, 1998;Winkler et al., 2002). The properties are investigated as a function of chain length N and interaction strength ξ. Starting with Eq. 5 and using k = 1 (cgs-system of units) the dimensionless interaction parameter can be introduced, where the Bjerrum length ξ B is given by: where k B is the Boltzmann constant, T is temperature, ǫ is the energy scale from the Lennard-Jones potential of Eq. 11, σ defines the length scale (size of one monomer) and e is the electronic charge.
using the first term of the bonded potential of Eq. 2. The interaction with the solvent is taken into account by a stochastic force ( Γ i ) and a friction force with a damping constant χ,a c t i n g on each mass point. The equations of motion of the system are thus given by the Langevin equations The force F i comprises the force due to the sum of the potentials of Eq. 11 with cutoff r cut = 1.5, Eq. 6 with k = 1, z i/j = ±1, and the first term on the right-hand side of the bonded potential in Eq. 2 with κ = 5000ε/σ and bond length l 0 = σ 0 = 1.0. The stochastic force Γ i is assumed to be stationary, random, and Gaussian (white noise). The electrically neutral system is placed in a cubic simulation box and periodic boundary conditions are applied for the intermolecular Lennard-Jones interaction according to Eq. 11, thereby keeping the density ρ = N/V = 2.1 × 10 −7 /σ 3 constant when changing the chain length N. The number of monomers N per chain was chosen as N = 10, 20, 40, 80 and 160 so as to cover at least one order of magnitude. For the Coulomb interaction a cutoff that is half the boxlength r cut = 1/2L B was chosen. This can be done as the eventually collapsed polyelectolyte complexes which are analyzed are confined to a small region in space which is much smaller than r cut . In the following we discuss exemplarily some scaling properties of charged linear macromolecules in the collapsed state. The simulations are started with two well separated and equilibrated chains in the simulation box. After turning on the Coulomb interactions with opposite charges z i/j = ±1 in the monomers of both chains, the chains start to attract each other. In a first step during the aggregation process the chains start to twist around each other and form helical like structures as presented in Fig. 5. In a second step, the chains start to form a compact globular structure because of the attractive interactions between dipoles formed by oppositely charged monomers, see the snapshots in Fig. 6(a).  (Srivastava & Muthukumar, 1994). The change of R g is connected with a change of the density ρ of the polyelectrolyte aggregate. However, in Fig. 6(b), which presents an example of ρ for ξ = 4, only a slight dependence of the density on the chain length N can be observed. ρ measures the radial monomer density with repsect to the center of mass of the total system. For longer chains, there is a plateau while for short chains there is a pronounced maximum of the density for small distances from the center of mass. While this maximum vanished with decreasing ξ it appears also at higher interaction strengths for longer chains. Monomers on the outer part of the polyelectrolyte complex experience a stronger attraction by the inner part of the cluster than the monomers inside of it, and for smaller ξ,c h a i n so f different lengths are deformed to different degrees which leads to a chain length dependence of the density profile.

Equilibrium and Non-Equilibrium Molecular Dynamics (NEMD)
An understanding of the behavior of fully flexible linear polymers in dilute solutions and of dense melts has been achieved decades ago by the fundamental works of Rouse and Zimm (Zimm, 1956), as well as of Doi and Edwards (Doi & Edwards, 1986) and deGennes (de Gennes, 1979). In contrast to the well understood behavior of fully flexible linear polymers in terms of the Rouse (Prince E. Rouse, 1953) and the reptation model (de Gennes, 1979;Doi & Edwards, 1986), semiflexible polymer models have received increasing attention recently, as on the one hand they can be applied to many biopolymers like actin filaments, proteins or DNA (Käs et al., 1996;Ober, 2000) and on the other hand even for polymers considered flexible, like Polyethylene, the Rouse model fails as soon as the local chemical structure can no longer be neglected. One of the effects of this structure is a certain stiffness in the polymer chain due to the valence angles of the polymer backbone (Paul et al., 1997). Thus, semiflexible polymers are considerably more difficult to treat theoretically, as they require the fulfillment of additional constraints, such as keeping the total chain length fixed, which render these models more complex and often nonlinear. A deeper understanding of the rheological, dynamical and structural properties of semiflexible or even rod-like polymers is thus of great practical and fundamental interest.

Theoretical framework
The Kratky-Porod chain model (or worm-like chain model, WLC) (Kratky & Porod, 1949) provides a simple description of inextensible semiflexible polymers with positional fluctuations that are not purely entropic but governed by their bending energy Φ bend and characterized e.g. by their persistence length L p . The corresponding elastic energy of the inextensible chain of length L depends on the local curvature of the chain contour s,w h e r e r(s) is the position vector of a mass point (a monomer) on the chain and κ is a constant (Doi & Edwards, 1986).

Harris and Hearst formulated an equation of motion for the WLC model by applying
Hamilton's principle with the constraint that the second moment of the total chain length be fixed and obtained the following expressions for the bending F bend and tension forces F tens Applying this result to elastic light scattering, this model yields correct results in the flexible coil limit (Harris & Hearst, 1966), but it fails at high stiffness, where it deviates from the solution obtained for rigid rods (Harris & Hearst, 1967).
A different model was proposed by Soda (Soda, 1973), where the segmental tension forces are modeled by stiff harmonic springs. This approach avoids large fluctuations in the contour length but has the disadvantage that an analytic treatment of the model is possible only for few limiting cases. Under the assumption that the longitudinal tension relaxes quickly, the bending dynamics can be investigated using a normal mode analysis (Aragón & Pecora, 1985;Soda, 1991). However, this approach cannot account for the flexible chain behavior which is observed on large length scales in the case L p ≪ L.
Winkler, Harnau and Reineker (Harnau et al., 1996;R.G. Winkler, 1994) considered a Gaussian chain model and used a Langevin equation similar to the equation employed by Harris and Hearst, but introduced separate Lagrangian multipliers for the end points of the chain, thus avoiding the problems of the Harris and Hearst equation in the rod-like limit. Thus, the equation used in (Harris & Hearst, 1966) is contained in the model used by Winkler, Harnau and Reineker and can be regained by setting all Lagrangian multipliers equal along the chain contour and at the end points. The expansion of the position vector r(s) in normal coordinates in the approach used in (R.G. Winkler, 1994) and resolving the obtained equations for the relaxation times τ p or the normal mode amplitudes X p (t) leads to a set of transcendental equations, the solution of which cannot be given in closed form. For some limiting cases, Harnau, Winkler and Reineker showed the agreement of the approximate solution of the transcendental equations with atomistic simulation results of a n-C 100 H 202 polymer melt (Harnau et al., 1999) that were performed by Paul, Yoon and Smith (Paul et al., 1997).
In contrast to fully flexible polymers, the modeling of semiflexible and stiff macromolecules has received recent attention, because such models can be successfully applied to biopolymers such as proteins, DNA, actin filaments or rodlike viruses (Bustamante et al., 1994;Ober, 2000). Biopolymers are wormlike chains with persistence lengths l p (or Kuhn segment lengths l K ) comparable to or larger than their contour length L and their rigidity and relaxation behavior are essential for their biological functions.

Modeling and simulation of semiflexible macromolecules
Molecular Dynamics simulations were performed using the MD simulation package "MD-Cube", which was originally developed by Steinhauser . A coarse-grained bead-spring model with excluded volume interactions as a model for dilute solutions of polymers in solvents of varying quality, respectively for polymer melts, is employed (Steinhauser, 2008;. A compiler switch allows for turning on and off the interaction between different chains. Thus, one can easily switch the type of simulation from single polymers in solvent to polymer melts. The excluded volume for each monomer is taken into account through the potential of Eq. 11.
Neighboring mass points along the chains are connected by harmonic bonds with the following potential for the bonded interactions which is often used in polymer simulations of charged, DNA-like biopolymers, see e.g. (Steinauser, 1998;Winkler et al., 2002). Note, that the potential in Eq. 33 corresponds to the first term on the right-hand side of Eq. 2. In order to keep fluctuations of the bond lengths and thus the fluctuations of the overall chain length L small (below 1%), a large value for the force constant K = 10000ε/σ 2 is chosen, where ε and σ are parameters of the truncated LJ-potential in Eq. 11. The average bondlength d 0 is taken to be 0.97σ which is the equilibrium distance of a potential that is composed of the FENE (Finite Extensible Non-Linear Elastic) potentialwhich is frequently employed in polymer simulations  -and the truncated LJ-potential of Eq. 11. In combination with the LJ-potential this particle distance keeps the chain segments from artificially crossing each other (Steinhauser, 2008). The FENE potential exhibits very large fluctuations of bond lengths which are unrealistic for the investigation of semiflexible or stiff polymers. This is the reason for choosing the simple harmonic potential in Eq. 33). It is noted, that in principle the exact analytic form of the bonded potential when using a coarse-grained polymer model is actually irrelevant as long as it ensures that the the basic properties of polymers are modeled correctly such as the specific connectivity of monomers in a chain, the non-crossability of monomer segments (topological constraints), or the flexibility/stiffness of a chain. Thus, very often, a simple potential that can be quickly calculated in a numerical approach is used.
The stiffness, i.e. the bending rigidity of the chains composed of N mass points, is introduced into the coarse-grained model by the following bending potential where u is the unit bond vector u i =( r i+1 − r i )/| r i+1 − r i | connecting consecutive monomers, and r i is the position vector to the i-th monomer. The total force acting on monomer i is thus given by Assuming a NVT ensemble at temperature k B T = 1ε, the trajectories of all particles i = (1, ..., N) are generated by integrating the stochastic equations of motion with a particle friction coefficient ξ and a Gaussian stochastic force F i S that satisfies and the fluctuation-dissipation-theorem The equations of motion are integrated using the Brownian dynamics algorithm proposed by van Gunsteren and Berendsen for a canonical ensemble (van Gunsteren & Berendsen, 1982) which -for vanishing particle friction ξ -changes into the velocity Verlet algorithm for a microcanonical ensemble. The algorithm is used with a constant timestep of Δt = 5 × 10 −3 τ, where τ is the time unit of the simulation.

Results
The crossover scaling from coil-like, flexible structures on large length scales to stretched conformations at smaller scales can be seen in the scaling of the structure function S(q) when performing simulations with different values of k θ (Steinhauser, Schneider & Blumen, 2009).
In Fig. 7(a) the structure functions of the simulated linear polymer chains of length N = 700 are displayed for different persistence lengths. The chains show a scaling according to q ν . The stiffest chains exhibit a q −1 -scaling which is characteristic for stiff rods. The dotted and dashed lines display the expected theoretical scaling behavior.
Thus, by varying parameter k θ , the whole range of bending stiffness of chains from fully flexible chains to stiff rods can be covered. The range of q-values for which the crossover from flexible to semiflexible and stiff occurs, shifts to smaller q-values with increasing stiffness k θ of the chains. The scaling plot in Fig. 7(b) shows that the transition occurs for q ≈ 1/l K ,i.e.ona length scale of the order of the statistical Kuhn length. In essence, only the fully flexible chains (red data points) exhibit a deviation from the master curve on large length scales (i.e. small q-values), which corresponds to their different global structure compared with semi-flexible macromolecules. Examples for snapshots of stiff and semiflexible chains are displayed in Fig. 8.
For a theoretical treatment, following Doi and Edwards (Doi & Edwards, 1986), we expand the position vector r(s, t) of a polymer chain, parameterized with time t and contour length s  in normal modes X p (t) as follows: Resolving Eq. 39 for the normal modes and inserting the result into the Langevin equations of motion for r(s, t) one obtains after some algebraic manipulations for the relaxation time τ p which can be interpreted physically as arising from contributions due to an entropic force term ∝ p 2 and a bending force term ∝ p 4 . In Fig. 9 we show that our simulation results for semiflexible chains scale according to Eq. 40. Figure 10 exhibits the results of a NEMD step-shear simulation, from which the shear modulus G(t) has been determined. The NEMD scheme produces the same results for the shear modulus G(t) as conventional methods at equilibrium, which are based on the Green-Kubo equation.
Finally, in Fig. 11 we illustrate our step-shear simulation scheme for polymer melts with two corresponding NEMD simulation snapshots.

Shock wave failure of granular materials
In the following we discuss a recently proposed concurrent multiscale approach for the simulation of failure and cracks in brittle materials which is based on mesoscopic particle dynamics, the Discrete Element Method (DEM), but which allows for simulating macroscopic properties of solids by fitting only a few model parameters (Steinhauser, Grass, Strassburger & Blumen, 2009

Introduction: Multiscale modeling of granular materials using MD
Instead of trying to reproduce the geometrical shape of grains on the microscale as seen in two-dimensional (2D) micrographs, in the proposed approach one models the macroscopic solid state with soft particles, which, in the initial configuration, are allowed to overlap, cf. Fig. 12(a). The overall system configuration, see Fig. 12(b), can be visualized as a network of links that connect the centers of overlapping particles, cf. Fig. 12(c).
The degree of particle overlap in the model is a measure of the force that is needed to detach particles from each other. The force is imposed on the particles by elastic springs. This simple model can easily be extended to incorporate irreversible changes of state such as plastic flow in metals on the macro scale. However, for brittle materials, where catastrophic failure occurs after a short elastic strain, in general, plastic flow behavior can be completely neglected. Additionally, a failure threshold is introduced for both, extension and compression of the springs that connect the initial particle network. By adjusting only two model parameters for the strain part of the potential, the correct stress-strain relationship of a specific brittle material as observed in (macroscopic) experiments can be obtained. The model is then applied to other types of external loading, e.g. shear and high-speed impact, with no further model adjustments, and the results are compared with experiments performed at EMI.
Fig. 12. The particle model as suggested in (Steinhauser, Grass, Strassburger & Blumen, 2009). (a) Overlapping particles with radii R 0 and the initial (randomly generated) degree of overlap indicated by d 0 ij . Here, only two particles are displayed. In the model the number of overlapping particles is unlimited and each individual particle pair contributes to the overall pressure and tensile strength of the solid. (b) Sample initial configuration of overlapping particles (N = 2500) with the color code displaying the coordination number: red (8), yellow (6), and green (4). (c) The same system displayed as an unordered network.

Model potentials
The main features of a coarse-grained model in the spirit of Occam's razor with only few parameters, are the repulsive forces which determine the materials resistance against pressure and the cohesive forces that keep the material together. A material resistance against pressure load is introduced by a simple Lennard Jones type repulsive potential Φ ij rep which acts on every pair of particles {ij} once the degree of overlap d ij decreases compared to the initial overlap d 0 ij : Parameter γ scales the energy density of the potential and prefactor R 0 3 ensures the correct scaling behavior of the calculated total stress Σ ij σ ij = Σ ij F ij /A which, as a result, is independent of N. Figure 13 shows that systems with all parameters kept constant, but only N varied, lead to the same slope (Young's modulus) in a stress-strain diagram. In Eq. 41 R 0 is the constant radius of the particles, d ij = d ij (t) is the instantaneous mutual distance of each interacting pairs {ij} of particles, and d 0 ij = d ij (t = 0) is the initial separation which the pair {ij} had in the starting configuration. Every single pair {ij} of overlapping particles is associated with a different initial separation d 0 ij and hence with a different force. The minimum of each individual particle pair {ij} is chosen such that the body is force-free at the start of the simulation. Here, only the 2D case is shown for simplicity. The original system (Model M a )withedgelengthL 0 and particle radii R 0 is downscaled by a factor of 1/a into the subsystem Q A of M A (shaded area) with edge length L, while the particle radii are upscaled by factor a.A sar e s u l t ,m o d e l M B of size aL = L 0 is obtained containing much fewer particles, but representing the same macroscopic solid, since the stress-strain relation (and hence, Young's modulus E)upon uni-axial tensile load is the same in both models. (b) Young's modulus E of systems with different number of particles N in a stress-strain (σ − ε) diagram. In essence, E is indeed independent of N.
When the material is put under a low tension the small deviations of particle positions from equilibrium will vanish as soon as the external force is released. Each individual pair of overlapping particles can thus be visualized as being connected by a spring, the equilibrium length of which equals the initial distance d 0 ij . This property is expressed in the cohesive potential by the following equation: In this equation, λ (which has dimension [energy/length]) determines the strength of the potential and prefactor R 0 again ensures a proper intrinsic scaling behavior of the material response. The total potential is the following sum: The repulsive part of Φ tot acts only on particle pairs that are closer together than their mutual initial distance d 0 ij , whereas the harmonic potential Φ coh either acts repulsively or cohesively, depending on the actual distance d ij . Failure is included in the model by introducing two breaking thresholds for the springs with respect to compressive and to tensile failure, respectively. If either of these thresholds is exceeded, the respective spring is considered to be broken and is removed from the system. A tensile failure criterium is reached when the overlap between two particles vanishes, i.e. when: Failure under pressure load occurs when the actual mutual particle distance is less by a factor α (with α ∈ (0, 1)) than the initial mutual particle distance, i.e. when Particle pairs without a spring after pressure or tensile failure still interact via the repulsive potential and cannot move through each other.
An appealing feature of this model, as opposed to many other material models used for the description of brittle materials, see e.g. (Cundall & Strack, 1979;Leszczynski, 2003;Walton & Braun, 1986), is its simplicity. The proposed model has a total of only three free parameters: γ and λ for the interaction potentials and α for failure. These model parameters can be adjusted to mimic the behavior of specific materials.

Shock wave simulations and comparison with experiments
Finally, in Fig. 14, non-equilibrium MD simulation (NEMD) results for systems with varying shock impact velocities are presented and compared with high-speed impact experiments performed at EMI with different ceramic materials (Al 2 O 3 and SiC) in the so-called edge-on-impact (EOI) configuration. These oxide and non-oxide ceramics represent two major classes of ceramics that have many important applications. The impactor hits the target at the left edge. This leads to a local compression of the particles in the impact area.
The top series of snapshots in Fig. 14(a) shows the propagation of a shock wave through the material. The shape of the shock front and also the distance traveled by it correspond very well to the high-speed photographs in the middle of Fig. 14(a). These snapshots were taken at comparable times after the impact had occurred in the experiment and in the simulation, respectively. In the experiments which are used for comparison, specimens of dimensions (100 × 100 × 10)mm were impacted by a cylindrical blunt steel projectile of length 23 mm, mass m = 126 g and a diameter of 29.96 mm (Steinhauser et al., 2006). After a reflection of the pressure wave at the free end of the material sample, and its propagation back into the material, the energy stored in the shock wave front finally disperses in the material. One can study in great detail the physics of shock waves traversing the material and easily Fig. 14. Results of a simulation of the edge-on-impact (EOI) geometry, except this time, the whole macroscopic geometry of the experiment can be simulated while still including a microscopic resolution of the system. The impactor is not modeled explicitly, but rather its energy is transformed into kinetic energy of the particle bonds at the impact site. (a) Top row: A pressure wave propagates through the material and is reflected at the free end as a tensile wave (not shown). Middle row: The actual EOI experiment with a SiC specimen. The time interval between the high-speed photographs is comparable with the simulation snapshots above. The red arrows indicate the propagating shock wave front. Bottom row: The same simulation run but this time only the occurring damage in the material with respect to the number of broken bonds is shown. (b) Number of broken bonds displayed for different system sizes N, showing the convergence of the numerical scheme. Simulation parameters (α, γ, λ) are chosen such that the correct stress-strain relations of two different materials (Al 2 O 3 and SiC) are recovered in the simulation of uniaxial tensile load. The insets show high-speed photographs of SiC and Al 2 O 3 , respectively, 4μs after impact.
identify strained or compressed regions by plotting the potential energies of the individual pair bonds. Also failure in the material can conveniently be visualized by plotting only the failed bonds as a function of time, cf. the bottom series of snapshots in Fig. 14(a). A simple measure of the degree of damage is the number of broken bonds with respect to the their total initial number. This quantity is calculated from impact simulations of Al 2 O 3 and SiC,a f t e r previously adjusting the simulation parameters γ, λ and α accordingly. Figure 14(b) exhibits the results of this analysis. For all impact speeds the damage in the SiC-model is consistently larger than in the one for Al 2 O 3 which is also seen in the experiments.
The impactor is not modeled explicitly, but rather its total kinetic energy is transformed into kinetic energy of the particles in the impact region. Irreversible deformations of the particles such as plasticity or heat are not considered in the model, i.e. energy is only removed from the system by broken bonds. Therefore, the development of damage in the material is slightly overestimated.

Conclusions
In summary, we presented an introduction into the molecular dynamics method, discussing the choice of potentials and fundamental algorithms for the implementation of the MD method. We then discussed proto-typical applications of MD, namely the collapse of two oppositely charged macromolecules (polyelectrolytes) and the simulation of semiflexible bio-macromolecules 6 . We demonstrated how semiflexibility, or stiffness of polymers can be included in the potentials describing the interactions of particles. We finally showed a somewhat unusual application of MD in the field of solid state physics where we modeled the brittle failure behavior of a typical ceramic and simulated explicitly the set-up of corresponding high-speed impact experiments. We showed that the discussed multiscale particle model reproduces the macroscopic physics of shock wave propagation in brittle materials very well while at the same time allowing for a resolution of the material on the micro scale and avoiding typical problems (large element distortions, element-size dependent results) of Finite Elements, which constitutes a different type of discretization for simulation problems that are closely connected with macroscopic experiments. The observed failure and crack pattern in impact MD simulations can be attributed to the random initial distribution of particle overlaps which generates differences in the local strength of the material. By generating many realizations of systems with different random initial overlap distributions of particles, the average values obtained from these many simulations lead to the presented fairly accurate results when compared with experimental high-speed impact studies.