**Meet the editor**

Lichang Wang, Ph.D., is Professor of Chemistry and Biochemistry at Southern Illinois University Carbondale and the founder and president of MASIS, Inc. Her research activities encompass the area of molecular modeling, such as developing molecular dynamics simulation tools for complex systems and performing dynamics simulations and electronic structure calculations

to study various systems of chemical and energy interests: (a) Transition metal nanoparticles, their toxicity to humans, and their catalytic activities for dehydrogenation of CH4, reduction of O2, direct conversion of coal to liquid fuel, and reactions of ethanol; and (b) Organic dyes for solar cell applications. Details of her research can be found on the website: http:// www.chem.siu.edu/wang. The ultimate goal of her research activities is to replace as many bench top R & D experiments as possible by molecular modeling in order to achieve a greener, as well as economically better, R & D in Chemical and Energy Industries.

## Contents

### **Preface XIII**

	- **Part 2 Dynamics of Biomolecules 83**

X Contents


Chapter 18 **High-Throughput Simulations of Protein Dynamics in Molecular Machines: The 'Link' Domain of RNA Polymerase 419**  Robert O. J. Weinzierl

VI Contents

Chapter 8 **Essential Dynamics on Different** 

Héctor A. Baldoni

Chapter 9 **MM-GB(PB)SA Calculations of** 

**Part 3 Dynamics of Plasmas 191**

Chapter 11 **Application of Molecular Dynamics** 

Chapter 12 **Molecular Dynamics Simulations** 

**Part 4 Dynamics at the Interface 273**

Chapter 14 **Simulations of Unusual Properties of** 

**Part 5 Dynamics of Nanomachines 339** 

Chapter 16 **Analysis of the Atomization Process**

**Water Inside Carbon Nanotubes 297**  Yoshimichi Nakamura and Takahisa Ohno

**Molecular Dynamics to Nanofluidics 319** 

**by Molecular Dynamics Simulation 341** 

Chapter 13 **Studies of Cardio Toxin Protein** 

Chapter 15 **Applications of All-Atom** 

Mauro Chinappi

Yeh Chun-Lang

Chapter 17 **Molecular Dynamics Simulation** 

Akinjide Oluwajobi

**of Nanoscale Machining 389** 

Koji Eriguchi

**Biological Systems: Fis Protein, tvMyb1** 

Lucas J. Gutiérrez, Ricardo D. Enriz and

Chapter 10 **Micro-Heterogeneity in Complex Liquids 193** 

**of Complex (Dusty) Plasmas 245** Céline Durniak and Dmitry Samsonov

**Transcriptional Factor and BACE1 Enzyme 151**

**Protein-Ligand Binding Free Energies 171**  Joseph M. Hayes and Georgios Archontis

Aurélien Perera, Bernarda Kežić, Franjo Sokolić and Larisa Zoranić

**Metal-Oxide-Semiconductor Field-Effect Transistors 221** 

**Simulations to Plasma Etch Damage in Advanced** 

**Adsorption on Mixed Self-Assembled Monolayers Using Molecular Dynamics Simulations 275** Shih-Wei Hung, Pai-Yi Hsiao and Ching-Chang Chieng

## Preface

Molecular dynamics (MD) simulations have played increasing roles in our understanding of physical and chemical processes of complex systems and in advancing science and technology. Over the past forty years, MD simulations have made great progress from developing sophisticated theories for treating complex systems to broadening applications to a wide range of scientific and technological fields. The chapters of *Molecular Dynamics* are a reflection of the most recent progress in the field of MD simulations.

This is the second book of *Molecular Dynamics* that focuses on the MD studies of synthetic and biological macromolecules. This book is divided into five parts. The first part deals with the molecular dynamics simulations of polymers. Steinhauser provides a general introduction of MD simulations, both equilibrium and nonequilibrium, on the studies of macromolecules of hard and soft matters in Chapter 1. In Chapter 2, Sandoval presents the MD results to understand various phenomena of synthetic polymers, three amphiphilic polymers at the air-water interface and two polymers in condensed phase. Ma & Hu discuss in Chapter 3 the MD simulations on the backbone connectivity of polymer chains and the aggregation process or phase separation. In Chapter 4, Islami & Mehdipour provide a summary of MD methods to investigate the solubility and diffusion of gaseous molecules, such as Ar, N2, CO2, CH4, and C3H8, in polymers.

The second part consists of five chapters that employ MD simulations to study biomolecules. Laaksonen et al. describe in Chapter 5 their MD simulation package M.Dyna*Mix* as well as its application to study lipid bilayers and the hydration and coordination of counterions around DNA. In Chapter 6, Tsurui & Takahashi give a description of dielectric and coarse-grained models and use TCR-pMHC complexes as examples to illustrate the application of these models. Chapter 7 provides a MD study of the conformational profile of Neuromedin B by Sharma, et al. Chapter 8 summarizes the MD studies of Fis protein, tvMyb1 transcriptional factor, and BACE1 enzyme by Gutiérrez et al. In Chapter 9, Hayes & Archontis present a review of calculations in the study of protein-ligand binding events.

Part III is about the MD studies of plasmas. In Chapter 10, Perera et al. provide a description of micro-structure and micro-heterogeneity of liquid mixtures. Eriguchi

#### XIV Preface

presents MD studies of plasma etch damage mechanism in Chapter 11. Chapter 12 by Durniak & Samsonov gives a description of the MD simulations of dusty plasmas on the structure, linear and nonlinear waves, shocks, and other related phenomena.

The fourth part is on the MD simulations of interfaces that play important roles in the development of biomaterials, implant biocompatibility, and biosensors. In Chapter 13, Hung et al. describe the MD studies of cardio toxin protein adsorption on selfassembled monolayers. In Chapter 14, Nakamura & Ohno provide the MD results of water behavior inside of carbon nanotubes. In Chapter 15, Chinappi provides a summary of MD results on the study of nanofluids.

In the last part of the book, MD studies of nanomachines are discussed. The study of nanomachines is important in our understanding of how biological systems work and in the technological development of nano-instruments. In Chapter 16, Chun-Lang describes the dynamics of fluid from nanojets and provides the analysis of the atomization process based on the MD simulation results. Oluwajobi gives a detailed description of MD simulations of nanomachines and presents the MD results in the studies of cutting processes in Chapter 17. The last chapter of this book is written by Weinzierl, which deals with protein dynamics in molecular machines.

With strenuous and continuing efforts, a greater impact of MD simulations will be made on understanding various processes and on advancing many scientific and technological areas in the foreseeable future.

In closing I would like to thank all the authors taking primary responsibility to ensure the accuracy of the contents covered in their respective chapters. I also want to thank my publishing process manager Ms. Daria Nahtigal for her diligent work and for keeping the book publishing progress in check.

> **Lichang Wang**  Department of Chemistry and Biochemistry Southern Illinois University Carbondale USA

**Part 1** 

**Dynamics of Polymers** 

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

Martin Oliver Steinhauser

*Research Group Shock Waves in Soft Biological Matter, Department of Composite Materials Fraunhofer Institute for High-Speed Dynamics, Ernst-Mach-Institut, EMI, Eckerstrasse 4, Freiburg Germany*

#### **1. Introduction**

Today, computer experiments play a very important role in science. In the past, physical sciences were characterized by an interplay between experiment and theory. In theory, a model of the system is constructed, usually in the form of a set of mathematical equations. This model is then validated by its ability to describe the system behavior in a few selected cases, simple enough to allow a solution to be computed from the equations. One might wonder why one does not simply derive all physical behavior of matter from an as small as possible set of fundamental equations, e.g. the Dirac equation of relativistic quantum theory. However, the quest for the fundamental principles of physics is not yet finished; thus, the appropriate starting point for such a strategy still remains unclear. But even if we knew all fundamental laws of nature, there is another reason, why this strategy does not work for ultimately predicting the behavior of matter on any length scale, and this reason is the growing complexity of fundamental theories – which are based on the dynamics of particles – when they are applied to systems of macroscopic (or even microscopic) dimensions. In almost all cases, even for academic problems involving only a few particles1, a strict analytical solution is not possible and solving the problem very often implies a considerable amount of simplification. In contrast to this, in experiments, a system is subject to measurements, and results are collected, very often in the form of large data sets of numbers from which one strives to find mathematical equations describing the data by generalization, imagination and by thorough investigation. Very rarely, normally based on symmetries which allow inherent simplifications of the original problem, is an analytical solution at hand which describes exactly the evidence of the experiment given by the obtained data sets. Unfortunately, many academic and practical physical problems of interest do not fall under this category of "simple" problems, e.g. disordered systems, where there is no symmetry which helps to simplify the treatment.

The advent of modern computers which basically arose from the Manhattan project in the United States in the 1950s added another element to (classical) experiment and theory, namely

<sup>1</sup> In fact, already the three-body problem which involves three coupled ordinary differential equations is not solvable analytically.

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" experiments3, too expensive to 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 mathematics4, 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 analytically5. 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

<sup>2</sup> In this context one has to realize that often, in experiments, too, considerable simplifications of the investigated "real system" are done, e.g. when preparing it in a particular state in terms of pressure, temperature or other degrees of freedom.

<sup>3</sup> For example, systems at a pressure comparable to that in the interior of the sun.

<sup>4</sup> The famous David Hilbert once commented the question of what became of one of his students: "He became a writer - he didn't have enough imagination."

<sup>5</sup> For example, when computing the phase diagram of a substance modeled by a certain force law.

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

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" experiments3, too expensive to 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 mathematics4, physics and other scientific areas, are very important qualities of a computer

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 analytically5. 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

<sup>2</sup> In this context one has to realize that often, in experiments, too, considerable simplifications of the investigated "real system" are done, e.g. when preparing it in a particular state in terms of pressure,

<sup>4</sup> The famous David Hilbert once commented the question of what became of one of his students: "He

<sup>5</sup> For example, when computing the phase diagram of a substance modeled by a certain force law.

<sup>3</sup> For example, systems at a pressure comparable to that in the interior of the sun.

realistic systems can be investigated.

temperature or other degrees of freedom.

became a writer - he didn't have enough imagination."

simulator!

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.

#### **2. 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.

#### **2.1 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

$$m\_i \vec{\vec{r}}\_i = \vec{F}\_i = -\frac{\partial}{\partial \vec{r}\_i} \Phi\_i \tag{1}$$

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 *Fi* acting on the atoms, which are usually derived from a potential energy function Φ(*rN*), where *r<sup>N</sup>* = (*r*1,*r*2, ...,*rN*) represents the complete set of 3*N* atomic coordinates.

#### **2.1.1 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 *li*, bond angles *θ* and torsion angles *φ*. When 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*<sup>0</sup> and angles *θ*0, *φ*<sup>0</sup> 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:

$$\Phi\_{\text{bonded}}(r,\theta,\phi) = \frac{\kappa}{2} \sum\_{i} (|\vec{r}\_i - \vec{r}\_{i-1} - l\_0|)^2 + \frac{k\_\theta}{2} \sum\_{k} (\theta\_k - \theta\_0)^2 + \frac{\beta}{2} \sum\_{m} (\phi\_m - \phi\_0)^2. \tag{2}$$

Here, the summation indices sum up the number of bonds *i* at positions �*ri*, 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, *θ*<sup>0</sup> and *φ*<sup>0</sup> 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*<sup>0</sup> (Steinhauser, 2005):

$$\Phi\_{FENE}(r) = \begin{cases} -\frac{1}{2}\kappa R\_0^2 \ln(1 - \frac{r^2}{R\_0^2}) & r < R\_{0\prime} \\ \infty & \text{otherwise.} \end{cases} \tag{3}$$

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).

Fig. 1. Illustration of the potential parameters used for modeling bonded interactions.

#### **2.1.2 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)

$$\Phi\_{nb}(\vec{r}) = \sum\_{i}^{N} \Phi\_1(\vec{r}\_i) + \sum\_{i}^{N} \sum\_{j>i}^{N} \Phi\_2(\vec{r}\_i \cdot \vec{r}\_j) + \sum\_{i}^{N} \sum\_{j>i}^{N} \sum\_{k>j>i}^{N} \Phi\_3(\vec{r}\_i \cdot \vec{r}\_j \cdot \vec{r}\_k) + \cdots \tag{4}$$

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 <sup>O</sup>(*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* = |�*ri* −�*rj*| is a potential of the following form:

$$\Phi\_{nb}(r) = \Phi\_{Coulomb}(r) + \left(\frac{C\_1}{r}\right)^{12} + \left(\frac{C\_2}{r}\right)^6. \tag{5}$$

Parameters *C*<sup>1</sup> and *C*<sup>2</sup> are parameters of the attractive and repulsive interaction and the electrostatic energy Φ*Coulomb*(*r*) between the particles with position vectors�*ri* and�*rj* is given by:

$$\Phi\_{\text{Coulomb}}(r) = \frac{1}{\varepsilon}k \cdot \sum\_{i} \sum\_{j>i} \frac{z\_i z\_j e^2}{|\vec{r}\_i - \vec{r}\_j|}. \tag{6}$$

The constant *k* = 1 in the cgs–system of units and *�* is the dielectric constant of the medium, for example *�air* = 1 for air, *�prot* = 4 for proteins or *�H*<sup>20</sup> = 82 for water. The variables *zi* and *zj* 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)

$$\Phi\_{a,b}(r) = a\varepsilon \left[ \left( \frac{\sigma\_0}{r} \right)^a + \left( \frac{\sigma\_0}{r} \right)^b \right],\tag{7}$$

where

4 Will-be-set-by-IN-TECH

considers the internal degrees of freedom of polymers and biomacromolecules by using generic potentials that describe bond lengths *li*, bond angles *θ* and torsion angles *φ*. When 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*<sup>0</sup> and angles *θ*0, *φ*<sup>0</sup> 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

Here, the summation indices sum up the number of bonds *i* at positions �*ri*, 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, *θ*<sup>0</sup> and *φ*<sup>0</sup> are the

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

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.,

Fig. 1. Illustration of the potential parameters used for modeling 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

<sup>0</sup> ln(<sup>1</sup> <sup>−</sup> *<sup>r</sup>*<sup>2</sup>

*R*2 0

*kθ* <sup>2</sup> ∑ *k*

(*θ<sup>k</sup>* − *θ*0)

) *r* < *R*0,

<sup>∞</sup> otherwise. (3)

<sup>2</sup> <sup>+</sup> *<sup>β</sup>* <sup>2</sup> ∑ *m*

(*φ<sup>m</sup>* − *φ*0)

2. (2)

can be treated by using some or all parts of the following potential term:

(|�*ri* <sup>−</sup>�*ri*−<sup>1</sup> <sup>−</sup> *<sup>l</sup>*<sup>0</sup> <sup>|</sup>)<sup>2</sup> <sup>+</sup>

<sup>2</sup> ∑ *i*

equilibrium distance, bond angle and torsion angle, respectively.

length of a polymer bond to a prefixed value *R*<sup>0</sup> (Steinhauser, 2005):

 −1 <sup>2</sup> *<sup>κ</sup>R*<sup>2</sup>

Φ*FENE*(*r*) =

1996; Siu et al., 2008; Steinhauser, 2008).

**2.1.2 Non-bonded interactions**

<sup>Φ</sup>*bonded*(*r*, *<sup>θ</sup>*, *<sup>φ</sup>*) = *<sup>κ</sup>*

$$\alpha = \frac{1}{a-b} \left( \frac{a^a}{b^b} \right)^{\frac{1}{a-b}}, \Phi\_{\text{min}} = \varepsilon \text{ and } \Phi(\sigma) = 0. \tag{8}$$

The most often used LJ-(6-12) potential for the interaction between two particles with a distance *r* = |�*ri* −�*rj*| then reads (cf. Eq. 5):

$$\Phi\_{LJ}(r) = 4\varepsilon \left[ \left(\frac{\sigma\_0}{r}\right)^{12} + \left(\frac{\sigma\_0}{r}\right)^6 \right]. \tag{9}$$

Parameter *ε* determines the energy scale and *σ*<sup>0</sup> 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 *kB* = 1.38 × 10−23J/K. In these reduced (simulation) units, one MD timestep is measured in units of *τ*ˆ = (*mσ*2/*ε*)1/2, where *m* is the mass of a particle and *ε* and *σ*<sup>0</sup> are often simply set to *σ*<sup>0</sup> = *ε* = *kBT* <sup>=</sup> 1. Applied to real molecules, for example to Argon with *<sup>m</sup>* <sup>=</sup> 6.63 <sup>×</sup> <sup>10</sup>−23*kg*, *<sup>σ</sup>*<sup>0</sup> <sup>≈</sup> 3.4 <sup>×</sup> <sup>10</sup>−10*<sup>m</sup>* and *<sup>ε</sup>*/*kB* <sup>≈</sup> <sup>120</sup>*<sup>K</sup>* one obtains a typical MD timestep of *<sup>τ</sup>*<sup>ˆ</sup> <sup>≈</sup> 3.1 <sup>×</sup> <sup>10</sup>−13*s*.

Using an exponential function instead of the repulsive *r*−<sup>12</sup> term, one obtains the *Buckingham potential* (Buckingham, 1938):

$$\Phi(r) = b \exp(-ar) - \frac{c}{r^6} - \frac{d}{r^8}. \tag{10}$$

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 *<sup>r</sup>* <sup>=</sup> *<sup>r</sup>*min <sup>=</sup> <sup>√</sup><sup>6</sup> <sup>2</sup> and shifted to positive values by *<sup>ε</sup>* such that it is purely repulsive and smooth at *<sup>r</sup>* <sup>=</sup> *<sup>r</sup>*cut <sup>=</sup> <sup>√</sup><sup>6</sup> 2:

$$\Phi\_{Lf}^{\text{cut}}(r) = \begin{cases} 4\ \varepsilon \left\{ \left(\frac{\sigma\_0}{r}\right)^{12} - \left(\frac{\sigma\_0}{r}\right)^6 \right\} + \varepsilon & r \le 2^{1/6}\sigma\_{0\prime} \\ 0 & \text{otherwise.} \end{cases} \tag{11}$$

Another extension of the potential in Eq. 9 is proposed in (Steinhauser, 2005) where a smooth attractive part is introduced again, in order to allow for including different solvent qualities of the solvent surrounding the polymer:

$$\Phi\_{\cos}(r) = \left[\frac{1}{2} \cdot \cos(\alpha r^2 + \beta) + \gamma\right] \epsilon. \tag{12}$$

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* = 21/6 and at *r* = *r*cut, while it is zero at *r* = *r*cut and has value *γ* at *r* = 21/6, where *γ* is the depth of the attractive part. Further details can be found in (Steinhauser, 2005). When setting *r*cut = 1.5 one sets *γ* = − and obtains *α* and *β* as solutions of the linear set of equations

$$2^{1/3}\mathfrak{a} + \mathfrak{B} \quad = \quad \mathfrak{m} \tag{13}$$

$$
\mathfrak{a}2.25\mathfrak{a} + \mathfrak{k} = 2\pi.\tag{14}
$$

The total unbounded potential can then be written as:

$$\Phi\_{\text{Total}}(r,\lambda) = \begin{cases} \Phi\_{LI}^{\text{cut}}(r) - \lambda \varepsilon & 0 < r < 2^{1/6} \sigma\_0 \\ \lambda \, \Phi\_{\text{cos}}(r) & 2^{1/6} \sigma\_0 \le r < r\_{\text{cut}} \\ \infty & \text{otherwise} \end{cases} \tag{15}$$

10−23J/K. In these reduced (simulation) units, one MD timestep is measured in units of *τ*ˆ = (*mσ*2/*ε*)1/2, where *m* is the mass of a particle and *ε* and *σ*<sup>0</sup> are often simply set to *σ*<sup>0</sup> = *ε* = *kBT* <sup>=</sup> 1. Applied to real molecules, for example to Argon with *<sup>m</sup>* <sup>=</sup> 6.63 <sup>×</sup> <sup>10</sup>−23*kg*, *<sup>σ</sup>*<sup>0</sup> <sup>≈</sup> 3.4 <sup>×</sup> <sup>10</sup>−10*<sup>m</sup>* and *<sup>ε</sup>*/*kB* <sup>≈</sup> <sup>120</sup>*<sup>K</sup>* one obtains a typical MD timestep of *<sup>τ</sup>*<sup>ˆ</sup> <sup>≈</sup> 3.1 <sup>×</sup> <sup>10</sup>−13*s*.

Using an exponential function instead of the repulsive *r*−<sup>12</sup> term, one obtains the *Buckingham*

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

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 *<sup>r</sup>* <sup>=</sup> *<sup>r</sup>*min <sup>=</sup> <sup>√</sup><sup>6</sup> <sup>2</sup> and shifted to positive values by *<sup>ε</sup>* such that it is purely repulsive and smooth at *<sup>r</sup>* <sup>=</sup> *<sup>r</sup>*cut <sup>=</sup> <sup>√</sup><sup>6</sup> 2:

Another extension of the potential in Eq. 9 is proposed in (Steinhauser, 2005) where a smooth attractive part is introduced again, in order to allow for including different solvent qualities

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* = 21/6 and at *r* = *r*cut, while it is zero at *r* = *r*cut and has value *γ* at *r* = 21/6, where *γ* is the depth of the attractive part. Further details can be found in (Steinhauser, 2005). When setting *r*cut = 1.5 one sets *γ* = − and obtains *α* and *β* as solutions of the linear set of

<sup>2</sup> · cos(*αr*<sup>2</sup> <sup>+</sup> *<sup>β</sup>*) + *<sup>γ</sup>*

*LJ* (*r*) <sup>−</sup> *λε* <sup>0</sup> <sup>&</sup>lt; *<sup>r</sup>* <sup>&</sup>lt; 21/6 *<sup>σ</sup>*0, *<sup>λ</sup>* <sup>Φ</sup>*cos*(*r*) 21/6 *<sup>σ</sup>*<sup>0</sup> <sup>≤</sup> *<sup>r</sup>* <sup>&</sup>lt; *<sup>r</sup>*cut,

∞ otherwise,

*<sup>r</sup>*<sup>6</sup> <sup>−</sup> *<sup>d</sup>*

<sup>+</sup> *<sup>ε</sup> <sup>r</sup>* <sup>≤</sup> 21/6*σ*0,

0 otherwise.

�

21/3*α* + *β* = *π*, (13) 2.25*α* + *β* = 2*π*. (14)

*<sup>r</sup>*<sup>8</sup> . (10)

*ε*. (12)

(11)

(15)

<sup>Φ</sup>(*r*) = *<sup>b</sup>* exp(−*ar*) <sup>−</sup> *<sup>c</sup>*

*potential* (Buckingham, 1938):

*r* where it has to be modified accordingly.

Φcut *LJ* (*r*) =

of the solvent surrounding the polymer:

equations

⎧ ⎨ ⎩

4 *ε*

Φ*cos*(*r*) =

The total unbounded potential can then be written as:

Φ*Total*(*r*, *λ*) =

⎧ ⎪⎨ Φcut

⎪⎩

�� *σ*<sup>0</sup> *r* �<sup>12</sup> − � *σ*0 *r* �6 �

> � 1

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.

#### **2.2 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 <sup>O</sup>(*N*2) run time to <sup>O</sup>(*N*) run time. Starting from the original LJ potential between two particles *i* and *j* with distance *r* = |*ri* −*rj*| of Eq. 7, one obtains the potential function for *N* interacting particles as the following double sum over all particles:

$$\Phi(\vec{r}\_1, \dots, \vec{r}\_N) = \sum\_{i=1}^N \sum\_{j=i+1}^N \Phi\_{LJ}(r) = 4\varepsilon \sum\_{i=1}^N \sum\_{j=i+1}^N \left(\frac{\sigma\_0}{r}\right)^6 \times \left(\left(\frac{\sigma\_0}{r}\right)^6 - 1\right). \tag{16}$$

The corresponding force *Fi* exerted on particle *i* by particle *j* is given by the gradient with respect to*ri* as:

$$\vec{F}\_{i} = -\nabla\_{\vec{r}\_{i}} \Phi\_{Lf}(\vec{r}\_{1}, \dots, \vec{r}\_{N}) = -24 \times \varepsilon \sum\_{j=1, j \neq i}^{N} \frac{1}{r^{2}} \times \left(\frac{\sigma\_{0}}{r}\right)^{6} \times \left(1 - 2 \times \left(\frac{\sigma\_{0}}{r}\right)^{6}\right) \vec{r}\_{\vec{i}j\prime} \tag{17}$$

where*rij* = (*ri* −*rj*) is the direction vector between particles *i* and *j* at positions*ri* and*rj*, and *<sup>r</sup>* <sup>=</sup> <sup>|</sup>*ri* <sup>−</sup>*rj*|. Hence, in general, the force *Fi* on particle *i* is the sum over all forces *Fij* := −∇*ri* Φ between particle *i* and all other particles *j*:

$$
\vec{F}\_{\vec{i}} = \sum\_{\substack{i=1, j \neq i}}^{N} \vec{F}\_{\vec{i}j}.\tag{18}
$$

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 <sup>1</sup> <sup>2</sup> *N*(*N* − 1) interactions with a *N*<sup>2</sup> efficiency. This algorithm becomes extremely inefficient for systems of more than a few thousand particles, cf. Fig. 3(a).

#### **2.3 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<sup>N</sup>* = (*r*1,*r*2, ...,*rN*), the potential energy Φ(*rN*) and the momenta *p<sup>N</sup>* = (*p*1,*p*2, ...,*pN*), in terms of which the kinetic energy may be written as *K*(*pN*), 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

$$
\vec{r}\_{\text{i}} = \vec{p}\_{\text{i}} / m\_{\text{i}} \quad \text{and} \quad \vec{p}\_{\text{i}} = \vec{F}\_{\text{i}}.\tag{19}
$$

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*ri*, *vi*,*ai* = *Ti*/*mi* are the position, velocity and acceleration of the *i*-th particle, respectively:

$$\text{Calculate } \vec{v}\_i(t + \frac{1}{2}\delta t) = \vec{v}\_i(t) + \frac{1}{2}\delta t \vec{F}\_i(t), \tag{20}$$

$$\text{Calculate}\vec{r}\_l(t+\delta t) = \vec{r}\_l(t) + \delta t \vec{v}\_l(t+\frac{1}{2}\delta t),\tag{21}$$

Derive*ai*(*t* + *δt*) from the interaction potential using*r*(*t* + *δt*), (22)

$$\text{Calculate } \vec{v}\_l(t + \delta t) = \vec{v}\_l(t + \frac{1}{2}\delta t) + +\frac{1}{2}\delta t \vec{a}\_l(t + \delta t). \tag{23}$$

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.

#### **2.3.1 Neighbor lists**

8 Will-be-set-by-IN-TECH

*N* ∑ *j*=1,*j*�=*i*

where*rij* = (*ri* −*rj*) is the direction vector between particles *i* and *j* at positions*ri* and*rj*, and

*N* ∑ *i*=1,*j*�=*i*

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 <sup>1</sup>

1) interactions with a *N*<sup>2</sup> efficiency. This algorithm becomes extremely inefficient for systems

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

When introducing an *N*-dimensional position vector*r<sup>N</sup>* = (*r*1,*r*2, ...,*rN*), the potential energy Φ(*rN*) and the momenta *p<sup>N</sup>* = (*p*1,*p*2, ...,*pN*), in terms of which the kinetic energy may be written as *K*(*pN*), 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

*ri* = *pi*/*mi and pi* =

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

> 2 *δt*

1 2 *<sup>δ</sup>t*)++<sup>1</sup> 2

1 2

Derive*ai*(*t* + *δt*) from the interaction potential using*r*(*t* + *δt*), (22)

 *Fi* =

1 *r*<sup>2</sup> ×

*Fi* exerted on particle *i* by particle *j* is given by the gradient with

*Fi* on particle *i* is the sum over all forces

1 − 2 ×

*Fij*. (18)

*Fi*. (19)

*Ti*/*mi* are the position, velocity and acceleration of the

*Fi*(*t*), (20)

*δt*), (21)

*δtai*(*t* + *δt*). (23)

 *σ*0 *r* 6 

*rij*, (17)

Φ

*Fij* := −∇*ri*

<sup>2</sup> *N*(*N* −

 *σ*0 *r* 6 × 

The corresponding force

*<sup>r</sup>* <sup>=</sup> <sup>|</sup>*ri* <sup>−</sup>*rj*|. Hence, in general, the force

between particle *i* and all other particles *j*:

of more than a few thousand particles, cf. Fig. 3(a).

following four steps, where*ri*, *vi*,*ai* =

Calculate*vi*(*t* +

1 2

Calculate*vi*(*t* + *δt*) = *vi*(*t* +

Calculate*ri*(*t* + *δt*) =*ri*(*t*) + *δtvi*(*t* +

*<sup>δ</sup>t*) = *vi*(*t*) + <sup>1</sup>

*i*-th particle, respectively:

Φ*LJ*(*r*1, ...,*rN*) = −24 × *ε*

respect to*ri* as:

 *Fi* = −∇*ri*

**2.3 The MD algorithm**

state of the field.

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 � *Fi* on particle *i* is approximated by

$$\vec{F}\_{\rm i} \approx -24 \times \varepsilon \sum\_{\substack{j=1, j \neq i \\ 0 < r \le \varepsilon\_{\rm cut}}}^{N} \frac{1}{r^2} \times \left(\frac{\sigma\_0}{r}\right)^6 \times \left(1 - 2 \times \left(\frac{\sigma\_0}{r}\right)^6\right) \vec{r}\_{\rm ij}. \tag{24}$$

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 <sup>O</sup>(*N*2) to <sup>O</sup>(*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 � *Fi* on particle *i* in cell number *n* one obtains a sum of the form

$$\vec{F}\_l = \sum\_{\substack{cell \ m \\ m \in \Omega(n)}} \sum\_{\substack{j \in \{\text{all particles in cell } m\} \\ j \neq i}} \vec{F}\_{ij\nu} \tag{25}$$

where Ω(*n*) denotes cell *n* itself together with all cells that are direct neighbors of cell *n*. The 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*<sup>3</sup> cut. The total number of subcells is *<sup>N</sup>*/ �*ρ*�*r*<sup>3</sup> cut and

Fig. 3. MD Optimization schemes for the search of potentially interacting particles. (a) The least efficient all particle "brute force" approach with run time <sup>O</sup>(*N*2) (b) The linked-cell algorithm which reduces the search effort to O(*N*). (c) 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 + *rskin* 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:

$$\frac{26}{2} \left( \left< \rho \right> r\_{\rm cut}^2 \right) \frac{N}{\left< \rho \right> r\_{\rm cut}^3} = 13 (\left< \rho \right> r\_{\rm cut}^3)^2 N. \tag{26}$$

For this method to function, the size of the simulation box has to be at least 3*r*cut, cf. 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.

#### **2.3.2 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. 3. MD Optimization schemes for the search of potentially interacting particles. (a) The least efficient all particle "brute force" approach with run time <sup>O</sup>(*N*2) (b) The linked-cell algorithm which reduces the search effort to O(*N*). (c) 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 + *rskin* and the inner

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 3*r*cut, cf. 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

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

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

<sup>=</sup> <sup>13</sup>(�*ρ*�*r*<sup>3</sup>

cut)2*N*. (26)

 *N* �*ρ*�*r*<sup>3</sup> cut

circle, which contains actually interacting particles, has radius *r*cut.

26 2 �*ρ*�*r* 2 cut

*r*cut.

**2.3.2 Boundary conditions**

two dimensional visualization.

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.

#### **2.3.3 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/2*LB* with boxlenth *LB*. This 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.

### **3. 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).

#### **3.1 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

$$\mathfrak{F} = \mathfrak{F}\_{\mathsf{B}} k\_{\mathsf{B}} \mathsf{T} / \mathfrak{e}\sigma \tag{27}$$

can be introduced, where the Bjerrum length *ξ<sup>B</sup>* is given by:

$$
\mathfrak{T}\_B = e^2 / \epsilon k\_B T \, , \tag{28}
$$

where *kB* 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.

The interaction parameter for the here presented study is chosen in the range of *ξ* = 0, ..., 100 which covers electrically neutral chains (*ξ* = 0) in good solvent as well as highly charged chain systems (*ξ* = 100). The monomers in the chains are connected by harmonic bonds

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

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

*ξ<sup>B</sup>* = *e*

where *kB* 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*

The interaction parameter for the here presented study is chosen in the range of *ξ* = 0, ..., 100 which covers electrically neutral chains (*ξ* = 0) in good solvent as well as highly charged chain systems (*ξ* = 100). The monomers in the chains are connected by harmonic bonds

*ξ* = *ξBkBT*/*�σ* (27)

2/*�kBT*, (28)

can be found in (Holm et al., 2004).

**3.1 Two oppositely charged macromolecules**

units) the dimensionless interaction parameter

is the electronic charge.

can be introduced, where the Bjerrum length *ξ<sup>B</sup>* is given by:

Fig. 5. Twisted, DNA-like polyelectrolyte complexes formed by electrostatic attraction of two oppositely charged linear macromolecules with *N* = 40 at different time intervals *τ* = 0 (a), *τ* = 10500 (b), *τ* = 60000 (c) and *τ* = 120000 (d), where *τ* is given in reduced Lennard-Jones units (Allen & Tildesly, 1991). The interaction strength is *ξ* = 8 (Steinauser, 1998; Winkler et al., 2002).

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 *χ*, acting on each mass point. The equations of motion of the system are thus given by the Langevin equations

$$
\hbar m \ddot{\vec{r}}\_i = \vec{F}\_i - \chi m \dot{\vec{r}}\_i + \vec{\Gamma}\_i. \tag{29}
$$

The force *Fi* comprises the force due to the sum of the potentials of Eq. 11 with cutoff*r*cut = 1.5, Eq. 6 with *k* = 1, *zi*/*<sup>j</sup>* = ±1, and the first term on the right-hand side of the bonded potential in Eq. 2 with *κ* = 5000*ε*/*σ* and bond length *l*<sup>0</sup> = *σ*<sup>0</sup> = 1.0. The stochastic force Γ*<sup>i</sup>* 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 *<sup>ρ</sup>* <sup>=</sup> *<sup>N</sup>*/*<sup>V</sup>* <sup>=</sup> 2.1 <sup>×</sup> <sup>10</sup>−7/*σ*<sup>3</sup> constant when changing the chain length *<sup>N</sup>*. 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/2*LB* 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 *zi*/*<sup>j</sup>* = ±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).

