Severin, , Nikolai: Molecular Dynamics Simulations of Polymers and Micelles at Interfaces 
3
Molecular Dynamics (MD) simulations solve Newton’s equation of motion for a system of i particles.
where is the force, acting on the ith particle, m_{i} is the mass and is the acceleration of the ith particle. The force can be computed directly from the derivative of the potential energy U:
With the known expression for the potential energy, known masses and known initial positions in coordinatemomentum 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 finitedifference 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 timestep 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 timestep to be used. The Verlet leapfrog algorithm function as follows:
4
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 timestep 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 timestep should be used. However, a too large timestep causes instability in the integration process.
The timestep used depends both on the system as well as on the integration algorithm. The main limitation to the timestep is a high frequency motion which should be taken into account. The vibrational period should be splitted into at least 810 segments to satisfy the Verlet integrator. In most molecular systems, the highest vibrational frequency is the one of CH bond vibrations, whose period is on the order of 10^{14} s. The integration timestep should therefor be about 1 fs, which was used during all simulations.
5
Forcefield
A forcefield is a set of potential energy functions of interatomic 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 forcefield 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 forcefield supplied with the InsightII package (MSI) was used in all MD simulations. The pcff forcefield is adopted to use for the simulation of polymers. The exact form of the pcff forcefield is presented below^{18}.
(1) 
(2) 
(3) 
(4) 
(5) 
(6) 
(7) 
(8) 
(9) 

(10) 

(11) 
6
Figure 1 presents graphic illustrations for the forcefield terms.
Figure 1 Graphic illustration of the pcff forcefiled 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
7
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.
Cutoff
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 nonbond cutoff, i.e. to neglect the nonbond interactions for pairs of atoms separated by distances greater than a cutoff value. The use of cutoff 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/r^{6}. At r=810 Å the energies and forces are quite small compared to k_{b}T 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 dipoledipole interaction dies off as 1/r^{3}. 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 monopolemonopole interaction at r=1 nm would be misleading, while cutting off dipoledipole 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.
8
Temperature
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 MaxwellBoltzman 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, E_{kin} is the total kinetic energy of the system and N_{f } is the number of degrees of freedom. For a nonperiodic isolated molecular system N_{f} is equal to 3N6 (six degrees of freedom are subtracted because both the translation and rotation of the centre of mass are ignored). For a periodic system N_{f} is 3N3 (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 T_{new } is the desired temperature and T_{old} is the current temperature^{18}.
9
Pressure
The pressure may be calculated using the virial theorem, which is written in the form of generalised equipartition:
where p_{k} and q_{k} 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 r_{i} and forces f_{i} on the atoms:
The symbol f_{i}^{tot} 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 f_{i} as the sum of forces f_{ij} 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 coordinates^{9}.
10
In general pressure is a tensor P_{ij} 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 Berendsen^{21} and the ParinelloRahman^{22} 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 P_{0 } is the target pressure. The cartesian components of the unit cell vectors are scaled by the same factor .
With the ParrinelloRahman 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.
© 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 HumboldtUniversität zu Berlin 
HTML  Version erstellt am: Thu Feb 24 14:14:32 2000 