Brownian Dynamics Simulations of Macromolecules: Algorithm Development and Polymers under Confinement

To my parents

List of Papers

This thesis is based on the following papers, which are referred to in the text by their Roman numerals.

I Carlsson, T., Sundberg, J., Arteca, G. A., Elvingson, C. Off-equilibrium response of grafted polymer chains subject to a variable rate of com- pression. Phys. Chem. Chem. Phys., 2011, (13), 11757-11765

II Carlsson, T., Kamerlin, N., Arteca, G. A., Elvingson, C. Brownian dynamics of a compressed polymer brush model. Off-equilibrium response as a function of surface coverage and compression rate. Phys.

Chem. Chem. Phys., 2011, (13), 16084-16094

III Carlsson, T., Kamerlin, N., Arteca, G. A., Elvingson, C. Structure and compression dynamics of two mutually approaching polymer brushes.

Manuscript

IV Kamerlin N., Ekholm, T., Carlsson T., Elvingson, C. Construction of a non-periodic closed network for computer simulations. Manuscript V Carlsson, T., Ekholm, T., Elvingson, C. Algorithm for generating a

Brownian motion on a sphere. J. Phys. A:Math. Theor. (43), 2010, 505001

Reprints were made with permission from the publishers.

Contents

1 Introduction . . . . 11

2 Polymeric systems . . . . 13

2.1 Description of polymer molecules . . . . 13

2.2 Ideal chains . . . . 15

2.3 Non-ideal chains . . . . 16

2.4 Polymer mushrooms and brushes . . . . 18

3 Computer simulations . . . . 19

3.1 Brownian motion . . . . 19

3.2 Boundary conditions . . . . 22

3.2.1 Periodic boundary conditions . . . . 22

3.2.2 Spherical boundary conditions . . . . 22

3.3 Brownian dynamics simulations . . . . 24

4 Results and discussion . . . . 27

4.1 Polymer mushrooms . . . . 27

4.1.1 Polymer mushrooms during compression . . . . 27

4.1.2 Relaxation from the compressed state . . . . 30

4.1.3 Compression with a finite sized object . . . . 32

4.2 A single polymer brush . . . . 33

4.2.1 Role of surface coverage during compression . . . . 33

4.2.2 Role of solvent quality during compression . . . . 35

4.2.3 Relaxation from the compressed state . . . . 37

4.3 Two approaching polymer brushes . . . . 38

4.3.1 Role of surface coverage during compression . . . . 38

4.4 Constructing a closed polymer network . . . . 40

4.5 Diffusion on a spherical surface . . . . 42

4.5.1 Computation of the probability density . . . . 42

4.5.2 Computation of the diffusion coefficient . . . . 43

5 Concluding remarks . . . . 45

6 Svensk sammanfattning . . . . 47

7 My contribution . . . . 51

8 Acknowledgements . . . . 53

Bibliography . . . . 55

Nomenclature

Latin alphabet

a Particle radius

b i Bond vector i

c(r,t) Concentration

d Particle diameter

D ₀ Diffusion coefficient

j(r,t) Particle flux

k _b Bending force constant

k _B Boltzmann constant

k _s Stretching force constant

L Inter-wall distance

L ^ Scaled inter-wall distance

L Contour length

n Number of bonds in a polymer chain

N Number of monomers in a polymer chain

N ¯ Mean overcrossing number (entanglement)

P Persistence length

P ^ Scaled persistence length

P( θ,t) Probability density function

q _i Point charge i

r S ² -radius

R S ³ -radius

r i Position vector of particle i

R ee End-to-end vector

R g Radius of gyration

R ² _{i j} Radius of gyration tensor

S Degree of oblateness/prolateness

T Temperature

U _b Bending potential

U _s Stretching potential

Greek alphabet

α Bond angle

β “Opening angle” of a paraboloid tip

ε Lennard-Jones energy parameter

ε ^ Scaled Lennard-Jones energy parameter ε r Relative permittivity

ε 0 Permittivity of vacuum

ζ Friction coefficient

η Viscosity

θ Angular distance

σ Surface coverage

σ LJ Lennard-Jones size parameter

τ Dimensionless time

ξ Side length of a subcube in I ³

ψ Angle between position vectors in R ⁴

1. Introduction

Macromolecules are large molecules which can have different architectures, like chains, combs, ladders etc. but they are usually comprised of a limited set of smaller building blocks (monomers) which are chemically connected to large molecules (polymers). Different types of polymers are made from differ- ent sets of monomers, e.g. biopolymers are constructed by monomers of bio- logical origin like the base pairs in DNA and the amino acids of proteins. The effective interactions between the monomers lead to different preferred con- formations, which will influence the average size and shape of the molecules and therefore, for biopolymers, also their biological function. ¹ Polymers are not only found as biomacromolecules in our bodies, but they are also used in a large number of industrial processes as well as in many applications en- countered in everyday life. Here we can mention the various use of plastics, fibres (examples being polyacryl and cotton) in the textile industry and cross- linked polymers (gels) in pharmaceutical applications. Certain polymers have end groups which can be grafted onto different types of surfaces and which can for instance be used to reduce surface friction ^2–6 , for steric stabilization of colloidal suspensions ^{7, 8} and as biocompatibilizers. ⁹ The abundance and large applicability of polymers make them an important class of molecules to understand in more detail.

For practical applications as well as from a more scientific point of view, it

is important to understand how the physical properties depend on e.g. molec-

ular weight and type of solvent. A conceptually appealing approach is to use

scaling arguments as introduced by Flory ¹⁰ to predict how the size of a free

polymer chain in solution will increase with the number of monomers. This

method was extended by de Gennes ¹¹ and Alexander ¹² to also include the

effects of interfaces. The grafting to an interface or the confinement of a poly-

mer molecule in a dramatic way change the average shape and size of the

polymer ¹³ as a chain may in this case adopt conformations unlikely for free

polymers in solution. This is due to an interplay between configurational en-

tropy, bending energy, excluded volume effects and other non-bonded inter-