Figure 6(a) exhibits the universal scaling regime of *Rg* obtained for intermediate interaction strengths *<sup>ξ</sup>* and scaled by (*<sup>N</sup>* <sup>−</sup> <sup>1</sup>)2/3. Here, the data of various chain lengths fall nicely on top of each other. This scaling corresponds to the scaling behavior of flexible chains in a bad solvent and is also in accordance with what was reported by Shrivastava and Muthukumar (Srivastava & Muthukumar, 1994). The change of *Rg* 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

Fig. 6. (a) Radii of gyration as a function of the interaction strength *ξ* for various chain lengths according to (Steinhauser, 2008; Winkler et al., 2002). The radius of gyration *Rg* 2 which is measure for the size of a polymer chain is scaled by (*<sup>N</sup>* <sup>−</sup> <sup>1</sup>)2/3, where (*<sup>N</sup>* <sup>−</sup> <sup>1</sup>) is the number of bonds of a single chain. Also displayed are sample snapshots of the collapsed globules with *N* = 40 and interaction strengths *ξ* = 0.4, 4, 40. (b) Radial density of monomers with respect to the center of mass of a globule and interaction strength *ξ* = 4 for different chain lengths, *N* = 20 (blue), *N* = 40 (red), *N* = 80 (green) and *N* = 160 (brown).

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 *ξ*, chains of different lengths are deformed to different degrees which leads to a chain length dependence of the density profile.

#### **4. 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.

#### **4.1 Theoretical framework**

14 Will-be-set-by-IN-TECH

Fig. 6. (a) Radii of gyration as a function of the interaction strength *ξ* for various chain lengths according to (Steinhauser, 2008; Winkler et al., 2002). The radius of gyration *Rg*

chain lengths, *N* = 20 (blue), *N* = 40 (red), *N* = 80 (green) and *N* = 160 (brown).

**4. Equilibrium and Non-Equilibrium Molecular Dynamics (NEMD)**

rod-like polymers is thus of great practical and fundamental interest.

of the density profile.

which is measure for the size of a polymer chain is scaled by (*<sup>N</sup>* <sup>−</sup> <sup>1</sup>)2/3, where (*<sup>N</sup>* <sup>−</sup> <sup>1</sup>) is the number of bonds of a single chain. Also displayed are sample snapshots of the collapsed globules with *N* = 40 and interaction strengths *ξ* = 0.4, 4, 40. (b) Radial density of monomers with respect to the center of mass of a globule and interaction strength *ξ* = 4 for different

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 *ξ*, chains of different lengths are deformed to different degrees which leads to a chain length dependence

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

2

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 *Lp*. The corresponding elastic energy

$$\Phi\_{\text{bend}} = \frac{\kappa}{2} \int\_0^L ds \left( \frac{\partial^2 \mathbf{r}}{\partial s^2} \right)^2 \tag{30}$$

of the inextensible chain of length *L* depends on the local curvature of the chain contour *s*, where *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

$$
\vec{F}\_{\text{bend}} = \mu \frac{\partial^4 \vec{r}}{\partial s^4} \tag{31}
$$

$$
\vec{F}\_{\text{lens}} = \omega \frac{\partial^2 \vec{r}}{\partial s^2}.\tag{32}
$$

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 *Lp 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 *τ<sup>p</sup>* or the normal mode amplitudes *Xp*(*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*-C100H202 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 *lp* (or Kuhn segment lengths *lK*) comparable to or larger than their contour length *L* and their rigidity and relaxation behavior are essential for their biological functions.

#### **4.2 Modeling and simulation of semiflexible macromolecules**

Molecular Dynamics simulations were performed using the MD simulation package "MD-Cube", which was originally developed by Steinhauser (Steinhauser, 2005). 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; Steinhauser & et. al., 2005). 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

$$\Phi\_{\text{bonded}}(r) = \frac{K}{2}(r - d\_0)^2 \,\text{,}\tag{33}$$

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*ε*/*σ*<sup>2</sup> is chosen, where *ε* and *σ* are parameters of the truncated LJ-potential in Eq. 11. The average bondlength *d*<sup>0</sup> 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) potential – which is frequently employed in polymer simulations (Steinhauser, 2005) – 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

$$\Phi\_{\text{bend}} = \frac{\kappa}{2} \sum\_{i=1}^{N-1} \left( \vec{u}\_{i+1} - \vec{u}\_i \right)^2 = \kappa \sum\_{i=1}^{N-1} \left( 1 - \vec{u}\_{i+1} \cdot \vec{u}\_i \right) \tag{34}$$

where �*u* is the unit bond vector �*ui* = (�*ri*+<sup>1</sup> −�*ri*)/|�*ri*+<sup>1</sup> −�*ri* | connecting consecutive monomers, and �*ri* is the position vector to the *i*-th monomer. The total force acting on monomer *i* is thus given by

$$\vec{F}\_{\text{i}} = -\partial / \partial \vec{r}\_{i} \Phi\_{\text{total}} = -\partial / \partial \vec{r}\_{i} \left(\Phi\_{\text{LJ}}^{\text{cut}} + \Phi\_{\text{bonded}} + \Phi\_{\text{bond}}\right) \tag{35}$$

Assuming a NVT ensemble at temperature *kBT* = 1*ε*, the trajectories of all particles *i* = (1, ..., *N*) are generated by integrating the stochastic equations of motion

$$\frac{d^2\vec{r}\_i}{dt^2} = -\xi \frac{d\vec{r}\_i}{dt} + \frac{1}{m}\vec{F}\_i + \frac{1}{m}\vec{F}\_i^S \,'\,. \tag{36}$$

with a particle friction coefficient *ξ* and a Gaussian stochastic force � *Fi <sup>S</sup>* that satisfies

< � *Fi <sup>S</sup>*(*t*) >= 0 , (37)

and the fluctuation-dissipation-theorem

$$<\langle F\_{l}^{S}(t)F\_{l}^{S}(t')> = 2kT\xi\delta\_{lj}\delta(t-t')\ . \tag{38}$$

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 <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>−3*τ*, where *τ* is the time unit of the simulation.

#### **4.2.1 Results**

16 Will-be-set-by-IN-TECH

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 *lp* (or Kuhn segment lengths *lK*) comparable to or larger than their contour length *L* and their rigidity and relaxation behavior

Molecular Dynamics simulations were performed using the MD simulation package "MD-Cube", which was originally developed by Steinhauser (Steinhauser, 2005). 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; Steinhauser & et. al., 2005). 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

Neighboring mass points along the chains are connected by harmonic bonds with the

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*ε*/*σ*<sup>2</sup> is chosen, where *ε* and *σ* are parameters of the truncated LJ-potential in Eq. 11. The average bondlength *d*<sup>0</sup> 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) potential – which is frequently employed in polymer simulations (Steinhauser, 2005) – 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

The stiffness, i.e. the bending rigidity of the chains composed of *N* mass points, is introduced

*N*−1 ∑ *i*=1

(*ui*+<sup>1</sup> <sup>−</sup>*ui*)<sup>2</sup> <sup>=</sup> *<sup>κ</sup>*

<sup>2</sup> (*<sup>r</sup>* <sup>−</sup> *<sup>d</sup>*0)

<sup>2</sup> , (33)

(1 −*ui*+<sup>1</sup> · *ui*) , (34)

<sup>Φ</sup>bonded(*r*) = *<sup>K</sup>*

are essential for their biological functions.

following potential for the bonded interactions

calculated in a numerical approach is used.

into the coarse-grained model by the following bending potential

*N*−1 ∑ *i*=1

2

<sup>Φ</sup>bend <sup>=</sup> *<sup>κ</sup>*

**4.2 Modeling and simulation of semiflexible macromolecules**

