Severin, , Nikolai: Molecular Dynamics Simulations of Polymers and Micelles at Interfaces


Chapter 2. Molecular dynamics (MD) Simulation

Molecular Dynamics (MD) simulations solve Newton’s equation of motion for a system of i particles.

where is the force, acting on the i-th particle, mi is the mass and is the acceleration of the i-th particle. The force can be computed directly from the derivative of the potential energy U:

(Eq. 1)

With the known expression for the potential energy, known masses and known initial positions in coordinate-momentum space it is possible to solve this equation, i.e. to find the trajectory of the system. Analytically only systems of 1 or 2 independent particles can be solved. The numerical solution of Newton’s equation using an empirical fit to the potential energy surface is called molecular dynamics simulations.

A standard method of solving an ordinary differential equation such as Eq.1 numerically is the finite-difference method. The general idea is as follows. Given the initial coordinates and velocities and other dynamic information at time t, the positions and velocities at time are calculated. The time-step depends on the integration method as well as the system itself.

The basic criteria for a good integrator are: it should be fast, it should require little computer memory, it should permit the use of a relatively long time step and it must conserve the total energy.

For the MD calculations two different integration algorithms were applied, which are described below.

Verlet leap frog integrator

The advantages of the Verlet algorithm is that it requires only one energy evaluation per time step, modest memory, and it also allows a relatively large time-step to be used. The Verlet leapfrog algorithm function as follows:


Where are coordinates, is the velocity and is an acceleration of an atom. The force is evaluated from at .

The Verlet leapfrog method has one major disadvantage: positions and velocities calculated are displaced by half a time-step out of synchrony.

Verlet velocity integrator

The algorithm is as follows:

A key parameter in the integration algorithms is the integration time step . In order to use computational time more efficiently, as large as possible time-step should be used. However, a too large time-step causes instability in the integration process.

The time-step used depends both on the system as well as on the integration algorithm. The main limitation to the time-step is a high frequency motion which should be taken into account. The vibrational period should be splitted into at least 8-10 segments to satisfy the Verlet integrator. In most molecular systems, the highest vibrational frequency is the one of C-H bond vibrations, whose period is on the order of 10-14 s. The integration time-step should therefor be about 1 fs, which was used during all simulations.



A forcefield is a set of potential energy functions of inter-atomic distances and of molecular internal coordinates, which represents the energy of molecules and of assemblies of molecules normally in the electronic ground state. From the knowledge of the force-field one can calculate all molecular properties which are determined by the energy of the system. The calculation of the potential energy, along with its first and second derivatives with respect to the atomic coordinates, yields the information necessary for minimisation, harmonic vibrational analysis and dynamics simulations.

The pcff force-field supplied with the InsightII package (MSI) was used in all MD simulations. The pcff force-field is adopted to use for the simulation of polymers. The exact form of the pcff force-field is presented below18.













Figure 1 presents graphic illustrations for the force-field terms.

Figure 1 Graphic illustration of the pcff force-filed terms.

Periodic boundary conditions

The size of the systems treated by MD simulation is limited to several tens of thousand particles. This limitation is governed by the limitation of computer power. Thus if one particle represents one atom, the size of a simulated cell is limited to several nanomiters. A simulation carried out on a isolated system is a poor approximation of what would happen in a true bulk, because of the influence of boundaries. To remedy this, the cell is replicated in three


dimensions, thus the simulated cell should have translational periodicity. This is a much better representation of a bulk system, because the molecules near the surface now interact with molecules in adjacent cells. The image atoms are used to calculate energies and forces on the real atoms.


The periodic boundary conditions method allows to simulate a virtually infinite system which consist of a base cell and its images. The energies and forces on the imaged atoms are not calculated because their motions are computed as symmetry operations on the real atoms. It is now important to define long range interactions in the system. The simplest way is to allow molecules to interact with the molecule or molecular image closest to it. This method is called minimum image model. Each molecule interacts only with molecules and images within a distance of half the cell size. The disadvantage of this method is that computational time has a cubic dependence on the number of atoms in the system. Therefore it is common to introduce a non-bond cut-off, i.e. to neglect the non-bond interactions for pairs of atoms separated by distances greater than a cut-off value. The use of cut-off with a van der Waals interaction potential is quite reasonable, since the van der Waals potential is relatively short range and dies out as 1/r6. At r=8-10 Å the energies and forces are quite small compared to kbT at room temperature. Thus, it is reasonable to set the cut off length to about 10 Å for van der Waals interactions. The Coulombic interactions on the other hand die off as 1/r, so at much longer distances Coulombic interactions are not negligible. However for most of cases polymeric systems are composed of neutral groups with electric dipoles or higher multipoles. The dipole-dipole interaction dies off as 1/r3. The interaction energy of two monopols of one elementary charge at r=1 nm is about 33 kcal mol-1, while the interaction energy of two dipoles formed from elementary monopols at the distance r=1 nm is 0.3 kcal mol-1. It is clear that cutting off monopole-monopole interaction at r=1 nm would be misleading, while cutting off dipole-dipole interaction at r=1 nm is only a modest approximation.