actions which are competing to minimize the free energy. One model which

captures many important features of the equilibrium behaviour of a polymer

in confinement is a chain in a slit between repulsive walls, and this model sys-

tem has been extensively used for studies of polymers in various solvents and

level of confinement. ^14–18

If we, in addition, are interested in the dynamics or the behaviour for short polymer chains, we no longer can rely on scaling theory in the same way. A valuable tool in this case is the use of computer simulations where we also can take into account more complex interactions. In this thesis we are pri- marily concerned with the change in average size, shape and entanglement of grafted polymers chains, when subjected to a dynamical confinement. Using Brownian dynamics simulations to capture also, the time dependence, we have studied the effect of how different surface coverages as well as different rates of compression will influence the conformation of semi-flexible polymers for different effective interactions. A coarse grained model makes it possible to perform simulations of larger molecules compared to atomistic simulations, as the molecular nature of the solvent is replaced by an effective friction and a stochastic force at this level of description. In the present study, we started by simulating a single polymer chain during confinement. A common method at the nanoscale to measure, manipulate and modify surfaces is by using an atomic force microscope. ^19–21 In addition to the behaviour of a single chain compressed by a planar surface, we also introduced a model to study the effect of a compressing paraboloid tip with variable sharpness. We further performed simulations of a polymer coated surface confined by a bare upper surface as well as of two polymer coated surfaces under compression for various surface coverages.

Having investigated the features of polymer brushes under confinement, a

natural step would be to extend the studies to a crosslinked polymer network,

in which the motion of the polymer chains is further restricted. We present in

this thesis a method to construct a closed network which can be used in com-

puter simulations. This is done by generating the network in S ³ , the surface of

a ball in R ⁴ . Investigating the motion of a Brownian particle or a grafted poly-

mer chain confined to move on a surface is of both theoretical interest ²² as

well as important in various applications. ^{23, 24} The thesis ends with describing

a novel method for generating a Brownian motion on a spherical surface.

2. Polymeric systems

2.1 Description of polymer molecules

The structure of a linear polymer chain can be described using different pa- rameters characterizing e.g. its size and shape. Due to the bending and tor- sional degrees of freedom, the polymer can be more or less flexible and has a large variety of available conformations. A certain value for one property can correspond to several values for another descriptor which makes it essential to use several parameters to discriminate between different conformations. One such descriptor is the end-to-end distance which, with the notations in figure 2.1, is defined as the vector sum of all bond vectors

R ee = ∑ ⁿ

i=1

b i (2.1)

where b i = r i+1 − r i is the bond vector between the center of monomer i and i + 1, n = N − 1 is the number of bonds and N denotes the number of chain units which we choose to be spherical beads.

1 2

3

4 5

i i + 1 i + 2

N − 1 N N − 2

r2

r1

b1 b2

b3

b4

bi

bN−1

bN−2

α12

αi−1i

bi+1

Figure 2.1: A schematic picture of a polymer modeled as N spherical beads connected by n = N − 1 bonds. The position vectors r

and r

defines the bond vector b

= r

−r

. For freely jointed chains the N −2 angles α

between bond vectors b

and b

are randomly distributed.

The maximum magnitude of the end-to-end distance R ee = |R ee | for a given

model is called the contour length L , which also gives an estimate of the size

of a linear chain and is defined by

L = ∑ ⁿ

i=1

|b i |. (2.2)

Another size descriptor, which in contrast to the end-to-end distance R _ee and the contour length, L , is well defined also for more complex molecular struc- tures such as branched or ring polymers, is the radius of gyration, R g , which for a homopolymer, where all monomers have the same mass is defined by

R ² _g = 1 N

∑ N i=1

(r i − r cm ) ² (2.3)

where r cm is the center of mass of the polymer molecule. This descriptor can also be written as a function of the monomer positions

R ² _g = 1 N ²

N−1 ∑

i=1

∑ N j=i+1

(r i − r j ) ² (2.4)

The radius of gyration can also be obtained from the radius of gyration ten- sor ²⁵

R ² _{i j} = 1 N

∑ N k=1

(r ⁱ _k − r ⁱ _cm )(r _k ^j − r _cm ^j ) (2.5)

where the indexes i and j refer to the three directions in space. Denoting the three eigenvalues of the corresponding matrix by λ i , R g can be calculated us- ing ²⁵

R ² _g = ∑ ³

i=1

λ i . (2.6)

These eigenvalues describe the extension of the polymer chain in the three principal directions. By combinations of these eigenvalues, it is possible to construct size independent shape descriptors. One such parameter is the de- gree of oblateness/prolateness, S, defined by ²⁶

S = 27 

Π ³ _i (λ i − ¯λ) 

 (∑ ³ i λ i ) ³  (2.7)

where ¯ λ is the average eigenvalue. This descriptor has a range between

−0.25 ≤ S ≤ 2 where a negative value corresponds to an oblate shape and a

positive value signifies a more prolate shape, with S = 0 for a sphere.

A polymer system is not exclusively characterized by the average size and shape of the individual molecules. It can also be important to describe the level of entanglement within and between polymer chains and one parameter which describes this property is the mean overcrossing number which, for a linear chain, is defined by ²⁷

N = ¯ 1 2π

N−3 ∑

i=1 N−1 ∑

j=i+2

 ₁ 0

 ₁ 0

|( ˙γ i (s) × ˙γ j (t)) · (γ i (s) − γ j (t))|

γ i (s) − γ j (t) ³ dsdt (2.8) where γ i is the parameterized bond-vector connecting monomers i and i+1 and γ ˙ i its parametric derivative, satisfying

γ i : [0,1] → R ³ , γ i (0) = r i γ i (1) = r i+1 . (2.9)

2.2 Ideal chains

For an ideal chain there are no interactions between two beads far apart along the backbone, even if they are close in distance. A chain is thus allowed to cross itself, which is not physically realistic. These models can, however, serve as good approximations, when the attractive and repulsive interactions between the monomers approximately cancel. By introducing restrictions for the available positions of consecutive beads we can construct different models.