for each monomer is taken into account through the potential of Eq. 11.

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<sup>θ</sup>* (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<sup>θ</sup>* of the chains. The scaling plot in Fig. 7(b) shows that the transition occurs for *q* ≈ 1/*lK*, i.e. on a 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*

Fig. 7. (a) Structure function *S*(*q*) of single linear chains with *N* = 700 and varying stiffness *k<sup>θ</sup>* . The scaling regimes (fully flexible and stiff rod) are indicated by a straight and dashed line, respectively. (b) Scaling plot of *S*(*q*)/*lK* versus *q* · *lK* using the statistical segment length *lK*, adapted from (Steinhauser, Schneider & Blumen, 2009).

Fig. 8. Simulation snapshots of (a) flexible chains (*k<sup>θ</sup>* = 0), (b) semiflexible chains (*k<sup>θ</sup>* = 20), (c) stiff, rod-like chains (*k<sup>θ</sup>* = 50).

in normal modes *Xp*(*t*) as follows:

$$\vec{r}(s,t) = \vec{X}\_0(t) + 2\sum\_{p=1}^{\infty}(t)\cos\left(\frac{p\pi}{L}s\right). \tag{39}$$

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 *τ<sup>p</sup>*

$$\tau\_p = \frac{3k\_B T \pi^2 p^2}{2L^2 \tilde{\zeta}} \left[ \frac{1}{L\_p} + \frac{L\_p \pi^2}{NL} p^2 \right],\tag{40}$$

which can be interpreted physically as arising from contributions due to an entropic force term ∝ *p*<sup>2</sup> 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.

Fig. 7. (a) Structure function *S*(*q*) of single linear chains with *N* = 700 and varying stiffness *k<sup>θ</sup>* . The scaling regimes (fully flexible and stiff rod) are indicated by a straight and dashed line, respectively. (b) Scaling plot of *S*(*q*)/*lK* versus *q* · *lK* using the statistical segment length

Fig. 8. Simulation snapshots of (a) flexible chains (*k<sup>θ</sup>* = 0), (b) semiflexible chains (*k<sup>θ</sup>* = 20),

∞ ∑ *p*=1

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 *τ<sup>p</sup>*

> 1 *Lp* + *Lpπ*<sup>2</sup> *NL <sup>p</sup>*<sup>2</sup>

which can be interpreted physically as arising from contributions due to an entropic force term ∝ *p*<sup>2</sup> and a bending force term ∝ *p*4. In Fig. 9 we show that our simulation results for

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

Finally, in Fig. 11 we illustrate our step-shear simulation scheme for polymer melts with two

(*t*) cos

 *pπ L s* 

. (39)

, (40)

*r*(*s*, *t*) = *X* <sup>0</sup>(*t*) + 2

*<sup>τ</sup><sup>p</sup>* <sup>=</sup> <sup>3</sup>*kBTπ*<sup>2</sup> *<sup>p</sup>*<sup>2</sup> 2*L*2*ξ*

*lK*, adapted from (Steinhauser, Schneider & Blumen, 2009).

*Xp*(*t*) as follows:

semiflexible chains scale according to Eq. 40.

corresponding NEMD simulation snapshots.

(c) stiff, rod-like chains (*k<sup>θ</sup>* = 50).

in normal modes

equation.

Fig. 9. Scaling of relaxation time *τ<sup>p</sup>* for semiflexible chains (*N* = 700) with different *Lp*. The dotted and dashed lines show the scaling behavior according to Eq. 40. The inset shows *τp* of our simulated flexible chains compared with the Rouse model.

Fig. 10. Shear modulus *G*(*t*) obtained from NEMD simulations compared with conventional (red line) equilibrium methods. The dashed line displays the expected Rouse scaling behavior and *γ*<sup>0</sup> displays the respective shear rates of the systems.

#### **5. 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).

Fig. 11. Non-equilibrium shear-simulation of *N* = 100 polymer chains with *N* = 100 particles each (a) before and (b) after shearing the system. For reasons of clarity only 30 different chains of the system are displayed.

#### **5.1 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*<sup>0</sup> and the initial (randomly generated) degree of overlap indicated by *d*<sup>0</sup> *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.

#### **5.2 Model potentials**

20 Will-be-set-by-IN-TECH

Fig. 11. Non-equilibrium shear-simulation of *N* = 100 polymer chains with *N* = 100 particles each (a) before and (b) after shearing the system. For reasons of clarity only 30 different

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

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*<sup>0</sup> and the initial (randomly generated) degree of

(6), and green (4). (c) The same system displayed as an unordered network.

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

*ij*. Here, only two particles are displayed. In the model the number of

**5.1 Introduction: Multiscale modeling of granular materials using MD**

links that connect the centers of overlapping particles, cf. Fig. 12(c).

chains of the system are displayed.

overlap indicated by *d*<sup>0</sup>

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 *dij* decreases compared to the initial overlap *d*0 *ij*:

$$\phi\_{rep}^{ij}\left(\gamma, d^{ij}\right) = \begin{cases} \gamma \mathcal{R}\_0^{-3} \left( \left(\frac{d\_0^{\vee}}{d^{\vee}}\right)^{12} - 2\left(\frac{d\_0^{\vee}}{d^{\vee}}\right)^6 + 1 \right) & : \quad 0 < d^{ij} < d\_0^{\cdot ij} \\ 0 & : \quad d^{ij} \ge d\_0^{\cdot ij} \end{cases} . \tag{41}$$

Parameter *γ* scales the energy density of the potential and prefactor *R*<sup>0</sup> <sup>3</sup> ensures the correct scaling behavior of the calculated total stress Σ*ijσij* = Σ*ijFij*/*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*<sup>0</sup> is the constant radius of the particles, *dij* = *dij* (*t*) is the instantaneous mutual distance of each interacting pairs {*ij*} of particles, and *d*<sup>0</sup> *ij* = *dij* (*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*<sup>0</sup> *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.

Fig. 13. (a) Schematic of the intrinsic scaling property of the proposed material model. Here, only the 2D case is shown for simplicity. The original system (Model *Ma*) with edge length *L*<sup>0</sup> and particle radii *R*<sup>0</sup> is downscaled by a factor of 1/*a* into the subsystem *QA* of *MA* (shaded area) with edge length *L*, while the particle radii are upscaled by factor *a*. As a result, model *MB* of size *aL* = *L*<sup>0</sup> 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*<sup>0</sup> *ij*. This property is expressed in the cohesive potential by the following equation:

$$
\Phi\_{coh}^{\dot{ij}}\left(\lambda, d^{\dot{ij}}\right) = \lambda R\_0 \left(d^{\dot{ij}} - d\_0^{\dot{ij}}\right)^2, \ d^{\dot{ij}} > 0. \tag{42}
$$

In this equation, *λ* (which has dimension [energy/length]) determines the strength of the potential and prefactor *R*<sup>0</sup> again ensures a proper intrinsic scaling behavior of the material response. The total potential is the following sum:

$$
\Phi\_{\rm tot} = \Sigma\_{\rm ij} \left( \phi\_{\rm rep}^{\rm ij} + \phi\_{\rm coh}^{\rm ij} \right). \tag{43}
$$

The repulsive part of Φ*tot* acts only on particle pairs that are closer together than their mutual initial distance *d*<sup>0</sup> *ij*, whereas the harmonic potential <sup>Φ</sup>*coh* either acts repulsively or cohesively, depending on the actual distance *dij*. 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:

$$d^{\dagger\dagger} > (2\mathcal{R}\_0). \tag{44}$$

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

$$d^{ij} < \alpha \cdot d\_0^{ij}.\tag{45}$$

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.

#### **5.3 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 (Al2O3 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

In this equation, *λ* (which has dimension [energy/length]) determines the strength of the potential and prefactor *R*<sup>0</sup> again ensures a proper intrinsic scaling behavior of the material

> *φij rep* <sup>+</sup> *<sup>φ</sup>ij coh*

The repulsive part of Φ*tot* acts only on particle pairs that are closer together than their mutual

depending on the actual distance *dij*. 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

Failure under pressure load occurs when the actual mutual particle distance is less by a factor

*<sup>d</sup>ij* <sup>&</sup>lt; *<sup>α</sup>* · *<sup>d</sup>*<sup>0</sup>

Particle pairs without a spring after pressure or tensile failure still interact via the repulsive

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

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 (Al2O3 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

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

*ij*, whereas the harmonic potential <sup>Φ</sup>*coh* either acts repulsively or cohesively,

. (43)

*dij* > (2*R*0). (44)

*ij*. (45)

Φ*tot* = Σ*ij*

response. The total potential is the following sum:

overlap between two particles vanishes, i.e. when:

potential and cannot move through each other.

can be adjusted to mimic the behavior of specific materials.

**5.3 Shock wave simulations and comparison with experiments**

left edge. This leads to a local compression of the particles in the impact area.

*α* (with *α* ∈ (0, 1)) than the initial mutual particle distance, i.e. when

initial distance *d*<sup>0</sup>

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 (Al2O3 and *SiC*) are recovered in the simulation of uniaxial tensile load. The insets show high-speed photographs of *SiC* and Al2O3, 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 Al2O3 and *SiC*, after 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 Al2O3 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.

### **6. 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-macromolecules6. 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.

#### **7. References**

Allen, M. & Tildesly, D. (1991). *Computer Simulation of Liquids*, Oxford University Press.

Aragón, S. & Pecora, R. (1985). Dynamics of wormlike chains, *Macromolecules* 18: 1868–1875.

Barrat, J.-L. & Joanny, F. (2007). Theory of Polyelectrolyte Solutions, *Adv. Chem. Phys.* 94: 1–66.


<sup>6</sup> Semiflexibility is a characteristic feature of bio-macromolecules such as DNA, F-actin, intermediate filaments or microtubuli which are present in the cytoplasm of cells.

bio-macromolecules6. 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.

Allen, M. & Tildesly, D. (1991). *Computer Simulation of Liquids*, Oxford University Press. Aragón, S. & Pecora, R. (1985). Dynamics of wormlike chains, *Macromolecules* 18: 1868–1875. Barrat, J.-L. & Joanny, F. (2007). Theory of Polyelectrolyte Solutions, *Adv. Chem. Phys.* 94: 1–66. Buckingham, R. (1938). The Classical Equation of State of Gaseous Helium, Neon and Argon,

Bustamante, C., Marko, J., Siggia, E. & Smith, S. (1994). Entropic Elasticity of Lambda-Phage

Cundall, P. & Strack, O. (1979). A Discrete Numerical Model for Granular Assemblies,

de Gennes, P. G. (1979). *Scaling Concepts in Polymer Physics*, Cornell University Press, Ithaca,

Debye, P. & Hückel, E. (1923). The Theory of Electrolytes. I. Lowering of Freezing Point and

Gössl, I., Shu, L., Schlüter, A. & Rabe, J. (2002). Complexes of DNA and positively charged

Haberland, R., Fritsche, S., Peinel, G. & Heinzinger, K. (1995). *Molekulardynamik -*

Harnau, L., Winkler, R. & Reineker, P. (1996). Dynamic structure factor of semiflexible macromolecules in dilute solution, *J. Chem. Phys.* 104(16): 6255–6258.

<sup>6</sup> Semiflexibility is a characteristic feature of bio-macromolecules such as DNA, F-actin, intermediate

Haile, J. (1992). *Molecular dynamics simulation: Elementary Methods*, Wiley, New York.

*Grundlagen und Anwendungen*, Friedrich VIeweg & Sohn Verlagsgesellschaft mbH,

Doi, M. & Edwards, S. (1986). *The Theory of Polymer Dynamics*, Clarendon Press, Oxford. Feller, S. (2000). An Improved Empirical Potential Energy Function for Molecular Simulations

Related Phenomena, *Physikalische Zeitschrift* 24: 185–206.

of Phospholipids, *J. Phys. Chem. B* 104: 7510–7515. Flory, P. (1969). *Statistical Mechanics of Chain Molecules*, Wiley, New York.

dendronized polymers, *J. Am. Chem. Soc.* 124: 6860.

filaments or microtubuli which are present in the cytoplasm of cells.

**7. References**

*Proc. Roy. Soc* A106: 264.

*Geotechnique* 29: 47–65.

Braunschweig, Wiesbaden.

London.

DNA, *Science* 265: 1599–1600.


## **Molecular Dynamics Simulation of Synthetic Polymers**

Claudia Sandoval

*Center for Bioinformatics and Integrative Biology (CBIB), Facultad de Ciencias Biológicas Universidad Andres Bello, Santiago Fraunhofer Chile Research Foundation - Center for Systems Biotechnology, FCR-CSB Chile* 

### **1. Introduction**

26 Will-be-set-by-IN-TECH

28 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

Steinhauser, M. O. (2008). *Computational Multiscale Modeling of Solids and Fluids – Theory and*

Steinhauser, M. O. (2012). *Computer Simulation in Physics and Engineering*, de Gruyter. To be

Steinhauser, M. O. & et. al. (2005). MAVO MMM-Tools: Entwicklung von

Steinhauser, M. O., Grass, K., Strassburger, E. & Blumen, A. (2009). Impact Failure of Granular

Steinhauser, M. O., Grass, K., Thoma, K. & Blumen, A. (2006). A Nonequilibrium Molecular

Steinhauser, M. O., Schneider, J. & Blumen, A. (2009). Simulating Dynamic Crossover

van Gunsteren, W. & Berendsen, H. (1982). Algorithms for Brownian Dynamics, *Molec. Phys.*

Walton, O. R. & Braun, R. L. (1986). Viscosity, Granular-Temperature, and Stress Calculation for Shearing Assemblies of Inelastic, Frictional Disks, *J. Rheol.* 39(5): 949. Winkler, R. G., Steinhauser, M. O. & Reineker, P. (2002). Complex Formation in Systems of

Yamasaki, Y., Teramoto, Y. & Yoshikawa, K. (2001). Disappearance Of The Negative Charge in

Zimm, B. H. (1956). Dynamics of Polymer Molecules in Dilute Solution: Viscoelasticity, Flow

Giant DNA With A Folding transition, *Biophys. J.* 80: 2823U2832.

Birefringence and Dielectric Loss, *J. Chem. Phys.* 24: 269–278.

durchgängigen Multiskalen Material Modellierungen, *Technical Report 17/05*, Fraunhofer Ernst-Mach-Institute (EMI), Freiburg, Eckerstrasse 4, 79104 Freiburg i.

Materials – Non-Equilibrium Multiscale Simulations and High-Speed Experiments,

Behavior of Semiflexible Linear Polymers in Solution and in the Melt, *J. Chem. Phys.*

Oppositely Charged Polyelectrolytes: A Molecular Dyamics Simulation Study, *Phys.*

˝

*Applications*, Springer, Berlin, Heidelberg, New York.

published in October 2012, ISBN 978-3-11-025590-4.

*International Journal of Plasticity* 25: 161–182.

Dynamics Study on Shock Waves, *Europhys. Lett.* 73: 62.

Br., Germany.

130: 164902.

*Rev. E* 66: 021802.

45: 637.

The polymeric systems are great interesting both academic and industrial level, due to interesting applications in several areas. Nowadays, the trends for the development of new materials are focused on the study of macromolecular complex systems. A macromolecular complex system is formed by the interaction between the macromolecules of the same system. There are many macromolecular complex systems; this chapter will refer to complex systems of synthetic polymer blends(Donald R. P., 2000) and synthetic polymer at the airwater interface(Gargallo et al., 2011). The most synthetic polymers are characterized by flexible structure able to adopt multiple and variable forms. The study about molecular behavior in polymers has attracted the attention of many researchers over the years. Studies in dilute solutions of macromolecules have contributed to solving the problem relating to the molecular characterization. Despite, the emergence of characterization methods and theories that have allowed for interpretation of experimental behavior, researchers still have not fully resolved the problem with interactions at the atomic level in polymer systems. The most polymers are amorphous so in some cases are impossible to obtain information about atomic structure through x-ray techniques or atomic force microscopy. At this point, the computer simulations have played an important role in microscopic understanding of dynamical properties of some polymeric systems.

In this chapter a review of molecular simulation methodologies focused on polymer systems will be discussed. This review is based on two case studies. The first one corresponds to a monolayer of amphiphilic polymers and their study at the air-water interface(Gargallo et al., 2009). Three polymers were selected Poly(ethylene oxide) (PEO), Poly(tetrahydrofuran) (PTHF) and Poly(-caprolactone) (PEC) and their cyclodextrin inclusion complexes. The second system corresponds to study the motion of side chains of polymers in condensed phase. For this reason, a diblock copolymers of Polystyrene-block-poly(t-butylacrylate) has been considered(Encinar et al., 2008). The methodology to build the models in both systems will be provided, as well as, the methodology to obtain the atomic parameters, necessary to molecular simulations.

To carry out molecular dynamics studies of polymer systems is necessary to know some important physicochemical properties or thermodynamic parameters of polymers. For the other hand, it is important to know the different methodologies employed to carry out molecular dynamic simulation of polymers.

### **2. Physicochemical properties of polymers**

In this chapter we concéntrate efforts on giving the reader a brief description of some important physicochemical properties of the polymers for molecular simulation analysis analysis. Macromolecular conformation is an important point to consider. The macromolecular chains in dilute solution have variables internal rotational angles, Figure 1. The internal rotational angle can take values that fluctuate within a range so that different minimum energy can be obtained for various values of the angles . The energy difference associated with the variation of this angle is minimal. This leads to a large number of stable conformations and significant mobility to oscillate around these conformations. Both factors cause that the macromolecule has a significant flexibility. Experimentally the dimensional parameters that characterize a macromolecule can only be characterized in terms of statistical parameters that represents an average over all conformations. The *end-to-end distance* and *radius of gyration* are statistical parameters that characterize the spatial distribution of a chain macromolecule(Flory, 1969, 1995).

Fig. 1. Three successive bonds for a chain. Where and represents the bond angle and internal rotation angle, respectively.

The *end-to-end distance* ൏ ݎ<sup>ଶ</sup> represents the distance between end chains. This parameter can vary from a maximum value, where the polymer chain is fully extended to a minimum value corresponding to the sum of the van der Waals radii of the terminal groups, Figure 2.

The end-to-end distance is a statistical parameter that characterizes the spatial distribution of the chains relating the microscopic with macroscopic properties, also allows to relate the hydrodynamic parameters with thermodynamic parameter, and finally allows to obtain information respect to flexibility chain.

other hand, it is important to know the different methodologies employed to carry out

In this chapter we concéntrate efforts on giving the reader a brief description of some important physicochemical properties of the polymers for molecular simulation analysis analysis. Macromolecular conformation is an important point to consider. The macromolecular chains in dilute solution have variables internal rotational angles, Figure 1. The internal rotational angle can take values that fluctuate within a range so that different minimum energy can be obtained for various values of the angles . The energy difference associated with the variation of this angle is minimal. This leads to a large number of stable conformations and significant mobility to oscillate around these conformations. Both factors cause that the macromolecule has a significant flexibility. Experimentally the dimensional parameters that characterize a macromolecule can only be characterized in terms of statistical parameters that represents an average over all conformations. The *end-to-end distance* and *radius of gyration* are statistical parameters that

characterize the spatial distribution of a chain macromolecule(Flory, 1969, 1995).

Fig. 1. Three successive bonds for a chain. Where and represents the bond angle and

The *end-to-end distance* ൏ ݎ<sup>ଶ</sup> represents the distance between end chains. This parameter can vary from a maximum value, where the polymer chain is fully extended to a minimum value corresponding to the sum of the van der Waals radii of the terminal groups, Figure 2. The end-to-end distance is a statistical parameter that characterizes the spatial distribution of the chains relating the microscopic with macroscopic properties, also allows to relate the hydrodynamic parameters with thermodynamic parameter, and finally allows to obtain

molecular dynamic simulation of polymers.

internal rotation angle, respectively.

information respect to flexibility chain.

**2. Physicochemical properties of polymers** 

Fig. 2. End-to-end distance � �� � representation to a) fully extended chain model, where � �� � takes a maximum value and b) � �� � takes a minimum value and correspond to the sum of the van der Waals radii.

� �� � is defined by the equation 1:

$$ = \int\_0^\infty r^2 W(r) dr\tag{1}$$

where W(r) is probability distribution function, defined by

$$\int\_{0}^{\infty} W(r) dr = 1$$

But the chains may adopt different conformational states depending on the characteristics of the chain, features such as substituents and the solvent. For the different conformational states, there are different chains models. Is not intended to provide a thorough review of statistical mechanics of polymer chains(Flory, 1969), but is important to note that the value of � �� � depends on two main factors, "short-range" and "long-range" factors. For these two factors, different chain models were raised in the past. Short-range factors involved characteristics of the bonds and interaction between neighboring atoms or groups along a chain, this chain correspond to a rotation free model. Hence, � �� � can be calculated knowing the values of bond angles and bond distances for all the conformations, the value obtained is called dimension mean square of end-to-end distance 〈��〉��. Take into account factors such as bond angle restrictions and steric hindrance between substituent groups, we can obtain unperturbed dimension 〈��〉�, and may be expressed by equation 2:

$$
\langle r^2 \rangle\_0 = \sigma^2 \langle r^2 \rangle\_{of} \tag{2}
$$

where is the stiffness coefficient or conformational factor. The 〈��〉� is so-called "theta condition".

In long-range factors are involved the interactions between distant segments along the chain, as well as interaction with the solvent. Taking into account that, two distant segments of a chain, cannot occupy the same place at the same time, arise the excluded volume concept. Under effect of these interactions the end-to-end distance to linear chains in dilute solution, may be described by:

$$^{1/2} = ^{1/2}\_0 a \tag{3}$$

where linear expansion chain factor. In equation 3, the root-mean-square end-to-end distance � �� �� �� depends of number of segments, temperature and solvent. When >1 them the interaction between solvent and chains are favorable, this means that excluded volume effect is important. Instead, if the interactions between segments and solvent are unfavorable them tend to decrease; this causes folding of the chains to minimize the contact with the solvent. These effects can be balance and them the polymer chain adopts an unperturbed state.

The *radius of gyration* is the parameter more appropriate to describe the dimensions of a macromolecule. The mean-square radius of gyration � �� � � represents the mean square of the distance between each atom in the chain and center of mass of the chain and may be defined to following equation 4:

$$1 < R\_g^2 > = (1/n) \,\Sigma\_1^n < r\_l^2 > \tag{4}$$

where n is the segments numbers of the chain and � �� � � represents the mean-square distance of the *i*-th segment from the center of mass of the chains. Through the random flight model we can express the radius of gyration in terms to mean-square end-to-end distance, indicating their interchangeability. For a freely rotating chain of *n* segments of length *l*, the following relationship is established:

$$ = \frac{1}{6}  = \frac{nl^2}{6} \tag{5}$$

A complete derivation of these equations can be reviewed in various texts of statistical mechanics of polymers(Flory, 1995) (Yamakawa, 2001). We have said that the dimensional parameters are important to know because of their contribution to the characterization of polymers. In addition, these parameters can be calculated both experimentally and theoretically. For example, radius of gyration has been calculated, in several works, from molecular simulation trajectories and compared with the experimental values(Lee et al., 2006; Prosa et al., 2001, 1997).

#### **3. Molecular dynamic simulation of supramolecular structures at the air-water interface**

A great interesting in the last year toward to study theoretical of self-assembly of macromolecular systems have been report(Duncan et al., 2011; Korchowiec et al., 2011; Saldias et al., 2009). The self-assembly of macromolecular systems are great interesting due to applications in several areas, as medicine, drug delivery and smart materials(Wick and Cummings, 2011). The polymers at the air-water interface are considering self-assembly systems. The studies of polymers at the air-water interface considering several molecular simulation methodologies, such as atomistic(Gargallo et al., 2009) and coarsegrained(Duncan et al., 2011) models. In this chapter we will focus on atomistic molecular simulation methodology. The aim of studies of polymer at the air-water interface, using molecular simulation, is characterized the orientation of polymer chains at the interface, as well as, calculate the different interaction between amphiphilic polymer and the water subphase.

### **3.1 Build of the molecular models**

32 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

where linear expansion chain factor. In equation 3, the root-mean-square end-to-end distance � �� �� �� depends of number of segments, temperature and solvent. When >1 them the interaction between solvent and chains are favorable, this means that excluded volume effect is important. Instead, if the interactions between segments and solvent are unfavorable them tend to decrease; this causes folding of the chains to minimize the contact with the solvent. These effects can be balance and them the polymer chain adopts an

The *radius of gyration* is the parameter more appropriate to describe the dimensions of a

the distance between each atom in the chain and center of mass of the chain and may be

� �� �� �� ∑ � ��

distance of the *i*-th segment from the center of mass of the chains. Through the random flight model we can express the radius of gyration in terms to mean-square end-to-end distance, indicating their interchangeability. For a freely rotating chain of *n* segments of

A complete derivation of these equations can be reviewed in various texts of statistical mechanics of polymers(Flory, 1995) (Yamakawa, 2001). We have said that the dimensional parameters are important to know because of their contribution to the characterization of polymers. In addition, these parameters can be calculated both experimentally and theoretically. For example, radius of gyration has been calculated, in several works, from molecular simulation trajectories and compared with the experimental values(Lee et al.,

A great interesting in the last year toward to study theoretical of self-assembly of macromolecular systems have been report(Duncan et al., 2011; Korchowiec et al., 2011; Saldias et al., 2009). The self-assembly of macromolecular systems are great interesting due to applications in several areas, as medicine, drug delivery and smart materials(Wick and Cummings, 2011). The polymers at the air-water interface are considering self-assembly systems. The studies of polymers at the air-water interface considering several molecular simulation methodologies, such as atomistic(Gargallo et al., 2009) and coarsegrained(Duncan et al., 2011) models. In this chapter we will focus on atomistic molecular simulation methodology. The aim of studies of polymer at the air-water interface, using molecular simulation, is characterized the orientation of polymer chains at the interface, as well as, calculate the different interaction between amphiphilic polymer and the water

� �� �� ���

macromolecule. The mean-square radius of gyration � ��

where n is the segments numbers of the chain and � ��

length *l*, the following relationship is established:

� ��

� ��

� �� � �

**3. Molecular dynamic simulation of supramolecular structures at the** 

� �� � (3)

� � represents the mean square of

� � represents the mean-square

� (5)

� � � ⁄ � (4)

� �� �� �� �� �� ��

unperturbed state.

defined to following equation 4:

2006; Prosa et al., 2001, 1997).

**air-water interface** 

subphase.

In general, the models used to reproduce several polymer chains at air-water interface, consist of the build a bilayer of polymers with an enough thick water slab between of the chains. All the systems were simulated taking into account three experimental surface areas 10, 20 and 30Å2 per molecules. The water slab should be larger in one axis (i.e. z-axis), to avoid the contact between the polymer chains, see Figure 3. A vacuum space larger enough to avoid the contact between chains of different side of the bilayer, should be considered. To the constructor of different molecular structure (repeating unit of the polymers) any molecular builder is required, the important point is generated the pdb coordinate file. To generate the polymer chains a script was made in Tcl language and them executed in VMD(Humphrey et al., 1996) program and this way, the polymer chain of ten or twelve repeating unit were obtained. The build of water slab was made to using VMD program. For this purpose, was necessary obtained the coordinate of each polymer chain and them located these chains in the correct coordinates in order to build the water slap in the center of the coordinate and this way, obtained the model that represent several polymer chains oriented at the air-water interface, as Figure 3.

Fig. 3. Drawing in two dimension of polymeric monolayer at the air-water interface. The dimension in x and y-axis are smaller in relation with z-axis.

To the build of the inclusion complexes of cyclodextrin-polymers, and -cyclodextrin were considering. The cyclodextrin structures were obtained from Protein Data Bank, PDB entry 1CXF for -CD and 1D3C for -CD. These structure were fully minimized, using CHARMM force field and a good agreement with the experimental geometry was obtained. To carry out the minimization and molecular dynamic simulation of synthetic polymer was necessary to use a force field, which content the atomic parameter appropriate to a particular polymer. To obtained the atomic parameters to Poly(ethylene oxide) (PEO), Poly(tetrahydrofuran) (PTHF) and Poly(-caprolactone) (PEC) we used the CHARMM force field, this force field is open source and it is possible aggregate the parameter because the polymers implicate in the study was necessary calculate the atomic parameter, through quantum mechanics. For this purpose, Paratool package implemented in VMD(Humphrey et al., 1996) software was used. Six models were built, three corresponding to the inclusion complexes: PEC-CD, PEO-CD and PTHF-CD and the polymer respectively at the air water interface, see Figure 4.

#### **3.2 Considerations to molecular dynamic simulations**

All Molecular dynamic (MD) simulations were carried out using NAMD (Kalé et al., 1999) program and the CHARMM27(Kuttel et al., 2002) force field. Simulations were performed in

Fig. 4. Snapshot for the CD-PTHF inclusion complex, cyclodextrins are drawn with line style and PTHF with van der Waals style. The drawing method for water molecules is solid surface. Carbon, hydrogen and oxygen atoms are green, white and red respectively.

NPT ensemble; the non-bonded van der Waals and electrostatic interactions were calculated with a cut-off using a switching function starting at a distance of 10Å and reaching zero at 12Å. The TIP3P(Jorgensen et al., 1983) water model was employed for the solvent. Periodic boundary conditions were applied with a rigid cell. The particle mesh Ewald(Darden et al., 1993) (PME) method was employed for computation of electrostatic forces. An integration time step of 1fs was assumed, permitting a multiple time stepping algorithm(Grubmüller et al., 1991) to be employed in which interactions involving covalent bonds were updated every time step. In all simulations, Langevin dynamics were utilized to keep a constant temperature of 298K, likewise the hybrid Nosé-Hoover Langevin piston method(Izrailev et al., 1997) was used to control a constant pressure of 1 atm. The system was first equilibrated for 1ns. The backbone carbons of the polymers in the inclusion complexes were restrained with a 0.5kcal/mol/Å2 spring constant during relaxation. Then, simulations of 5ns were performed to all the systems. A novel alternative to obtain the charges and potentials to each atoms or groups is through CHARMM General Force Field (CGenFF)(Vanommeslaeghe et al., 2010). This force field is an extension of the CHARMM force field to small molecules and included a wide range of chemical groups. The molecule, in our case the monomer, is required in mol2 format to be uploading in the web site of Paramchem and this way, obtain the charges and potentials atomic. The server provides the topology and parameter files, necessaries to run the MD simulation. Topology file contains the monomer residue, atom types and the charges. Then, is mandatory to make changes to the topology file because we need create a topology for a polymer containing many repeat units. At the moment the topology file only have one residue that corresponding to the monomeric structure, so is require to build three types of residues. A "head" residue that defined the initial monomer chain, the "repeat" residue that is corresponding to the residue in-between chain and "tail" residue that defined the final monomer chain. Another way is to define only one residue corresponding to de monomer and build a path union, where is specified the hydrogen atoms that most be eliminated and the atoms that most be linked. This way is possible obtaining the coordinate and connectivity file to a polymeric chain.

#### **3.3 Properties polymeric systems, calculations and analysis**

From MD simulations trajectories several interesting properties cab be obtained and give insight about the behavior of inclusion complexes and polymer chains at the air-water

Fig. 4. Snapshot for the CD-PTHF inclusion complex, cyclodextrins are drawn with line style and PTHF with van der Waals style. The drawing method for water molecules is solid surface. Carbon, hydrogen and oxygen atoms are green, white and red respectively.

NPT ensemble; the non-bonded van der Waals and electrostatic interactions were calculated with a cut-off using a switching function starting at a distance of 10Å and reaching zero at 12Å. The TIP3P(Jorgensen et al., 1983) water model was employed for the solvent. Periodic boundary conditions were applied with a rigid cell. The particle mesh Ewald(Darden et al., 1993) (PME) method was employed for computation of electrostatic forces. An integration time step of 1fs was assumed, permitting a multiple time stepping algorithm(Grubmüller et al., 1991) to be employed in which interactions involving covalent bonds were updated every time step. In all simulations, Langevin dynamics were utilized to keep a constant temperature of 298K, likewise the hybrid Nosé-Hoover Langevin piston method(Izrailev et al., 1997) was used to control a constant pressure of 1 atm. The system was first equilibrated for 1ns. The backbone carbons of the polymers in the inclusion complexes were restrained with a 0.5kcal/mol/Å2 spring constant during relaxation. Then, simulations of 5ns were performed to all the systems. A novel alternative to obtain the charges and potentials to each atoms or groups is through CHARMM General Force Field (CGenFF)(Vanommeslaeghe et al., 2010). This force field is an extension of the CHARMM force field to small molecules and included a wide range of chemical groups. The molecule, in our case the monomer, is required in mol2 format to be uploading in the web site of Paramchem and this way, obtain the charges and potentials atomic. The server provides the topology and parameter files, necessaries to run the MD simulation. Topology file contains the monomer residue, atom types and the charges. Then, is mandatory to make changes to the topology file because we need create a topology for a polymer containing many repeat units. At the moment the topology file only have one residue that corresponding to the monomeric structure, so is require to build three types of residues. A "head" residue that defined the initial monomer chain, the "repeat" residue that is corresponding to the residue in-between chain and "tail" residue that defined the final monomer chain. Another way is to define only one residue corresponding to de monomer and build a path union, where is specified the hydrogen atoms that most be eliminated and the atoms that most be linked. This way is possible

obtaining the coordinate and connectivity file to a polymeric chain.

**3.3 Properties polymeric systems, calculations and analysis** 

From MD simulations trajectories several interesting properties cab be obtained and give insight about the behavior of inclusion complexes and polymer chains at the air-water interface. The most important point to know about the supramolecular systems(Komarov et al., 2011; Krishnamoorthy et al., 2011) at the interface is the orientational correlation about the order of monolayer. In the case of study of the inclusion complexes(Almeida et al., 2011; Boonyarattanakalin et al., 2011; Pan et al., 2011; Shi et al., 2011), we calculate the *order parameter* (Sr). This property is defined in terms of second-order Legendre polynomials and is the average of an angle made by a section of the polymer chain in a specified direction. This orientation correlation function S(r) is defined by equation:

$$S(r) = 0.5 \{ 3 \cos^2 \theta\_{lj}(r) - 1 \} \tag{6}$$

where is the angle between *i* and *j* vectors, when the center of mass of the associated chains of the vectors *i* and *j* are separated by the distance *r*. For the case of the inclusion complexes, we selected the atoms carbon C-C unit vectors of the different chains from PTHF, PEC and PEO as the vectors *i* and *j*. This parameter is indicative of the tendency towards polymers ordering in systems containing flexible chains. A value of 1.0 for S(r) mean that the chains are parallel to each other and cero value for a disordered distribution.

In the Figure 5, the S(r) variation as a function of time for these polymers is shown. It is clear that PEC is more organized at the air-water interface than PTHF and PEO. These results were in good agreement with the experimental data. Comparing the order parameter for PEC, PTHF and PEO and their respective ICs, can observe that the ICs are more ordered at air-water interface than their precursor polymer, Figure 6 a, b and c.

The order parameter shown differences evident in lose conformational freedom for cyclodextrins-polymers systems. The inclusion complexes are stiffer than single chains explaining the lower values of S(r) of the chains compared to those included in cyclodextrins.

The *radial distribution function*, g(r), is another important property that describe the probability to find an atom in a distance r from another atom and can be represent by equation 7, where Ri is the distance of atoms *i* from the center of mass and N the total number of atoms. By this way, the estimation of g(r) between atoms belonging to hydrophilic groups of the precursor polymer and water, using computational methods, give a view of the orientation of polymer chains at the air-water interface.

$$\langle g\_{\{r\}}^2 = \langle \frac{1}{N} \Sigma\_{l=1}^N \{ R\_l^2 \} \rangle \tag{7}$$

Figure 7 and 8 presents g(r) between the hydrophobic groups of the different polymer chains and also g(r) between hydrophilic groups of the different polymer chains and water molecules. The electrostatic interactions between hydrophilic groups of the precursor polymers and water molecules are established. Nevertheless, some differences were seen in g(r) plots to PEC, PEO and PTHF. In the case of, PEC and PTHF (Figure 7a,b) the distance between the hydrophilic groups and water molecules are slightly lower in compare to the distances in PEC.

The � �� � ��� � �� � � parameters can be estimate from molecular dynamics simulation trajectories, in the same way that the radial distribution function and order parameter. Scripts for each of the parameter were written in Tcl language and were ran in VMD program. The scripts can be to provide the reader through a mail to the author.

Fig. 5. Order parameter S(r) for precursor polymers from MD simulation.

Fig. 6. Order parameter S(r) for a) PEC and IC-CD-PEC, b) PTHF and IC-CD-PTHF, c) PEO and IC-CD-PEO

Fig. 5. Order parameter S(r) for precursor polymers from MD simulation.

Fig. 6. Order parameter S(r) for a) PEC and IC-CD-PEC, b) PTHF and IC-CD-PTHF, c) PEO

and IC-CD-PEO

Fig. 7. Radial distribution function (RDF) between oxygen atoms of the precursor polymers and water molecules for: a) PEC, b) PEO and c) PTHF.

Fig. 8. Radial distribution function between hydrophobic groups of different polymer chains for: a) PEC, b) PEO and c) PTHF.

#### **4. Free energy conformational calculation to the copolymer systems at condensed phase**

The di-block copolymers in melt state have been studied by divers researcher groups, due to great interesting from the point view of the dielectric and dynamic-mechanics behavior(Alig et al., 1992; Encinar et al., 2008; Karatasos et al., 1996; Kotaka and Adachi, 1997; Kyritsis et al., 2000; Ma et al., 2003; Moreno and Rubio, 2002; Vogt et al., 1992). Dynamic-mechanical spectroscopy is an important technique to characterize the deformation of polymeric materials applying, for example, a tension to material. The viscoelastic behavior of the polymers can be studied through the determination of properties such as the complex modulus and glass transition temperature. On the other hand, the dielectric relaxation explore the response of a sample submitted an electric field. The permanent dipolar moments of a flexible polymer can be reoriented by an electric field, giving the relaxation mechanics of the fluctuations of dipolar moieties of the sample. For this reason is interesting to know at atomic level the conformational motions that have a polymer or copolymer. The conformational free energy difference (G) for a polymer chain, can be calculated through molecular dynamics simulation at different temperatures, with the aims of obtain a conformational wide sampling. For this purpose, Boltzmann equation(Gilson et al., 1997) was used to obtain G, equation 8.

$$
\Delta G = -RT \ln \frac{f \phi\_1}{f \phi\_2} \tag{8}
$$

Where *f* is the frequency distribution of the conformational change of the dihedral angle of side chain of PtBa, R is the gas constant and T is the temperature.

The case study was carried out for two systems of diblock copolymer of Polystyrene (PS) and Poly(t-butylacrylate) (PtBa) with different Poly(t-butylacrylate) chain lengths. The two blocks of PS and PtBa are amorphous and immiscible, also these blocks have well separates their glass transitions temperature (Tg). The PS-b-PtBa copolymer is dielectrically active due to PtBa part only. Nevertheless, the PtBa block has no present dipole moment to along the chain, thus the segmental mode is not influenced by other dynamical contributions. Besides, both blocks are active in dynamic-mechanical spectroscopy. The results obtained from the two techniques were compared with the PtBa homopolymer. The experimental results shown that the shape of the relaxation curves has different dependences with the temperature in the copolymers and the PtBa homopolymer. However, the relaxation maps to -mode were almost equivalent in the copolymer than in the PtBa homopolymer, but the relaxation time to -mode were different in the copolymer respect to homopolymer. This result was unexpected, because the -relaxation has a local character due to the PtBa side chains relaxation. The above behavior could be owing to the influence of the PS on the PtBa side chains. For this reason, different environments were simulated in condensed phase and three models were built Figure 9 with the goal of study, if the PS blocks influenced the motion of side chains of PtBa. The parametrization of PtBa and PS monomers were done in the Paratool module of the VMD program. Paratool is a tool that used quantum mechanics methodology to the optimized structures and vibrational frequency to obtain the atomic parameter, both unbounded and bounded generating topology and parameter files, in CHARMM format. The MD simulations were carried out in NAMD program using CHARMM force field.

Fig. 9. Different structure models to study the side chain conformational movement of the PtBa containing different environments, frame in dashed line represent the monomer considered in measurements of frequency to the dihedral angle. a) PtBa monomer surrounded by different neighborhood, both PS and PtBa, labeled P(305) b) PtBa monomer with PS neighbors, labeled P(307) and c) PtBa monomer surrounded by PtBa units, labeled PtBa.

οܩ ൌ െܴ݈ܶ݊ థభ

Where *f* is the frequency distribution of the conformational change of the dihedral angle of

The case study was carried out for two systems of diblock copolymer of Polystyrene (PS) and Poly(t-butylacrylate) (PtBa) with different Poly(t-butylacrylate) chain lengths. The two blocks of PS and PtBa are amorphous and immiscible, also these blocks have well separates their glass transitions temperature (Tg). The PS-b-PtBa copolymer is dielectrically active due to PtBa part only. Nevertheless, the PtBa block has no present dipole moment to along the chain, thus the segmental mode is not influenced by other dynamical contributions. Besides, both blocks are active in dynamic-mechanical spectroscopy. The results obtained from the two techniques were compared with the PtBa homopolymer. The experimental results shown that the shape of the relaxation curves has different dependences with the temperature in the copolymers and the PtBa homopolymer. However, the relaxation maps to -mode were almost equivalent in the copolymer than in the PtBa homopolymer, but the relaxation time to -mode were different in the copolymer respect to homopolymer. This result was unexpected, because the -relaxation has a local character due to the PtBa side chains relaxation. The above behavior could be owing to the influence of the PS on the PtBa side chains. For this reason, different environments were simulated in condensed phase and three models were built Figure 9 with the goal of study, if the PS blocks influenced the motion of side chains of PtBa. The parametrization of PtBa and PS monomers were done in the Paratool module of the VMD program. Paratool is a tool that used quantum mechanics methodology to the optimized structures and vibrational frequency to obtain the atomic parameter, both unbounded and bounded generating topology and parameter files, in CHARMM format. The MD simulations

side chain of PtBa, R is the gas constant and T is the temperature.

were carried out in NAMD program using CHARMM force field.

by PtBa units, labeled PtBa.

Fig. 9. Different structure models to study the side chain conformational movement of the PtBa containing different environments, frame in dashed line represent the monomer considered in measurements of frequency to the dihedral angle. a) PtBa monomer surrounded by different neighborhood, both PS and PtBa, labeled P(305) b) PtBa monomer with PS neighbors, labeled P(307) and c) PtBa monomer surrounded

థమ

(8)

The polymeric fragments containing five repeating units of each sequence were placed in a box under periodic boundary condition Figure 10; ten chains were placed in three different boxes for each model. Molecular simulation of 1ns of long time were carried out using Valbond Force Field (CVFF)(Maple et al., 1998) and CHARMM forcé field. The three models were run at different temperatures between 1000K and 1700K with the purpose of obtain a wide conformational sample. A conformational analysis was performed to each model in order to estimate the free energy difference between two conformations for dihedral angles *1* and *<sup>2</sup>*, see Figure 9.

The frequency average of the dihedral angles was calculate in the three models at the different temperatures from molecular simulation trajectories Figure 11, and applied to the equation 8. The G values were compared with activation energy EA from dielectric experiments to the -relaxation, Table 1.

Fig. 10. Box with periodic boundary condition containing ten polymer chains of P(305).

Fig. 11. Frequency distribution versus degree of the dihedral angles, 1 and 2 of obtained at 1700K, the numbers in the top peak correspond to area under the curve of these angles. The area values were considered in Boltzmann equation.


Table 1. Activation free energy (EA) to -relaxation obtained from dielectric experiments and variations.

The experimental data showed that the characteristic time to -transition is faster in the case of homopolymer compared to the copolymers. This behavior provides evidence that the PS blocks are influencing to the PtBa blocks. From MD simulations was possible to obtain the frequencies distribution to the conformational changes of 1 and 2 dihedral angles, Figure 11. Distribution frequencies calculated for P(305), P(307) and PtBa for all temperatures were applied to equation 8 and them the G was obtain from the plots of ln(*f*1*/f*2) against 1/T. Conformational free energy had the following trend PtBa>P(305)>P(307). Table 1 shown that the G is unfavorable in the case to PtBa, indeed this mean that there is less mobility of ester group, then the relaxation processes should be more hindered than in the other fragments. These results explain the activation energies values obtain by experimental techniques, where the EA to PtBa is higher in compare to P(305) and P(307). Therefore, the results obtain by MD simulation of different fragment of the copolymer are in agree with the experimental data. In addition, the model containing the monomer surrounded by PS and PtBa (P305) have a higher mobility, so the activation energy to the relaxation process is less.

We have seen that the molecular simulation methods can be effective to study behavior of the system in solution as well as in condensed phase. The considerations to carried out molecular simulation in polymeric systems are:


Shake algorithm(Andersen, 1983; Ryckaert et al., 1977) to do cyclic computation of all constraints, in an iterative fashion until each constraint be satisfied.


### **5. Conclusions**

40 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

Polymer 10-3G /R 10-3EA/R PtBa -1.60.1 1.50.04 P(307) -1.90.1 1.30.05 P(305) -2.20.2 1.10.1 Table 1. Activation free energy (EA) to -relaxation obtained from dielectric experiments and

The experimental data showed that the characteristic time to -transition is faster in the case of homopolymer compared to the copolymers. This behavior provides evidence that the PS blocks are influencing to the PtBa blocks. From MD simulations was possible to obtain the frequencies distribution to the conformational changes of 1 and 2 dihedral angles, Figure 11. Distribution frequencies calculated for P(305), P(307) and PtBa for all temperatures were applied to equation 8 and them the G was obtain from the plots of ln(*f*1*/f*2) against 1/T. Conformational free energy had the following trend PtBa>P(305)>P(307). Table 1 shown that the G is unfavorable in the case to PtBa, indeed this mean that there is less mobility of ester group, then the relaxation processes should be more hindered than in the other fragments. These results explain the activation energies values obtain by experimental techniques, where the EA to PtBa is higher in compare to P(305) and P(307). Therefore, the results obtain by MD simulation of different fragment of the copolymer are in agree with the experimental data. In addition, the model containing the monomer surrounded by PS and PtBa (P305)

have a higher mobility, so the activation energy to the relaxation process is less.

1. To know if the process of study occur in solution or condensed phase

3. Establish that force field is more appropriate to our polymeric system.

in VMD program or using the web platform of Paramchem.

6. Building models should represent the issue to be solved. 7. Select to appropriate program to carried out the MD simulation

add to the force field, in the case of those of open source, as CHARMM.

molecular simulation in polymeric systems are:

coarse-grained methodology.

We have seen that the molecular simulation methods can be effective to study behavior of the system in solution as well as in condensed phase. The considerations to carried out

2. The time scale that take this process, so with that information is possible decide which methodology can be better to establish an equilibrium between simulation time and accuracy of results. Some process, as the mesoscopic structures, should be studied with

4. In the absence of atomic parameter is necessary get through quantum mechanics and

5. The atomic parameter can be to get of an easy way through Paratool module available

8. Depending of the issue of our system, may be necessary apply to the polymer system constrain. Oftentimes, is required the study of movement of the side chain and not necessary to considered the backbone polymeric movement, in this case can be considerer apply restriction to the movement of the backbone. The restraints to the motion have an important advantage that allows simulation time economy. The important is emphasizing de differences between restraints and constraints. The constraints eliminate the hard degrees of freedom of the systems that correspond to high frequencies of vibration. The constraints can be applied on the hydrogen bond distances and hydrogen bond angles, for this purpose the shake algorithm is used.

variations.

Molecular dynamic simulations can became a powerful tool to study of the behavior of polymer system both in solution and condensed phase. In this chapter, two different works were approached, in both cases the intermolecular interactions were calculated for the purpose of established a correlation between theoretical and experimental results. In the case of inclusion complexes, molecular simulation was carried out to understand the behavior at the air-water interface. The results obtained from molecular simulation trajectories to each system were used to calculate the order parameter and correlated with - A isotherms from Langmuir experiments. Order parameter established that both PEC and IC-CD-PEC were more ordered systems related to PTHF and IC-CD-PTHF and PEO and IC-CD-PEO. These results were agreement with experimental results. However, radial distribution function was calculated with the purpose of obtain information about orientation of polymeric systems at the interface. From radial distribution function could be established that the interaction between polar groups of the polymers and water molecules took place around 2Å distance, then was possible estimate that there is a diffusion of the polymer to the water. Also, radial distribution function between hydrophobic groups of the different chains was calculated. The distances between polymeric chains were considerably short so the chains adopted a conformation rather widespread; with the purpose of interact among them. In the second case of study, molecular dynamics of Polystyrene-bock-Poly(tbutylacrylate) (PS-b-PtBa) copolymers were carried out at different Polystyrene composition for the purpose of observed the influence of PS in the experimental dielectric relaxation of PtBa. From molecular simulation trajectories the conformational free energy was calculated to the different PS composition in the copolymers. Conformational free energies calculate from Boltzmann equation were correlated with activation free energies calculated from experimental dielectric relaxation. Activation free energy to PtBa was higher compared to P(305) and P(307), this mean that a more hindered change conformational was observed. A less favored change conformational was obtained to PtBa compared to P(305). Therefore, good agree was obtained between the activation free energy results and conformational free energy from molecular simulation. Finally, it is possible conclude that molecular dynamics simulation is a important and challenging tool that can be used to know properties and processes that occurs in polymer systems both dissolution and condensed phase. A important point to take in account is dispose of appropriate atomic parameters, that mean an adequate force field and the other hand, a model that represent the polymeric system.

### **6. Acknowledgment**

This work was supported in part by "Proyecto de Inserción en la Academia CONICYT Nº79090038" and InnovaChile CORFO Code FCR-CSB 09CEII-6991.

#### **7. References**


Alig, I., Kremer, F., Fytas, G., and Roovers, J., 1992. DIELECTRIC-RELAXATION IN

Almeida, E.W.C., Anconi, C.P.A., Novato, W.T.G., De Oliveira, M.A.L., De Almeida, W.B.,

Andersen, H., 1983. Rattle: A Velocity Version of the Shake Algorithm for Molecular

Boonyarattanakalin, K.S., Wolschann, P., and Lawtrakul, L., 2011. Molecular dynamics of

Duncan, S.L., Dalal, I.S., and Larson, R.G., 2011. Molecular dynamics simulation of phase

Encinar, M., Guzman, E., Prolongo, M.G., Rubio, R.G., Sandoval, C., Gonzalez-Nilo, F.,

Gargallo, L., Vargas, D., Becerra, N., Sandoval, C., Saldias, C., Leiva, A., and Radic, D., 2009.

Gargallo, L., Becerra, N., Sandoval, C., Pitsikalis, M., Hadjichristidis, N., Leiva, A., and

Gilson, M.K., Given, J.A., Bush, B.L., and McCammon, J.A., 1997. The statistical-

Izrailev, S., Stepaniants, S., Balcera, M., Oono, Y., and Schulten, K., 1997. Biophys. J. 72, 1568. Jorgensen, W., Chandrasekhar, J., Madura, J., Impey, R., and Klein, M., 1983. J. Chem. Phys.

Kalé, L., Skeel, R., Bhadarkar, M., Brunner, R., Gursoy, A., Krawetz, N., and al., e., 1999. J.

Karatasos, K., Anastasiadis, S.H., Floudas, G., Fytas, G., Pispas, S., Hadjichristidis, N., and

Pakula, T., 1996. Composition fluctuation effects on dielectric normal-mode

Grubmüller, H., Heller, H., Windemuth, A., and Schulten, K., 1991. Mol. Simul. 6, 121. Humphrey, W., Dalke, A., and Schulten, K., 1996. VMD: Visual molecular dynamics. Journal

Dynamics Calculations. J. Comput. Phys. 52, 24-34.

Darden, T., York, D., and Pedersen, L., 1993. J. Chem. Phys. 98, 10089.

poly(t-butylacrylate). Polymer 49, 5650-5658. Flory, J.P., 1969. Statistical mechanics of chain molecules, NY.

Applied Polymer Science 122, 1395-1404.

Biophysical Journal 72, 1047-1069.

of Molecular Graphics 14, 33-&.

Comput. Phys. 151, 283.

79, 926.

Flory, J.P., 1995. Principles of polymer chemistry, United state of america.

DISORDERED POLY(ISOPRENE STYRENE) DIBLOCK COPOLYMERS NEAR THE MICROPHASE-SEPARATION TRANSITION. Macromolecules 25, 5277-5282.

and Dos Santos, H.F., 2011. Box-Behnken design for studying inclusion complexes of triglycerides and alpha-cyclodextrin: application to the heating protocol in molecular-dynamics simulations. Journal of Inclusion Phenomena and Macrocyclic

beta-CD in water/co-solvent mixtures. Journal of Inclusion Phenomena and

transitions in model lung surfactant monolayers. Biochimica Et Biophysica Acta-

Gargallo, L., and Radic, D., 2008. Dielectric and dynamic-mechanical study of the mobility of poly(t-butylacrylate) chains in diblock copolymers: Polystyrene-b-

Supramolecular Structures, Organization and Surface Behavior at Interfaces.

Radic, D., 2011. Amphiphilic Diblock Copolymers Containing Poly(Nhexylisocyanate): Monolayer Behavior at the Air-Water Interface. Journal of

thermodynamic basis for computation of binding affinities: A critical review.

**7. References** 

Chemistry 71, 103-111.

Macrocyclic Chemistry 70, 279-290.

Donald R. P., B., C.B., 2000. Polymer Blends.

Biomembranes 1808, 2450-2465.

Macromol. Symp. 278, 80-88.

relaxation in diblock copolymers .2. Disordered state in proximity to the ODT and ordered state. Macromolecules 29, 1326-1336.


Ryckaert, J.P., Ciccotti, G., and Berendsen, H.J.C., 1977. J. Comput. Phys. 23, 327.

Saldias, C., Gargallo, L., Sandoval, C., Leiva, A., Radic, D., Caballero, J., Saavedra, M., and Gonzalez-Nilo, F.D., 2009. Inclusion complexes containing poly(epsiloncaprolactone)diol and cyclodextrins. Experimental and theoretical studies. Polymer 50, 2926-2932.


## **Backbone Connectivity and Collective Aggregation Phenomena in Polymer Systems**

Wen-Jong Ma1,2 and Chin-Kun Hu<sup>2</sup>

<sup>1</sup>*Graduate Institute of Applied Physics, National Chengchi University, Taipei 11605* <sup>2</sup>*Institute of Physics, Academia Sinica, Nankang, Taipei 11529 Taiwan*

#### **1. Introduction**

44 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

Shi, J.H., Hu, Y., and Ding, Z.J., 2011. Theoretical study on chiral recognition mechanism of

Vanommeslaeghe, K., Hatcher, E., Acharya, C., Kundu, S., Zhong, S., Shim, J., Darian, E.,

Vogt, S., Gerharz, B., Fischer, E.W., and Fytas, G., 1992. LOCAL MOTIONS IN A

Wick, C.D., and Cummings, O.T., 2011. Understanding the factors that contribute to ion

DIBLOCK COPOLYMER MELT. Macromolecules 25, 5986-5990.

interfacial behavior. Chemical Physics Letters 513, 161-166.

Yamakawa, H., 2001. Modern Theory of Polymer Solutions, NY.

and Theoretical Chemistry 973, 62-68.

Chemistry 31, 671-690.

ethyl-3-hydroxybutyrate with permethylated beta-cyclodextrin. Computational

Guvench, O., Lopes, P., Vorobyov, I., and MacKerell, A.D., 2010. CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. Journal of Computational

MICROPHASE-SEPARATED POLY(STYRENE-B-METHYLPHENYLSILOXANE)

A convenient theoretical entry to tackle the protein aggregation problems is to track on the origin of the aggregation phenomena from the viewpoint of polymer physics, enquiring how the underlying factors, such as the geometric packing, the topology and their mutual interplays in presence of solvent, play their roles in the processes. Such a view has been inspired by recent experiments [1, 2] and molecular dynamics simulations [3, 4] which suggest the ubiquitous presence of fibril formation [5] in various natural and laboratory prepared proteins or peptides. While the variety of amino acid sequences interferes with the occurrence of long-range structural ordering, a material-insensitive tendency of aggregation is observed [2, 3, 5]. By coarsening the sequence-sensitive details, the aggregation problem can be formulated in its minimal form as the clustering process of polymer chains [6, 7]. With such simplified models, the approach focuses more on the entropic effect caused by the constraint of chain connectivity [6–8], rather than on the material-dependent characteristics. In this chapter, we summarize our molecular dynamics simulation studies [6, 7] that reveal the relationship between the backbone connectivity of polymer chains and some benchmark features displayed in the aggregation processes of the model polymer chains.

In the model, the backbone connectivity of a polymer chain is realized by assigning a string of monomers with specific monomer-monomer two-body forces, perturbed with three-body and four-body angle dependent interactions [6], with their strengths measured by the parameters, *k*nn, *K*<sup>b</sup> and *K*t, for the nearest neighbor (n.n.) interaction, the bending angle and the torsion angle potentials, respectively. In a collection of polymer chains, all the non n.n. pairs along a chain and those pairs on different chains are subject to Lennard-Jones (L-J) pair interactions. The presence of angle potentials, with significantly nonzero *K*<sup>b</sup> or *K*<sup>t</sup> values, breaks the isotropy surround the backbone. In such a system, the local inter-chain hindrance prevents the chains from clustering into ordered domains. The situation can, however, be reverted by reducing the values of *K*<sup>b</sup> and *K*t. We find that the formation of bundle-like domains (Fig. 1) is robust as soon as *K*<sup>b</sup> and *K*<sup>t</sup> are small enough. The observations are assured even if we introduce a small amount of dispelling background fluid molecules or impose a small fraction of impurity monomer sites.

Fig. 1. Bundled domains formed in systems of homopolymer chains, with zero (left plot) or tiny (right plot) values in *K*<sup>b</sup> and *K*t. They are snapshots for the systems with rigid bonds, labeled by I-∞ and III, respectively, in Table 1.

In quantifying the overall structural changes over the relaxation processes, we reveal the shared scenario for the cluster forming procedure in a range of systems (see Table 1) . In this scenario, the ordering of the clustered (bundled) domains can occur spontaneously as soon as the energy barriers contributed by the local structural hindrance are overcome. Under such circumstances, the completion of aggregation always occurs at a temperature well above the characteristic temperature determined by the strength of the site-site L-J pair potential. This is in contrast to the formation of ordered clusters in the two phase coexistence region for a L-J fluid of the same number of monomers. (The system is equivalently obtained by removing those n.n. bonding along the polymer chains.) The clustering can only occur below the critical point temperature, which is about 1.3 times of the characteristic temperature /*kB* (*kB*: Boltzmann constant). To have solid-like structured cluster, the temperature has to go below the triple point temperature, which is about 0.6/*kB*. The result suggests that the key role played by the chain connectivity is via the lowering of the entropy of the system, so that the formation of ordered domains, becomes feasible at a temperature not necessarily very low for the system to overcome the local structural hindrance.

The bonding strength, parameterized by *k*nn, of the nearest neighboring monomers along the chains affects the detailed procedures in the formation of clusters during the aggregation processes. More than one stages are present during the processes, with the local alignment of segments in the individual chains followed by the coalescence of these short patches. The second stage is found further to differentiate into several stages with the increased *k*nn. We allow *k*nn to have infinity value, which is realized in simulation as a (rigid) bond with fixed length maintained by the constraint force (see Section 2.1). In manipulating the parameters *k*nn, *K*<sup>b</sup> and *K*<sup>t</sup> in a range of systems, we can obtain more precise information about how the backbone anisotropy prevent the polymer chains from the formation of bundled structure. The scenario unveiled in this study may be useful for the analysis of protein aggregation phenomena in realistic complex biological systems.

Fig. 1. Bundled domains formed in systems of homopolymer chains, with zero (left plot) or tiny (right plot) values in *K*<sup>b</sup> and *K*t. They are snapshots for the systems with rigid bonds,

In quantifying the overall structural changes over the relaxation processes, we reveal the shared scenario for the cluster forming procedure in a range of systems (see Table 1) . In this scenario, the ordering of the clustered (bundled) domains can occur spontaneously as soon as the energy barriers contributed by the local structural hindrance are overcome. Under such circumstances, the completion of aggregation always occurs at a temperature well above the characteristic temperature determined by the strength of the site-site L-J pair potential. This is in contrast to the formation of ordered clusters in the two phase coexistence region for a L-J fluid of the same number of monomers. (The system is equivalently obtained by removing those n.n. bonding along the polymer chains.) The clustering can only occur below the critical point temperature, which is about 1.3 times of the characteristic temperature /*kB* (*kB*: Boltzmann constant). To have solid-like structured cluster, the temperature has to go below the triple point temperature, which is about 0.6/*kB*. The result suggests that the key role played by the chain connectivity is via the lowering of the entropy of the system, so that the formation of ordered domains, becomes feasible at a temperature not necessarily very low

The bonding strength, parameterized by *k*nn, of the nearest neighboring monomers along the chains affects the detailed procedures in the formation of clusters during the aggregation processes. More than one stages are present during the processes, with the local alignment of segments in the individual chains followed by the coalescence of these short patches. The second stage is found further to differentiate into several stages with the increased *k*nn. We allow *k*nn to have infinity value, which is realized in simulation as a (rigid) bond with fixed length maintained by the constraint force (see Section 2.1). In manipulating the parameters *k*nn, *K*<sup>b</sup> and *K*<sup>t</sup> in a range of systems, we can obtain more precise information about how the backbone anisotropy prevent the polymer chains from the formation of bundled structure. The scenario unveiled in this study may be useful for the analysis of protein aggregation

labeled by I-∞ and III, respectively, in Table 1.

for the system to overcome the local structural hindrance.

phenomena in realistic complex biological systems.


Table 1. List of systems

‡ systems of pure isotropic (*K*<sup>b</sup> = *K*<sup>t</sup> = 0) homopolymer chains (see Ref.[6, 7]); § systems of chains with tiny backbone anisotropy (*K*<sup>b</sup> <sup>=</sup> *<sup>K</sup>*<sup>t</sup> <sup>=</sup> <sup>10</sup>−4) and with <sup>≈</sup> 5% (199 out of *nNP* =4000) impurity monomers (see Ref.[6, 7]), mixed with fluid atoms (*NF* > 0) ; � a pure system of rigidly bonded homopolymer chains with tiny backbone anisotropy and without impurity monomers (see Ref.[7]);

† systems of anisotropic (*K*<sup>b</sup> <sup>=</sup> *<sup>K</sup>*<sup>t</sup> <sup>=</sup> 0.1) chains that fail to aggregate; there being <sup>≈</sup> 5% (199 out of *nNP* =4000) impurity monomers (see Ref.[6, 7]) and mixing with fluid atoms (*NF* > 0).

#### **2. Model and analysis**

#### **2.1 Model systems**

#### **nearest neighbor bonding**

We consider systems containing *NP* = 40 polymer chains and each chain has *n* = 100 monomers. Two major sets of systems are studied. The first set (set I) are pure (*NF* = 0) systems of homopolymer chains with no built-in anisotropy with respect to the backbone (*K*<sup>b</sup> = *K*<sup>t</sup> = 0), which can be classified phenomenologically as a bead-spring (if *k*nn < ∞) or a bead-rod (if *k*nn = ∞) model [9, 10]. The nearest neighbors *i* and *j* at a distance *rij* in our bead-spring chains are connected by the potential

$$V\_{\mathbf{s}}(r\_{i\bar{j}}) = \frac{1}{2}k\_{\mathbf{s}}(r\_{i\bar{j}} - r\_0)^2$$

Fig. 2. Bending angle *θi*,*j*,*<sup>k</sup>* and torsion angle *φijkl* extended by the consecutive monomers *i*, *j* = *i* + 1, *k* = *i* + 2 and *l* = *i* + 3 along a chain.

with its strength constant *ks*, the *k*nn = 10*<sup>s</sup>* multiple of the value *k*<sup>0</sup> (i.e. *ks* = *k*nn*k*0) and balance length *r*0, for the five systems with *s* = 1, 2, 3, 4 (*k*<sup>0</sup> and *r*<sup>0</sup> are given below). They are labelled as I-*s* (I-1, ..,I-4, etc) systems in Table 1. The bead-rod chains in our model are chains with strict constant nearest neighbor bond lengths and are realized numerically by Lagrange's constraint forces using RATTLE numerical scheme [11]. The system is labelled as system I-∞ in Table 1.

In all these model systems, the pair interactions between the non-neighboring monomers along a chain or between pairs in different molecules, *i* and *j* at a separation *rij*, are L-J potentials

$$V\_{L\boldsymbol{l}}(r\_{i\boldsymbol{j}}) = 4\boldsymbol{\epsilon} ((r\_{i\boldsymbol{j}}/\sigma)^{-12} - (r\_{i\boldsymbol{j}}/\sigma)^{-6})\boldsymbol{\cdot}$$

In terms of the distance parameter *σ* and the strength parameter *�* of L-J potential, we choose *<sup>k</sup>*<sup>0</sup> <sup>=</sup> 1.5552 <sup>×</sup> 105*�σ*−<sup>2</sup> and *<sup>r</sup>*<sup>0</sup> <sup>=</sup> 0.357*<sup>σ</sup>* [6, 12].

#### **backbone anisotropy along the chains**

In the second set (set II) of systems, the degrees of backbone anisotropy along the chains are introduced by perturbing bending potential,

$$\Theta(\theta\_{i,j,k}) = \mathcal{K}\_{\mathsf{b}}\mathbf{c}\_{\mathsf{b}} \left(\cos\theta\_{ijk} - \cos\theta\_{0}\right)^{2}\mathsf{b}$$

determined by the bending angle *θi*,*j*,*<sup>k</sup>* of three consecutive monomers *i*, *j* and *k* along a chain; and the torsion potential [13]

$$\Phi(\phi\_{ijkl}) = K\_{\text{ft}} \sum\_{\iota=0}^{3} a\_{\iota} (\cos \phi\_{i,j,k,l})^{\iota}$$

as a function of the torsion angle *φijkl* extended by four consecutive monomers *i*, *j*, *k* and *l* (see Fig. 2). The constants *c*b, *θ*<sup>0</sup> [12]; and *a*1, *a*<sup>2</sup> and *a*<sup>3</sup> [13] are given so that a model polyethelyne

Fig. 2. Bending angle *θi*,*j*,*<sup>k</sup>* and torsion angle *φijkl* extended by the consecutive monomers *i*,

with its strength constant *ks*, the *k*nn = 10*<sup>s</sup>* multiple of the value *k*<sup>0</sup> (i.e. *ks* = *k*nn*k*0) and balance length *r*0, for the five systems with *s* = 1, 2, 3, 4 (*k*<sup>0</sup> and *r*<sup>0</sup> are given below). They are labelled as I-*s* (I-1, ..,I-4, etc) systems in Table 1. The bead-rod chains in our model are chains with strict constant nearest neighbor bond lengths and are realized numerically by Lagrange's constraint forces using RATTLE numerical scheme [11]. The system is labelled as system I-∞

In all these model systems, the pair interactions between the non-neighboring monomers along a chain or between pairs in different molecules, *i* and *j* at a separation *rij*, are L-J

*VLJ*(*rij*) = <sup>4</sup>*�*((*rij*/*σ*)−<sup>12</sup> <sup>−</sup> (*rij*/*σ*)−6). In terms of the distance parameter *σ* and the strength parameter *�* of L-J potential, we choose

In the second set (set II) of systems, the degrees of backbone anisotropy along the chains are

<sup>Θ</sup>(*θi*,*j*,*k*) = *<sup>K</sup>*b*c*b(*cosθijk* <sup>−</sup> *cosθ*0)2,

determined by the bending angle *θi*,*j*,*<sup>k</sup>* of three consecutive monomers *i*, *j* and *k* along a chain;

3 ∑ *ι*=0

as a function of the torsion angle *φijkl* extended by four consecutive monomers *i*, *j*, *k* and *l* (see Fig. 2). The constants *c*b, *θ*<sup>0</sup> [12]; and *a*1, *a*<sup>2</sup> and *a*<sup>3</sup> [13] are given so that a model polyethelyne

*aι*(*cosφi*,*j*,*k*,*l*)*<sup>ι</sup>*

Φ(*φijkl*) = *K*<sup>t</sup>

*j* = *i* + 1, *k* = *i* + 2 and *l* = *i* + 3 along a chain.

*<sup>k</sup>*<sup>0</sup> <sup>=</sup> 1.5552 <sup>×</sup> 105*�σ*−<sup>2</sup> and *<sup>r</sup>*<sup>0</sup> <sup>=</sup> 0.357*<sup>σ</sup>* [6, 12].

introduced by perturbing bending potential,

**backbone anisotropy along the chains**

and the torsion potential [13]

in Table 1.

potentials

Fig. 3. *<sup>R</sup>*<sup>0</sup> (black), *<sup>R</sup>*<sup>1</sup> (blue) and *<sup>R</sup>*<sup>2</sup> (red) versus −*U*inter−chain, for systems set I: I-1, I-2, I-3, I-4 and set II: II-0, II-1, II-2, II-3 and II-4, which contain chains with their nearest neighbor connected by spring forces (see Table 1).

chain is obtained by choosing *r*<sup>0</sup> = 0.357*σ*, and let the strength control parameters *K*<sup>b</sup> and *K*<sup>t</sup> be unity [12, 14]. With a range of values for *K*<sup>b</sup> and *K*<sup>t</sup> and a choice of bending angle parameter *θ*0, we are able to prepare chains with different degrees of local conformational hindrance [6, 15]. In our simulations, we prepare several systems, each composed of identical chains with *K*<sup>b</sup> = *K*t, which are allowed to have the value 0, 10−<sup>4</sup> or 0.1 (see Table 1). The model mimics qualitatively the dispersed local conformational degrees of freedom that is present in those structured monomers, such as the amino-acids in protein molecules.

For systems in set II, the chains is further perturbed by randomly imposing a fraction of 5% monomers to have smaller sizes. In merging with monomer-repelling Lennard-Jones fluid atoms, the heterogeneity is further enhanced by rendering these 5% monomers to interact, in addition to the repulsive force, with an attractive part with the fluid atoms [6, 7]. Under the convention mentioned above, we label those systems as II-*s* (*s* = 0, 1, 2, 3, 4 or ∞) in Table 1, according to the exponent of the multiple (*k*nn = 10*<sup>s</sup>* or ∞) of their nearest bonding strength *ks* = *k*nn*k*0.

Three additional systems of pure (*NF* = 0) homogeneous rigidly-bonded chains are also studied. Two of them that contain the same rigidly bonded chains, of shorter chain lengths

Fig. 4. Time evolutions of (a) instantaneous temperature *T*∗; (b) the parameters *R*<sup>0</sup> (solid line), *R*<sup>1</sup> (dash dot line) and *R*<sup>2</sup> (dashed line); and (c) the inter-chain potential energy, for system I-4 (adapted from Ref. [6]).

(*n* = 25 and *n* = 50, respectively), are included in set I. They are used to underscore the effect of chain length. The other system consists of homopolymer chains to have the same backbone anisotropy (*K*<sup>b</sup> > 0 and *K*<sup>t</sup> > 0) as those in set II, but in absence of the 5% impurity monomers. We list this system as system III. Table I lists all the systems that are analyzed in this chapter. Each of them contains the same number *nNp* = 4000 of monomers.

To show that the increased backbone anisotropy indeed hinders the ability to aggregate, we list the systems set IV in Table 1, which are different from their counterparts in set II, only in the 1000 times larger values in *K*<sup>b</sup> and *K*<sup>t</sup> (=0.1).

#### **quenching**

According to the scenario of clustering, the polymer chains aggregate when the temperature and the density of the systems fall within the coexistence region of the phase diagram. In the systems of polymer chains with rigid bonds, a quenching can be achieved by the numerical effect for the cases *<sup>K</sup>*<sup>b</sup> <sup>=</sup> *<sup>K</sup>*<sup>t</sup> <sup>≤</sup> <sup>10</sup>−4, that a convergence to satisfy holonomic constraints on bond lengths lead the dynamic system toward the low temperature attractor [7]. The system is, therefore, quenched spontaneously. With the bond constraints replaced by the confining of soft interaction potentials, on the other hand, the incoherence among the springs in different bonds does not favor a global numerical convergence. In order to lower the temperature of the polymer chains in this case, we have to control the simulated system manually by using a thermostat.

Since the distributions of the monomer velocities are known to deviate systematically from the standard Maxwell-Boltzmann with the increased strength of nearest neighbor bonding [16, 17] and are described by the Tsallis *q*-statistics [18], if the system reaches a "quasi-steady state", it is inappropriate to use the existing thermostating methods that would force the simulated system to converge to a state described by standard Maxwell-Boltzmann statistics. A naive choice without imposing presumed statistics is simply controlling *only* the mean value of the kinetic distributions. In the simplest version of the "velocity scaling" thermostat [19], we

Fig. 4. Time evolutions of (a) instantaneous temperature *T*∗; (b) the parameters *R*<sup>0</sup> (solid line), *R*<sup>1</sup> (dash dot line) and *R*<sup>2</sup> (dashed line); and (c) the inter-chain potential energy, for

Each of them contains the same number *nNp* = 4000 of monomers.

the 1000 times larger values in *K*<sup>b</sup> and *K*<sup>t</sup> (=0.1).

(*n* = 25 and *n* = 50, respectively), are included in set I. They are used to underscore the effect of chain length. The other system consists of homopolymer chains to have the same backbone anisotropy (*K*<sup>b</sup> > 0 and *K*<sup>t</sup> > 0) as those in set II, but in absence of the 5% impurity monomers. We list this system as system III. Table I lists all the systems that are analyzed in this chapter.

To show that the increased backbone anisotropy indeed hinders the ability to aggregate, we list the systems set IV in Table 1, which are different from their counterparts in set II, only in

According to the scenario of clustering, the polymer chains aggregate when the temperature and the density of the systems fall within the coexistence region of the phase diagram. In the systems of polymer chains with rigid bonds, a quenching can be achieved by the numerical effect for the cases *<sup>K</sup>*<sup>b</sup> <sup>=</sup> *<sup>K</sup>*<sup>t</sup> <sup>≤</sup> <sup>10</sup>−4, that a convergence to satisfy holonomic constraints on bond lengths lead the dynamic system toward the low temperature attractor [7]. The system is, therefore, quenched spontaneously. With the bond constraints replaced by the confining of soft interaction potentials, on the other hand, the incoherence among the springs in different bonds does not favor a global numerical convergence. In order to lower the temperature of the polymer chains in this case, we have to control the simulated system manually by using a

Since the distributions of the monomer velocities are known to deviate systematically from the standard Maxwell-Boltzmann with the increased strength of nearest neighbor bonding [16, 17] and are described by the Tsallis *q*-statistics [18], if the system reaches a "quasi-steady state", it is inappropriate to use the existing thermostating methods that would force the simulated system to converge to a state described by standard Maxwell-Boltzmann statistics. A naive choice without imposing presumed statistics is simply controlling *only* the mean value of the kinetic distributions. In the simplest version of the "velocity scaling" thermostat [19], we

system I-4 (adapted from Ref. [6]).

**quenching**

thermostat.

Fig. 5. Snapshots of the configurations of system I-4, for the process described by Fig. 4, in the pure system (*K*<sup>b</sup> = *K*<sup>t</sup> = 0) with *k*spring = *k*4, (a) before (*t*=-47, the starting time in Fig. 4) and (b) on the initiation of aggregation (*t* = 0); (c), (d), and (e), respectively, at the onset step of the first, the second and the third fast decaying for *R*<sup>2</sup> (the time spots of which are marked by vertical lines in Fig. 4(b)); (f) the last step in Fig. 4 (adapted from Ref. [6]).

Fig. 6. Same as Fig. 4, for system II-4, which has tiny values for *K*<sup>b</sup> = *K*<sup>t</sup> = 10−4. (adapted from Ref. [6]).

re-scale the velocities of the monomers and the L-J molecules at *every* step by a factor forcing the mean kinetic energy for the *whole* system to the value of the target temperature *T*∗ 0 .

#### **2.2 Parameters to identify clustering**

To quantify the overall structural change of the simulated systems over their non-equilibrium relaxation processes, we define the following parameters. For a system of *NP* polymer chains each with *n* monomers in a simulation cubic box of size *L* × *L* × *L* subject to periodic boundary conditions, we use the parameters [7]

$$R\_k \equiv \frac{\int\_0^{\frac{L}{2}} g\_{\text{interchain}}(r) r^k dr}{\int\_0^{\frac{L}{2}} \Xi(g\_{\text{interchain}}(r)) r^k dr} \tag{1}$$

for *k* = 0, 1, 2, to monitor the degree of clustering, where Ξ(*x*) = 1 if *x* �= 0 and Ξ(*x*) = 0 if *x* = 0. The parameters are evaluated by finding the interchain pair distribution function *g*interchain(*r*), which is defined as

$$\begin{aligned} \mathcal{S}\_{\text{interchain}}(r) 4\pi r^2 dr &= \frac{1}{N\_m} \sum\_{i=1}^{N\_m} \{ \frac{V}{N\_m - n} \\ &\times \sum\_{\substack{j(i \text{ and } j \text{ in different chains}) \\ r < r' \le r + dr}} \delta(|\vec{r}\_{ij}| - r') \} \end{aligned} \tag{2}$$

for *r* < *L*/2, where *Nm* = *nNP* is the total number of monomers. *R*0, *R*<sup>1</sup> and *R*<sup>2</sup> provide information about the conformational features of the clusters at different stages of the growth process. In particular, the condition for uniformly distributed situation, *R*<sup>0</sup> = *R*<sup>1</sup> = *R*<sup>2</sup> = 1 can be used as a criteria to locate the time spot of initiation of the aggregation process when such a condition starts to become untrue.

Fig. 6. Same as Fig. 4, for system II-4, which has tiny values for *K*<sup>b</sup> = *K*<sup>t</sup> = 10−4. (adapted

re-scale the velocities of the monomers and the L-J molecules at *every* step by a factor forcing the mean kinetic energy for the *whole* system to the value of the target temperature *T*∗

To quantify the overall structural change of the simulated systems over their non-equilibrium relaxation processes, we define the following parameters. For a system of *NP* polymer chains each with *n* monomers in a simulation cubic box of size *L* × *L* × *L* subject to periodic boundary

<sup>0</sup> *<sup>g</sup>*interchain(*r*)*rkdr*

<sup>0</sup> Ξ(*g*interchain(*r*))*rkdr*

*Nm* ∑ *i*=1

{ *<sup>V</sup> Nm* − *n*

*δ*(|�*rij*| − *r*�

for *k* = 0, 1, 2, to monitor the degree of clustering, where Ξ(*x*) = 1 if *x* �= 0 and Ξ(*x*) = 0 if *x* = 0. The parameters are evaluated by finding the interchain pair distribution function

> <sup>2</sup>*dr* <sup>=</sup> <sup>1</sup> *Nm*

for *r* < *L*/2, where *Nm* = *nNP* is the total number of monomers. *R*0, *R*<sup>1</sup> and *R*<sup>2</sup> provide information about the conformational features of the clusters at different stages of the growth process. In particular, the condition for uniformly distributed situation, *R*<sup>0</sup> = *R*<sup>1</sup> = *R*<sup>2</sup> = 1 can be used as a criteria to locate the time spot of initiation of the aggregation process when

*j*(*i* and *j* in different chains) *r*<*r*� ≤*r*+*dr*

*Rk* ≡

*g*interchain(*r*)4*πr*

 *L* 2

 *L* 2

× ∑

0 .

)} (2)

(1)

from Ref. [6]).

**2.2 Parameters to identify clustering**

conditions, we use the parameters [7]

*g*interchain(*r*), which is defined as

such a condition starts to become untrue.

Fig. 7. Snapshots of the configurations for the aggregation process in system II-4, i.e. the mixed system with *<sup>K</sup>*<sup>b</sup> <sup>=</sup> *<sup>K</sup>*<sup>t</sup> <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> and *<sup>k</sup>*spring <sup>=</sup> *<sup>k</sup>*4, at (a) *<sup>t</sup>*=0; (b) the initiation of aggregation (*R*<sup>0</sup> ≈ *R*<sup>1</sup> ≈ *R*<sup>2</sup> ≈ 1); the onset step of fast decaying of (c) *R*0, (d) *R*<sup>1</sup> and (e) *R*2; (f) the last step of Fig. 6. The time spots for (c), (d) and (e) are marked by vertical lines in Fig. 6(b) (adapted from Ref. [6]).

The formation of clusters can be identified by the trend to have the lowered total interchain (L-J) potential energy *U*interchain and a decreasing virial *υ*interchain (see Ref. [7]). Since both *υ*interchain and *U*interchain follow the same trend [6, 7] and the latter is subject to much less fluctuations than the former is, we adapt here *U*interchain as the major quantity to monitor the process.

#### **3. Aggregation of polymer chains**

#### **3.1 Stages of the clustering process**

The formation of clusters can be identified by the tendency to have lowered total interchain (L-J) potential energy *U*interchain. Figure 3 shows the parameters *R*0, *R*<sup>1</sup> and *R*<sup>2</sup> decrease as the functions of *U*interchain, over the dynamic relaxation processes for the nine systems with soft n.n. bonding (*k*nn < ∞), for either pure homopolymer chains (*NF*=0 and *K*<sup>b</sup> = *K*<sup>t</sup> = 0, i.e. systems I-1, I-2, I-3 and I-4) or mixed heterogeneous chains with backbone anisotropy (*NF* > 0 and *K*<sup>b</sup> > 0 or *K*<sup>t</sup> > 0, i.e. systems II-0, II-1, II-2, II-3 and II-4) listed in Table 1. The data are collected in equal time intervals [6, 7]. The increasing of the quantity '−*U*inter−chain', which is used as the variable for the horizontal axis, corresponds to the direction for the growth of clusters. The parameter *R*<sup>0</sup> decreases, more or less, smoothly and is concave upward for all systems before reaching the tail regimes with densely accumulated data points, when the size of the clusters virtually saturate under the given temperature and the space constraint by the simulation box. The parameters *R*<sup>1</sup> and *R*2, on the other hand, show the reflected s-shaped decreasing curves before entering the same stalled growth regimes. The reflections in the curves of *R*<sup>2</sup> are so strong that overall the shapes are close to the letter "z".

The parameter *R*<sup>0</sup> reflects the increasing contact between the pairs of monomers in different chains during the relaxation processes, over which the line geometry of the chains governs the values. The parameters *R*<sup>1</sup> and *R*2, on the other hand, contain information about the clustering of higher order structures. The three parameters meet at *R*<sup>0</sup> = *R*<sup>1</sup> = *R*<sup>2</sup> = 1, corresponding to a spatially uniformly distributed configuration. In simulation, it can be taken as a criteria to find the time spot of initiation of an aggregation process. By the same token, the reflection points in *R*<sup>1</sup> and *R*<sup>2</sup> can be identified to divide the processes into stages of growth.

Figure 4 shows the time evolution of a typical process of aggregation in a pure homopolymer system (I-4). The temperature of the system is controlled with a target temperature *T*∗ <sup>0</sup> = 6, starting from *t* = 0. We can identify the presence of multi-stage plateaus in the subsequent process.

To refine an aggregation process into a number of stages of relaxation, we identify each time spot at which the value of one of the parameters *R*0, *R*<sup>1</sup> or *R*2, starts to fall rapidly at the end of each plateau. In general, while the time spot of the triggered rapid decay in *R*<sup>0</sup> may not be distinguishable from the initiation time (at which *R*0=*R*1=*R*2=1), the onset of the fast falling in the value of *R*<sup>1</sup> and, especially, that of *R*<sup>2</sup> can easily be identified. The latter signals the growth transverse to the backbone. The chains in the system I-4 described in Fig. 4 have the stiffest springs among the spring bonded systems in set I (Table 1). The clustering process would become stuck and the system would need to wait until the initiation of the next major growth before there is another rapid fall in the value of *R*2. In Fig. 4(b), the fast decaying for both *R*<sup>0</sup> and *R*<sup>1</sup> start at *t* = 0. Such time spots for *R*<sup>2</sup> can be identified (marked by dashed lines) for three sections of plateaus, at *t* = 36.75, 91.75 and 133.75, respectively. The snapshots at these time spots, together with those at the beginning and the ending of the time span plotted in Fig. 4, are shown in Fig. 5.

The formation of clusters can be identified by the trend to have the lowered total interchain (L-J) potential energy *U*interchain and a decreasing virial *υ*interchain (see Ref. [7]). Since both *υ*interchain and *U*interchain follow the same trend [6, 7] and the latter is subject to much less fluctuations than the former is, we adapt here *U*interchain as the major quantity to monitor the

The formation of clusters can be identified by the tendency to have lowered total interchain (L-J) potential energy *U*interchain. Figure 3 shows the parameters *R*0, *R*<sup>1</sup> and *R*<sup>2</sup> decrease as the functions of *U*interchain, over the dynamic relaxation processes for the nine systems with soft n.n. bonding (*k*nn < ∞), for either pure homopolymer chains (*NF*=0 and *K*<sup>b</sup> = *K*<sup>t</sup> = 0, i.e. systems I-1, I-2, I-3 and I-4) or mixed heterogeneous chains with backbone anisotropy (*NF* > 0 and *K*<sup>b</sup> > 0 or *K*<sup>t</sup> > 0, i.e. systems II-0, II-1, II-2, II-3 and II-4) listed in Table 1. The data are collected in equal time intervals [6, 7]. The increasing of the quantity '−*U*inter−chain', which is used as the variable for the horizontal axis, corresponds to the direction for the growth of clusters. The parameter *R*<sup>0</sup> decreases, more or less, smoothly and is concave upward for all systems before reaching the tail regimes with densely accumulated data points, when the size of the clusters virtually saturate under the given temperature and the space constraint by the simulation box. The parameters *R*<sup>1</sup> and *R*2, on the other hand, show the reflected s-shaped decreasing curves before entering the same stalled growth regimes. The reflections in the

The parameter *R*<sup>0</sup> reflects the increasing contact between the pairs of monomers in different chains during the relaxation processes, over which the line geometry of the chains governs the values. The parameters *R*<sup>1</sup> and *R*2, on the other hand, contain information about the clustering of higher order structures. The three parameters meet at *R*<sup>0</sup> = *R*<sup>1</sup> = *R*<sup>2</sup> = 1, corresponding to a spatially uniformly distributed configuration. In simulation, it can be taken as a criteria to find the time spot of initiation of an aggregation process. By the same token, the reflection

Figure 4 shows the time evolution of a typical process of aggregation in a pure homopolymer system (I-4). The temperature of the system is controlled with a target temperature *T*∗

starting from *t* = 0. We can identify the presence of multi-stage plateaus in the subsequent

To refine an aggregation process into a number of stages of relaxation, we identify each time spot at which the value of one of the parameters *R*0, *R*<sup>1</sup> or *R*2, starts to fall rapidly at the end of each plateau. In general, while the time spot of the triggered rapid decay in *R*<sup>0</sup> may not be distinguishable from the initiation time (at which *R*0=*R*1=*R*2=1), the onset of the fast falling in the value of *R*<sup>1</sup> and, especially, that of *R*<sup>2</sup> can easily be identified. The latter signals the growth transverse to the backbone. The chains in the system I-4 described in Fig. 4 have the stiffest springs among the spring bonded systems in set I (Table 1). The clustering process would become stuck and the system would need to wait until the initiation of the next major growth before there is another rapid fall in the value of *R*2. In Fig. 4(b), the fast decaying for both *R*<sup>0</sup> and *R*<sup>1</sup> start at *t* = 0. Such time spots for *R*<sup>2</sup> can be identified (marked by dashed lines) for three sections of plateaus, at *t* = 36.75, 91.75 and 133.75, respectively. The snapshots at these time spots, together with those at the beginning and the ending of the time span plotted in

<sup>0</sup> = 6,

curves of *R*<sup>2</sup> are so strong that overall the shapes are close to the letter "z".

points in *R*<sup>1</sup> and *R*<sup>2</sup> can be identified to divide the processes into stages of growth.

process.

process.

Fig. 4, are shown in Fig. 5.

**3. Aggregation of polymer chains 3.1 Stages of the clustering process**

#### **3.2 Backbone anisotropy and spring stiffness**

The process for a mixed system (system II-4) of polymer chains, with perturbed backbone anisotropy (*K*<sup>b</sup> = *K*<sup>t</sup> = 10−4) and impurity monomers, shows similar features ( Fig. 6) as those found in system I-4 (Fig. 4). While the chains possess the same strengths in the monomer-monomer bonding springs as those in system I-4, the time evolution in Fig. 6(b) seems to show less tendency to become stuck over the relaxation process, as compared with its counterpart in Fig. 4(b). The corresponding snapshots of those time spots at the onsets of the fast decaying in *R*0, *R*<sup>1</sup> or *R*<sup>2</sup> are shown in Fig. 7.

A quantitative comparison between the configurations at the onsets of the rapid decaying in *R*<sup>2</sup> for system I-4 (Fig. 5(c)) and II-4 (Fig. 7(e)) is carried out by plotting the full monomer-monomer pair distribution functions *g*(*r*) in Fig. 8. To underscore the effect of the strength of the springs, we also put the data for the configuration at onset of *R*<sup>2</sup> for system II-0, which is a mixture system having the softest spring bonding in set II listed in Table 1. Each function *g*(*r*) shows a power-law like decaying curve in increasing *r* (see the log-log plot in the inset of Fig. 8), until reaching a minimum, which corresponds to the void regions next to the crowded domains of clusters. The thickness of the segments of the aggregated bundles can be estimated by locating *r* for such a minimum in *g*(*r*). While the two mixed systems (II-0 and II-4, marked by "A" and "B" respectively) have about the same bundle thickness, the packing in each system has its dependence on the strength of the monomer-monomer connecting springs. The plot shows the packing for chains in system II-0 has a crossover (the curve in the inset of Fig. 8), in contrast to the single power-law like behavior for those connected by strong springs (system II-4). The comparison between the data for system I-4 and II-4 indicates that the bundles formed in the pure polymer system (I-4, marked by "C")) is more compact than its counterpart in the mixed system (II-4, marked by "B"). The onset of the (first) fast growing stage requires larger size subunits in the mixed system with backbone anisotropy than those in the isotropic pure system.