If periodic boundary conditions are applied to a molecular system, volume and pressure can be defined. Different statistical ensembles can be generated, depending on which thermodynamic variables are kept fixed. In the present simulations the following ensembles were used: 1. NVT: number of particles, volume and temperature are constant, and 2. NPT: number of particles, pressure and temperature are constant.



The temperature is a macroscopic quantity, which is related to the microscopic description of molecular simulations through the kinetic energy, which is calculated from the atomic velocities. The temperature and the distribution of atomic velocities in a system are related by the Maxwell-Boltzman equation:


In the initial stage the velocities are generated according to Eq. 2, and the gaussian distribution is generated by a random number generator.

During MD simulation the temperature is related to the average kinetic energy of the system through the equipartition principle. This principle states that every degree of freedom has an average energy of associated with it:

where T is the instantaneous temperature, Ekin is the total kinetic energy of the system and Nf is the number of degrees of freedom. For a nonperiodic isolated molecular system Nf is equal to 3N-6 (six degrees of freedom are subtracted because both the translation and rotation of the centre of mass are ignored). For a periodic system Nf is 3N-3 (only the three degrees of freedom corresponding to translational motion can be ignored, since the rotation of a central cell imposes a torque on its neighbouring cells).

The direct Velocity scaling method was used for temperature scaling. Velocities of all atoms were uniformly scaled as follows:

where Tnew is the desired temperature and Told is the current temperature18.



The pressure may be calculated using the virial theorem, which is written in the form of generalised equipartition:

where pk and qk are generalised momentums and coordinates and H is the hamiltonian of the system. If we choose cartesian coordinates the second equation can be rewritten in terms of cartesians coordinates of atoms ri and forces fi on the atoms:

The symbol fitot is used because it represents the sum of intermolecular forces and external forces. The external forces are related to the pressure:

If we define the internal virial W

where we now restrict attention to intermolecular forces, then

For pairwise interactions it is more convenient to express W in a form, which is independent of the origin of coordinates:

The first equality follows from writing fi as the sum of forces fij on atom i due to atom j. The second equality follows because the indices i and j are equivalent. Using Newton’s third law this expression may be rewritten in the form:

where . This form is independent on the origin of coordinates9.


In general pressure is a tensor Pij where i and j subscripts run over x,y and z coordinates. The first subscript i refers to the direction normal to the plane on which the force acts and the second refers to the direction of the applied force.

Two methods of pressure control were used during MD simulation: the isotropic Berendsen21 and the Parinello-Rahman22 method. With the Berendsen method there is no change in the shape of the system, and only the pressure is controlled. At each time step, the x, y and z coordinates of each atom are scaled by the factor:

where is the time step, P is the instantaneous pressure, and P0 is the target pressure. The cartesian components of the unit cell vectors are scaled by the same factor .

With the Parrinello-Rahman method, the system’s shape can change, and therefore both pressure and stress can be controlled.

Energy minimisation

Energy minimisation is the important method for potential energy surface investigations. The minimisation finds configurations that are stable points on the potential surface. This means finding a point in the configuration space where the total force acting on each atom vanishes. By simply minimising the energy, stable conformations can be identified. The advantage of the minimisation technique is that it is very fast in comparison with MD simulations, but it minimises the system with respect to potential energy and does not take into account entropy effects. Second it does not give any information about the evolution of the system in time. It is very useful to combine both minimisation and MD simulations techniques. In the initial structure some atoms may experience strong forces from neighbouring atoms which would lead to a high acceleration of these atoms during the first steps of MD simulations. That is why it is necessary to run a minimisation before a MD simulation.

[Titelseite] [1] [2] [3] [4] [5] [6] [Bibliographie] [Lebenslauf] [Selbständigkeitserklärung]

© Die inhaltliche Zusammenstellung und Aufmachung dieser Publikation sowie die elektronische Verarbeitung sind urheberrechtlich geschützt. Jede Verwertung, die nicht ausdrücklich vom Urheberrechtsgesetz zugelassen ist, bedarf der vorherigen Zustimmung. Das gilt insbesondere für die Vervielfältigung, die Bearbeitung und Einspeicherung und Verarbeitung in elektronische Systeme.

DiML DTD Version 2.0
Zertifizierter Dokumentenserver
der Humboldt-Universität zu Berlin
HTML - Version erstellt am:
Thu Feb 24 14:14:32 2000