The simplest model without any non-bonded interactions is called the freely jointed chain, where the bond length, b = |b i |, is constant and consecutive bonds vectors are uncorrelated. With the notations in figure 2.1 the average square end-to-end vector in this model can be calculated according to

 R ² _ee 

=

 n

∑ i

b i · ∑ ⁿ

j

b j



= ∑ ⁿ

i

∑ n j

 b i · b j

 (2.10)

= b ² ∑ ⁿ

i

∑ n j

 cos α i j

 = b ² n (2.11)

where α i j denotes the angle between bond vectors i and j. This shows that the average end-to-end distance for a random walk may be written, in a way commonly used in polymer physics, as 

R ² _ee  _1/2

= bn ^ν , where the scaling exponent is ν = 1/2 in this particular case.

The bond angles in a polymer are not completely random, due to e.g. steric

hindrance, and by introducing bond angle correlations, one can take a step

towards a more generally applicable model. This is done in the freely rotating

model where each bond has a constant length and the bond angles α ii+1 are

all equal to a constant α. The bond vector is, however still free to rotate in the azimuthal angle as depicted in figure 2.2.

b _i b i+1

b _i+2 b i+3

α α

α

Figure 2.2: Figure showing bond vectors in the freely rotating model. All bond angles are the same ( α

= α) but the azimuthal angle is random.

In this way the bond angles and therefore the positions of the beads will depend on the positions of the previous bonds. Given a position along the backbone of a polymer chain, the average bond angle correlation decreases exponentially with the distance s along the chain according to ²⁸

cosα i j (s) = e ^−s/P . (2.12) Here P is called the persistence length and is a measure of how far from a given position in the chain we have to travel, in order for the bond angle correlations to become negligible. In this way P will quantify the stiffness of the chain, with a large value corresponding to stiff chains and a small value to more flexible polymers.

For the freely rotating model the average square end-to-end distance can be calculated analytically and for large n it is given by ¹⁰

R ² _ee  = nb ² 1 + cosα

1 − cosα . (2.13)

2.3 Non-ideal chains

If we take into consideration the non-bonded (long range) interactions bet-

ween the monomers as well as monomer-solvent interactions, the available

space for the chain is in a first approximation reduced by the introduction of an

excluded volume. The chain will now become a self-avoiding walk obtaining

a more extended conformation relative a freely jointed chain. This is reflected

in a larger scaling exponent which now becomes ν = 3/5 (a more exact nu-

merical calculation ²⁹ will give ν = 0.588). The excluded volume effect can be

incorporated into a model by a potential, U, which describes the interactions

between non-connected units along the chain and preventing them from over-

lapping. The effective interaction between two monomers can, in addition to

this short range repulsion originating from the overlap of two close monomers,

also include a long range attraction. If we denote the distance between two in- teracting monomers by r, this contribution to the interaction potential can be shown ³⁰ to be ∼ 1/r ⁶ . A widely used potential which models the combination of these two effects is the Lennard-Jones potential given by ³¹

U _LJ = 4ε ∑ ^N

i< j

 σ LJ

r _{i j 12}

−

 σ LJ

r _{i j 6}

. (2.14)

Here r _{i j} = |r i − r j | is the distance between units i and j, σ LJ is a measure of the effective size of a monomer and ε is the Lennard-Jones parameter describing the strength of the interaction. A large ε for the inter-bead potential means that the polymer effectively tries to minimize its interface with the solvent and is referred to as a bad or poor solvent condition and conversely for a small ε, referred to as a good solvent condition.

In the simulations of polymeric systems in this thesis we have used a dis- crete wormlike chain model, where a harmonic stretching force allows for a small variation in the bond length and a bending force makes it possible to vary the stiffness of the chains. The bending potential is given by

U _b = k _b 2

N−2 ∑

i=1

(cosα ii+1 − cosα 0 ) ² (2.15)

where α 0 is the bond angle of zero potential energy and k b is the bending force constant, while the bond stretching is modeled as a harmonic potential given by

U _s = k s

2

N −1 i=1 ∑

(b _i − 2a) ² (2.16)

where a is the monomer radius and k _s is the stretching force constant. The persistence length for this model can be analytically calculated and is given by ³²

P = b 1 − cosα ^N−1

1 − cosα (2.17)

where b and cosα is the average bond length and average cosine bond an-

gle respectively. By providing P and b as input parameters, the bending force

constant can be calculated in a self-consistent way using equations (2.15) and

(2.17). ^{33, 34}

2.4 Polymer mushrooms and brushes

The models above refer primarily to a free polymer molecule in solution, but in many situations of practical interest, the polymer can have an endgroup grafted to a surface. A single grafted molecule or a polymer chain in a dilute layer in a good solvent at equilibrium, have a structure often called a “mush- room”. For polymers in a denser layer, a mushroom conformation can not be accomodated and the favourable conformation will be a “brush”, describing the more stretched out conformations due to excluded volume forces. Accord- ing to the mean field theory by de Gennes, ¹¹ the equilibrium height, h, of a polymer brush in a good solvent is, for a large number of monomers, propor- tional to N

h ∝ Nσ ^1/3 (2.18)

where σ is the surface coverage. The polymer layer has in this case the pos-

sibility to extend beyond the range of the attractive van der Waals forces and

e.g. induce steric stabilization between colloidal particles in a dispersion.

3. Computer simulations

As we have seen in the previous section we can, for example, obtain the scal- ing behaviour of polymer size with the number of monomers for an ideal chain at equilibrium. When we build more detailed models, the equations often be- come increasingly difficult to solve and computer simulations is an invaluable tool. In a computer simulation we have to follow only a fraction of the parti- cles in the macroscopic sample we want to model, and a way to avoid surface effects is to impose periodic boundary conditions as will be explained below.