Figure 9 shows the *Ri* parameters versus *U*interchain for the three cases with monomers connected by rigid bonds, *k*nn = ∞ (system I-∞, system II-∞ and system III), listed in Table 1. The curves share the same qualitative features as those in Fig. 3. A quantitative comparison of the curves for *R*<sup>0</sup> shows the larger degree of curving for systems with rigid bonding than the cases with soft bonding, as a result of less flexible patching among packed chains.

#### **3.3 Chain length dependence**

For the two pure homopolymer systems (*NF* = 0) with isotropic (*K*<sup>b</sup> = *K*<sup>t</sup> = 0) rigidly bonded (*k*nn = ∞) chains (systems I-A, *n* = 50 and I-B *n* = 25, in Table 1), each containing the same total number of monomers (*nNp*=4000), there are two major aggregated clusters in each case. The clusters have bundled structures and are tilted from the axes of simulation boxes, which suggest the independence of the boundary conditions in the formation of bundles. On a close examination of the snapshots of the structures (Fig. 10), we find that the effects of close sphere-packing lead to the tendency to deviate from a parallel bundled structure and lead to the formation of spiral-like local ordering (left plot of Fig. 10). Such orderings cannot, however, extend to the whole chain as the length of the chain becomes larger (middle and right plots of Fig. 10). In Fig. 11, we compare the monomer-monomer pair distributions among the three cases *n* = 25, *n* = 50 and *n* = 100. The humps behind the sharp, delta-function-like peaks, signal the ordering of sphere-packing. They are apparently larger for *n* = 25 and *n* = 50 than those for *n* = 100. The new ordering, other than that imposed by chain connectivity, is stronger for short chains.

Fig. 8. (Color online) Monomer-monomer pair distribution functions for configurations in the systems II-0 (dots labeled as "A"), II-4 (solid line labeled as "B") and I-4 (solid line labeled as "C") at the onset time for fast falling *R*2. The inset shows the data in log-log plot (adapted from Ref. [6]).

Fig. 9. Same as Fig. 3 for systems of rigidly-bonded chains.

An examination on the time evolutions (Fig 12 and Fig. 13) in these two systems, I-A and I-B, show also corresponding dynamic effects. Since the parameters *R*0, *R*<sup>1</sup> and *R*<sup>2</sup> feature the line-like, plane-like, and volume-like packed spatial structures, respectively, the returning of *R*<sup>2</sup> to a larger value in the later stage of evolution indicates the prevailing of three dimensional packing. The latter trend results in the larger values in volume. Figure 14 shows the *R*<sup>0</sup> and *R*<sup>1</sup> for *n* = 25 and *n* = 50 follow the decreasing trend in the growth processes, corresponding to reducing *<sup>U</sup>*inter−chain (or increasing −*U*inter−chain) which is the same as those for systems of longer chains (the cases of *n* = 100 in Figs. 3 and 9), in contrast to the returning of *R*<sup>2</sup> to a larger value in the later stage. The latter trend is stronger for the shorter chains. The result strongly suggests the relevance of chain length on the structural formation in aggregation phenomena.