We are often not only interested in the equilibrium properties but we also want to follow the time dependent dynamic properties. The time evolution of a polymeric system can often be described by Newton’s equations of motion. A molecular dynamics simulation solving these equations will involve the calcu- lation of a large number of interactions due to the solvent molecules. This will ultimately lead to a large computational cost which nessesarily will reduce the time during which we can follow our system using an atomistic model. A way to reduce the CPU-time and still be able to mimic the important interactions is by coarse graining, which means reducing the number of degrees of freedoms in a system. The effect of this process is shown in the next section where the solvent interactions are replaced by a friction force together with a stochastic force. In this way we reduce the computational cost and obtain results on a time scale relevant for the structural features of a polymeric system.

3.1 Brownian motion

The random motion of dispersed colloid particles in a solvent, called Brow- nian motion, has its name after the botanist Robert Brown who studied the motion of pollen grains immersed in water in the beginning of the 19th cen- tury. ³⁵ For that particular case the irregular motion is due to the constant bom- bardment of the water molecules experienced by the pollen grains, but can always be observed when sufficiently small macromolecules is immersed in a fluid and experience thermal motion. For the theoretical treatment we can use the Einstein-Smoluchowski approach ³⁶ which follows the time evolution of a probability distribution or the Langevin approach, ³⁷ in which a stochastic differential equation is governing the motion of the particles.

If we consider the diffusion of, for example, colloidal particles in water,

a nonuniform concentration c(r,t) will create a flux until a uniform particle

distribution is obtained when equilibrium is reached. As described by Ficks’

first law, ³⁸ this flux, j(r,t), is proportional to the gradient of the concentration

j(r,t) = −D 0 ∇c(r,t) (3.1)

where D ₀ is the diffusion coefficient. Combining the above equation with the continuum equation

∂c(r,t)

∂t = −∇ · j(r,t) (3.2)

leads to the diffusion equation

∂c(r,t)

∂t = D 0 ∇ ² c( r,t). (3.3)

If the particles are under the influence of an external force field, the above equation 3.3 must be modified. ³⁹ The velocity due to this force, F = −∇U(r) can be written as

v = − 1

ζ ∇U(r) (3.4)

where 1/ζ is the mobility and ζ is the friction coefficient, and this additional force leads to a modified flux

j(r,t) = −D 0 ∇c(r,t) − 1

ζ ∇U(r)c(r,t). (3.5)

At equilibrium the net flux is zero and assuming the concentration to be Boltz- mann distributed

c _eq ∝ e ^−U(r)/k

^T (3.6)

it follows that the diffusion coefficient is related to the friction constant by D ₀ = k B T

ζ (3.7)

where k _B is the Boltzmann constant and T is the temperature. The friction co- efficient is throughout this thesis assumed to obey the Stokes-Einstein equa- tion ³⁸ which relates the friction coefficient to the viscosity and size of the Brownian particle by

ζ = 6πηa (3.8)

where η is the viscosity of the solvent and a being the radius of the particle.

Using the continuum equation again but now for the flux of particles in an external force field, the equation for the concentration becomes ³⁹

∂c(r,t)

∂t = ∇ 1

ζ (k B T ∇c(r,t) + c(r,t)∇U(r)). (3.9) This equation is called the Smoluchowski equation. The concentration is pro- portional to the probability to find a Brownian particle at position r at time t so this equation is also satisfied by the probability distribution function from which statistical averages of a dynamic variable can be calculated.

The Langevin approach starts instead from the Newton’s equation for the Brownian particle and the effect of the solvent is taken into account by a fric- tion force F _i ^f and a stochastic force F ^s _i

m ¨r i = F ^tot _i = −∇U(r) + F _i ^f + F ^s _i . (3.10) The friction force is assumed to be proportional to the velocity of the particle relative to the quiescent fluid

F _i ^f = −ζ ˙r i (3.11)

and the random stochastic force F ^s _i , describing the diffusive motion is assumed to be Gaussian, defined by the moments ³⁹

F ^s _i  = 0 (3.12)

and F ^s _i (t) · F ^s _i (t ^ ) = 6ζk B T δ(t −t ^ ). (3.13)

The dynamics of macromolecules is usually well described in the overdamped limit in which the acceleration is set to zero. With this approximation the equation of motion 3.10 becomes

ζ ˙r i = −∇U(r) + F ^s _i . (3.14) An algorithm describing the time evolution of the position of a Brownian par- ticle from time t to a later time t + Δt, may be written ⁴⁰

r i (t + Δt) = r i (t) + D ₀ Δt

k B T F i + B i . (3.15) Here F i = −∇U(r) and B i is a random Gaussian variable modeling the col- lisions with the solvent medium having zero mean and with a variance given

by B i · B j  = 6D 0 Δtδ i j . (3.16)

In the computer simulations in this thesis we neglected hydrodynamic inter- actions and the diffusion coefficient is position independent in which case the Langevin equation (3.14) and the Smoluchowski equation (3.9) are equiva- lent. ³⁹ By following the time evolution of the positions of the Brownian parti- cles using equation 3.15 for a number of initial configurations, often referred to as trajectories, we can calculate the average value for any position depen- dent dynamical variable. In this case one should create many trajectories and calculate the average at different time steps.

3.2 Boundary conditions

3.2.1 Periodic boundary conditions

In a simulation it is not computationally possible to follow the time evolution of all the atoms in a macroscopic sample and we are forced to follow only a fraction of the corresponding number of particles. This number is far from Avogadros number and one way to mimic the bulk behaviour is to employ periodic boundary conditions. Exact copies of our simulation cell with the particles, also called the primary cell, is replicated in all directions (see figure 3.1). A particle leaving the simulation box is entering at the opposite side of the box. In this way the particles in the primary cell should represent the bulk behaviour and when calculating the forces on these particles one should take into account any interactions from particles in the primary cell, as well as from the infinite number of replicas. A truncation of the force at a certain distance, r _cut , will introduce an error, but the forces involved are often of a short ranged nature making the error negligable in this case. In the so called minimum image convention, as indicated in figure 3.1, the nearest replica within r cut is considered in the force evaluation.

In simulations where also long-ranged forces are included, e.g. electrostatic interactions, which decay slowly due to the Coulumb law, it is essential to include the contribution from several replicas.

3.2.2 Spherical boundary conditions