Fig. 8. (Color online) Monomer-monomer pair distribution functions for configurations in the systems II-0 (dots labeled as "A"), II-4 (solid line labeled as "B") and I-4 (solid line labeled as "C") at the onset time for fast falling *R*2. The inset shows the data in log-log plot (adapted

An examination on the time evolutions (Fig 12 and Fig. 13) in these two systems, I-A and I-B, show also corresponding dynamic effects. Since the parameters *R*0, *R*<sup>1</sup> and *R*<sup>2</sup> feature the line-like, plane-like, and volume-like packed spatial structures, respectively, the returning of *R*<sup>2</sup> to a larger value in the later stage of evolution indicates the prevailing of three dimensional packing. The latter trend results in the larger values in volume. Figure 14 shows the *R*<sup>0</sup> and *R*<sup>1</sup> for *n* = 25 and *n* = 50 follow the decreasing trend in the growth processes, corresponding to reducing *<sup>U</sup>*inter−chain (or increasing −*U*inter−chain) which is the same as those for systems of longer chains (the cases of *n* = 100 in Figs. 3 and 9), in contrast to the returning of *R*<sup>2</sup> to a larger value in the later stage. The latter trend is stronger for the shorter chains. The result strongly suggests the relevance of chain length on the structural formation in aggregation

Fig. 9. Same as Fig. 3 for systems of rigidly-bonded chains.

from Ref. [6]).

phenomena.

Fig. 10. Snapshots of the largest aggregated clusters for (left plot) *n* = 25, consisting of 81 chains, for (middle plot) *n* = 50), consisting of 51 chains, and for (right plot) *n* = 100, consisting of all 40 chains of the systems (adapted from Ref. [7]).

Fig. 11. (Color online) Comparison of monomer-monomer pair distributions of the aggregated configuration for *n* = 25 (blue line), *n* = 50 (red line) and *n* = 100 (black line) (adapted from Ref. [7]).

#### **4. Growth hindered by backbone anisotropy**

We have compared the aggregated structures between the systems with zero (system I-4) and with non-zero, but tiny angle potentials (system II-4),. While they both possess bundled clusters (Fig. 1), the non-zero tiny angle potentials break the symmetry to have a branched structure (right plot of Fig. 1). One more sensible question to solve is to ask if the bundled cluster can survive under even larger strengths in angle potentials. The simulations with chains having *K*<sup>b</sup> = *K*<sup>t</sup> = 0.1 (set IV in Table 1), which is 1000 times larger than those of the systems in set II have shown no cluster domains formed over time spans equivalent to or larger than those considered for systems in set II [6, 7]. The systems are brought to states at lower temperatures in a stepwise manner, with decreasing target temperature *T*∗ <sup>0</sup> at 8, 6 and 3, respectively, in three stages. We do observe the formation of local parallel conformations. But this local ordering can hardly extend to either the larger segment or to the gathering of more chains. The results indicate that the effect of frustration between the straight bundled packing and the local curvature imposed by the angle potentials prevent the aggregation. Such an observation is suggestive in considering more complicated material systems.

Fig. 12. Same as Fig. 4, for a simulated system with rigidly bonded homogeneous chains (*k*nn = ∞) of length *n* = 50 each (adapted from Ref. [7]).

Fig. 13. Same as Fig. 12 for a simulated system with rigidly bonded homogeneous chains (*k*nn = ∞) of length *n* = 25 each (adapted from Ref. [7]).

The typical time evolution of such processes is shown in Fig. 15, for system IV-4. We observe that the values of the parameter *R*<sup>0</sup> fluctuate in the larger amplitude as compared with those for *R*<sup>1</sup> and *R*2, indicating the occurrence of local alignment between chains fails to survive for triggering the formation of stable structures. In lowering the temperature, the fluctuations change to have the larger amplitude and the longer periods, which seems closer to a stage with monotonic decaying *R*0. Entering the regime of *T*∗ *<sup>P</sup>* ≈ 7, the stronger tendency of alignment between chains signalled by the larger decaying in the value of *R*<sup>0</sup> still fails to prevent the reverting of the trend, that the enhanced changes in *R*<sup>0</sup> can only stay with the status as the (larger amplitude) fluctuations. The local mechanical conflicts caused by the angle potentials hinder the pile-up of the packed local segments to form stable higher order structures. We

Fig. 12. Same as Fig. 4, for a simulated system with rigidly bonded homogeneous chains

Fig. 13. Same as Fig. 12 for a simulated system with rigidly bonded homogeneous chains

The typical time evolution of such processes is shown in Fig. 15, for system IV-4. We observe that the values of the parameter *R*<sup>0</sup> fluctuate in the larger amplitude as compared with those for *R*<sup>1</sup> and *R*2, indicating the occurrence of local alignment between chains fails to survive for triggering the formation of stable structures. In lowering the temperature, the fluctuations change to have the larger amplitude and the longer periods, which seems closer to a stage with

between chains signalled by the larger decaying in the value of *R*<sup>0</sup> still fails to prevent the reverting of the trend, that the enhanced changes in *R*<sup>0</sup> can only stay with the status as the (larger amplitude) fluctuations. The local mechanical conflicts caused by the angle potentials hinder the pile-up of the packed local segments to form stable higher order structures. We

*<sup>P</sup>* ≈ 7, the stronger tendency of alignment

(*k*nn = ∞) of length *n* = 50 each (adapted from Ref. [7]).

(*k*nn = ∞) of length *n* = 25 each (adapted from Ref. [7]).

monotonic decaying *R*0. Entering the regime of *T*∗

Fig. 14. Same as Fig. 3 for systems with chains of lengths *n* = 50 (circle) and *n* = 25 (diamond), respectively.

Fig. 15. (Color online) Same as Fig. 4, for the mixed system of polymer and fluid, system IV-4, where *Kb* = *Kt* = 0.1. (adapted from Ref. [6]).)

conclude that the values of *K*<sup>b</sup> and *K*<sup>t</sup> have to be small for the presence of aggregation process for the polymer chains.

#### **5. Perspectives**

The study of aggregations of peptides or proteins and their roles in the occurrence of neurodegenerative diseases, such as Alzheimer's disease [20–22], and the transmissible Creutzfeldt-Jakob disease [23], has become an important cross-disciplined issue for scientists. The important factors for polymer aggregation problems revealed by our analysis, such as the backbone anisotropy and the chain length are suggestive. The suspected origin of Alzheimer's disease at the molecular level is related to the aggregation of proteins A*β*<sup>40</sup> and A*β*<sup>42</sup> [20–22] consisting of 40 and 42 residues, respectively. The chain of A*β*<sup>40</sup> or A*β*<sup>42</sup> forms two *β*-strands and such *β*-strands form *β*-sheets [21, 22]. A scenario that describes their aggregation is to consider each *β*-sheet as an effective polymer chain and has a small degree of effective backbone anisotropy. Thus such effective polymer chains has the trend to aggregate into parallel-domain conformation. It will be interesting in the future work to include more refined backbone anisotropy [24, 25] and quantify its effect on the aggregation dynamics.

#### **6. Acknowledgement**

This work was supported by Grants NSC100-2112-M-004-001, NSC100-2112-M-001-003-MY2 and NCTS (North).

#### **7. References**


## **Molecular Dynamics Simulation of Permeation in Polymers**

Hossein Eslami and Nargess Mehdipour *Persian Gulf University Iran* 

### **1. Introduction**

16 Will-be-set-by-IN-TECH

60 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

and such *β*-strands form *β*-sheets [21, 22]. A scenario that describes their aggregation is to consider each *β*-sheet as an effective polymer chain and has a small degree of effective backbone anisotropy. Thus such effective polymer chains has the trend to aggregate into parallel-domain conformation. It will be interesting in the future work to include more refined

This work was supported by Grants NSC100-2112-M-004-001, NSC100-2112-M-001-003-MY2

[1] T. P. J. Knowles, J. F. Smith, A. Craig, C. M. Dobson and M. E. Welland: *Phys. Rev. Lett.*

[2] T. P. J. Knowles, A. W. Fitzpatrick, S. Meehan, H. R. Mott, M. Vendruscolo, C. M. Dobson

[4] B. Urbanic, L. Cruz, S. Yun, S. V. Buldyrev, G. Bitan, D. B. Teplow and H. E. Stanley:

[8] R. V. Pappu, X. Wang, A. Vitalis and S. L. Crick: *Arch. Biochem. Biophys.* 469 132 (2008). [9] M. Doi and S. F. Edwards: *The Theory of Polymer Dynamics* (Oxford University Press,

[10] R. B. Bird, C. F. Curtiss, R. C. Armstrong and O. Hassager: *Dynamics of Polymer Liquids*

[15] D. D. Voet and J. G. Voet: *BioChemistry* (John Weily and Sons, 1995) 2nd ed. Chap. 7.

[20] A. Goate, M. C. Chartierharlin, M. Mullan, J. Brown, F. Crawford, L. Fidani, L. Giuffra, A. Haynes, N. Irving, L. James, R. Mant, P. Newton, K. Rooke, P. Roques, C. Talbot, M. Pericakvance, A. Roses, R. Williamson, M. Rossor, M. Owen and J. Hardy: *Nature* 349

[21] A. T. Petkova, Y. Ishii, J. J. Balbach, O. N. Antzutkin, R. D. Leapman, F. Delaglio and R

[22] T. Luhrs, C. Ritter, M. Adrian, D. Riek-Loher, B. Bohrmann, H. Doeli, D. Schubert and

[23] N. J. Cobb, F. D. Sonnichsen, H. Mchaourab and W. K. Surewicz: *Proc. Natl. Acad. Sci.*

[3] H. D. Nguyen and C. K. Hall: *Proc. Natl. Acad. Sci.* U. S. A. 101 (2004)16180 .

backbone anisotropy [24, 25] and quantify its effect on the aggregation dynamics.

**6. Acknowledgement**

96 (2006) 238301.

New York, 1986).

(1991) 704.

USA 104 (2007) 18946.

and M. E. Welland: *Science* 318 (2007)1900.

*Proc. Natl. Acad. Sci.* U. S. A. 101 (2004) 17345.

Vol.2, 2nd ed. (John Wiley and Sons 1987). [11] H. C. Andersen: *J. Comput. Phys.* 52, 24 (1983).

[19] L. V. Woodcock: *Chem. Phys. Lett.* 10 (1971) 257.

[12] J. T. Padding, W. J. Briels: *J. Chem. Phys.* 114 (2001) 8685. [13] D. Steele: *J. Chem. Soc. Faraday Trans.* 2 81 (1985) 1077.

[16] W.-J. Ma and C.-K. Hu: *J. Phys. Soc. Jpn.* 79 (2010) 024005. [17] W.-J. Ma and C.-K. Hu: *J. Phys. Soc. Jpn.* 79 (2010) 024006.

Tycko: *Proc. Natl. Acad. Sci.* USA 99 (2002) 16742.

R. Riek: *Proc. Natl. Acad. Sci.* USA 102 (2005) 17342.

[24] R. Pellarin and A. Caflisch: *J. Mol. Biol.* 360 (2006) 882.

[25] N. Y. Chen, Z. Y. Su, and C. Y. Mou: *Phys. Rev. Lett.* 96 (2006) 078103.

[14] D. Brown and J. H. R. Clarke: *Comput. Phys. Commun.* 62, (1991) 360.

[6] W. J. Ma and C.-K. Hu: *J. Phys. Soc. Japan* 79 (2010) 100402. [7] W. J. Ma and C.-K. Hu: *J. Phys. Soc. Japan* 79 (2010) 54001.

[5] C. M. Dobson: *Nature* 426 (2003)884.

[18] C. Tsallis: *J. Stat. Phys.* 52 (1988) 479.

and NCTS (North).

**7. References**

Sorption and diffusion of small molecules through polymers is a topic of broad range of applications in many industrially important phenomena. There is a rapidly growing demand for polymers of specified permeabilities, such as selective memberanes for separation technologies, barrier membranes for packaging applications, foaming, and plasticization. Polymers with a high degree of permeability and permselectivity have been widely used for gas or liquid separation systems based on membranes (Kesting & Fritzsche, 1993), while those with a low degree of permeability have been used in barrier packaging films as containers (Koros, 1990). Therefore, understanding the underlying mechanism of solubility and/or diffusion of peneterants in macromolecular substances is very useful in obtaining a clear picture of the molecular level mechanism of polymer permeability and in design of new membranes.

The first study of gas permeation through a polymer goes back to 1829 by Thomas Graham (Stannett, 1978). According to Graham's postulate, the permeation process involves the dissolution of penetrant followed by transmission of the dissolved species through the membrane. In the late 1870's, Stefan and Exner demonstrated that gas permeation through a membrane was proportional to the product of solubility and diffusion coefficients. Based on the findings of Stefan and Exner, von Wroblewski (1879) constructed a quantitative solution to the Graham's solution-diffusion model. The dissolution of gas was based on Henry's law of solubility, where the concentration of the gas in the membrane was directly proportional to the applied gas pressure. Later Daynes (1920) showed that it was impossible to evaluate both diffusion and solubility coefficients by steady-state permeation experiments. He presented a mathematical solution, using Fick's second law of diffusion, for calculating the diffusion coefficient, which was assumed to be independent of concentration.

Microscopic theories of the diffusion and permeation of penetrants in polymers, on the other hand, were developed during 1950-1970's. In these theories, to explain the mechanism behind the transport of gas molecules through the free volume present in the polymer matrix, the gas-polymer system was defined in terms of liquid molecules traveling through a liquid membrane. The Cohen and Turnbull (1959) theory considered transport through liquid molecules as hard spheres. In this model the gas diffusion in polymers occurs through the redistribution of free volume. Another theory, proposed by Dibenedetto (1963), used the same concept as that of Cohen and Turnbull (1959) regarding free volume distribution, but with a different chain packing theory at the molecular level. An elaborate theory was developed by Pace and Datyner (1979), assuming that the penetrant molecule is moving along the polymer chain bundles and being stopped only by chain entanglements or a crystallite. The penetrnat molecules were further assumed to jump into the adjacent bundles, similar to the Cohen and Turnbull's model. Such a jump event was considered to be the rate controlling step, with the diffusion along the bundle being three times faster than the perpendicular jump of the molecule. None of these theories, however, provides a molecular lever understanding of the permeability in polymers.

Significant advances in the understanding of gas permeation in polymers have been made only in the last few decades with the advent of powerful computers. Monte Carlo (MC) and molecular dynamics (MD) simulations on small-molecule permeation of amorphous polymers have become feasible in recent decades. The simulations now cover a range of different polymers with varieties of chemical complexities ranging from flexible polymers with simple chemical structure, like polyethylene (Boyd & Pant, 1991), to stiff polymers with detailed chemical structure, like polyamides (Eslami & Mehdipour, 2011). In the following we give a detailed discussion of the calculation methods as well as the polymeric samples employed to study the sorption/diffusion mechanism. The main difficulty in the calculation of gas permeabilities stems from the calculation of solubility coefficients. In fact the calculation of gas solubilities necessitates the condition of equilibrium between the permeant (sorbate) in the gas phase and that dissolved in polymer. The equality of the temperature, pressure, and chemical potentials of all species is the necessary condition to establish such an equilibrium situation. The chemical potential is, however, coupled to the number of particles and cannot be easily calculated using molecular simulation methods. Therefore, while there exist excellent reviews on the pearmation of gases in polymers, with the main emphasis on the diffusion process (Müller-Plathe, 1994), in this review we explain the polymer permeation putting more emphasis on the solubility process.

#### **2. Sorption of gases in polymers**

Historically, the dual-mode sorption theory (Stannett et al., 1978; Frederickson and Helfand, 1985) describes the concentration *C* (solubility) of the gas inside a polymer at equilibrium with the gas at a partial pressure *p*, i.e.,

$$\mathbf{C} = k\_H p + \mathbf{C}\_{\infty} \frac{bp}{1 + bp} \tag{1}$$

where *kH* is Henry's law solubility coefficient, *C* is the saturation concentration of the gas, and *b* is the affinity coefficient. This model assumes that there are two distinct modes by which a glassy polymer can sorb gas molecules; Henry's law and a Langmuir mechanism which corresponds to the sorption of gases into specific sorption sites in the polymer. Henry's constant has the same physical meaning for glassy polymers as it does for rubbery polymers and liquids, whereas the Langmuir-type term is believed to account for gas sorption into interstitial spaces and microvoids, which are consequences of local heterogeneities and are intimately related to the slow relaxation processes associated with the glassy state of the polymer. Local equilibrium is assumed between the two modes. Figure 1 shows a schematic representation of Henry, Langmuir, and dualmode sorptions.

Fig. 1. Schematic representation of Henry, Langmuir, and dual-mode sorptions.

In the low pressure region, Eq. 1 provides the following linear relationship against the pressure:

$$\mathbf{C} = \left(k\_H + \mathbf{C}\_\alpha b\right) P = \mathbf{S}\_0 P \tag{2}$$

where *S0* is called the apparent solubility coefficient in the zero-pressure limit in glassy polymers.

The determination of the free energy change in sorption of gases in polymer is a problem of primary importance, since the free energy is the main thermodynamic parameter governing the equilibrium between the gas in the pure state and the dissolved one in the polymer. For a system of *N* particles located at *r1, r2, ...rN*, the statistical mechanical expression for the Helmholtz free energy, *A*, reads as (McQuarrie, 1976):

$$A = -k\_B T \ln(Q) = \frac{1}{N! \Lambda^{3N}} \int\_0^\infty \dots \int\_0^\infty \exp\left(-lI\_N\left(r^N\right)/k\_B T\right) \, dr^N \tag{3}$$

with

62 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

distribution, but with a different chain packing theory at the molecular level. An elaborate theory was developed by Pace and Datyner (1979), assuming that the penetrant molecule is moving along the polymer chain bundles and being stopped only by chain entanglements or a crystallite. The penetrnat molecules were further assumed to jump into the adjacent bundles, similar to the Cohen and Turnbull's model. Such a jump event was considered to be the rate controlling step, with the diffusion along the bundle being three times faster than the perpendicular jump of the molecule. None of these theories, however, provides a

Significant advances in the understanding of gas permeation in polymers have been made only in the last few decades with the advent of powerful computers. Monte Carlo (MC) and molecular dynamics (MD) simulations on small-molecule permeation of amorphous polymers have become feasible in recent decades. The simulations now cover a range of different polymers with varieties of chemical complexities ranging from flexible polymers with simple chemical structure, like polyethylene (Boyd & Pant, 1991), to stiff polymers with detailed chemical structure, like polyamides (Eslami & Mehdipour, 2011). In the following we give a detailed discussion of the calculation methods as well as the polymeric samples employed to study the sorption/diffusion mechanism. The main difficulty in the calculation of gas permeabilities stems from the calculation of solubility coefficients. In fact the calculation of gas solubilities necessitates the condition of equilibrium between the permeant (sorbate) in the gas phase and that dissolved in polymer. The equality of the temperature, pressure, and chemical potentials of all species is the necessary condition to establish such an equilibrium situation. The chemical potential is, however, coupled to the number of particles and cannot be easily calculated using molecular simulation methods. Therefore, while there exist excellent reviews on the pearmation of gases in polymers, with the main emphasis on the diffusion process (Müller-Plathe, 1994), in this review we explain the polymer permeation putting more

Historically, the dual-mode sorption theory (Stannett et al., 1978; Frederickson and Helfand, 1985) describes the concentration *C* (solubility) of the gas inside a polymer at equilibrium

> 1 *<sup>H</sup> bp C kp C bp*

where *kH* is Henry's law solubility coefficient, *C* is the saturation concentration of the gas, and *b* is the affinity coefficient. This model assumes that there are two distinct modes by which a glassy polymer can sorb gas molecules; Henry's law and a Langmuir mechanism which corresponds to the sorption of gases into specific sorption sites in the polymer. Henry's constant has the same physical meaning for glassy polymers as it does for rubbery polymers and liquids, whereas the Langmuir-type term is believed to account for gas sorption into interstitial spaces and microvoids, which are consequences of local heterogeneities and are intimately related to the slow relaxation processes associated with the glassy state of the polymer. Local equilibrium is assumed between the

(1)

molecular lever understanding of the permeability in polymers.

emphasis on the solubility process.

**2. Sorption of gases in polymers** 

with the gas at a partial pressure *p*, i.e.,

$$
\Lambda = \left(\frac{h}{2\pi mk\_BT}\right)^{1/2} \tag{4}
$$

Here, *Q* is the canonical partition function, *UN* is the potential energy of the system, *rN* stands for the whole set of coordinates, *r1, r2, …rN*, is the de Broglie wavelength, *h* is the Planck's constant, *m* is the molecular mass, *kB* is the Boltzmann's constant, and *T* is the temperature. Assuming pair-wise additiveity of the potential energy between particles, we have:

$$\mathcal{LI}\_N = \sum\_{i}^{N} \sum\_{j>i}^{N} \mu\_{ij}(r\_{ij}) \tag{5}$$

where *uij* is the pair potential interacting between particles *i* and *j* and *rij* is the interparticle distance. The expression for the chemical potential is obtained by taking the logarithm of the ratios of partition functions for a system composed of *N* particles, with the sorbate (penetrant) density of *<sup>s</sup>*, and a system composed of *N+1* particles, which is obtained by adding one sorbate particle to the previous *N*-particle system, i.e,

$$\mu\_s = A\_{N+1} - A\_N = k\_B T \ln\left(\rho\_s \Lambda^3\right) - k\_B T \ln\left(\left\{\exp\left(-\Delta L \mid k\_B T\right)\right\}\right) \tag{6}$$

with

$$
\Delta \mathcal{U} = \mathcal{U}\_{N+1} - \mathcal{U}\_N \tag{7}
$$

Here *s* and *<sup>s</sup>* are the chemical potential and the density of sorbate, respectively, and the brackets indicate the ensemble average. Comparing *<sup>s</sup>* with the chemical potential of the ideal gas at the same state, i.e.,

$$
\mu\_s = \mu\_s^{\epsilon x} + k\_B T \ln \left( \rho\_s \Lambda^3 \right) \tag{8}
$$

various contributions to the chemical potential can be interpreted according to Ben-Naim (1987). In Eq. (8) *ex <sup>s</sup>* is called the excess chemical potential of sorbate. The second term on the right hand side of Eq. (8) is the chemical potential of the ideal gas. Now, the chemical potential of the sorbate is split into two terms; a term arising from putting the sorbate molecule at a fixed position in the polymer matrix, *ex <sup>s</sup>* , and a term arising from releasing the constraint, i.e., letting the solute molecule move freely, which results in the contribution <sup>3</sup> *k TB s* ln to the chemical potential. In fact this term is the contribution of translational motion to the chemical potential.

Looking at the sorption process, we consider the pure sorbate gas *s* in equilibrium with the dissolved gas in the polymer. Again we connect the thermodynamic states of the system to their corresponding ideal gas states (at the same temperature, and the same density of sorbate *s*), as it is depicted in Scheme 1.

Scheme 1. Description of equilibrium between the gas *s* in the gaseous and polymeric phases.

Here, and indicate the densities of gas *s* in the polymeric and gaseous phases, respectively. The condition of equilibrium for gas *s*, between the gas phase and the polymer phase is the equality of the chemical potentials, i.e.,

*N N N ij ij i ji U ur* 

where *uij* is the pair potential interacting between particles *i* and *j* and *rij* is the interparticle distance. The expression for the chemical potential is obtained by taking the logarithm of the ratios of partition functions for a system composed of *N* particles, with the sorbate

> 

(penetrant) density of

with

Here *s* and  ideal gas at the same state, i.e.,

motion to the chemical potential.

ideal gas *s* (*T*,

sorbate *s*), as it is depicted in Scheme 1.

gas *s* sorbed in polymer (*T*,

phase is the equality of the chemical potentials, i.e.,

(1987). In Eq. (8) *ex*

 <sup>3</sup> *k TB s* ln 

phases. Here, and  1

brackets indicate the ensemble average. Comparing

molecule at a fixed position in the polymer matrix, *ex*

adding one sorbate particle to the previous *N*-particle system, i.e,

( )

*<sup>s</sup>*, and a system composed of *N+1* particles, which is obtained by

<sup>3</sup>

*s N NB s B A A kT kT U kT* ln ln exp( / *<sup>B</sup>* (6)

*<sup>s</sup>* are the chemical potential and the density of sorbate, respectively, and the

 <sup>3</sup> ln *ex ss B s*

various contributions to the chemical potential can be interpreted according to Ben-Naim

the right hand side of Eq. (8) is the chemical potential of the ideal gas. Now, the chemical potential of the sorbate is split into two terms; a term arising from putting the sorbate

the constraint, i.e., letting the solute molecule move freely, which results in the contribution

Looking at the sorption process, we consider the pure sorbate gas *s* in equilibrium with the dissolved gas in the polymer. Again we connect the thermodynamic states of the system to their corresponding ideal gas states (at the same temperature, and the same density of

Scheme 1. Description of equilibrium between the gas *s* in the gaseous and polymeric

respectively. The condition of equilibrium for gas *s*, between the gas phase and the polymer

to the chemical potential. In fact this term is the contribution of translational

) ideal gas *s* (*T*,

) gas *s* (*T*,

indicate the densities of gas *s* in the polymeric and gaseous phases,

  *<sup>s</sup>* is called the excess chemical potential of sorbate. The second term on

(5)

*UU U N N* 1 (7)

*k T* (8)

*<sup>s</sup>* with the chemical potential of the

*<sup>s</sup>* , and a term arising from releasing

)

)

$$
\mu\_s(T,\rho) = \mu\_s'(T,\rho') \tag{9}
$$

where and *µ* are the chemical potentials of the sorbed gas in polymer and in gas phase, respectively. According to the thermodynamic cycle shown in Scheme 1, the equilibrium condition between the pure gas and the gas sorbed in the polymer is formulated as:

$$
\mu\_s^{\rm ex} = \mu\_s^{\prime \rm ex} - k\_B T \ln \frac{\rho}{\rho^{\prime}} \tag{10}
$$

Equation 10 is the main equation governing the equilibrium between the sorbed gas molecules in the polymer and those in the gaseous phase. According to this equation, the main parameter to be determined for the calculation of solubilities is the excess chemical potential. In the following section we study the computational approach, with the main emphasis on the MD simulation methods, for calculation of chemical potentials.

#### **3. Calculation of excess chemical potentials**

From the historical point of view, traditional approaches such as the equation-of-state models (Lacombe and Sanchez, 1976 ; Sanchez and Rodgers, 1990 ; von Solms et al., 2005) and the activity coefficient models (Ozkan and Teja, 2005) have been applied to calculate the phase equilibria and sorption of penetrant molecules in polymers. Molecular simulations are the other attractive method for this type of calculation. These methods do not invoke any approximations and predictions are based on well-defined molecular characteristics. In the following sections we describe the application of molecular simulation methods in the case of sorption of gases in polymers. As stated above, the main problem is to calculate the excess chemical potentials, which is given address to in the following sections.

#### **3.1 Widom's test particle insertion method**