To avoid the time consuming summations, involved in the electrostatic force calculations associated with periodic boundary conditions using e. g.

Ewald summation, ⁴¹ one can instead use spherical boundary condition. This

method was proposed by Kratky, ⁴² simulating the equilibrium properties of a

Lennard-Jones fluid and have been further developed to include electrostatic

interactions for this geometry in a series of paper by Caillol ^43–45 as well

as for the dynamics of polymeric systems. ^{46, 47} In this method the system

is embedded on the three dimensional surface of a ball in R ⁴ , denoted

by S ³ and particles are interacting along geodesics in this curved space.

rcut

Figure 3.1: Periodic boundary conditions in two dimensions. The simulation cell is in the middle with exact replicas in all directions.

Proper attention has now to be given to the functional form of the equations governing the interactions and the equation of motion. For example the (non normalized) solution of the diffusion equation in S ³ , from which random numbers describing the Brownian motion can be generated, is no longer a simple Gaussian and will instead be given by ⁴⁶

P(θ,t) = sinθ ∑ ^∞

k=−∞ (θ + 2πk)e ^−R

^(θ+2πk)

^/4D

^t (3.17) where θ is an angle characterizing the size of the random step along the geodesic and R is the radius of the hypersphere. The electrostatic force be- tween two electric charges q i and q j is derived from the solution of the Poisson equation in this geometry and may be written as ⁴³

F i = q i q j

4π ² R ² ε r ε 0

cot ψ + π − ψ sin ² ψ

t i j . (3.18)

Here ψ is the angle between the position vectors r i ,r j ∈ R ⁴ and t i j is the tangent vector to the geodesic connecting the two charges.

When using spherical boundary conditions the equations describing possi-

ble constraints of a particular system must also be changed accordingly. In

R ³ , a confining surface can e.g. be described by z =constant, in contrast to

in S ³ where the corresponding surfaces confining the molecules are ordinary spheres, S ² .

3.3 Brownian dynamics simulations

We can now describe the main steps during the course of a Brownian dynamics simulation which are shown in the flow chart in figure 3.2.

• We first specify input parameters needed in the simulation, such as the time step Δt, the temperature T and the viscosity η of the solvent.

• Then we generate an initial configuration which preferably should be close to an equilibrium structure.

• We apply boundary conditions and calculate the total force on each particle, which also include calculating the forces due to any geometrical constraint for the particles, as well as generating the random numbers describing the stochastic motion.

• Move the particles according to an integrated Langevin equation of motion.

• Steps 3 and 4 are repeated a given number of times during which we store the positions of the particles at regular intervals. If we are at the end of the trajectory, we generate a new initial configuration.

• This procedure is performed for a predetermined number of initial

configurations, leading to different trajectories from which statistical

averages of position dependent dynamical variables can be calculated.

Move the particles

Store particle positions

New

No

Yes parameters

trajectory Specify input

Generate initial

Calculate forces

Calculate averages configuration

Figure 3.2: Flow chart describing the main steps during a Brownian dynamics simu-

lation.

4. Results and discussion

In this section we present a summary of the findings in the manuscripts in- cluded at the end of this thesis. We first discuss the size, shape and entangle- ment during compresssion of polymers grafted to a solid surface and also the subsequent relaxation process both during a continued confinement but also after removing the upper surface. We start the analysis for a single chain and continue with the behaviour of a polymer brush for variable surface coverages.

We then consider the situation when two polymer coated surfaces approach each other. Having considered the structure and dynamics of polymer coated surfaces, we continue by presenting an alghorithm for generating a crosslinked polymer network, without having to include periodic boundary conditions.

In the last part we present a novel algorithm for generating a Brownian mo- tion on a spherical surface. For all simulations the temperature, T , is equal to 293.15 K, the bead diameter, d, is 0.4 nm and the viscosity, η, of the solvent is 1.002 · 10 ⁻³ Pa s which is equal to that for water at the corresponding tem- perature. Part of the results are presented using the following dimensionless parameters: the persistence length P ^ = P/d, the inter-wall distance L ^ = L/d, and the Lennard-Jones parameter ε ^ = ε/k B T .

4.1 Polymer mushrooms

4.1.1 Polymer mushrooms during compression

As a starting point for understanding the behaviour of a grafted polymer layer during compression, and for characterizing the behaviour of a polymer chain under confinement in general, we first consider a single grafted polymer mush- room. This situation can be regarded as a limiting case for a very low surface coverage. The polymer is modeled as a linear chain of 50 connected spherical beads with the first unit grafted to the lower surface. The scaled Lennard-Jones parameter, ε ^ , reflecting the magnitude of non-connected inter-bead interac- tions, is in the range ε ^ = 0.1 to 1.25 modeling a polymer in good to poor solvent conditions. The equilibrium structures of these chains are usually re- ferred to as a soft and hard polymer mushroom repectively (see figure 4.1).

The polymer is subjected to compression by an upper planar surface and we

follow the change in the radius of gyration, R g , the asphericity, S, and the

overcrossing number, ¯ N, at different rates of compression. The values of the

descriptors are ensemble averages over 50 independent trajectories and for

convenience we introduce the notation R _g ≡ R ² _g  ^1/2 . A slow compression rate corresponds to one compressing move of the upper surface every hundred time steps, whereas a fast compression rate corresponds to moving down the upper surface every second time step. The incremental change in wall-wall distance, ΔL, is always set to 0.01a, with a being the radius of a bead. The results below are for chains with a persistence length P ^ = 2, modeling a rather flexible poly- mer. Before starting the compression, we created an ensemble of equilibrated chains for given P ^ and ε ^ . The upper surface was then introduced just above z = z _max , with z max being the maximum z-coordinate over the total ensemble of 50 trajectories for every set of parameters. The polymer was then in a first instance compressed to z = z max  by moving the upper surface downwards, after which the chains were equilibrated once again monitoring both R g and S.

In this way we started the actual compression phase from a state where most chains in the ensemble were only slightly effected by the upper surface. The compression was then continued until the inter-wall distance was L ^ = 1.5, which corresponds to a highly compressed state.

Figure 4.1: Snapshots from two different simulations showing typical configurations before the compression. The left figure shows a soft polymer mushroom (P

= 2 and ε

= 0.1) and the right picture shows a hard polymer mushroom (P

= 2 and ε

= 1.25).

In figure 4.2, the radius of gyration and the asphericity are shown as a function of distance between two hard walls for soft mushrooms.The polymer chains have a persistence length P ^ = 2 and two compression rates have been used. During the compression, R g goes through a shallow minimum and, in the case of slow compression, finally reaches a value close to the initial value, whereas a faster compression will lead to an almost monotonic decrease in the average molecular size. At the same time the asphericity shows that a soft mushroom becomes more prolate during the compression compared to the ini- tial value. For a slow compression rate, the chains go through a maximum in S while a fast compression results in a monotonically increasing asphericity.

In figure 4.3 the radius of gyration and asphericity for hard mushrooms with

a larger ε ^ are presented. These chains have a smaller initial size compared to a

softer mushroom, due to more compact structures modeling a poor solvent. A

fast compression will also in this case results in the chain going through a min-

imum in R g . For a slow compression rate the size of the chain increases dur-

ing the compression, after an initial resistance to the confinement. A different

behaviour for hard mushrooms is also mirrored in the observed asphericity.

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2

0.5 1 1.5 2 2.5 3 3.5 4

Rg / nm

L / nm

ε ε

ε

´=0.1

´=0.25

´=0.5

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.5 1 1.5 2 2.5 3 3.5 4

S

L / nm

ε´=0.1 ε´^=0.25

ε´^=0.5

Figure 4.2: Radius of gyration and asphericity vs. the inter-wall distance L during the compression for soft mushrooms with P

= 2. The strength of the interactions is given by the dimensionless parameter ε

= ε/k

T . Solid lines are used for a slow compression (one compression step every 100 time step) and the dashed lines are the result for fast compression (one compression step every second time step).

A compression leads here to more oblate structures compared to the initial shapes. For a fast compression, a hard mushroom becomes more oblate go- ing through a small maximum in the asphericity, while a slower compression, when the chains have more time to relax between every compression step, results in a faster decrease in S.

0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4

Rg / nm

L / nm

ε´^=0.75

ε´=1.0 ε´^=1.25

−0.2

−0.15

−0.1

−0.05 0 0.05 0.1 0.15

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4

S

L / nm

ε´^=0.75

ε´=1.25

Figure 4.3: Radius of gyration and asphericity vs. the inter-wall distance L during the compression. The results are for hard mushrooms with P

= 2 with the notation and parameters as in figure 4 .2. The results for ε

= 1.0 is omitted for clarity.

In figure 4.4 the entanglement is shown as a function of inter-wall distance

for soft- and hard mushrooms respectively. For soft mushrooms (figure 4.4

left) a slow compression leads to an increased entanglement in the initial stage

of the compression, while a faster compression rate results in an almost con-

stant value or only a slight increase before, also in this case, finally decrease

at the end of the compression. For chains with a larger ε ^ (figure 4.4 right) this

descriptor is almost unchanged in the initial stage after which a decrease is

observed for both compression rates.

4 4.5 5 5.5 6 6.5

0.5 1 1.5 2 2.5 3 3.5 4

Entanglement

L / nm

ε´=0.1

ε´=0.25

10 11 12 13 14 15 16 17 18 19

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4

Entanglement

L / nm

ε´=1.25 ε´=1.0

ε´=0.75

Figure 4.4: Entanglement vs. the inter-wall distance L during the compression. The figures show the results for for soft mushrooms (left) and hard mushrooms (right).

The results for a slow compression are shown with solid lines while the results using a fast compression are shown with dotted lines.

4.1.2 Relaxation from the compressed state

After the compression phase, the chains were allowed to relax following two different routes: (i) instantly removing the upper wall after the compression and (ii) by first letting the chain relax during 5 · 10 ⁶ time steps in the com- pressed state, before also in this case removing the upper wall. In figure 4.5 the time dependence for R _g is shown for a slow compression followed by the two different relaxation processes. An immediate removal of the upper surface (the lower curve in every pair) leads to a fast relaxation to the equi- librium value. A relaxation in the compressed state (the upper curve in every pair) instead leads to a drift in the radius of gyration as the chains have the possibility to extend in the slit and this drift, more pronounced in the case of a smaller ε ^ , indicates a nonequilibrium state at the end of the compression phase. A release of the upper surface after the relaxation in the compressed state, will also now be followed by a relaxation to the equilibrium value for an uncompressed chain.The equilibrium value for the radius of gyration is close to but not coinciding with the initial value, since we started the compression from a slightly compressed state.

In figure 4.6 the equivalent processes for hard polymer mushrooms are pre-

sented. An immediate removal of the upper surface (route (i)) is followed by a

fast relaxation back to the equilibrium value for the radius of gyration. When

letting the chains relax in the compressed state (route (ii)), they find a new

stable state before, also in this case, relaxing back to the equilibrium value for

unconfined chains when removing the upper surface.

1 1.2 1.4 1.6 1.8 2 2.2 2.4

0 10 20 30 40 50 60 70

Rg / nm

Time / ns

ε=0.1´ ε=0.25´

ε=0.5´

Figure 4.5: The relaxation dynamics for a soft polymer mushroom. The figure shows the results from simulations of chains interacting with three different Lennard-Jones energy parameters: ε

= 0.1 (solid), ε

= 0.25 (dashed) and ε

= 0.5 (dotted). The lower curve for every ε

follows the relaxation protocol (i), while the upper curve are the results from relaxation protocol (ii).

0.7 0.8 0.9 1 1.1 1.2 1.3 1.4

0 5 10 15 20 25

Rg / nm

Time / ns

ε=0.75´

ε=1.0´ ε=1.25´

Figure 4.6: The relaxation dynamics for a hard polymer mushroom given as the time

dependent behaviour for the radius of gyration. The figure shows the results from