The Widom's test particle insertion method (Widom, 1963) is an elegant method for calculating the excess chemical potentials. In this method a test particle is momentarily inserted at random points in the simulation box and the interaction energy between the test particle and the host particles are averaged to yield the excess chemical potential. Since the method relies on the statistical accuracy of sampling of configurations that permit the insertion of molecules with low values of the binding energies, its application is more feasible for small-sized solutes. In a dense fluid there is a small probability of finding a cavity to insert a particle in. Therefore, particle insertion becomes a rare event and long simulation times are needed to obtain a good statistics. Similarly, at high densities, the probability of particle deletion is also reduced. The deletion probability depends on the probability of generation of a high density configuration, from which to delete a particle. Due to several limitations, such as poor sampling at high densities or inadequacy of the method in the case of chemical potential calculation for polar systems, the utility of this method to many important applications is precluded. The problem of poor sampling has been solved to some extent by developing such techniques as the umbrella sampling (Shing and K. E. Gubbins, 1981) and map sampling (Deitrick et al., 1989), as the extensions of original Widom's method. The Widom method or its modifications can be used in both MC and MD simulations, but in this method the chemical potential will be calculated at the end of the simulation. This method has widely been applied in calculating the excess chemical potentials of small penetrant molecules in polymers including polypropylene (Müller-Plathe, 1991), polydimethyl siloxane (Sok et al., 1992), polyisobutylene (Müller-Plathe et al., 1993), polystyrene (Eslami & Müller-Plathe, 2007a), and polyimides (Neyertz & Brown, 2004; Pandiyan et al. 2010). However, the method is not applicable in the case of excess chemical potential calculation for polar penetrant in polar polymers. Moreover, the sampling probability gets poor in dense polymeric systems and in the case of big penetrants.

#### **3.2 Thermodynamic integration and fast growth thermodynamic integration methods**

Thermodynamic integration method is another useful method for the calculation of the excess chemical potentials. In this method, a parameter λ is used to couple the sorbate with the rest of the system and the excess chemical potential, or coupling work, is obtained by integrating along λ as:

$$
\mu\_s^{ex} = \int \left\langle \left\langle \frac{\partial \mathcal{U}}{\partial \mathcal{Z}} \right\rangle \right\rangle \, d\mathcal{Z} \tag{11}
$$

In this method, several simulations are required to obtain a series of points for the integration. To improve the sampling, recently Hess and van der Vegt (2008) have developed a fast-growth thermodynamic integration method. The approach is based on the non-equilibrium work theorem of Jarzynski (1997), which relates the work that is being performed on a system when going from an initial state to a final state along a coordinate λ with the free energy change. Through the applications of this method, polymer microstructure can be considered as a potential landscape as in Figure 2.

In this context, the free energy can be calculated by performing fast growth thermodynamic integration in order to calculate the work required to insert the sorbate in a solvent matrix from sampling the possible coordinates where interaction potentials determine the potential landscape. This method has been applied successfully to calculate the excess chemical potential of relatively large molecules in polymer matrices (Hess & van der Vegt 2008; Hess et al. 2008; Fritz et al., 2009).

#### **3.3 Open system simulations**

Since the chemical potential is conjugated with the number of particles, the statistical mechanical ensemble representative for its calculation is the grand canonical ensemble. There are several open-system simulation techniques in the literature, allowing to set the target chemical potential as an independent thermodynamic variable. These methods are preferential to the methods explained above for the calculation of phase equilibria, in that the chemical potential in each phase is set from the beginning as an independent variable, while in closed system simulations the chemical potential is calculated at the end of the run. In the following two subsections we describe the grand canonical ensemble simulations using MC and MD techniques.

method or its modifications can be used in both MC and MD simulations, but in this method the chemical potential will be calculated at the end of the simulation. This method has widely been applied in calculating the excess chemical potentials of small penetrant molecules in polymers including polypropylene (Müller-Plathe, 1991), polydimethyl siloxane (Sok et al., 1992), polyisobutylene (Müller-Plathe et al., 1993), polystyrene (Eslami & Müller-Plathe, 2007a), and polyimides (Neyertz & Brown, 2004; Pandiyan et al. 2010). However, the method is not applicable in the case of excess chemical potential calculation for polar penetrant in polar polymers. Moreover, the sampling probability gets poor in dense polymeric systems and in the case of big

**3.2 Thermodynamic integration and fast growth thermodynamic integration methods**  Thermodynamic integration method is another useful method for the calculation of the excess chemical potentials. In this method, a parameter λ is used to couple the sorbate with the rest of the system and the excess chemical potential, or coupling work, is obtained by

> 1 0

*d*

In this method, several simulations are required to obtain a series of points for the integration. To improve the sampling, recently Hess and van der Vegt (2008) have developed a fast-growth thermodynamic integration method. The approach is based on the non-equilibrium work theorem of Jarzynski (1997), which relates the work that is being performed on a system when going from an initial state to a final state along a coordinate λ with the free energy change. Through the applications of this method, polymer

In this context, the free energy can be calculated by performing fast growth thermodynamic integration in order to calculate the work required to insert the sorbate in a solvent matrix from sampling the possible coordinates where interaction potentials determine the potential landscape. This method has been applied successfully to calculate the excess chemical potential of relatively large molecules in polymer matrices (Hess & van der Vegt 2008; Hess

Since the chemical potential is conjugated with the number of particles, the statistical mechanical ensemble representative for its calculation is the grand canonical ensemble. There are several open-system simulation techniques in the literature, allowing to set the target chemical potential as an independent thermodynamic variable. These methods are preferential to the methods explained above for the calculation of phase equilibria, in that the chemical potential in each phase is set from the beginning as an independent variable, while in closed system simulations the chemical potential is calculated at the end of the run. In the following two subsections we describe the grand canonical ensemble simulations

*U*

(11)

*ex s*

microstructure can be considered as a potential landscape as in Figure 2.

penetrants.

integrating along λ as:

et al. 2008; Fritz et al., 2009).

**3.3 Open system simulations** 

using MC and MD techniques.

Fig. 2. A schematic one-dimensional representation of a free-energy landscape for an additive molecule (*x* is a component of its center of mass) in a polymer matrix as a function of the coupling parameter . During a fast-growth simulation, a molecule is forced along from state A to B, while sampling in *x*; three example paths are given. With fast-growth thermodynamic integration, it is possible to explore the multiple minima in the free energy landscape by multiple trajectories starting from a flat free energy landscape (a solute not interacting with the matrix). Figure is taken from Hess and van der Vegt (2008) with permission.

### **3.3.1 Grand canonical ensemble Monte Carlo simulation**

The first grand canonical simulation method (Norman & Filinov, 1969 ) was based on a MC simulation technique. This method was then developed by Adams (1974, 1975). In these methods full particles are inserted into or removed from the simulation box, in a Markov process with a probability, which depends on the target chemical potential. These kinds of movements are compatible with the stochastic nature of moves in MC technique. The main difficulty with the MC simulation methods in the grand canonical ensemble is their inefficiency when applied to the condensed phases. Besides, the insertion or deletion of full particles highly perturbs the simulation box, and it takes some time for the fluid to relax. These methods envisage similar sampling problems, as addressed above in the case of Widom's method, at high densities. However, the sampling problem has been overcome somewhat by searching in the simulation box for suitable cavities and biasing the insertion event (Escobedo & de Pablo, 1995; Shelley & Patey, 1995; Boda et al., 1997). Another problem with the MC simulation methods in the grand canonical ensemble is that these methods do not directly reflect the dynamic information about the system. Therefore, attention has been attracted toward development of the grand canonical ensemble simulations using MD technique.

#### **3.3.2 Grand canonical ensemble molecular dynamics simulation**

The discontinuity of the number of particles, introduced through particle insertion/deletion, makes it difficult to do MD simulation in the grand canonical ensemble. Therefore, in some of the existing MD-based simulations in the grand canonical ensemble, the sampling is performed through using the MC-based stochastic procedures (Boinepalli & Attard, 2003). Another alternative, is the introduction of the particle number as a continuous dynamic variable through the MD simulation (Eslami & Müller-Plathe, 2007b; Eslami et al., 2011). This approach is based on the so-called extended system method, in which a Hamiltonian including the kinetic and potential energy terms for the extension variable is constructed and the appropriate equations of motion are solved.

Here we assume that our physical system is composed of *N* real particles plus one "fractional" particle and is coupled to a heat reservoir and a particle reservoir. The fractional particle is a particle whose potential energy is weighted by a variable which varies from zero to one. The inclusion of the fractional particle in the system provides a variable number of particles in a dynamical way. The Hamiltonian of such a combined system; the real *N* particle system, the fractional particle, and the additional degrees of freedom, i.e., the material reservoir and the heat reservoir, is written as (Eslami & Müller-Plathe, 2007b):

$$H = \sum\_{i=1}^{N} \frac{p\_{x\_i}^2 + p\_{y\_i}^2 + p\_{z\_i}^2}{2m\_i s^2} + \sum\_{i=1}^{N-1} \sum\_{j>i}^{N} \mathcal{U}\_{\vec{\eta}} + \frac{p\_{x\_f}^2 + p\_{y\_f}^2 + p\_{z\_f}^2}{2m\_f s^2} + \sum\_{i=1}^{N} \mathcal{U}\_{\vec{\mathcal{H}}} + \frac{p\_{\vec{\lambda}}^2}{2\mathcal{W}\_{\vec{\lambda}}} + \mathcal{U}\_{\nu} + \frac{p\_s^2}{2\mathcal{W}\_s} + \mathcal{U}\_s \tag{12}$$

where *H* is the Hamiltonian, subscripts *i* and *j* refer to the real particles, subscript *f* refers to the fractional particle, *p*i is the momentum of the *i*th particle, *s* is the velocity scaling variable which couples the system to the thermostat, and is the particle-number scaling variable which couples the fractional particle to the rest of the system and varies from zero to one. In this equation, the first four terms on the right hand side are the kinetic energy of the real particles, potential energy of interaction between real particles, kinetic energy of the fractional particle, and the potential energy of interaction between real particles and the fractional particle. The number extension variable, , and the temperature extension variable, s, have their own kinetic and potential energy terms, shown as the last four terms on the right hand side of Eq. (12). In fact *U* is the work required to carry a fractional particle from a medium of zero potential to the medium of interest with the potential energy *U*; i.e., it is the so-called excess Helmholtz free energy.

Adopting one of the real particles (molecules) in the simulation box as the fractional particle, the fractional particle is grown to a full (real) particle or deleted from the system in a dynamical way, by solving the equations of motion. The relevant equation of motion, governing the dynamics of fractional particle (the magnitude of coupling to the ideal gas reservoir) is (Eslami & Müller-Plathe, 2007b):

$$\mathcal{W}\ddot{\mathcal{X}} = -\sum\_{i=1}^{N} \frac{\partial \mathcal{U}\_{\dot{\mathcal{Y}}}}{\partial \mathcal{X}} + \mu - \mu^{\dot{\mathcal{U}}} \tag{13}$$

where d *λ* =d2/dt2 and *id* is the chemical potential of the ideal gas. In fact, combination of the last two terms on the right hand side of Eq. (13) is the target excess chemical potential. By adding/deleting the fractional particle into/from the simulation box, decision is made either to choose another real particle in the box as the new fractional particle, or to add a new fractional particle to the box. This is done according to the algorithm by Eslami and Müller-Plathe (2007b). On adding the new fractional particle to the box, the coupling parameter, , is set close to zero. Similarly, when a real particle is chosen as the new fractional particle, is set close to 1. Repeating such a procedure, equilibrium is achieved, in which the density (number) of particles in the system oscillates about an average value, corresponding to the target chemical potential, temperature, and volume.

### **4. Phase equilibrium calculation**

68 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

The discontinuity of the number of particles, introduced through particle insertion/deletion, makes it difficult to do MD simulation in the grand canonical ensemble. Therefore, in some of the existing MD-based simulations in the grand canonical ensemble, the sampling is performed through using the MC-based stochastic procedures (Boinepalli & Attard, 2003). Another alternative, is the introduction of the particle number as a continuous dynamic variable through the MD simulation (Eslami & Müller-Plathe, 2007b; Eslami et al., 2011). This approach is based on the so-called extended system method, in which a Hamiltonian including the kinetic and potential energy terms for the extension variable is constructed

Here we assume that our physical system is composed of *N* real particles plus one "fractional" particle and is coupled to a heat reservoir and a particle reservoir. The fractional

zero to one. The inclusion of the fractional particle in the system provides a variable number of particles in a dynamical way. The Hamiltonian of such a combined system; the real *N* particle system, the fractional particle, and the additional degrees of freedom, i.e., the material reservoir and the heat reservoir, is written as (Eslami & Müller-Plathe, 2007b):

<sup>222</sup> <sup>222</sup> <sup>1</sup> 2 2

*xyz xyz <sup>s</sup> ij if <sup>s</sup>*

(12)

1 1 <sup>1</sup> 2 2 2 2

*i i f i ji i s ppp ppp p p <sup>H</sup> <sup>U</sup> U UU m s m s W W*

where *H* is the Hamiltonian, subscripts *i* and *j* refer to the real particles, subscript *f* refers to the fractional particle, *p*i is the momentum of the *i*th particle, *s* is the velocity scaling variable

which couples the fractional particle to the rest of the system and varies from zero to one. In this equation, the first four terms on the right hand side are the kinetic energy of the real particles, potential energy of interaction between real particles, kinetic energy of the fractional particle, and the potential energy of interaction between real particles and the

variable, s, have their own kinetic and potential energy terms, shown as the last four terms

from a medium of zero potential to the medium of interest with the potential energy *U*; i.e.,

Adopting one of the real particles (molecules) in the simulation box as the fractional particle, the fractional particle is grown to a full (real) particle or deleted from the system in a dynamical way, by solving the equations of motion. The relevant equation of motion, governing the dynamics of fractional particle (the magnitude of coupling to the ideal gas

*<sup>N</sup> if id*

 

*U*

the last two terms on the right hand side of Eq. (13) is the target excess chemical potential. By adding/deleting the fractional particle into/from the simulation box, decision is made

 

1

*i*

 

*W* is the particle-number scaling variable

, and the temperature extension

is the work required to carry a fractional particle

(13)

*id* is the chemical potential of the ideal gas. In fact, combination of

which varies from

**3.3.2 Grand canonical ensemble molecular dynamics simulation** 

particle is a particle whose potential energy is weighted by a variable

2 2

*f ff iii N N N N*

which couples the system to the thermostat, and

fractional particle. The number extension variable,

on the right hand side of Eq. (12). In fact *U*

it is the so-called excess Helmholtz free energy.

reservoir) is (Eslami & Müller-Plathe, 2007b):

/dt2 and

where d *λ*

 =d2

and the appropriate equations of motion are solved.

There are several methods for the calculation of phase coexistence points using molecular simulations, such as thermodynamic scaling method (Valleau, 1991; Valleau and Graham, 1990), histogram reweighting method (McDonald and Singer, 1967; Potoff and Panagiotopoulos, 1998), the Gibbs-Duhem integration method (Kofke, 1993a,b), NpT plus test particle method (Boda et al., 1995; 1996), various extensions of it to other ensembles (Moller and Fischer, 1990; Boda et al., 2001), and the Gibbs eensemble MC method (Panagiotopoulos,1987). In the Gibbs eensemble MC method, which is one of the most popular methods for the calculation of phase equilibria, the simulation box is divided into two subsystems and simultaneous simulation of both subsystems are done. To establish the equilibrium (the equality of chemical potentials), particles are exchanged between the two subsystems. This technique has been applied to coexistence properties of simple systems, such as fluids of spherical Lennard-Jones or Yukawa particles (Panagiotopoulos et al., 1998; Rudisill and Cummings, 1989), as well as more complex systems, such as polyatomic hydrocarbons (de Pablo and Prausnitz, 1989; de Pablo et al., 1992) and chain molecules (de Pablo, 1995). There are also reports on the mixed methods in which the molecular simulation approaches have been utilized, to calculate the chemical potentials in the condensed phase and the results from equations-of-state predictions are used to calculate the phase coexistence point (Lim et al., 2003), or to calculate the interaction energy parameters of solvent and polymer, in combination with statistical-mechanical theories for the study of phase equilibria of polymer-solvent mixtures (Jang and Bae, 2002).

Many computational studies of the permeation of small gas molecules through polymers, using afore-cited techniques, have appeared, which were designed to analyze, on an atomic scale, diffusion mechanisms or to calculate the diffusion coefficient and the solubility parameters. Most of these studies have dealt with flexible polymer chains of relatively simple structure such as polyethylene, polypropylene, and poly-(isobutylene) (Müller-Plathe, 1991; Pant & Boyd, 1993; Tamai et al., 1995; Pricl and Fermeglia, 2003; Abu-Shargh, 2004). There are however some reports on polymers consisting of stiff chains. Of these we may address to the works by Mooney and MacElroy (1999) on the diffusion of small molecules in semicrystalline aromatic polymers, by Cuthbert *et al.* (1997) on the calculation of the Henry's law constant for a number of small molecules in polystyrene and the effect of box size on the calculated Henry's law constants, by Lim *et al.* (2003) on the sorption of methane and carbon dioxide in amorphous polyetherimide. In some of these reports the authors have used Gibbs ensemble Monte Carlo method (Vrabec and Fischer, 1995; Panagiotopoulos, 1987) to do the calculations, and in some other cases alternative methods, like the equation-of-state approaches are employed to describe the gas phase. In fact, as explained above, one needs to satisfy the equality of chemical potentials in both phases. There are, however, some recent methods employing the grand canonical ensemble simulation formalism (to set the target chemical potentials in equilibrating phases), without the necessity to do simultaneous simulations on two boxes (unlike the Gibbs ensemble MC method) to calculate the gas solubilities in polymers. In the following section we explain these methods.

#### **4.1 Grand equilibrium method**

In the Gibbs ensemble simulation method one specifies the thermodynamic variables temperature, global composition, and global pressure for the simulation of both phases in separate volumes. Practically, this set of thermodynamic variables is in many cases not convenient and simultaneous simulation of both phases has the disadvantage that fluctuations occurring in one phase influence the other one. Recently a new method, grand equilibrium method, has been developed by Vrabec and Hasse (2002). This method circumvents the afore-cited problems for the study of phase equilibria. The specified thermodynamic variables are temperature and composition and two independent simulations are performed for the two phases without the need to exchange particles in the condensed phase. According to this method for a mixture composed of several components, it is possible to do a simulation in the isothermal-isobaric (*NpT*) ensemble at constant temperature, a constant composition of the condensed phase, and at an arbitrary constant pressure, preferably close to the pressure at the phase coexistence point, to obtain the density of the condensed phase.

In the grand equilibrium method, a simulation of the condensed phase is done to calculate the excess chemical potentials, *s ex*, and the partial molar volumes, *Vs*, of all components. One may use the test-particle insertion method (Widom, 1963) to calculate the excess chemical potentials and the partial molar volumes as:

$$
\mu\_s^{ex} = -k\_B T \ln \left\langle \frac{pV}{Nk\_B T} \exp \left( -\Delta U \,/ \, k\_B T \right) \right\rangle \tag{14}
$$

and

$$V\_s = \frac{\left\langle V^2 \exp\left(-\Delta U \,/\, k\_B T\right) \right\rangle}{\left\langle V \exp\left(-\Delta U \,/\, k\_B T\right) \right\rangle} - \left\langle V \right\rangle \tag{15}$$

Knowing the parameters *Vs* and *s ex* from this simulation, the desired excess chemical potentials as functions of pressure are obtained from a first order Taylor series expansion, i.e.,

$$
\mu\_s^{ex}(p) = \mu\_i^{ex}(p^\*) + \frac{V\_s}{kT}(p - p^\*) \tag{16}
$$

where *p\** is the target pressure at which the *NpT* ensemble simulation is done. Once the *s ex(p)* is determined from Eq. (16) by one *NpT* ensemble simulation of the condensed phase, one vapour/gas simulation has to be performed in the pseudo-grand canonical ensemble (pseudo*-VT*). In a common grand canonical ensemble (Eslami and MüllerPlathe, 2007b) the parameters temperature, volume, and the chemical potential of all species are fixed, while in this pseudo*-VT* ensemble simulation, the parameters *T* and *V* are fixed in the common way, but instead of fixing the chemical potentials, they are set as a function of the pressure in the gas phase. This procedure ensures that equilibrium between the condensed phase and the gas phase is imposed. In a common *VT* ensemble simulation (Eslami and Müller-Plathe, 2007b), the chemical potentials are set through insertion and deletion of particles by a comparison between the resulting potential energy change and the desired fixed chemical potential. Here, starting from a low density state point, the gas phase simulation is forced to change its state to the corresponding phase equilibrium state point.

#### **5. Applications**

70 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

There are, however, some recent methods employing the grand canonical ensemble simulation formalism (to set the target chemical potentials in equilibrating phases), without the necessity to do simultaneous simulations on two boxes (unlike the Gibbs ensemble MC method) to calculate the gas solubilities in polymers. In the following section we explain

In the Gibbs ensemble simulation method one specifies the thermodynamic variables temperature, global composition, and global pressure for the simulation of both phases in separate volumes. Practically, this set of thermodynamic variables is in many cases not convenient and simultaneous simulation of both phases has the disadvantage that fluctuations occurring in one phase influence the other one. Recently a new method, grand equilibrium method, has been developed by Vrabec and Hasse (2002). This method circumvents the afore-cited problems for the study of phase equilibria. The specified thermodynamic variables are temperature and composition and two independent simulations are performed for the two phases without the need to exchange particles in the condensed phase. According to this method for a mixture composed of several components, it is possible to do a simulation in the isothermal-isobaric (*NpT*) ensemble at constant temperature, a constant composition of the condensed phase, and at an arbitrary constant pressure, preferably close to the pressure at the phase coexistence point, to obtain the

In the grand equilibrium method, a simulation of the condensed phase is done to calculate

One may use the test-particle insertion method (Widom, 1963) to calculate the excess

ln exp / *ex s B B B pV k T U kT Nk T*

> <sup>2</sup> exp / exp /

*V U kT V V V U kT*

potentials as functions of pressure are obtained from a first order Taylor series expansion,

 \* \* () ( ) *ex ex <sup>s</sup> s i V <sup>p</sup> p pp kT*

where *p\** is the target pressure at which the *NpT* ensemble simulation is done. Once the

*ex(p)* is determined from Eq. (16) by one *NpT* ensemble simulation of the condensed phase, one vapour/gas simulation has to be performed in the pseudo-grand canonical