simulations with chains interacting with three different Lennard-Jones energy param-

eters: ε

= 0.75 (solid), ε

= 1.0 (dashed) and ε

= 1.25 (dotted). The lower curve for

every ε

follows the relaxation protocol (i), while the upper curve are the results from

relaxation protocol (ii).

4.1.3 Compression with a finite sized object

To qualitatively investigate the behaviour of chains subjected to compres- sion by a finite-sized object, we performed a series of simulations using a paraboloid shaped object. The opening angle, β, was changed by choosing different radii of the circular base and simulations were performed with β = 15 ^◦ ,30 ^◦ ,45 ^◦ ,60 ^◦ and 75 ^◦ . For a large opening angle, β = 75 ^◦ , the results al- most coincided with the results for an infinite upper wall for corresponding values of (P ^ ,ε ^ ). In the case of a sharper tip, a soft mushroom has the pos- sibility to perform an escape transition to avoid the main part of the obstacle, resulting in a almost unchanged radius of gyration (see figure 4.7). For a hard mushroom a sharper tip will have an increased influence on the size and shape during the compression, due to the more compact conformations which will not as easily avoid the tip.

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1

0.5 1 1.5 2 2.5 3 3.5 4

Rg / nm

L / nm

ε^´=0.1 ε^´=0.25

ε^´=0.5

Figure 4.7: Radius of gyration vs. inter-wall distance for soft mushrooms compressed

by a paraboloid shaped tip. The solid lines correspond to results using a slow com-

pression rate and the dashed lines are results from a fast compression. The upper two

curves for every ε

are results from simulations with a sharp tip ( β = 15

) and the

lower two curves for every ε

are from simulations with β = 75

.

4.2 A single polymer brush

In this section we study the dynamic behaviour during compression and re- laxation of a polymer brush for different surface coverages. The non-bonded intra- and inter-chain interactions were modeled by an excluded volume in- teraction using the repulsive part of a Lennard-Jones potential with ε ^ = 1.0.

In addition, we also studied the effect of including the full Lennard-Jones po- tential for a selected set of chain parameters. The dynamics of the polymer layer was studied following 15 grafted chains and applying periodic boundary conditions in the two directions parallel to the confining surfaces. The results below are from simulations of chains with N = 20 and for different surface coverages, σ, which were obtained by varying the length of the simulation box in the periodic directions. The upper surface was introduced at a posi- tion determined, in the same way as for the simulations of a single polymer mushroom, from the equilibrium maximum z-coordinate in the ensemble for a given set of parameters (P ^ ,ε ^ ) and the chains were compressed until an inter- wall distance comparable with the persistence length. In figure 4.8 we can see a typical initial structure of a polymer brush before the compression.

L´

z x y

Figure 4.8: Snapshot from one of the simulations of a typical polymer brush before the compression.

4.2.1 Role of surface coverage during compression

We first consider the results for polymer brushes compressed by an upper

wall, where the interactions between the beads are modeled by the repulsive

part of a Lennard-Jones potential. In figure 4.9 the radius of gyration and the

asphericity are shown as a function of the inter-wall distance L ^ for chains

with P ^ = 2 and for three different σ. We note that the uncompressed size of a chain and the prolateness increase with a higher σ as expected due to more stretched equilibrium conformations. In the initial stage of the compression, the radius of gyration decreases in the same way as for a polymer mushroom and for a low surface coverage a chain can increase its size in the lateral direc- tions at the end of the compression. For a more dense layer the size decreases monotonically due to the more limited space to swell parallel to the confining surfaces in this case. At the same time the asphericity goes through a min- imum for all surface coverages during the compression, with the minimum appearing somewhat earlier for smaller σ.

2.2 2.4 2.6 2.8 3 3.2 3.4 3.6

0 2 4 6 8 10 12 14 16

Radius of gyration, Rg

Inter-wall separation, L′ P=2′ (with excluded volume

interactions) σ=0.025 σ=0.05 σ=0.10

0.4 0.5 0.6 0.7 0.8 0.9 1 1.1

0 2 4 6 8 10 12 14 16

Asphericity, S

Inter-wall separation, L′

P′=2 (with excluded volume

interactions) σ=0.025 σ=0.05 σ=0.10

Figure 4.9: Radius of gyration (given in terms of bead diameter) and asphericity vs.

scaled inter-wall distance during compression for a polymer brush for three different surface coverages with P

= 2 and the beads interacting by excluded volume forces.

In figure 4.10 the intra- and inter-chain entanglement are shown as a function of inter-wall distance. The equilibrium inter-chain entanglement is smaller compared to the intra-chain entanglement for chains with corresponding surface coverage, again showing the stretching of chains in a polymer brush. Now, following the compression, the inter-chain entanglement goes through a small maximum for all σ which indicates that the chains have only a limited interpenetration during the compression.

The increase in  ¯N inter should mainly be due to the distance dependence (proportional to the inverse distance squared) between chain segments belonging to different chains and mirrors the decreased space between the beads.

The equilibrium intra-chain entanglement decreases with higher surface

coverage in line with the results above. During the compression, this order

will eventually reverse, as the chains in a layer with low σ can somewhat

increase their size laterally at the end of the compression.

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

0 2 4 6 8 10 12 14 16

Entanglement,〈¯N〉

Inter-wall separation, L^′

⎧⎨

⎩⎧

⎨⎩ P=2′ (with excluded volume

interactions)

〈^N¯〉intra

〈N^¯〉inter

σ=0.025 σ^=0.05 σ^=0.10 σ^=0.025 σ^=0.05 σ=0.10

Figure 4.10: The intra- and inter-chain entanglement  ¯N vs. scaled distance L

be- tween the surfaces.

We have also looked at the monomer concentration profile for the equi- librium and compressed states. In the left part of figure 4.11 we can see the monomer concentration profile for polymer brushes at equilibrium and we no- tice the increased streching away from the grafting surface for higher σ. The figure to the right shows instead the concentration profile in the final stage of the compression, where we can see a well ordered state with distinct layers at both surfaces.