*VT*). In a common grand canonical ensemble (Eslami and Müller-

 

 

*B*

*B*

*ex*, and the partial molar volumes, *Vs*, of all components.

(14)

(15)

*ex* from this simulation, the desired excess chemical

(16)

*s*

*s*

*s*

chemical potentials and the partial molar volumes as:

these methods.

**4.1 Grand equilibrium method** 

density of the condensed phase.

the excess chemical potentials,

Knowing the parameters *Vs* and

and

i.e.,

*s*

ensemble (pseudo*-*

Recently, we have applied the grand equilibrium method for the calculation of solubilities of gases in polymers. The method has been well tested for the case of gas solubilities in polystyrene over a wide range of temperatures and pressures (Eslami and Müller-Plathe, 2007a). To calculate the phase equilibeium points, MD simulations were performed at a specified temperature and concentrations of sorbed gases in polymer. After equilibration, several configurations are extracted at different times from the dynamic simulation and used to insert the test particles, to calculate the excess chemical potentials. Calculating the excess chemical potentials of the sorbed gases in polymer as a function of pressure, a simulation of the gas phase is done in the grand canonical ensemble, to find the phase coexistence point. We have shown in Figure 3 typical results for the zero-pressure limit solubility coefficients of methane in polystyrene (Eslami and Müller-Plathe, 2007a).

Fig. 3. Comparison of calculated () and experimental ( Sada et al. 1987; Vieth et al., 1966a; Vieth et al., 1966b; Barrie et al. 1980; Raymond et al. 1990; Lundberg et al. 1969) zero-pressure limit solubility coefficients of methane in polystyrene. Figure is taken from Eslami and Müller-Plathe (2007a).

Another more interesting application of the method is in the case of gas solubility calculation for polar penterants in polar polymers. This is a severe test, as most of the existing methods, like the Widom's test particle insertion method or its modifications, are not applicable in the case of polar fluids. We have recently applied this method to the case of water solubilities of poly(ethylene terephthalate) (Eslami & Müller-Plathe, 2009) and polyamide-6,6 (Eslami & Mehdipour, 2011). Preparing an initially relaxed configuration of the polymer, we have performed a long simulation, in the grand canonical ensemble, of the polymer phase to get the density of water, which corresponds to the target values of , *V*, and *T*. The number of molecules (density) during the simulation changes until achieving equilibrium, at which the density fluctuates about the average value. We have shown in Figure 4 the variation of number of water molecules in the simulation box at two different temperatures (one above and one below the glass transition temperature).

In Figure 4, the results at *T*= 280 K belong to *exs* = -18.5 kJmole-1 in a cubic box with a box length, *L*, of 3.082 nm. Calculating the partial molar volume of water from the results of this simulation, we expanded the chemical potential of water in the polymer phase as a function of pressure, according to Eq. (16). Then we have performed another simulation of the gas phase in the grand canonical ensemble, as explained above. Here the target chemical potential is varied during the simulation as a function of running-average pressure in the gas phase (see Eq. (16)). The result for time evolution of the density of water in the gas phase (*L*=50 nm) is also shown in Figure 4. The results in Figure 4 show that the gas phase quickly adjusts its state to the phase coexistence point. Also shown in the same figure are the results of similar calculations at 450 K (*exs* = -9.0 kJmole-1). Here the values of *L* for simulations of polymer and gas phases are 3.128 nm and 7.0 nm, respectively.

Fig. 4. Time evolution of the number of water molecules in the grand canonical ensemble simulation at 280 K and 450 K. The solid and dotted curves represent the results for water molecules in poly(ethylene terephthalate) and in gaseous phases, respectively.

Another more interesting application of the method is in the case of gas solubility calculation for polar penterants in polar polymers. This is a severe test, as most of the existing methods, like the Widom's test particle insertion method or its modifications, are not applicable in the case of polar fluids. We have recently applied this method to the case of water solubilities of poly(ethylene terephthalate) (Eslami & Müller-Plathe, 2009) and polyamide-6,6 (Eslami & Mehdipour, 2011). Preparing an initially relaxed configuration of the polymer, we have performed a long simulation, in the grand canonical ensemble, of the polymer phase to get the density of water, which corresponds to the target values of

and *T*. The number of molecules (density) during the simulation changes until achieving equilibrium, at which the density fluctuates about the average value. We have shown in Figure 4 the variation of number of water molecules in the simulation box at two different

length, *L*, of 3.082 nm. Calculating the partial molar volume of water from the results of this simulation, we expanded the chemical potential of water in the polymer phase as a function of pressure, according to Eq. (16). Then we have performed another simulation of the gas phase in the grand canonical ensemble, as explained above. Here the target chemical potential is varied during the simulation as a function of running-average pressure in the gas phase (see Eq. (16)). The result for time evolution of the density of water in the gas phase (*L*=50 nm) is also shown in Figure 4. The results in Figure 4 show that the gas phase quickly adjusts its state to the phase coexistence point. Also shown in the same figure are the results

Fig. 4. Time evolution of the number of water molecules in the grand canonical ensemble simulation at 280 K and 450 K. The solid and dotted curves represent the results for water

molecules in poly(ethylene terephthalate) and in gaseous phases, respectively.

temperatures (one above and one below the glass transition temperature).

polymer and gas phases are 3.128 nm and 7.0 nm, respectively.

In Figure 4, the results at *T*= 280 K belong to

of similar calculations at 450 K (

, *V*,

*exs* = -18.5 kJmole-1 in a cubic box with a box

*exs* = -9.0 kJmole-1). Here the values of *L* for simulations of

Similar calculations are done in the case of water solubility of polyamide-6,6 (Eslami & Mehdipour, 2011) over a wide range of temperatures and pressures. This is a particularly interesting example of the applicability of the method, since polyamide-6,6 is a highly polar polymer matrix with the ability of forming a strong hydrogen bond network and can dissolve water to about 10 wt %. We have shown in Figure 5 a typical calculated sorption isotherm for water in polyamide-6,6 at 300 K and its comparison with experiment.

Fig. 5. The sorption isotherm of water in PA-6,6 at 300 K. The symbols represent calculations () and experimental data at 298 K () (Lim et al., 1999) and 300 K () (Watt, 1980). Figure is taken from (Eslami & Mehdipour, 2011).

### **6. Diffusion in polymers**

Most theories describing the mechanism of diffusion in polymeric materials are based on the free-volume approximation (Vrentas & Duda, 1979). In the free-volume theories, there is a volume which is directly occupied by polymeric molecules, and there is the remainder of the volume, which is called the free volume. A part of the free volume is assumed to be uniformly distributed among the molecules and is identified as the interstitial free volume, which requires a large energy for redistribution and is not affected by random thermal fluctuations. The other part, which is called the hole free volume, is assumed to require negligible energy for its redistribution. Therefore, the hole free volume is being continuously redistributed due to random fluctuations, and is assumed to be occupied by penetrant molecules. This redistribution of hole free volume will move the penetrant molecule with it. According to this model, by movement of segments of the polymer chain, a void will be created adjacent to the penetrant molecule. If the size of this hole is sufficient to host a penetrant molecule and if the penetrant have sufficient energy to jump into the hole, a successful jump of the penetrant molecule is made into the hole. Although the free volume model has been used extensively to describe the mechanism of transport through molten or glassy polymers (Mauritz & Storey, 1990), this model does not show a microscopic view point of penetrant transport in polymers, since it just connects bulk transport properties, like diffusion coefficient, into bulk properties, like molecular volume or thermal expansion coefficient.

The transition-state theory (TST), introduced by Gusev et al. (1993), is another useful method for the calculation of diffusion coefficient of a low-molar-mass substance through the polymer matrix. In the TST it is assumed that the movement of the penetrant from an initial cavity to the saddle point and to a neighboring cavity is a unimolecular rearrangement. For such a transition the reaction trajectory in configuration space is tracked and the transition rate constant is evaluated. In the first studies on the application of TST to study the dynamics of light gases dissolved in rigid microstructures of glassy polycarbonate and rubbery polyisobutylene Gusev et al. (1993), the method was shown to be only capable to study just the dynamics of light gases like He. The method was developed by Gusev and Suter (1993), by Gusev et al. (1994), and by Karayiannis et al. (2004), to calculate the diffusion coefficient of bigger penetrants in glassy polymers.

Computer modeling of molecular systems at a detailed atomistic level has become a standard tool in investigation of sorption and diffusion of small molecules in polymeric media (Müller-Plathe, 1991, 1994, Milano et al., 2002, Mozaffari et al. 2010). MD simulation is a useful tool for exploring the structure and properties of bulk amorphous polymers. The length of the trajectories that can be generated in practice presently is on the order of many nanoseconds. Thus the range of properties that can be studied directly is limited to those that evolve over this time scale. One of the phenomena that appear to be suitable for investigation is the diffusion of small penetrant molecules in an amorphous polymeric matrix. That is, the diffusion coefficients of small penetrants in many rubbery or liquid polymers are such that, at temperatures close to room temperature and above, the average displacement of the diffusant is large enough in a nanosecond interval to be determined via molecular dynamics simulation. Performing such simulations is of practical importance in predicting diffusion coefficients and also in understanding the mechanism of diffusion.

Although there has been significant progress in the use of MD methods in the simulations of diffusion coefficients, early studies were focused on the simulation of gas diffusion in rubbery polymers which could be investigated using full atomistic or united-atom simulations in reasonable computational times (Boyd, 1991; Sok et al., 1992, Pant & Boyd, 1993; Gee & Boyd, 1995). Due to the recent development of improved force fields and the wider availability of sophisticated commercial softwares and high-speed computing facilities, attention is shifting to direct to the more challenging task of simulating the slower diffusional processes occurring in glassy polymers (Han & Boyd, 1996, Milano et al. , 2002, Lime at al., 2003, Mozaffari et al., 2010). In the following section we present the MD simulation results done recently by the authors on the diffusion of gases in polystyrene over a wide range of temperatures and pressures.

#### **6.1 MD simulation of diffusion in polymers (case study)**

Recently we have studied the diffusion mechanism of some gases in polystyrene over a wide range of temperatures, ranging from below the glass transition temperature to temperatures well above the glass transition temperature (Mozaffari et al., 2010). Center-ofmass mean-square displacements have been measured during the simulation to calculate the diffusion coefficients. In the limit of long times, which the penetrant molecules perform random walks in the polymer matrix, the mean-square-displacement becomes linear in time, and the diffusion coefficients can be calculated using the Einstein relation:

The transition-state theory (TST), introduced by Gusev et al. (1993), is another useful method for the calculation of diffusion coefficient of a low-molar-mass substance through the polymer matrix. In the TST it is assumed that the movement of the penetrant from an initial cavity to the saddle point and to a neighboring cavity is a unimolecular rearrangement. For such a transition the reaction trajectory in configuration space is tracked and the transition rate constant is evaluated. In the first studies on the application of TST to study the dynamics of light gases dissolved in rigid microstructures of glassy polycarbonate and rubbery polyisobutylene Gusev et al. (1993), the method was shown to be only capable to study just the dynamics of light gases like He. The method was developed by Gusev and Suter (1993), by Gusev et al. (1994), and by Karayiannis et al. (2004), to calculate the

Computer modeling of molecular systems at a detailed atomistic level has become a standard tool in investigation of sorption and diffusion of small molecules in polymeric media (Müller-Plathe, 1991, 1994, Milano et al., 2002, Mozaffari et al. 2010). MD simulation is a useful tool for exploring the structure and properties of bulk amorphous polymers. The length of the trajectories that can be generated in practice presently is on the order of many nanoseconds. Thus the range of properties that can be studied directly is limited to those that evolve over this time scale. One of the phenomena that appear to be suitable for investigation is the diffusion of small penetrant molecules in an amorphous polymeric matrix. That is, the diffusion coefficients of small penetrants in many rubbery or liquid polymers are such that, at temperatures close to room temperature and above, the average displacement of the diffusant is large enough in a nanosecond interval to be determined via molecular dynamics simulation. Performing such simulations is of practical importance in predicting diffusion coefficients and also in understanding the

Although there has been significant progress in the use of MD methods in the simulations of diffusion coefficients, early studies were focused on the simulation of gas diffusion in rubbery polymers which could be investigated using full atomistic or united-atom simulations in reasonable computational times (Boyd, 1991; Sok et al., 1992, Pant & Boyd, 1993; Gee & Boyd, 1995). Due to the recent development of improved force fields and the wider availability of sophisticated commercial softwares and high-speed computing facilities, attention is shifting to direct to the more challenging task of simulating the slower diffusional processes occurring in glassy polymers (Han & Boyd, 1996, Milano et al. , 2002, Lime at al., 2003, Mozaffari et al., 2010). In the following section we present the MD simulation results done recently by the authors on the diffusion of gases in polystyrene over

Recently we have studied the diffusion mechanism of some gases in polystyrene over a wide range of temperatures, ranging from below the glass transition temperature to temperatures well above the glass transition temperature (Mozaffari et al., 2010). Center-ofmass mean-square displacements have been measured during the simulation to calculate the diffusion coefficients. In the limit of long times, which the penetrant molecules perform random walks in the polymer matrix, the mean-square-displacement becomes linear in time,

diffusion coefficient of bigger penetrants in glassy polymers.

mechanism of diffusion.

a wide range of temperatures and pressures.

**6.1 MD simulation of diffusion in polymers (case study)** 

and the diffusion coefficients can be calculated using the Einstein relation:

$$D = \frac{1}{6N} \lim\_{t \to \infty} \frac{d}{dt} \left\langle \sum\_{i}^{N} [\mathbf{r}\_i(t) - \mathbf{r}\_i(0)]^2 \right\rangle \tag{17}$$

The diffusion coefficients of penetrant gases in polystyrene have been calculated over a wide range of temperatures, 300 K-500 K. The motion pattern of penetrant gases in host polymer can be qualitatively studied by monitoring the penetrant's displacement |*r*(t)-*r*(0)|, from its initial position. We have shown in Figure 6 the displacements of argon and propane in polystyrene at 300 K.

Fig. 6. Displacement of argon (upper curve) and propane (lower curve) molecules from their initial positions at 300 K. In order to avoid overlapping between the curves, the displacements of argon are offset by 0.6 nm. Figure is taken from Mozaffari et al. (2010).

The curves are representative of a common hopping mechanism, showing that for a considerable time interval the penetrants dwell in existing voids in the polymer and occasionally do a jump into the neighboring voids. When dwelling in the voids, the penetrants just perform oscillatory motions around their equilibrium positions, therefore, no net motion of a penetrant molecule occurs with these positional fluctuations. The amplitude of oscillations varies according to the size of voids. From time to time, the penetrants can do a quick jump into their neighboring voids, see Figure 6. The jump frequency depends on penetrant's size, therefore the bigger penetrants (see Figure 6 in the case of propane) can rarely jump between the voids.

Particularly useful information can be obtained by analyzing the trajectories of penetrants in the polymer. The two-dimensional center-of-mass *x-y* trajectories of nitrogen, as a typical example, in polystyrene at temperatures below (300 K) and above (500 K) the glass transition temperature are indicated in Figure 7.

The results in Figure 7 are indicative of faster movement of penetrants at higher temperatures, as represented by the broadened range of displacements at higher temperatures. This indicates that at higher temperatures the hole free volume redistributes faster and penetrant molecules have higher energy to overcome the activation energy required to jump into new voids.

Fig. 7. Typical trajectories of nitrogen molecules in polystyrene at 300 K (upper curve) and 500 K (lower curve). Figure is taken from Mozaffari et al. (2010).

As a measure of diffusion coefficients, the center-of-mass mean-square displacements of carbon dioxide at temperatures below (340 K) and above (480 K) the glass transition temperature, are shown in Figure 8.

Fig. 8. Center-of-mass mean-square displacement for carbon dioxide at temperatures below and above the glass transition temperature. The curve at 480 K is scaled down by a factor of 15 for the sake of clarity. The dashed lines show the least-squares fits to the linear parts of the curves. Figure is taken from Mozaffari et al. (2010).

The linearity of mean-square displacements versus time, as indicated in Figure 8, confirms Einstein diffusion. In the glassy polymer at times below 500 ps, the penetrant motion is highly anomalous and the diffusion regime begins at longer times. Intercavity jumps rarely occur at this time scale. However, as it can be seen from the results in Figure 8 at high temperatures, the diffusion regime sets in a shorter time. These findings clearly indicate the difference between the diffusion mechanisms in the glassy and melt polymers.

#### **7. Permeation in polymers**

76 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

faster and penetrant molecules have higher energy to overcome the activation energy


Fig. 7. Typical trajectories of nitrogen molecules in polystyrene at 300 K (upper curve) and

As a measure of diffusion coefficients, the center-of-mass mean-square displacements of carbon dioxide at temperatures below (340 K) and above (480 K) the glass transition

y (nm)

0 500 1000 1500 2000 2500 3000 3500

480 K

340 K

t (ps)

Fig. 8. Center-of-mass mean-square displacement for carbon dioxide at temperatures below and above the glass transition temperature. The curve at 480 K is scaled down by a factor of 15 for the sake of clarity. The dashed lines show the least-squares fits to the linear parts of

required to jump into new voids.

0

0

the curves. Figure is taken from Mozaffari et al. (2010).

0.2

0.4

0.6

0.8

1

1.2

temperature, are shown in Figure 8.

<(r(t)-r(0))2> (nm2

)

500 K (lower curve). Figure is taken from Mozaffari et al. (2010).

1

2

x (nm)

3

For a permeant gas subject to a pressure gradient at two different sides of a membrane, the permeability, P, is defined as

$$\text{P} = \frac{\text{(quantity of } \text{ permeant)} \times \text{(film } \text{ thickness)}}{\text{(area)} \times \text{(time)} \times \text{(pressure drop across film)}} \tag{18}$$

The quantity of permeant gas is often expressed as the volume of gas at STP condition. In fact, the pressure difference between up- and down-stream sides of polymer membrane is the driving force of permeation phenomenon. As stated above, penetrant transport through membrane is described by solution-diffusion mechanism. According to this model the permeation occurs in three steps: 1-sorption of gas molecules at upstream side of polymer 2 diffusion of the gas molecules into the polymer matrix 3-desorption of gas molecules at downstream side of polymer. This is indicated schematically in Figure 9.

Fig. 9. The schematic of gas permeation in a membrane

The permeability coefficient of gas molecules across a polymer film of thickness *L* is described by the flux *J*:

$$J = P \frac{\Delta \mathcal{L}}{L} \tag{19}$$

where *Δc* is the concentration gradient of penetrant molecules through the membrane. Accordingly, the permeation coefficient can be expressed by the product of the diffusion coefficient *D*, and the solubility coefficient *S*, as:

$$P = DS \tag{20}$$

Since direct simulation of the permeation process is not possible due to the system size, Eq. (20) is a particularly useful expression connecting the permeability coefficient with the diffusion and solubility coefficients. Therefore, diffusion and solubility coefficients, calculated in MD simulation methods, are used to calculate the permeability coefficients according to Eq. (20).

We have shown in Figure 10 our recently calculated results on the permeability of water in poly(ethylene terephthalate) (Eslami & Müller-Plathe, 2009). The results, indicating the permeability coefficients over a wide range of temperatures, are shown in Figure (10).

Fig. 10. Comparison of the calculated () and experimental permeability coefficients of water in poly(ethylene terephthalate). The experimental data are from ( Launay et al., 1999); ( Shigetomi et al., 2000); ( Changjing & Junying, 1987); ( Rueda & Varkalis, 1995); and ( Kloppers et al., 1993). Figure is taken from Eslami & Müller-Plathe (2009).

The Arrhenius plot of ln(*P*) *vs.* 1/*T* shows a break at the glass transition temperature. As a result of calculating a higher solubility coefficient, our calculated permeation coefficients are also higher than the experimental values. Solubility coefficients are higher than the experimental values, especially at lower temperatures. This is trivial for calculations of this type. Similar differences between experimentally and computed solubility coefficients have been observed in previous studies (Cuthbert et al., 1997; Lim et al., 2003, Eslami & Müller-Plathe, 2007a, 2009; Eslami & Mehdipour, 2011). According to Gusev and Suter (1993), errors of the order of (2-4)*k*B*T* are common for the calculation of Helmholtz energies by molecular simulations. It has been speculated that the main contribution to the solubility comes from single holes in the simulated polymer structure, which might not be present in similar proportion in real polymers (Müller-Plathe et al., 1993). Errors of this type become even more serious in cases where the solubility is very low.

We have also compared our calculated permeability coefficients of a number of gases in polystyrene (Mozaffari et al., 2010) with experimental values (Csernica et al., 1987) and with simulation results by Kucukpinar and Doruker (2003). The results in Table 1 indicate that the calculated values are in good agreement with experiment. The results indicate that polystyrene is much permeable to CO2, compared to other gases studied by Mozaffari et al. (2010), because of the higher solubility coefficient of CO2 in polystyrene and its relatively higher diffusion coefficient. On the other hand, polystyrene is less permeable to the biggest penetrant molecule studied, propane, because of its very small diffusion coefficient.

We have shown in Figure 10 our recently calculated results on the permeability of water in poly(ethylene terephthalate) (Eslami & Müller-Plathe, 2009). The results, indicating the permeability coefficients over a wide range of temperatures, are shown in Figure (10).

2.0 2.5 3.0 3.5 4.0 4.5

1000/*T* (K-1

Fig. 10. Comparison of the calculated () and experimental permeability coefficients of water in poly(ethylene terephthalate). The experimental data are from ( Launay et al., 1999); ( Shigetomi et al., 2000); ( Changjing & Junying, 1987); ( Rueda & Varkalis, 1995); and ( Kloppers et al., 1993). Figure is taken from Eslami & Müller-Plathe (2009).

The Arrhenius plot of ln(*P*) *vs.* 1/*T* shows a break at the glass transition temperature. As a result of calculating a higher solubility coefficient, our calculated permeation coefficients are also higher than the experimental values. Solubility coefficients are higher than the experimental values, especially at lower temperatures. This is trivial for calculations of this type. Similar differences between experimentally and computed solubility coefficients have been observed in previous studies (Cuthbert et al., 1997; Lim et al., 2003, Eslami & Müller-Plathe, 2007a, 2009; Eslami & Mehdipour, 2011). According to Gusev and Suter (1993), errors of the order of (2-4)*k*B*T* are common for the calculation of Helmholtz energies by molecular simulations. It has been speculated that the main contribution to the solubility comes from single holes in the simulated polymer structure, which might not be present in similar proportion in real polymers (Müller-Plathe et al., 1993). Errors of this type become even

We have also compared our calculated permeability coefficients of a number of gases in polystyrene (Mozaffari et al., 2010) with experimental values (Csernica et al., 1987) and with simulation results by Kucukpinar and Doruker (2003). The results in Table 1 indicate that the calculated values are in good agreement with experiment. The results indicate that polystyrene is much permeable to CO2, compared to other gases studied by Mozaffari et al. (2010), because of the higher solubility coefficient of CO2 in polystyrene and its relatively higher diffusion coefficient. On the other hand, polystyrene is less permeable to the biggest penetrant molecule studied, propane, because of its very small diffusion

)


more serious in cases where the solubility is very low.

coefficient.

ln (*p* (g cm/cm2Pa s))

An interesting test is to compare the calculated permeability coefficient ratios (selectivities) in the zero pressure limit with the corresponding experimental values (Csernica et al., 1987). This is done in Table 2, in which we have listed the ratios of permeability coefficients at 300 K to that of nitrogen's permeability and compared the results with experimental measurements (Csernica et al., 1987) and with the calculations by Kucukpinar and Doruker (2003). The results in Table 2 show that the calculated ratios are quite close to the experimental ratios. This shows that our calculated permeability coefficients are higher that the experimental values by nearly the same factor. The same conclusion was made in our work on the calculated solubility coefficients (Eslami & Müller-Plathe, 2007a). This indicates that there is an excessive free volume in the polymer sample compared to the experimental polymer samples. Better agreement of the results by Mozaffari et al. (2010) with experiment (Csernica et al., 1987) compared to those of Kucukpinar and Doruker (2003) has been attributed to using a united atom models for both polymer and diffusant gases, in the later case.


Table 1. Comparison of calculated gas permeability coefficients at 300 K (Mozaffari et al., 2010) with experimental values (Csernica et al., 1987) and with simulation results by Kucukpinar and Doruker (2003).


Table 2. Comparison of calculated ratios of permeability coefficients at 300 K (Mozaffari et al., 2010) with experimental values (Csernica et al., 1987) and with simulation results by Kucukpinar and Doruker (2003).

### **8. Summary**

In this chapter the polymer permeability is discussed and reviewed from a computational point of view. Although it is impossible to directly simulate the permeability process, the simulation techniques are powerful tools to simulate and to give a molecular level insight to the solubility and diffusion mechanisms of gases in polymers. While there are lots of simulations in the literature, giving address to the diffusion mechanism of gases in glassy and/or rubbery polymers, the computations are less straightforward in the case of gas solubilities. Therefore, the main problem in calculation of gas perameability of polymers, using molecular simulation methods, stems from the calculation of solubility coefficients. In fact the calculation of gas solubilities necessitates the condition of equilibrium between the permeant in the gas/liquid phase and the permeant dissolved in polymer. The equality of the temperature, pressure, and chemical potentials of all species is the necessary condition to establish such an equilibrium situation. The chemical potential is, however, coupled to the number of particles and cannot easily be calculated employing molecular simulation methods.

The applicability and feasibility of different techniques for the calculation of chemical potentials is evaluated and discussed. It is explained that in some straightforward cases, like the thermodynamic integration method, several simulations are needed to calculate the chemical potential. Widom's test particle insertion method (Widom, 1963) is of practical importance in many cases, but the method encounters sampling problems at high densities and in the case of big solutes. The method is also not applicable in the case of chemical potential calculation for polar species. It is shown that recent grand canonical MD simulation methods (Eslami and Müller-Plathe, 2007b; Eslami et al., 2010). are trustable methods for the sake of phase equilibrium calculation. In these methods the chemical potential is set as an independent thermodynamic quantity, and the number of particles in the box is changed in a gradual and dynamical way. The applicability of the method to the challenging cases like gas solubilities in polymers with stiff chemical structures (Eslami and Müller-Plathe, 2007a) and solubility of water (as a polar penetrant) in polar polymers, like polyamide-6,6 (Eslami & Mehdipour, 2011) is reviewed and discussed.

The process of gas diffusion in polymers and different approaches for studying the gas diffusion are addressed to. It is explained that MD simulation well describes the jumping mechanism of diffusions in glassy polymers. According to this mechanism, in a glassy polymer the penetrants dwell in the polymer holes and occasionally do a jump into the neighboring voids (Mozaffari et al., 2010). Therefore, no net motion of a penetrant molecule occurs while the penterant just performs random oscillations inside a hole. From time to time, the penetrants can do a quick jump into their neighboring voids, with a jump frequency depending on penetrant's size, temperature, and chemical structure of the polymer. The permeation coefficients of gases in polymers are, therefore, reviewed in this chapter in terms of gas solubilities and diffusions.

### **9. References**

Abu-Shargh, B. F. (2004) *Polymer* 45, 6383. Adams, D. J. (1974) *Mol. Phys.* 28, 1241. Adams, D. J. (1975) *Mol. Phys.* 29, 1307. Barrie, J. A.; Williams, M. J. L.; Munday, K. (1980) *Polym. Eng. Sci.*, 20, 21. Ben-Naim, A. (1987) *Solvation thermodynamics*, Plenum, New York.

Boda, D.; Chan, K.-Y.; Szalai, I. (1997) *Mol. Phys.*, 92, 1067.

80 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

In this chapter the polymer permeability is discussed and reviewed from a computational point of view. Although it is impossible to directly simulate the permeability process, the simulation techniques are powerful tools to simulate and to give a molecular level insight to the solubility and diffusion mechanisms of gases in polymers. While there are lots of simulations in the literature, giving address to the diffusion mechanism of gases in glassy and/or rubbery polymers, the computations are less straightforward in the case of gas solubilities. Therefore, the main problem in calculation of gas perameability of polymers, using molecular simulation methods, stems from the calculation of solubility coefficients. In fact the calculation of gas solubilities necessitates the condition of equilibrium between the permeant in the gas/liquid phase and the permeant dissolved in polymer. The equality of the temperature, pressure, and chemical potentials of all species is the necessary condition to establish such an equilibrium situation. The chemical potential is, however, coupled to the number of particles and cannot easily be calculated employing molecular simulation methods. The applicability and feasibility of different techniques for the calculation of chemical potentials is evaluated and discussed. It is explained that in some straightforward cases, like the thermodynamic integration method, several simulations are needed to calculate the chemical potential. Widom's test particle insertion method (Widom, 1963) is of practical importance in many cases, but the method encounters sampling problems at high densities and in the case of big solutes. The method is also not applicable in the case of chemical potential calculation for polar species. It is shown that recent grand canonical MD simulation methods (Eslami and Müller-Plathe, 2007b; Eslami et al., 2010). are trustable methods for the sake of phase equilibrium calculation. In these methods the chemical potential is set as an independent thermodynamic quantity, and the number of particles in the box is changed in a gradual and dynamical way. The applicability of the method to the challenging cases like gas solubilities in polymers with stiff chemical structures (Eslami and Müller-Plathe, 2007a) and solubility of water (as a polar penetrant) in polar polymers, like

polyamide-6,6 (Eslami & Mehdipour, 2011) is reviewed and discussed.

Barrie, J. A.; Williams, M. J. L.; Munday, K. (1980) *Polym. Eng. Sci.*, 20, 21. Ben-Naim, A. (1987) *Solvation thermodynamics*, Plenum, New York.

chapter in terms of gas solubilities and diffusions.

Abu-Shargh, B. F. (2004) *Polymer* 45, 6383. Adams, D. J. (1974) *Mol. Phys.* 28, 1241. Adams, D. J. (1975) *Mol. Phys.* 29, 1307.

**9. References** 

The process of gas diffusion in polymers and different approaches for studying the gas diffusion are addressed to. It is explained that MD simulation well describes the jumping mechanism of diffusions in glassy polymers. According to this mechanism, in a glassy polymer the penetrants dwell in the polymer holes and occasionally do a jump into the neighboring voids (Mozaffari et al., 2010). Therefore, no net motion of a penetrant molecule occurs while the penterant just performs random oscillations inside a hole. From time to time, the penetrants can do a quick jump into their neighboring voids, with a jump frequency depending on penetrant's size, temperature, and chemical structure of the polymer. The permeation coefficients of gases in polymers are, therefore, reviewed in this

**8. Summary** 

Boda, D.; Liszi, J.; Szalai, I .(1995) *Chem. Phys. Lett.* 235, 140. Boda, D.; Winkelmann, J.; Liszi, J.; Szalai, I. (1996) *Mol. Phys.* 87, 601. Boda, D.; Kristof, T.; Liszi J.; Szalai, I. (2001) *Mol. Phys.* 99, 2011. Boinepalli S.; Attard, P. (2003) *J. Chem. Phys.* 119, 12769. Boyd, R. H.; Pant, P. V. K. (1991) *Macromolecules*, 24, 6325. Changjing, L.; Junying, G. (1987) *Chin. J. Polym. Sci.*, *5*, 261. Cohn, M. H.; Turnbull, D. (1959) *J. Chem. Phys.*, 31, 1164. Csernica, J.; Baddour, R. F.; Cohen, R. E. (1987) *Macromolecules*, 20, 2468. Cuthbert, T. R.; Wagner, N. J.; Paulaitis, M. E. (1997) *Macromolecules* 30, 3058 Daynes, H. A. (1920) *Proc. R. Soc. London Ser. A*, 97, 286. Deitrick, G. L.; Scriven, L. E.; Davis, H. T. (1989) *J. Chem. Phys.* 90, 2370. de Pablo, J. J. (1995) *Fluid Phase Equilib.* 104, 195. de Pablo, J. J.; Prausnitz, J. M. (1989) *Fluid Phase Equilib.* 53, 177. de Pablo, J. J.; Bonnin, M.; Prausnitz, J. M. (1992) *Fluid Phase Equilib.* 73, 187. Dibenedetto, A. T. (1963) *J. Polym. Sci. A*, 1, 3477. Escobedo F. A.; de Pablo, J. J. (1995) *J. Chem. Phys.* 103, 2703. Eslami, H.; Mehdipour, N. (2011) *Phys. Chem. Chem. Phys.* 13, 669-673. Eslami, H.; Mojahedi, F.; Moghadasi, J. (2010) *J. Chem. Phys.* 133, 084105.