0 0.5 1 1.5 2 2.5

0 2 4 6 8 10 12 14 16

Concentration

Inter-wall separation, L′

P=2, N=20′ (with excluded volume

interactions) σ=0.025 σ=0.05 σ=0.10

0 1 2 3 4 5 6 7 8

0 0.5 1 1.5 2 2.5

Concentration

Inter-wall separation, L′

P=2, N=20′ (with excluded volume

interactions) σ=0.025 σ=0.05 σ=0.10

Figure 4.11: The initial-(left figure) and final(right figure) monomer concentration profile for chains with N = 20 and P

= 2 for different σ.

4.2.2 Role of solvent quality during compression

To model the effect of solvent quality, we also included the full Lennard-

Jones interaction for a number of simulations. In figure 4.12 the results for a

low ε ^ are shown for chains with P ^ = 2 and N = 20. We can notice the effect

of including a small attraction by a decreased equilibrium radius of gyration together with less prolate structures compared with the results in figure 4.9.

The results during the compression resembles those obtained for excluded volume chains. For low surface coverages the radius of gyration goes through a minimum as there is enough space at the end of the compression to allow for the chains to swell. A larger initial equilibrium radius of gyration associated with a larger surface coverage, will also in this case result in a smaller size in the most compressed state.

2.4 2.5 2.6 2.7 2.8 2.9 3 3.1 3.2

0 2 4 6 8 10 12 14 16

Radius of gyration, Rg

Inter-wall separation, L′ P=2′ ε=0.25′ σ=0.025 σ=0.05 σ=0.10

0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8

0 2 4 6 8 10 12 14 16

Asphericity, S

Inter-wall separation, L′ P=2′ ε=0.25′ σ=0.025 σ=0.05 σ=0.10

Figure 4.12: Radius of gyration (given in terms of bead diameter) and asphericity vs.

scaled inter-wall distance for chains with P

= 2 and ε

= 0.25 for three different σ.

The behaviour for chains with N = 20, P ^ = 2 and ε ^ = 0.75 is shown in figure 4.13. Starting from an even more compact initial structure with smaller R g and S, the chains have in this case a larger resistance to confinement, and there is also a larger difference in shape for different σ.

2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1

1 2 3 4 5 6 7 8 9 10

Radius of gyration, Rg

Inter-wall separation, L′ P=2′ ε=0.75′ σ=0.025 σ=0.05 σ=0.10

0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75

1 2 3 4 5 6 7 8 9 10

Asphericity, S

Inter-wall separation, L′ P=2′ ε=0.75′ σ=0.025 σ=0.05 σ=0.10

Figure 4.13: Radius of gyration (given in terms of bead diameter) and asphericity vs.

scaled inter-wall distance for chains with P

= 2 and ε

= 0.75 and for three different

σ.

4.2.3 Relaxation from the compressed state

After compression, the relaxation was followed using the same two protocols as for a single grafted molecule: (i) the upper surface was instantly removed after the compression and (ii) the chains were first allowed to relax in the com- pressed state before removing the upper surface. In figure 4.14 the time depen- dence of R _g in terms of the dimensionless parameter τ = D o t/a ² is shown on a logarithmic scale, for the compression and immediate relaxation for chains interacting by excluded volume interactions and having a persistence length P ^ = 2. The chains return rapidly to equilibrium upon removing the upper surface for all surface coverages. The relaxation towards equilibrium is faster for chains within a layer with a high surface coverage, reflecting the larger number of interactions, as well as indicating only a small contribution from entanglement effects during the relaxation process.

2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

0.1 1 10 100 1000

Radius of gyration, Rg

Dimensionless time, τ

P=2, N=20′ (with excluded volume

interactions) σ^=0.025 σ^=0.05 σ=0.10

Figure 4.14: Radius of gyration (given in terms of bead diameter) vs. dimensionless time τ = D

t /a

on logaritmic scale following route (i) where the upper wall was immediately removed after the compression for chains with P

= 2 and N = 20.

When letting the chains relax in the compressed state before removing the

upper wall (route (ii), figure 4.15), we can see a drift in R g for the lowest

surface coverage in contrast to for the same process for higher surface cov-

erages, where the chains have a limited possibilty to move due to a restricted

space. The removal of the upper surface will also now trigger a fast relaxation

towards the equilibrium value.

2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

0.1 1 10 100 1000

Radius of gyration, Rg

Dimensionless time, τ P=2, N=20′

(with excluded volume interactions) σ^=0.025 σ=0.05 σ^=0.10

Figure 4.15: Radius of gyration (given in terms of bead diameter) vs. dimensionless time τ = D

t /a

on logarithmic scale following route (ii), where the chains were allowed to relax in the compressed state before the removal of the upper surface, for chains with P

= 2 and N = 20.

4.3 Two approaching polymer brushes

In this section we follow the dynamic behaviour of two approaching polymer brushes, each having the same surface coverage, using the same model for the polymer chains as above. We performed the simulations for two different per- sistence lengths using excluded volume forces modeled by the repulsive part of a Lennard-Jones force. Periodic boundary conditions were applied in the x- and y directions and the surface coverage, σ, was varied using different sizes of the simulation box. In addition to σ = 0.05 and σ = 0.10 we also performed simulations with σ = 0.15, to investigate the effect of a large crowding. The compression started from an inter-wall distance where the equilibrated chains were almost at a touching distance, which was confirmed from the monomer concentration profile as well as by inspection using a graphics program. Each compression step, performed by moving the upper surface downwards in the z-direction, was followed by hundred relaxation steps and this process contin- ued, until a predetermined final distance between the surfaces was reached.

4.3.1 Role of surface coverage during compression

In figure 4.16 we can see the radius of gyration and asphericiy for the chains

on the lower surface as a function of the inter-wall distance, L ^ . Qualitatively

the results are in line with the results above for a compression with a bare

upper surface and we can see that the effect from the polymer chains grafted

to the upper wall is to bring forward the effect of compression. The initial

size increases with surface coverage, and during the compression, this relation