Lecture Notes on Polymer Physics

Angelo Rosa∗
SISSA - Scuola Internazionale Superiore di Studi Avanzati, Via Bonomea 265, 34136 Trieste (Italy)
(Dated: February 9, 2015)
In these notes, we present some amongst the most general and popular models used to describe the
physics of polymer equilibrium and dynamics. First, we analyse in detail static (equilibrium) models
describing single-chain properties. We start from those models (FJC, GC, WLC) which admit an
exact (or almost exact!) solution, and then we continue by discussing the properties of self-avoiding
polymers which – being very difficult to treat analytically – provide a nice example of the power of
Monte Carlo numerical simulations. Second, we move to polymer dynamics by discussing in depth
the Rouse model and its analytical solution. Finally, we present the basics of Molecular Dynamics
which constitutes a powerful computational tool to explore polymer dynamics (but, equilibrium
too!) in those cases where – again! – polymer models fail any simple analytical attempt. The whole
material presented here was assembled from a representative (but clearly, not exhaustive) selection
of the available literature (please, see the bibliography at the end for reference).

PACS numbers:

Contents I. SOME DEFINITIONS TO START WITH

I. Some definitions to start with 1 Let us start by summarising some standard notation,
see also Fig. 1.
II. Equilibrium. Part I: Ideal chains 3
A. The Freely-Jointed Chain (FJC) model 3 • Contour length, L. Also called curvilinear length,
B. The Gaussian Chain (GC) model 5 is the distance measured by an imaginary walker
1. The Gaussian chain and the analogy to who moves along the chain from one of its ends to
heat equation 6 the other. It is expressed in units of [Length].
C. Semiflexibile polymers: The Kratky-Porod
Chain (KPC) model 8 • Kuhn length, lK . The Kuhn length is a measure of
D. Semiflexible polymers: The Worm-like Chain chain rigidity. It can be described as the shortest
(WLC) model 10 segment along the polymer backbone above which
thermal fluctuations start bending the polymer sig-
III. Equilibrium. Part II: Real chains 12 nificantly. As a matter of fact, on large-scales, a
A. Flory theory for self-avoiding polymers 13 polymer can be described as a sequence of N sta-
B. Monte Carlo simulations for self-avoiding tistically independent contour length units, each
polymers 14 unit having length = lK . Mathematically, it can
2
1. General principles 14 defined as limL→∞ hR L(L)i , where hR2 (L)i is the
2. Choice of MC moves 15 mean-square end-to-end distance of the chain, see
3. Model interactions 15 Eq. 1.
4. Averages & system equilibration 16
• Persistence length, lp . The persistence length is
IV. Dynamics 17 a measure of the rigidity of a semi-flexible poly-
A. The Rouse model 17 mer chain. It corresponds to the length-scale along
1. General solution 17 the polymer contour above which the chain looses
2. Time correlations for X~ 0 (t) 20 memory of its initial direction. Mathematically (see
~
3. Time correlations for Xp (t) 20 Sec. II C, Eq. 30), it can be defined through the
4. Time correlations for ~rn (t) 20 average value hcos θ(s)i ≡
exp(−s/lp ) [1], where
B. Molecular Dynamics simulations 21 θ(s) = cos−1 t̂(s0 + s) · t̂(s0 ) is the angle between
1. General principles & simple algorithms 21 tangent vectors to the polymer chain taken at fixed
2. Force-field parametrization 23 contour length separation, s (see Fig. 1). For suf-
ficiently long chains (L  lp ), the concepts of per-
A. The Langevin equation 24 sistence length and Kuhn length are equivalent to
each other, since it can be proven (see Sec. II C,
References 26 Eq. 32) that lK = 2lp .

• Mean-square end-to-end distance, hR2 (L)i. This

quantity (its square root, in fact) is commonly used
∗ Electronic address: anrosa@[Link] as a measure of a polymer chain average size and
 2

-$ !")$ !"+,%$ is defined as the square distance between the chain

!"'$ )&# ends (see Fig. 1), averaged over all possible spatial
!"&$ !"*$
conformations that the chain may assume owing to
)%# !"+$
!"($ random thermal fluctuations:
!"%$ (#
)$#
!%# !&'$#
!&#
!"#$ !$#

!"#
hR2 (L)i = h(~rL − ~r0 )2 i . (1)

0$
)%# ."/$
)$# .".&$
.".%$
(# )+#
)"#
!*%# !+#
."#$ !*$#

!"#
• Mean-square gyration radius, hRg2 (L)i. Along with
the mean-square end-to-end distance, this is an-
other quantity routinely employed to characterize
FIG. 1: Alternative representations of a model polymer chain. the average polymer size. It is defined as the av-
(A) A polymer chain can be described as a sequence of N + 1
erage square distance between each monomer and
beads, indexed as i = 0, 1, ..., N − 1, N and connected by rigid RL
joints of constant length b. By construction then, the total the center of mass of the chain, ~rcm ≡ L1 0 ds ~rs :
contour length is L = N b. The spatial vector pointing at bead
i is ~ri , with end-to-end vector R ~ ≡ ~rN − ~r0 . For convenience,
we define the oriented bond vector between ~ri−1 and ~ri as
~ti ≡ ~ri − ~ri−1 and the normalised oriented bond vector t̂i ≡
~
ti
b
. (B) Alternatively, a polymer chain can be described as
a linear, continuous string of total contour length, L. This
model can be considered as the continuum version of model 1 L
Z
2
(A) in the limits b → 0, N → ∞ provided L = N b = constant. hRg2 (L)i ≡ ds h(~rs − ~rcm ) i, (2)
Monomer position along the chain is given by the curvilinear
L 0
As for the mean-square end-to-end distance (Eq. 1),
coordinate s ∈ [0, L], with corresponding spatial coordinate ~rs
~ ≡ ~rL −~r0 . The (normalised) tangent hRg2 (L)i is averaged over all possible conformations
and end-to-end vector R
vector t̂s is defined as t̂s = |∂~rs1/∂s| ∂~
rs
.
of the chain. Eq. 2 can be recast in a different form,
∂s
which is sometimes more useful:
 3

!2
L
1 L 1 L 0
Z Z Z
1 2
hRg2 (L)i = ds h(~rs − ~rcm ) i = ds h ~rs − ds ~rs0 i
L 0 L 0 L 0
!
1 L
Z Z L Z L Z L
2 2 0 1 00 0
= h ds ~rs − ~rs · ds ~rs0 + 2 ds ~rs00 · ds ~rs0 i
L 0 L 0 L 0 0

1 L
Z Z L Z L Z L Z L
2 1
= h ds ~rs2 − 2 ds ~rs · ds0 ~rs0 + 2 ds00~rs00 · ds0 ~rs0 i
L 0 L 0 0 L 0 0
1 L
Z Z L Z L
1
= h ds ~rs2 − 2 ds ~rs · ds0 ~rs0 i
L 0 L 0 0
Z L Z L Z L Z L
1 1 1
= h ds ~rs2 − 2 ds ~rs · ds0 ~rs0 + ds0 ~rs20 i
2L 0 L 0 0 2L 0
Z L Z L Z L Z L Z L Z L
1 0 0
= h ds ds ~rs
2
− 2 ds ~
rs · ds ~
rs 0 + ds ds0 ~rs20 i
2L2 0 0 0 0 0 0
Z L Z L
1
ds0 ~rs2 − 2~rs · ~rs0 + ~rs20 i
 
= 2
h ds
2L 0 0
Z L Z L Z L Z L
1 0 2 1 2
= ds ds h(~rs − ~rs0 ) i = 2 ds ds0 h(~rs − ~rs0 ) i. (3)
2L2 0 0 L s=0 s

II. EQUILIBRIUM. PART I: IDEAL CHAINS which corresponds to the classical result for a random-
walk (RW). This is not surprising, since the FJC as-
A. The Freely-Jointed Chain (FJC) model sumes that cross-correlations between different bond vec-
tors are zero. An important consequence of Eq. 5 is
The FJC model is one of the basic models for describ- that the Kuhn length, lK , of a FJC is (Sec. I) lK =
2
ing the properties of ideal (i.e., with no excluded vol- limL→∞ hR L(L)i = b.
ume) polymers at equilibrium. A typical FJC consists
of N + 1 beads linearly connected one after the other by
rigid joints of equal length = b, see Fig. 1. The FJC
model assumes, that the spatial orientation of each bond
is independent from the orientations of the other bonds.
Within this assumption, the configurational space of a
FJC is described by the following statistical distribution:
N N
Y Y 1
ΨF JC ({~t}) ≡ ψF JC (~ti ) ≡ 2
δ(|~ti | − b), (4)
i=1 i=1
4πb

where: ~ti = ~ri+1 − ~ri and δ(x) is the Dirac δ-function.

The typical size of a FJC is described by the correspond-
ing mean-square end-to-end distance (Eq. 1), hR2 (L =
N b)i ≡ h(~rN − ~r0 )2 i given by:
N
!2
X
2 2 ~ti i
hR (L = N b)i = h(~rN − ~r0 ) i = h
i=1
N X
X N N X
X N
= h ~ti · ~tj i = h~ti · ~tj i
i=1 j=1 i=1 j=1
N
A second quantity describing the typical size
=
X
h~t2i i + 2
X
h~ti · ~tj i of a polymer is the average square gyration ra-
PN 2
i=1 i<j
dius hRg2 (L = N b)i ≡ N1+1 i=0 h(~ri − ~rcm ) i =
1
P N −1 P N 2
N
X (N +1)2 i=0 j=i+1 h(~
ri − ~rj ) i (see Eqs. 2 and 3). By
= b2 + 0 2
using Eq. 5, h(~ri − ~rj ) i = b2 |j − i|, and
i=1
= N b2 = L b, (5)
 4

N −1 XN N −1 X
N N −1 N −i
b2 X 2 b2 X b2 X X
hRg2 (L = N b)i = h(~
ri − ~
rj ) i = (j − i) = j
(N + 1)2 i=0 j=i+1 (N + 1)2 i=0 j=i+1 (N + 1)2 i=0 j=1
−1 −1 −1
N N N
!
b2 X (N − i + 1)(N − i) b2 X
2
X
= = (N − i) + (N − i)
(N + 1)2 i=0 2 2(N + 1)2 i=0 i=0
N N
!
b2 b2
 
X
2
X N (N + 1)(2N + 1) N (N + 1)
= i + i = +
2(N + 1)2 i=1 i=1
2(N + 1)2 6 2

N (N + 2)  N Lb
= b2 ' N 1 b2 = . (6)
6(N + 1) 6 6

Quantities as the mean-square end-to-end distance or where the second δ-function constraints the end-to-end
the gyration radius define the typical size of the poly- ~ Eq. 7 is easy to calculate, upon substi-
vector to be = R.
mer. Thermal fluctuations of the typical polymer size tution of the second δ-function by its Fourier representa-
and other higher-order momenta can be captured by the tion:
~ defined as:
end-to-end distribution function, PN (R),
N N
!
~t
Z Y
~ ≡ d X
~ ,
PN (R) δ(|~ti | − b)δ ~ti − R (7)
i=1
4πb2 i=1

N
" N
!#
Z Y
d~t
Z
d ~k X
~ =
PN (R) δ(|~ti | − b) exp ik · ~ ~
~ti − R
i=1
4πb2 (2π)3 i=1
N N
!
d~k ~t
Z Z Y
 d X
= exp −i~k · R ~ δ(|~ti | − b) exp i~k · ~ti
(2π)3 i=1
4πb 2
i=1
~ N Z   sin(kb) N
~ d~k
Z Z
dk 
~k · R
~
 dt 
~t| − b) exp i~k · ~t
 
~k · R
~
= exp −i δ(| = exp −i
(2π)3 4πb2 (2π)3 kb

N
d~k 2
d~k 2
Z   Z  
~ 1 − (b) k 2 ~ exp − N b k 2
    
≈ exp −i~k · R ≈ exp −i~k · R

N 1 3 3
 (2π) 6 (2π) 6
 3/2 Z  1/2 !
3 3
= d~k exp −i 2 ~k · R~ exp(−k 2 )
2π 2 N b2 2N b2
 
3/2 1/2 !2
3R2
  Z 
3 3
= exp − d~k exp − ~k + i ~ 
R
2π 2 N b2 2N b2 2N b2
3/2 3/2
3R2 3R2
     
3 3/2 3
= exp − ×π = exp − . (8)
2π N b2
2 2N b2 2πN b2 2N b2

Eq. 8 is one of the central result of these notes. It says bility). Of course, this is a result of the approximations
that, the distribution function of a flexible chain is a leading to Eq. 8.
Gaussian. Yet, we should keep in mind that this result is
only an approximation (although a very good one, in the Elastic response of a FJC – We now study the elastic
“N  1”-limit): in particular, Eq. 8 predicts that large response of a FJC to an external force f~ = f (0, 0, 1) ap-
(R > N b), unphysical end-to-end distances might occur plied at chain ends and oriented along the z-axis. Under
(although, with exponentially-damped negligible proba- these conditions, the statistical distribution of the system
 5

is given by (see Eq. 4): !")$ !"+,%$

!"'$ )&#
N
Y 1 ~ ~ !"&$
ΨF JC ({~t}, f~) ≡ 2
δ(|~ti | − b) ef ·R , (9) )%#
4πb !"*$ !"+$
i=1 !"($
where R~ = PN ~ti . The partition function, ZF JC (f ), of !"%$ (#
i=1
the system is given by: )$#
!%# !&'$#
N !&#
!$#
Z Y
ZF JC (f ) = d~ti ΨF JC ({~t}, f~) !"#$
i=1
N
Z Y
1 ~ ~ !"#
= d~ti δ(|~ti | − b) ef ·R
i=1
4πb2
N
Z Y
1 ~~ FIG. 2: Schematic picture of a Gaussian chain. Each bead in
= d~ti 2
δ(|~ti | − b) ef ·ti
i=1
4πb the linear sequence of N +1 beads (indexed as i = 0, 1, ..., N −
Z N 1, N ) is connected to its neighbours by springs with elastic
~ 1 ~ f~·~
t constant = 3kbB2 T , see Eq. 12.
= dt δ(|t| − b) e
4πb2
 Z +1 N
1 We notice the peculiarity of Eq. 12, which contains an
= d cos θ ef b cos θ
2 −1 explicit dependence on the temperature, and hence it

sinh f b
N should be viewed as an effective free-energy rather than
= . (10) a true energetic term. Since bond orientation is uncorre-
fb
lated from the orientation of the other bonds, the average
The average elongation of the chain hRz i along the di- square end-to-end distance (and gyration radius) of a GC
rection of the applied force f~ is given by: is given by the same expressions of a FJC (Eqs. 5 and
6). Moreover, by repeating the calculation for the end-
∂ log ZF JC (f ) to-end distribution function (Eq. 8) it turns out that the
hRz i = −
∂f resulting Gaussian function provides the exact answer for
a Gaussian chain. There are two more interesting prop-
 
1
= N b coth f b − erties concerning the GC:
fb
 
1 1. Given two beads at contour length position
= L coth f lK − m and n, the probability distribution function,
f lK
 L lK Pm,n (~rm , ~rn ), of the spatial distance between these

 3 f, f lK  1 two beads is given by:
≈   . (11) 3/2
 L 1 − 1 , f lK  1 3(~rm − ~rn )2
  
 3
f lK Pm,n (~rm , ~rn ) = exp − ,
2π|m − n| b2 2|m − n|b2
In particular, the small-force regime is equivalent to the (13)
linear response of a Hookean spring with elastic constant which is even more general than Eq. 8;
= L 3lK .
2. Given 3 monomers l, m, and n:
Z
B. The Gaussian Chain (GC) model Pm,n (~rm , ~rn ) = d~rl Pm,l (~rm , ~rl ) Pl,n (~rl , ~rm ), (14)

Because of the fixed-bond constraint, the FJC model which is a straightforward application of Gaussian
is of little practical interest. A more advantageous integrals.
alternative is given by the GC model. In the GC An important application of the GC formalism is given
model, a polymer chain is represented as a linear ar- by the following problem (adapted from problem 2.38 of
ray of beads connected by harmonic springs (Fig. 2). book [2]):
Its mathematical formulation ~2 is similar to Eq. 4, with What is the (3-dimensional) probability, P (~r0 , ~rN ; ~rn ),

3/2 3ti
~
ψGC (ti ) = 2πb23
exp − 2b2 . So, the bond interac- of finding a monomer of index n along an ideal chain
tion between consecutive beads of coordinates ~ri and ~ri+1 whose ends are fixed positions ~r0 and ~rN ?
(i = 0, ..., N − 1) is given by the harmonic function: By adopting the GC model, P (~r0 , ~rN ; ~rn ) is propor-
tional to:
3 kB T ~2 3 kB T
U (~ti ) = 2
ti = (~ri+1 − ~ri )2 . (12) P (~r0 , ~rN ; ~rn ) ∝ Pn (~rn − ~r0 ) PN −n (~rn − ~rN ), (15)
2b 2b2
 6

where Pn (~rn − ~r0 ) (resp., PN −n (~rn − ~rN )) is given by the expression for P (~r0 , ~rN ; ~rn ) is given by:
corresponding GC expression R (see Eq. 8). By imposing
standard normalisation d~rn P (~r0 , ~rN ; ~rn ) = 1, the final

   3/2  

3/2 r0 )2
n −~ rN )2
rn −~
3
2πn b2
exp − 3(~r2nb 2
3
× 2π(N −n) b2
exp − 3(~ 2(N −n)b2
P (~r0 , ~rN ; ~rn ) =  

3/2 N −~r0 )2
3
2πN b2
exp − 3(~r2N b2
!3/2  2
~rn + ~r02 − 2~rn · ~r0 ~rn2 + ~rN
2
~r02 + ~rN
2
 
3 3 − 2~rn · ~rN − 2~r0 · ~rN
= exp − + −
2π n(NN−n) b2 2 b2 n N −n N
!3/2
3 (N − n)2~r02 + N 2~rn2 + n2~rN 2
 
3 − 2(N − n)N~rn · ~r0 − 2nN~rn · ~rN + 2n(N − n)~r0 · ~rN
= exp −
2π n(NN−n) b2 2 b2 n (N − n) N
!3/2
3 N 2~rn2 − 2((N − n)N~r0 + nN~rN ) · ~rn + ((N − n)~r0 + n~rN )2
 
3
= exp − 2
2π n(NN−n) b2 2b n (N − n) N
3/2
3 ~rn2 − 2((1 − n/N )~r0 + n/N ~rN ) · ~rn + ((1 − n/N )~r0 + n/N ~rN )2
  
3
= exp −
2π n(1 − n/N ) b2 2 b2 n (1 − n/N )
3/2
3 (~rn − ((1 − n/N )~r0 + n/N ~rN ))2
  
3
= exp − , (16)
2π n(1 − n/N ) b2 2 b2 n (1 − n/N )

which is equivalent to the situation
where monomer n is PN,0 (~rN , ~r0 )|boundary = 0 [3]. One can immediately
n
the end of a GC made
n

of nn 1 − N and the other end recognise the traditional form of the HE by substitut-
2
is at position 1 − N ~r0 + N ~rN . In the special case of ing “N ” with “time” and “ b6 ” with the “diffusion coeffi-
~r0 ≡ ~rN (i.e. a circular or ring polymer), Eq. 16 reduces cient”.
to:
!
Polymer chain inside a regular box – As an applica-
r0 )2 tion of the HE formalism, we consider the instructive ex-
3/2
rn − ~

3 3 (~
Pring (~
r0 ; ~
rn ) = exp − .
2π n(1 − n/N ) b2 2 b2 n (1 − n/N ) ample of a GC inside a box with (non-equivalent) sides
(17) “Lx 6= Ly 6= Lz ”. Eq. 19 can be solved by employing the
So, internal distances for Gaussian ring polymers are well-known method of “separation of variables”, i.e. by
given by: assuming that a possible solution is given by:
 n
h(~rn − ~r0 )2 i = b2 n 1 − , (18)
N PN,0 (~rN , ~r0 ) ∝ PN,0 (xN , x0 ) PN,0 (yN , y0 ) PN,0 (zN , z0 ),
and, by using Eq. 3, the corresponding average-square (20)
N
gyration radius hRg2 i = b2 12 = bL with initial and boundary conditions given by
12 , i.e. the radius of
gyration of a Gaussian ring polymer is the same of a PN →0,0 (xN , x0 ) = δ(xN − x0 ) and PN,0 (xN , x0 )|xN =0 =
Gaussian chain of half the contour length. PN,0 (xN , x0 )|xN =Lx = 0, and analogous expressions in
the y- and z-directions. The corresponding solution is
given by:
1. The Gaussian chain and the analogy to heat equation
∞    
2 X pπxN pπx0
Given the close analogy between polymer conforma- PN,0 (xN , x0 ) = sin sin ×
Lx p=1 Lx Lx
tions and random-walks, it is no surprise that the end-  2 2
p π N b2

to-end distribution function, PN,0 (~rN , ~r0 ), of a Gaus-
× exp − , (21)
sian chain with monomer “0” at spatial position ~r0 and 6L2x
monomer “N ” at spatial position ~rN satisfies the (three-
dimensional) heat equation (HE): with analogous solutions for PN,0 (yN , y0 ) and

∂ b2
 PN,0 (zN , z0 ). The partition function of the system
− ∇~rN PN,0 (~rN , ~r0 ) = δ(~rN − ~r0 )δ(N ), (19)
2 is then given by
∂N 6
Z Z
with: (1) “initial” condition PN →0,0 (~rN , ~r0 ) =
δ(~rN − ~r0 ), and (2) “absorbing” boundary conditions Z= d~rN d~r0 PN,0 (~rN , ~r0 ) = Zx Zy Zz , (22)
 7

8Lx P∞ 1
where Zx ≈ π2 p=0 (2p+1)2 = Lx and Πx = Πy =
kB T kB T
Z Z Πz = Lx Ly Lz = V , which is the standard law
Zx = dxN dx0 PN,0 (xN , x0 ) of ideal gas.
∞
(2p + 1)2 π 2 N b2
 
8Lx X 1
= exp − . √
π 2 p=0 (2p + 1)2 6L2x 2. In the opposite limit Lx , Ly , Lz  N b2 ,
(23) only the “p = 0”-term  2 in  Eq. 23 mat-
π N b2
ters, Zx ≈ 8L π2
x
exp − 6L2x , and Πx ≈
The free energy of the system is calculated by F = kB T
 2 2
 2 2

−kB T log Z = Fx + Fy + Fz , and the pressure of the Ly Lz

1 π Nb
Lx + 3L3x ≈ kBV T π 3LNb
2 . Hence, for an
x
polymer chain acting on, say, the x-wall is given by anisotropic box with Lx 6= Ly 6= Lz , the force
acting on each wall is anisotropic as well. This
1 ∂F 1 ∂Fx kB T 1 ∂Zx anisotropy in the pressure due to the anisotropic
Πx = − =− = .
Ly Lz ∂Lx Ly Lz ∂Lx Ly Lz Zx ∂Lx confinement is responsible for many unusual me-
(24) chanical properties of polymers [3].

Two noteworthy limits deserve attention:

√ Polymer chain inside a sphere of radius R – Alterna-
1. Lx , Ly , Lz  N b2 , i.e. the size of the box is tively, Eq. 19 can also be solved in the case of a Gaussian
much bigger than the bulk chain size. In this case, chain inside a rigid sphere of radius R:

h 2 i
Pl=∞ Pm=+l P∞ xli 2
2 1
jl xli rRN jl xli rR0 Ylm (θN , φN ) Ylm
∗
(θ0 , φ0 ) exp − b 6N




R3 l=0 m=−l i=1 [jl+1 (xli )]2 R
PN,0 (~rN , ~r0 ) = h
i ,
b2 N x0i 2
8R3
P∞ 1
π i=1 i2 exp − 6 R
(25)

where: of the distribution functions can be calculated.

• xli is the i-th zero of the l-th spherical Bessel func- In particular, one is interested in the typical square
tion, jl (x): jl (xli ) = 0; end-to-end distance, hR2 (N )i ≡ h(~rN − ~r0 )2 i =
2


• Ylm (θ, φ) is the spherical harmonic with “quantum” 2 hrN i − h~rN · ~r0 i :
∗
numbers l = 0, ..., ∞ and m = −l, +l, and Ylm (θ, φ)
denotes its complex conjugate.
Notice, that Eq. 25 is normalized according to
2
R
d~rN d~r0 PN,0 (~rN , ~r0 ) = 1. From Eq. 25 all momenta 1. The first term, hrN i, gives

h 2 i
x0i 2
Pi=∞ R1 R1
1
dρ j0 (x0i ρ) ρ4 0 dρ0 j0 (x0i ρ0 ) ρ02 exp − b 6N


4πR5 i=1 [j1 (x0i )]2 0 R
2
hrN i = h
i
b2 N x0i 2
4R3
Pi=∞ 1
π i=1 i2 exp − 6 R
h i h 2
i
b N x0i 2
Pi=∞ 1 j2 (x0i )
i=1 (x0i )2 1 − 2 x0i j1 (x0i ) exp − 6 R
= π 2 R2 h
i . (26)
b2 N x0i 2
Pi=∞ 1
i=1 i2 exp − 6 R

Pm=+1 ∗
2. The second term, h~rN · ~r0 i ≡ hrN r0 r̂N · r̂N i = hrN r0 4π
3 m=−1 Y1m (θ, φ) Y1m (θ0 , φ0 )i, gives
 8

hR i2 h 2 i
Pi=∞ 1 1 x1i 2
− b 6N


4πR5 i=1 dρ
[j2 (x1i )]2
0
j 1 (x 1i ρ) ρ 3
exp R
h~rN · ~r0 i = h i
b2 N x0i 2
4R3
Pi=∞ 1

π i=1 i2 exp − 6 R
h i2 h 2
i
b N x1i 2
Pi=∞ 1 1
i=1 [j2 (x1i )] 2 x1i j 2 (x 1i ) exp − 6 R
= π 2 R2 h
i
b2 N x0i 2
P∞ 1
i=1 i2 exp − 6 R
h 2 i
Pi=∞ 1 b N x1i

2
i=1 (x1i )2 exp − 6 R
= π 2 R2 P h
i . (27)
∞ 1 b N x0i 2
2
i=1 i2 exp − 6 R

<R2 HNL>  lK 2
50 polymer is the so called KPC model, which “completes”
the FJC model (see Eq. 4) by adding an energy penalty
40
<R2 HNL> = lK 2 N for bending:
30
radius=2.5 lK
−1
N N
!
radius=5.0 lK
Y 1 X
20
ΨKP C ({~t}) ≡ δ(|~ti | − b) exp Kb t̂i · t̂i+1
radius=7.5 lK
i=1
4πb2 i=1
10

−1
N N
!
Y 1 ~
X
10 20 30 40 50
N
= δ(|ti | − b) exp Kb cos θi
i=1
4πb2 i=1
FIG. 3: Average square end-to-end distance, hR2 (N )i, of a (28)
Gaussian chain made of N Kuhn segments confined inside ~
where t̂ = bt is the normalised bond vector, θi = cos−1 (t̂i ·
rigid spheres of different radii (given in units of Kuhn length,
b), compared to the typical random-walk behavior of a free,
t̂i+1 ) is the angle between consecutive bonds, and Kb is
unconfined Gaussian chain. Notice the plateauing in the case a measure of the bending stiffness of the chain. Notice
of extreme confinement. also, that the temperature dependence of the Boltzmann
factor has been implicitly taken into account through Kb .
KP C
The partition function of the KPC model, ZN , is
Fig. 3 shows the result of numerical summation of Eqs. then given by:
26 and 27 compared to the standard solution hR2 (N )i =
−1 −1
Z NY N
!
b2 N of an unconfined Gaussian chain. X
ZKP C = d cos θi exp Kb cos θi
i=1 i=1
C. Semiflexibile polymers: The Kratky-Porod Z +1 N −1
Chain (KPC) model = d cos θ eKb cos θ
−1
 N −1
Semiflexible polymers form a more general class of 2
= sinh Kb . (29)
polymers. Roughly speaking, a semiflexible polymer at- Kb
tains a rigid, rod-like conformation on short (L  b)
length-scales while is Gaussian-like on large (L  b) Similarly, the correlation function h~tl · ~tm i = b2 ht̂l · t̂m i
scales. A simple mathematical model for a semiflexible (m > l) reads:
 9

−1 −1
Z NY N
!
2 1 X
h~tl · ~tm i = b KP C
d cos θi exp Kb cos θi t̂l · t̂m
ZN i=1 i=1
R Qm−1  P 
m−1
i=l d cos θ i exp K b i=l cos θ i t̂l · t̂m
= b2  m−l
2
Kb sinh Kb
R Qm−1  P 
m−1
i=l d cos θi exp Kb i=l cos θi (t̂l · t̂l+1 )(t̂l+1 · t̂l+2 ) ... (t̂m−2 · t̂m−1 )(t̂m−1 · t̂m )
= b2  m−l
2
Kb sinh K b
R Qm−1  P 
m−1
i=l d cos θi exp Kb i=l cos θi (t̂l · t̂l+1 ) (t̂l+1 · t̂l+2 ) ... (t̂m−2 · t̂m−1 ) (t̂m−1 · t̂m )
= b2  m−l
2
Kb sinh Kb
R Qm−1  P 
m−1
i=l d cos θ i exp K b i=l cos θ i cos θl cos θl+1 ... cos θm−2 cos θm−1
= b2  m−l
2
Kb sinh Kb
R +1 !m−l
2 −1
d cos θ eKb cos θ cos θ
= b 2
Kb sinh Kb
R +1 !m−l
Kb cos θ
−1
d cos θ e cos θ
= b2 2
Kb sinh Kb
 m−l  
2 1 2 b
= b coth Kb − ≡ b exp − (m − l) . (30)
Kb lp

2
By using the identity hR b(N )i PN
We have solved Eq. 30 by noticing that, for each bond = i=1 ht̂2i i + 2 i<j ht̂i ·
P
2
vector, only its projection along the preceding bond is t̂j i (see Eq. 5) and Eq. 30, the mean-square end-to-end
transmitted to the next one. In Eq. 30, we have also in- distance of the KPC model reads:
troduced an important quantity, the so called persistence
length, lp , defined as:

b
lp ≡ −   ≈ |Kb →∞ b Kb . (31)
1
log coth Kb − Kb
 10

N −1 XN
hR2 (N )i
 
X b
= N +2 exp − (j − i)
b2 i=1 j=i+1
lp
N −1 N −i  
X X b
= N +2 exp − j
i=1 j=1
lp
   
−1 exp − b b
lp − exp − lp (N − i + 1)
N
X
= N +2  
i=1 1 − exp − lbp
   
exp − lbp N
X −1 exp − b (N − i + 1)
lp
= N +2   (N − 1) − 2  
b
1 − exp − lp i=1 1 − exp − lbp
   
exp − lbp X N exp − lbp i
= N +2   (N − 1) − 2  
1 − exp − lbp i=2 1 − exp − lp
b

     
exp − lbp exp −2 lbp − exp − lbp (N + 1)
= N +2   (N − 1) − 2   2
1 − exp − lbp 1 − exp − lbp
 
1 + exp − lbp
≈ |N 1  N .
1 − exp − lbp
(32)

We thus conclude, that the long-chain behaviour of and N → ∞ by keeping L = N b constant. From Eq. 31,
l
a Kratky-Porod chain is still random-walk-like
  with a we deduce that the bending stiffness Kb ≈ bp → ∞ as
1+exp − lbp  well, since being a “macroscopic” quantity lp does not
Kuhn length given by: lK =   b ≈ l →∞ 2 lp ,
b
1−exp − lp
p
depend on the microscopic description of the model. We
i.e. the persistence length lp and the Kuhn length lK are obtain the following quantities:
of the same order of magnitude. (1) The mean-square end-to-end distance of the WLC
modelcan be derived from Eq. 32 which gives hR2 (L)i =
l
D. Semiflexible polymers: The Worm-like Chain
2 lp L 1 − Lp (1 − e−L/lp ) . In particular, we notice that
(WLC) model this expression interpolates between the stiff, short-chain
(hR2 (L  lp )i ≈ L2 ) and the Gaussian, long-chain be-
For most practical applications, the KPC model is haviour (hR2 (L  lp )i ≈ 2 lp L).
translated into its continuum version, the so called worm- (2) The Boltzmann weight of the WLC model can be
like chain (WLC) model. The WLC model is obtained adapted from Eq. 28 and is given by the following ex-
from the KPC model by letting simultaneously b → 0 pression:
 11

−1
N N
!
Y X
ΨW LC ({~t}) ∝ δ(|t̂i | − 1) exp Kb cos θi
i=1 i=1
∞
!
Y lp X
= | b → 0 δ(|t̂s | − 1) exp t̂s · t̂s+1
N → ∞ b s=1
Nb = L t̂s
∞
!
Y lp X (t̂s+1 − t̂s )2
∝ δ(|t̂s | − 1) exp − b
2 s=1 b2
t̂s
Z L  2 !
Y lp ∂ t̂s
≈ δ(|t̂s | − 1) exp − ds .
2 0 ∂s
t̂s
(33)

 2
We notice, that the argument of the integral in the weight ∂ t̂(s)
= ρ12 . Hence, the minimum energy configura-
∂s
of the exponential function is directly related to the local
tion for the Hamiltonian in the Boltzmann weight Eq. 33
curvature of the chain.
corresponds to ρ = ∞, i.e. to a straight rod.
To see this, let us assume that the chain forms a
perfect circle of radius ρ and L = 2πρ. In this case,
the spatial coordinates of the point on the chain of Elastic response of a WLC – We now employ Eq. 33
curvilinear coordinate s ∈ [0, 2πρ] are given by ~r(s) = in order to calculate the elastic response of a WLC to a
ρ (cos(s/ρ), sin(s/ρ), 0), while the tangent to the curve
constant external force f~ = f (0, 0, 1) applied at its ex-
is given by t̂(s) ≡ ∂~r∂s(s) = (− sin(s/ρ), cos(s/ρ), 0). tremes [4], and oriented, say, along the z-direction. The
∂ t̂(s) 1
Consequently, ∂s = ρ (− cos(s/ρ), − sin(s/ρ), 0) or corresponding Boltzmann weight is given by:

Z L  2 !
Y lp ∂ t̂s
ΨW LC ({~t}, f~) = δ(|t̂s | − 1) exp − ds + f~ · R
~
2 0 ∂s
t̂s
Z L  2 Z L
!
Y lp ∂ t̂s
= δ(|t̂s | − 1) exp − ds +f · ds t̂s,z
2 0 ∂s 0
t̂s
Z L  2 Z L  2 Z L
!
Y lp ∂ t̂s,⊥ lp ∂ t̂s,z
= δ(|t̂s | − 1) exp − ds − ds +f · ds t̂s,z , (34)
2 0 ∂s 2 0 ∂s 0
t̂s

where t̂s,⊥ and t̂s,z are the components of the tangential we get (not surprisingly!) the typical elastic response of
vector t̂s respectively orthogonal and longitudinal to the a linear spring (see Eq. 11 for f lK  1):
direction of the applied force f~. In particular, we want
to calculate the average elongation hRz i along the force 2 lp L
hRz i = f. (35)
direction. Unfortunately, because of the constraint t̂ = 1 3
imposed in Eq. 34 the problem admits no general exact f lp  1 – In the opposite, large-force limit the elas-
solution. tic response is much less trivial. First, we start by
We study then the problem in the two opposite limits, simplifying Eq. 34 in the following manner. We notice
f lp  1 and f lp  1: that when f → ∞ the norm of the longitudinal com-
f lp  1 – In the small-force limit, the polymer is only ponentqt̂s,z should be very close to 1, so we can expand
slightly perturbed from its typical, Gaussian-like config- t̂s,z = 1 − t̂2s,⊥ ≈ 1− 21 t̂2s,⊥ . At the same time, the term
uration, and the total free-energy of the system is given  2  2
3R2 ∂ t̂s,⊥
by F = 4l −f Rz . By minimising F with respect to Rz , in Eq. 34 ∂∂st̂s
≈ −t̂s,⊥ · ∂s ∼ O(t̂4s,⊥ ), and then
pL
 12

it can be neglected. To summarise, in the large-f limit Eq. 34 simplifies to:

Z L  2 Z L
!
Y lp ∂ t̂s,⊥ f
ΨW LC ({~t}, f~) ≈ exp − ds − · ds t̂2s,⊥ , (36)
2 0 ∂s 2 0
t̂s,⊥

which resembles a standard Gaussian distribution. • The “2” factor, taking into account the two com-
By employing the Fourier representation of t̂s,⊥ = ponents of the vector t̂q,⊥ .
R dq −iqs
2π e t̂q,⊥ , we get:
• The “L” factor, following from the “discrete” rep-
L L resentation of the wave vector qi = 2πL i which is
dq 0 +iq0 s ∗
Z Z Z Z
dq −iqs of practical use when dealing with functional inte-
(1) ds t̂2s,⊥ = dse t̂q,⊥ e t̂q0 ,⊥
0 0 2π 2π grals.
Z L
dq 0
Z Z
dq
We remark, that the functional relationship 1 − hRLz i ∝
0
= ds e−i(q−q )s t̂q,⊥ t̂∗q0 ,⊥
2π 2π 0 1
, is markedly different from the elastic response of a
f 1/2
dq 0
Z Z
dq
= 2π δ(q − q 0 ) t̂q,⊥ t̂∗q0 ,⊥ FJC (1 − hRLz i ∝ f1 , see Eq. 11). To conclude, our final
2π 2π
Z result, Eq. 41, can be put together [4] with Eq. 35 so
dq 2
= t̂ , (37) to “design” the following semi-empirical formula valid at
2π q,⊥ any applied force:
and hRz i 1 1

hRz i
−2
lp f = − + 1− . (42)
Z L  2 Z L  Z 2 L 4 4 L
∂ t̂s,⊥ ∂ dq −iqs
(2) ds = ds e t̂q,⊥
0 ∂s 0 ∂s 2π In particular, Eq. 42 was used to fit the elastic response
Z L Z 2 of DNA filaments to pulling played by atomic force mi-
dq −iqs
= ds e (−iq)t̂q,⊥ croscope [4]. Since the contour length of DNA is a known
2π experimental parameter, Eq. 42 is left with just one un-
Z0
dq 2 2 known fitting parameter, the DNA persistence length lp .
= q t̂q,⊥ . (38)
2π Best fit to experimental data gave ≈ 50 nm ≈ 150 base-
pairs.
By substitution of Eqs. 37 and 38 into Eq. 36, we get:
 Z   
Y dq lp 2 f III. EQUILIBRIUM. PART II: REAL CHAINS
ΨW LC ({~t}, f~) ≈ exp − q + t̂2q,⊥ .
2π 2 2
t̂q,⊥
(39) For simplicity, the models described above (FJC, GC,
The average chain elongation along the z-direction is WLC) neglect an important physical ingredient: the in-
given by: teraction experienced by two distinct monomers which –
even when they are located far along the chain – come
Z L
1
Z L close to each other in space. Clearly, these two monomers
hRz i = h ds t̂s,z i ≈ L − h ds t̂2s,⊥ i can not occupy the same place in space: this excluded
0 2 0
Z Z volume effect is very important, especially in dilute so-
1 dq 2 1 dq 2 lutions [2, 3], where it leads to overall polymer swelling
= L− h t̂q,⊥ i = L − h t̂ i
2 2π 2 2π q,⊥ in comparison to the corresponding ideal situation. More
Z
1 dq 1 importantly, excluded volume effects change the scaling
= L− 2×L× (40)
2 2π lp q 2 + f of polymer size (hR2 (N )i or hRg2 (N )i) with the degree of
! polymerisation (N ): in particular, in 3 dimensions the
1 relation hR2 (N )i ∼ hRg2 (N )i ∼ N 2ν with ν > 1/2 holds.
= L 1− 1/2
. (41)
(4 lp f ) In a more quantitative fashion, the total interaction
energy between chain monomers can be described by the
In order to derive the final Eq. 40, we considered: following functional form [3]:
R −ax2 π 1/2 N X N
• The identities e dx = a1/2
and 1 X
R −ax2 2 1 π 1/2
Vrep = kB T v(~rn − ~rm ) , (43)
e x dx = 2 a3/2 . 2 n=0 m=0
 13

where kB T v(~rn − ~rm ) is the two-body-like, short-range ing spatial coordinates ~rm and ~rn . Eq. 43 can be approx-
repulsive interaction between the pair of monomers hav- imated as:

N X N Z Z
1 X
Vrep = kB T d~r d~r 0 v(~r − ~r 0 ) δ(~r − ~rn ) δ(~r 0 − ~rm )
2 n=0 m=0
Z Z N X N
1 X
= kB T d~r d~r 0 v(~r − ~r 0 ) δ(~r − ~rn ) δ(~r 0 − ~rm )
2 n=0 m=0
Z Z
1
= kB T d~r d~r 0 v(~r − ~r 0 ) ρ(~r) ρ(~r 0 )
2
Z
1
≈ kB T v0 d~r ρ(~r)2 , (44)
2

becomes:
Z ∞  
u(r)
v0 = 1 − exp − d~r
0 kB T
Z ∞  
u(r)
= 1 − exp − 4πr2 dr
0 kB T
3 Z rB   
4πrA u(r)
= + 1 − exp − 4πr2 dr
3 rA k B T
3 Z rB
4πrA |u(r)|
≈ − 4πr2 dr
3 rA kB T
3
 
4πrA Tθ
≡ 1− , (45)
3 T
FIG. 4: Generic shape for the pairwise excluded volume in- Rr
teraction between chain monomers: u(r) = +∞ for r < rA , where Tθ ≡ kB3r3 rAB |u(r)| r2 dr, known as the θ-
A
u(r) ≈ −kB T for rA < r < rB (= 1.2 rA ), u(r) = 0 for r > rB .
temperature [3], matches the temperature where two-
body interactions are virtually canceled out. At T = Tθ
chain conformations are almost ideal [3].

A. Flory theory for self-avoiding polymers

A simple approach to estimate the typical size of self-

where we have introduced the monomer density func- avoiding polymer chains (T > Tθ ) under the influence
PN of excluded volume interactions was provided approxi-
tion ρ(~r) ≡ r − ~rn ). The final expression in
n=0 δ(~ mately around 60 years ago by Flory [5], see also Ref. [6]
Eq. 44 has been obtained by using the assumption that
for a recent review. In the Flory theory, the typical size
v(~r) is short-range: contributions to the integral are
R(N ) of a polymer chain made of N monomers and in d
non-negligible for monomer pairs close to each other in
dimensions can be obtained by minimisation of the fol-
space. The excluded-volume parameter v0 = v(0) has
lowing free energy FF lory :
units of [volume] and
 it is related to the virial coefficient
R∞h u(r)
i
v0 = 0 1 − exp − kB T d~r [3], where u(~r) = u(r) FF lory N2 R2
≈ v0 d + 2 , (46)
is the two-monomer interaction which depends only the kB T R N lK
pair distance, r, between the monomers. The typical
shape of u(r) is summarised in Fig. 4: it is repulsive where we (and we will systematically) neglect numerical
(u(r) = +∞) for r < rA , followed by an attractive hole prefactors of order 1. The first term in Eq. 46 is of energy
(u(r) ≈ −kB T ) for rA < r < rB , and u(r) = 0 for origin for it follows directly from our result Eq. 45, while
r > rB . Accordingly, the excluded volume parameter v0 the second term is entropic being it simply the entropy
 14

of a Gaussian chain (= kB T log PN (R), see Eq. 8). Upon or rejected according to some rule (which needs to be
∂FF lory specified), implying that St+1 = Snew or St+1 = Sold
minimisation of Eq. 46 with respect to R, ∂R (R =
RF lory ) = 0, we get: respectively. The process continues then iteratively, un-
til a “fair” thermodynamic exploration of the system is
1

d
 d+2
3
achieved.
2
RF lory = v0 lK N d+2 . (47) The main goal of MC is to sample the state space of
2
the system according to the equilibrium Boltzmann dis-
−βH(S)
The Flory theory predicts then a dimension-dependent tribution, π(S) = e Z , where β = kB1T , H is the
critical exponent ν = νF lory (d) = d+2 3
≥ 12 for d ≤ 4. Hamiltonian of the system and Z is the partition func-
For d > 4, νF lory is < 1/2 which is clearly an absurd tion. The motion in the state space can be described
as it would imply that excluded volume effects tend to as a Markov chain, and its evolution obeys the standard
compact the polymer. Hence, we must conclude that master equation [7]:
d = 4 is the critical dimension of self-avoiding polymers
d π(S 0 , t) X
and that ν = 1/2 for any d > 4. The same conclu- = [Γ(S → S 0 )π(S, t) − Γ(S 0 → S)π(S 0 , t)] ,
sion becomes even more transparent by replacing the dt 0
S6=S
Flory radius RF lory (Eq. 47) inside Eq. 46 which gives: (48)
2
2
where Γ(S → S 0 ) is the transition probability to move
  d+2 4−d
FF lory 2

d
d+2 v0
(R = R F lory ) = + 1 N d+2 .
kB T d 2 d
lK from state S to state S 0 . The transition probability sat-
Again, this equation is consistent as far as d ≤ 4, while isfiesPthe following conditions: (1) Γ(S → S 0 ) ≥ 0 and
for d > 4 it would predict an unphysical free energy de- (2) S 0 Γ(S → S 0 ) = 1.
creasing with N . By proper choice of transition probabilities, the goal
To conclude, let us ask: How good the Flory theory is? of MC is to reach a stationary distribution, namely
Indeed, it turns out to be extremely accurate, at least in d π(S,t)
= 0. To this purpose, we ask then two things
dt
the prediction of the critical exponent ν. In fact, νF lory to a MC simulation:
is exact in 1, 2 and 4 dimensions while νF lory (d = 3) =
0.6 is very close to the accepted numerical estimate of 1. Ergodicity – Every state must be reachable starting
0.588... [3]. from any other state by employing a finite num-
ber of moves. This condition ensures, that no sub-
region of the configuration space of the system is
B. Monte Carlo simulations for self-avoiding artificially removed from sampling.
polymers
2. Detailed balance – This is not a necessary condi-
In Statistical Physics, numerical techniques based on tion, although it is of very practical implementa-
Monte Carlo (MC) simulations are typically employed tion for most situations. Detailed balance implies
in order to provide an approximate solution to prob- that transition probabilities must be chosen accord-
lems which do not admit an exact one (see Walter and ing to the simple rule: Γ(S → S 0 )π(S) = Γ(S 0 →
Γ(S→S 0 ) π(S 0 ) −β (H(S 0 )−H(S))
Barkema [7] for a pedagogical review about MC meth- S)π(S 0 ), or Γ(S 0 →S) = π(S) = e . In
0
ods). To give an example, the Ising model in d > 2 was this case, Eq. 48 implies that d π(S ,t)
= 0 at equi-
dt
(and still is!) extensively studied by MC techniques [7]. librium. To simplify the problem further, we de-
As discussed in previous section, self-avoiding polymers compose the transition probability as Γ(S → S 0 ) =
do not admit an exact solution either. For this reason, τ (S → S 0 ) α(S → S 0 ), where (1) τ (S → S 0 ) is
they also have been the object of extensive numerical in- the trial probability to move from S to S 0 and (2)
vestigations using MC. Here, we will review briefly some α(S → S 0 ) is the corresponding acceptance proba-
MC techniques developed for the specific realisation of a bility. Then, the detailed balance condition reads:
self-avoiding polymer (1) on lattice (also conventionally
known as a self-avoiding walk (SAW), see the beautiful τ (S → S 0 ) α(S → S 0 ) 0

review by Sokal [8]), (2) off lattice, i.e for polymers em- 0 0
= e−β (H(S )−H(S)) . (49)
τ (S → S) α(S → S)
bedded in the conventional 3d-space.
Detailed balance allows for rescaling acceptance
probabilities by a common factor so to to increase
1. General principles performance. Nonetheless, being probabilities they
can not exceed 1. By choosing then the largest
A conventional MC simulation consists of a trajec- of the two acceptances, say α(S 0 → S), equal
tory in the configurational space of the system: being to 1, the other must be fixed to α(S → S 0 ) =
St = Sold the state of the system at MC “time” t, the τ (S 0 →S) −β (H(S 0 )−H(S))
τ (S→S 0 ) e .
typical MC strategy consists in moving to a new state
Snew obtained from a slight modification of the previ- These conditions are summarised into the so called
ous state Sold . Then, the new state Snew is accepted Metropolis algorithm [7], where acceptance probabilities
 15

are given by:

τ (S 0 → S) −β (H(S 0 )−H(S))
 
α(S → S 0 ) = min 1, e .
τ (S → S 0 )
(50)
We remark the fact, that the Metropolis algorithm is
not the only possible one, as other MC schemes can be
devised, see Ref. 50. It is one of the most (if not, THE
most!) popular one though, in particular it is commonly
employed in computational polymer models. For these
reasons, in this notes we will concentrate uniquely about
it.
An illustration of the Metropolis algorithm is given in
the following. To fix the ideas, let us assume that the FIG. 5: Example of a 16-step self-avoiding walk (in red) em-
system is found in the configuration St = S at MC time bedded on a regular square lattice (in black). Self-avoidance
t. Then: means that distinct monomers of the walk occupy distinct
sites of the lattice. This features introduces effective long-
1. Switch from configuration S to S 0 . range correlations along the walk, which make the problem
impossible to be solved analytically.
2. Calculate the quantity q(S → S 0) =
τ (S 0 →S) −β (H(S 0 )−H(S))
τ (S→S 0 ) e .

3. Extract a real random number rn ∈ [0, 1].

4. If rn < q(S → S 0 ) update the system to the new !"#$%"&'()"$*+$&,(-).$(

configuration, i.e. St+1 = S 0 , otherwise St+1 = S.

5. Repeat the procedure.

accepted configurations
6. Check that the ratio is
number of trials /0&1"($",()"$*+$&,(-).$(
large enough, say ≈ 50%.

2. Choice of MC moves

In a MC simulation of a self-avoiding polymer, a careful

implementation of moves is crucial. In general, moves can !"#$%"&'(#2)*+$&,(-).$3(
be of two kinds: (1) local, if they tend to modify small
(compared to the rest of the chain) internal portions of
the chain, and (2) non-local, if they modify substantially
the typical shape of the chain.
On-lattice polymers – For the special case of a self-
avoiding polymer of the regular cubic lattice (in d- FIG. 6: Examples of one- and two-bead local moves for SAW
dimensions), a so-called self-avoiding walk (SAW, see on the square lattice. Generalizations to three- (or more) bead
Fig. 5 for an example of a 2-dimensional SAW on the moves and in higher dimensions are straightforward. Figure
square lattice), examples of local and non-local moves adapted from [8].
are provided in Figs. 6 and 7 respectively.
Off-lattice polymers – Examples of local (crankshaft)
and non local (pivot) moves for a typical, off-lattice poly- accepted (or, rejected) according to the rules of the
mer model made of N rigid joints are illustrated in Figs. 8 Metropolis algorithm (see previous Sec. III B 1).
and 9, respectively. In general, the total energy of a single configuration S,
H(S), can be expressed as the sum of contributions from
PN −1 PN
all monomer pairs as H(S) = i=0 j=i+1 Uint (di,j )
3. Model interactions where di,j is the spatial distance between monomer i and
monomer j. Which is (are) the possible choice(s) for
After performing a MC move the energy of the ob- modelling the interaction energy Uint (d) between pairs of
tained new configuration is compared to the energy of monomers at distance d? Several options have appeared
the old configuration, and the new configuration is then in the literature.
 16

!"+,%$
-$ !"'$
!")$
!"&$
!"($ !"*$ !"+$
!" !" !"%$

FIG. 7: Example of a non-local pivot move (here, a 90◦ rota- !"#$

tion) for a SAW on the square lattice. The pivot monomer is
indicated with an ×. Dashed lines indicate the proposed new !"+,%$
segment. Figure adapted from [8]. .$ !"'$
!")$
!"&$
!"($ !"*$ !"+$
!"+,%$ !"%$
-$ !"'$
!")$
!"&$
!"+$ !"#$
!"($ !"*$
!"%$ !"($
!"+,%$
/$ !"'$
!")$

!"#$ !"*$
!"+$
!"+,%$
.$ !"'$
!")$ !"&$
!"&$ !"%$
!"($ !"+$ !"#$
!"*$
!"%$
FIG. 9: Pivot moves for an off-lattice polymer model made
of N rigid joints. The move consists in: (A) Selecting ran-
!"#$ domly a monomer inside the chain. (B) Rotating the short-
est of the two subchains attached to this monomer by an
!"($ angle randomly chosen inside the interval [0, 2π] around a
!"+,%$ randomly-oriented vector. The obtained new configuration
/$ !"'$
!")$
(C) is then accepted according to the standard Metropolis
!"&$ !"*$ rule, see Sec. III B.
!"+$
!"%$
For generic, off-lattice self-avoiding polymers monomer
pair distances d are not restricted to integer values, con-
!"#$ sequently the choice for Uint (d) is typically less straight-
forward. In general, the two options shown in Fig. 10
FIG. 8: Crankshaft moves for an off-lattice polymer model
look like reasonable choices.
made of N rigid joints. The move consists in: (A) Selecting
a random, small portion of the chain made of m consecutive
joints (m up to 5 − 6 units are typically used). (B) Rotat-
ing it by an angle randomly chosen inside the interval [0, 2π]
around the vector connecting the two end monomers. The 4. Averages & system equilibration
obtained new configuration (C) is then accepted according to
the standard Metropolis rule, see Sec. III B.
Let us suppose that we are now in the ideal situation of
having performed a long enough MC simulations which
allowed us to obtain a sequence of M statistically inde-
To what it concerns the self-avoiding walk on lattice, pendent configurations, {S1 , S2 , ..., SM }. Then, the aver-
the choice is pretty straightforward: the new configura- age value hA(S)i for a generic observable A(S) is given
tion is rejected whenever a single overlap between two by:
monomers is produced, otherwise it is accepted [8]. This
corresponds to a simple point-like hard core repulsion, M
equivalent to an interaction energy: Uint (d = 0) = +∞ 1 X
hA(S)i = A(Sm ) , (51)
and Uint (d > 0) = 0. M m=1
 17

IV. DYNAMICS

A. The Rouse model

1. General solution

The Rouse model [2, 3, 9] can be described as the

dynamic version of the Gaussian model (see Sec. II B),
and for its conceptual importance in polymer physics
can be fairly called THE standard model of polymer
dynamics. As in the Gaussian model (see Eq. 12), in
the d-dimensional Rouse model the polymer chain is de-
scribed as a linear array of N + 1 beads with coordinates
~rn (n = 0, 1, ..., N ) and connected by Gaussian springs,
FIG. 10: Examples of model interactions for MC simulations
the full Hamiltonian H describing the system being then
of off-lattice self-avoiding polymers, as functions of the spa-
tial distance d (expressed in units of the bond distance, [b])
given by:
between pairs of monomer: (1) Uhard−core / kB T = +∞ if
d ≤ b and 0 otherwise; (2) ULJ / kB T = 4 [(b/d)12 − (b/d)6 + N −1
1/4], d ≤ 21/6 b. k X 2
H= (~rn+1 − ~rn ) , (54)
2 n=0

with variance:
d kB T
with k = b2 .
1
hA(S)2 i − hA(S)i2 .


varA = (52) By assuming [3, 9] the same friction coefficient, ζn ≡ ζ,
M
for all monomers, the time evolution of spatial positions
Of course, the next question now is: How can we build ~rn (t) is described by the following set of N + 1 coupled
statistically independent configurations? In general, this Langevin equations (see Appendix A for a brief introduc-
can be one of the most difficult aspect of Monte Carlo tion to the Langevin equation for a single particle):
simulations [8]. Nonetheless, a good rule-of-thumb is the
following: r˙0 (t) = −k[~
ζ~ r1 (t)] + f~0 (t)
r0 (t) − ~

1. Given a MC trajectory {S1 , S2 , ..., SM } and the r˙n (t) = −k[2 ~

ζ~ rn (t) − ~ rn−1 (t)] + f~n (t), 1 ≤ n ≤ N − 1
rn+1 (t) − ~
generic observable A(S), calculate the correlation r˙N (t) = −k[~
ζ~ rN −1 (t)] + f~N (t)
rN (t) − ~ (55)
function

where the Gaussian random forces f~n (t) satisfy hf~n (t)i =
M −µ
1 X
CA (µ) = (A(Sm+µ ) − hA(S)i) (A(Sm ) − hA(S)i) .
M − µ m=1
0 and hf~m (t) · f~n (t0 )i = 2 d ζ kB T δmn δ(t − t0 ), where δmn
(53)
and δ(t − t0 ) are the usual Kronecker delta and Dirac
2. Provided the MC trajectory is long enough, CA (µ) delta function, respectively.
is expected to behave as C A (µ)
CA (0) ≈ exp(−µ/µ0 ). µ0
In order to solve Eqs. 55, we look for a new set of
is a quantitative measure of the total number of functions X~ p (t) (p = 0, 1, ..., N ) which: (1) are linear
MC steps needed before the system starts losing combinations of ~rn (t), X ~ p (t) = PN φpn ~rn (t), and (2)
n=0
memory of its initial condition. Otherwise said, a satisfy a set of N + 1, mutually-independent Langevin
correct thermodynamic sample of the system can
equations, ζp X ~˙ p (t) = −kp X ~ p (t)+f~p (t), with correspond-
be obtained by selecting a sufficiently high num-
ing p-dependent frictions ζp and elastic constants kp .
ber of configurations one after the other along the
trajectory and separated by, at least, µ0 MC steps. ~ p (t) and Eqs. 55, we thus have:
From the definition of X
 18

N
~˙ p (t) = ζp
X
ζp X φpn ~r˙n (t)
n=0
−1
" N
# N
ζp X ζp X
= −k φp0 (~r0 (t) − ~r1 (t)) + φpn (2 ~rn (t) − ~rn+1 (t) − ~rn−1 (t)) + φpN (~rN (t) − ~rN −1 (t)) + φpn f~n (t)
ζ n=1
ζ n=0
−1 −2
" N N N
#
ζp X X X
= −k φp0 (~r0 (t) − ~r1 (t)) + 2 φpn ~rn (t) − φp,n−1 ~rn (t) − φp,n+1 ~rn (t) + φpN (~rN (t) − ~rN −1 (t))
ζ n=1 n=2 n=0
N
ζp X
+ φpn f~n (t)
ζ n=0
−1 −1
" N N
ζp X X
= −k (φp0 − φp1 ) ~r0 (t) − φp0 ~r1 (t) + 2 φpn ~rn (t) − φp,n−1 ~rn (t)
ζ n=1 n=2
−2
N
# N
X ζp X
− φp,n+1 ~rn (t) − φpN ~rN −1 (t) + (φpN − φp,N −1 ) ~rN (t) + φpn f~n (t)
n=1
ζ n=0
−1 −1 −1
" N N N
#
ζp X X X
= −k (φp0 − φp1 ) ~r0 (t) + 2 φpn ~rn (t) − φp,n−1 ~rn (t) − φp,n+1 ~rn (t) + (φpN − φp,N −1 ) ~rN (t)
ζ n=1 n=1 n=1
N
ζp X
+ φpn f~n (t)
ζ n=0
−1
" N
# N
ζp X ζp X
= −k (φp0 − φp1 ) ~r0 (t) + (2 φpn − φp,n−1 − φp,n+1 ) ~rn (t) + (φpN − φp,N −1 ) ~rN (t) + φpn f~n (t)
ζ n=1
ζ n=0
(56)
N
X
~ p (t) + f~p (t) = −kp
≡ −kp X φpn ~rn (t) + f~p (t) . (57)
n=0

  
ζ
By equating each term of Eq. 56 to the corresponding with kp = 2 k ζp 1 − cos Npπ +1 =
term of Eq. 57, we get: ζ
 
4 k ζp sin2 2(Npπ+1) . In order to complete the analysis
ζp of the system, we need to calculate next the following
k (φp0 − φp1 ) = kp φp0
ζ average values: (1) hf~p (t)i and (2) hf~p (t) · f~q (t0 )i where
ζp ζ PN
f~p (t) ≡ ζp ~
k (2 φpn − φp,n−1 − φp,n+1 ) = kp φpn , 1 ≤ n ≤ N − 1 n=0 φpn fn (t), see Eqs. 56 and 57. Being
ζ
just a linear combination of f~n (t), the former is = 0.
ζp The latter gives instead:
k (φpN − φp,N −1 ) = kp φpN
ζ
(58)

which corresponds to an eigenvalue problem. The corre-

sponding solution is given by:
 
1 pπ 1


φpn = cos n + /2 , (59)
N +1 N +1
 19

N N
ζp ζq X X
hf~p (t) · f~q (t0 )i = φpn φqm hf~n (t) · f~m (t)i
ζ 2 n=0 m=0
N N
ζp ζq X X
= φpn φqm 2 d ζ kB T δmn δ(t − t0 )
ζ 2 n=0 m=0
N
ζp ζq 0
X
= 2 d kB T δ(t − t ) φpn φqn
ζ n=0
N    
ζp ζq 1 X pπ qπ
δ(t − t0 ) 1 1



= 2 d kB T cos n + / 2 cos n + /2
ζ (N + 1)2 n=0 N +1 N +1
"N   X N  #
ζp ζq 0 1 X (p + q)π 1

(p − q)π 1


= d kB T δ(t − t ) cos n + /2 + cos n + /2
ζ (N + 1)2 n=0 N +1 n=0
N +1
N (p+q)π 1 (p+q)π 1 (p−q)π 1 (p−q)π 1
ζp ζq 0 1 X ei N +1 (n+ /2 ) + e−i N +1 (n+ /2 ) + ei N +1 (n+ /2 ) + e−i N +1 (n+ /2 )
= d kB T δ(t − t )
ζ (N + 1)2 n=0 2
"
ζp ζq 1 (p+q)π 1 − ei(p+q)π (p+q)π 1 − e−i(p+q)π
= d kB T δ(t − t0 ) 2
ei 2(N +1) (p+q)π
+ e−i 2(N +1) (p+q)π
ζ 2(N + 1) 1 − ei (N +1) 1 − e−i (N +1)
#
i(p−q)π −i(p−q)π
i 2(N +1) 1 − e −i 2(N +1) 1 − e
(p−q)π (p−q)π
+e (p−q)π
+e (p−q)π
1 − ei (N +1) 1 − e−i (N +1)
" #
−i(p+q)π −i(p−q)π
ζp ζq 0 1 i 2(N +1) e
(p+q)π − ei(p+q)π i 2(N +1) e
(p−q)π − ei(p−q)π
= d kB T δ(t − t ) e +e
ζ 2(N + 1)2 (p+q)π
1 − ei (N +1)
(p−q)π
1 − ei (N +1)
 

= d kB T
ζp ζq
δ(t − t0 )
1  sin (p + q)π + sin (p − q)π 
ζ 2(N + 1)2 sin 2(N (p+q)π (p−q)π
sin 2(N
+1) +1)

ζp ζq 1  4(N + 1), p = q = 0
= d kB T δ(t − t0 ) 2(N + 1), p = q 6= 0
ζ 2(N + 1)2  0, p 6= q

ζp ζq 1  2, p = q = 0
= d kB T δ(t − t0 ) 1, p = q 6= 0
ζ N + 1  0, p 6= q

(60)

On the other hand, the final expression in Eq. 60 has to be with

= 2 d kB T ζp δpq δ(t − t0 ) which, by comparison, implies:

 (N + 1)ζ, p = 0
ζp = (61)   
 2(N + 1)ζ, 1 ≤ p ≤ N ζp pπ
kp = 2k 1 − cos
ζ N +1
To summarise, by defining the eigenmodes X~ p (t) (0 ≤  
ζp 2 pπ
p ≤ N ) of the Rouse chain as = 4k sin (64)
ζ 2(N + 1)
N  
~ 1 X pπ 1


Xp (t) = cos n + /2 ~rn (t) , (62)
N + 1 n=0 N +1
Eqs. 55 reduce to N independent linear equations for the
modes given by: and corresponding statistical correlations hf~p (t)i = 0 and
hf~p (t) · f~q (t0 )i = 2 d ζp kB T δpq δ(t − t0 ), where ζp is given
~˙ p (t) = −kp X
ζp X ~ p (t) + f~p (t), 0 ≤ p ≤ N (63) by Eqs. 61.
 20

~ 0 (t) ζ
2. Time correlations for X and τ1 ≡ 2 k [1−cos( Nπ+1 )]
≡ τRouse is the so called Rouse
time of the chain [2, 3], corresponding to the longest re-
For p = 0, k0 = 0 (Eq. 64) and X ~ 0 (t) = laxation time of the chain.
1
PN
N +1 n=0 ~
rn (t) (Eq. 62) corresponds to the centre of The general solution to Eq. 67 is given by:
mass of the chain. In this case, Eq. 63 with p = 0 corre-
sponds to the Langevin equation for simple diffusion. Its
Z t
~ p (t) = 1
X
0
e−(t−t )/τp f~p (t0 ) dt0 + X
~ p (0) e−t/τp , (69)
solution is given by: ζp 0
Z t
1
~
X0 (t) = f~0 (t0 ) dt0 + X
~ 0 (0) , (65) where X~ p (0) is the initial condition of X
~ p (t). The cor-
ζ0 0
~
responding time correlation function hXp (t) · X ~ q (0)i is
where X~ 0 (0) is the initial condition of X ~ 0 (t). The p = 0 hence equal to:
correlation function δX02 (t) ≡ h(X ~ 0 (t)− X ~ 0 (0))2 i is hence
~ q (0)i = δpq d kB T e− τp
t
given by: ~ p (t) · X
hX (70)
kp
1 t ~ 0 0 1 t ~ 00 00
Z Z
δX02 (t) = h f0 (t ) dt · f0 (t ) dt i
ζ0 0 ζ0 0 where we have used the relation hf~p (t) · X ~ q (0)i = 0 and
Z t Z t ~ p (0) i =
the equipartition relation hX 2 d k T
kp . Accordingly,
B
1
= 2 dt0 dt00 hf~0 (t0 ) · f~0 (t00 )i the correlation function
ζ0 0 0
Z t Z t
1 ~ p (t) − X ~ p (0))2 i
= 2 2 d ζ0 kB T dt 0
dt00 δ(t0 − t00 ) δXp2 (t) ≡ h(X
ζ0 0 0
= hX~ p (t)2 i + hX ~ p (0)2 i − 2 hX ~ p (t) · X
~ p (0)i
2 d kB T
= t  
ζ0 = 2 hX ~ p (0)2 i − hX ~ p (t) · X
~ p (0)i
≡ 2 d Dcm t , (66)
2 d kB T  − t

= 1 − e τp . (71)
i.e. the motion of the centre of mass is pure diffusive, kp
with Dcm ≡ kζB0T = (NkB T
+1)ζ being the diffusion coeffi-
cient of the chain. Notice that Dcm = N1+1 Dmon where
Dmon = kBζ T is the diffusion coefficient of one single 4. Time correlations for ~rn (t)
monomer.
We complete the analysis of the Rouse model, by con-
sidering time correlations of monomer vectors ~rn (t). The
3. ~ p (t)
Time correlations for X latter can be expressed in terms of Rouse modes as:

For p > 0, kp > 0 (Eq. 64) and Eq. 63 can be written N

X 
pπ

~ 0 (t) + 2 n + 1/2 X ~ p (t) .


as: ~rn (t) = X cos
p=1
N +1
~˙ p (t) = − 1 X
X ~ p (t) + 1 f~p (t) , (67) (72)
τp ζp We consider here two important examples of correlation
functions:
where ~ ~
  (1) φ(t) = hR(t)·
~
R(0)i
hR(0) 2i
~
where R(t) ≡ ~rN (t) − ~r0 (t) is the
π end-to-end vector of the chain. It can be expressed in
ζp 1 − cos N +1
τp ≡ = τ1  , (68) terms of Rouse modes (see Eq. 72) as:
kp 1 − cos Npπ
+1

N      
~
X pπ 1 pπ ~ p (t)
R(t) = 2 cos N+ − cos X
p=1
N +1 2 2(N + 1)
N  pπ   
X pπN ~ p (t)
= −4 sin sinX
p=1
2 2(N + 1)
 
X pπN ~ p (t) .
= 4 (−1)p sin X (73)
p=1,3,5,...
2(N + 1)
 21

By employing correlation relations which have been pre- To conclude then, the monomer mean-square displace-
sented in previous sections, the time correlation function ment is characterised by a cross over from sub-diffusive /
φ(t) finally reads: N -independent (δr2 (t) ∼ t1/2 ) behaviour at short-times

pπ
 to pure diffusive / N -dependent (δr2 (t) ∼ t) behaviour
2 1 + cos N +1 t
 1/2
 e− τp . (74) at long-times, the crossover time being N1 τmon t
X
φ(t) =  ∼
N (N + 1) p=1,3,5,... 1 − cos pπ
N +1 1 → t ∼ τmon N 2 ∼ τ1 .
(2) The single monomer time mean-square displacement:
δrn2 (t) ≡ h(~rn (t) − ~rn (0))2 i. In particular, δrn2 (t) can
be simply written as linear combination of δXp2 (t) (see
Eqs. 66 and 71):
N  
X pπ
δrn2 (t) = δX02 (t) + 4 cos2 n + 1/2 δX 2p (t)


p=1
N +1
(75)
Eq. 75 can be simplified further by averaging over the
monomer index n. We define then
N
2 1 X 2
δr (t) ≡ δr (t)
N + 1 n=0 n
N
( N  )
1 pπ  
2
X X 2 1 2
= δX0 (t) + 4 cos n + /2 δX p (t)
p=1
N + 1 n=0 N +1 B. Molecular Dynamics simulations
N
2
X 2
= δX0 (t) + 2 δX p (t) , (76)
p=1

where we have  used the trigonometric identity 1. General principles & simple algorithms
1
PN pπ
= 12 .
2 1


N +1 n=0 cos N +1 n + /2
It is interesting to analyse in detail the continuum
(N → ∞) limit [2, 3] of the Rouse model. This is equiv-
alent to set (see Eq. 68): In general, MC techniques allow to explore feasibly
2π 2 k 2 2π 2 d kB T 2 the thermodynamical equilibrium of the systems under
kp → p = p (77) study. In many cases, though, we are not just interested
N N b2
ζ ζb 2
2 in equilibrium, we would also like to explore the dynam-
τ1 → 2 N2 = 2 N 2 = 2 τmon N 2 (78) ical properties of the system. In these cases, Molecular
π k π d kB T π
τ1 Dynamics (MD) represents the most natural tool.
τp → 2 (79)
p
2 2
b
where τmon = 2 d D mon
= 2 d kbB T/ is the “microscopic” The general principles of MD are very simple [11]. As
ζ
time-scale characterising diffusion over a length-scale usual, let us consider a system made of N + 1 particles.
equal to the size of a single monomer [10]. According to the principles of classical mechanics [12],
In this limit, φ(t) (Eq. 74) takes the functional form: the time evolution of particle n = 0, 1, ..., N is governed
∞ by the corresponding Newton’s equation of motion:
8 X 1 2 t
φ(t) =2 2
e−(2p+1) τ1 , (80)
π p=0 (2p + 1)
2
and δr (t) (Eq. 76) becomes:
∞  
X 2 d kB T −p2 τt
δr2 (t) = 2 d Dcm t + 2 1−e 1

p=1
kp
d2~rn (t)
∞
m = f~n (t) , (82)
 
2 d kB T N X 1 −p2 τt
= 2 d Dcm t +
π2 k p2
1 − e 1 dt2
p=1
Z ∞  
2 d kB T N 1 −p2 τt
≈ 2 d Dcm t + 1 − e 1 dp
π2 k 0 p2
 1/2 Z ∞
2 d kB T N t 1  2

= 2 d Dcm t + 1 − e−x dx
π2 k τ1 0 x 2

where f~n (t) if the force acting on particle n at time t.

 1/2
2 d kB T N t
= 2 d Dcm t + π 1/2
π2 k τ1 By Taylor expansion, the spatial positions ~rn (t + δt) and
2 d kB T N
 1/2
t ~rn (t − δt) for “small” time increments δt are given by:
= 2 d Dcm t + 3/2
π k τ1
1/2
21/2
  
1 t t
= b2 + b2 1/2
N τmon π τmon
(81)
 22

d~rn (t) 1 d2~rn (t) 2 1 d3~rn (t) 3

~rn (t + δt) ≈ ~rn (t) + δt + δt + δt + O(δt4 ) (83)
dt 2 dt2 6 dt3
d~rn (t) 1 d2~rn (t) 2 1 d3~rn (t) 3
~rn (t − δt) ≈ ~rn (t) − δt + δt − δt + O(δt4 ) (84)
dt 2 dt2 6 dt3

Adding these two expressions together and taking into ac- with ensemble averages: (1) h~ηn (t)i = 0 and (2) h~ηm (t) ·
count Eq. 82 we obtain the following approximate expres- ~ηn (t0 )i = 6 kB T ζ δmn δ(t − t0 ), where δmn and δ(t − t0 )
sion for the “time-forward” particle position ~rn (t + δt): are the Kronecker- and Dirac-δ functions, respectively.
1 ~ Eqs. 89 can be written as:
~rn (t+δt) ≈ 2~rn (t)−~rn (t−δt)+ fn (t)δt2 +O(δt4 ) (85)
m
d~rn (t) m d~vn (t) 1 ~ 1
which is accurate up to O(δt4 ). Instead, by subtracting = − + fn (t) + ~ηn (t) (90)
dt ζ dt ζ ζ
the two equations, we get an approximate expression for
the velocity ~vn (t) ≡ d~rdt
n (t)
: d~vn (t) ζ 1 ~ 1
= − ~vn (t) + fn (t) + ~ηn (t) (91)
~rn (t + δt) − ~rn (t − δt) dt m m m
~vn (t) ≈ + O(δt2 ) , (86)
2 δt Eq. 91 has the following solution:
which is accurate up to O(δt2 ). Numerical implementa-
tion of Eq. 85 and Eq. 86 leads to the so called “Verlet ζ 1
Z t+δt
ζ 0
algorithm” [11]. However, since the numerical solution ~vn (t + δt) = ~vn (t) e− m δt + f~n (t0 ) e− m (t+δt−t ) dt0
of Eqs. 85 and 86 requires the simultaneous knowledge m t
Z t+δt
of the initial position ~r(t = 0) (which is quite natural), 1 ζ 0
+ ~ηn (t0 ) e− m (t+δt−t ) dt0
and of ~r(t = −δt) (which looks indeed less natural), the m t
Verlet algorithm is scarcely employed.
A possible alternative consists in the so called ζ 1 ~  ζ

≈ ~vn (t) e− m δt + fn (t) 1 − e− m δt
“velocity-Verlet” (vV) algorithm, which is described by ζ
the following equations: 1
Z t+δt
ζ 0
+ ~ηn (t0 ) e− m (t+δt−t ) dt0
1 ~ m
~rn (t + δt) ≈ ~rn (t) + ~vn (t)δt + fn (t)δt2 + O(δt3 ) t
2m
δt ~
(87) ≡ ~vn (t) c0 + fn (t) c1 + δ~vnG , (92)
m

f~n (t) + f~n (t + δt) where:

~vn (t + δt) ≈ ~vn (t) + δt + O(δt2 )
2m
(88) ζ ζ
 
1 ζ
2
c0 ≡ e− m δt ≈ 1 − δt + δt + ...
where only the knowledge of initial (i.e. at t = 0) po- m 2 m
sitions and velocities is required. Yet, this scheme is far 1 − c0 1 ζ
  
1 ζ
2
from being perfect either. In fact, the updating of veloc- c1 ≡ ζ ≈1− δt + δt + ...
m δt
2 m 6 m
ity (Eq. 88) requires the knowledge of the “time-forward”
force f~n (t+δt). This does not represent a problem, unless
Z t+δt
1 ζ 0
δ~vnG ≡ ~ηn (t0 ) e− m (t+δt−t ) dt0 (93)
the force itself depends on the velocity. Unfortunately, m t
this is precisely the case for particles coupled to a viscous
environment (where the dissipative force is linear in the The solution for Eq. 90 is instead:
velocity), which is the typical ensemble used in simula-
tions of polymer systems (see Appendix A for the tutorial m
case of one particle embedded in a viscous environment). ~rn (t + δt) = ~rn (t) − (~vn (t + δt) − ~vn (t))
ζ
In this latter case, time evolution of the system is de-
1 t+δt ~ 0 0 1 t+δt
Z Z
scribed by the Langevin dynamics [13]: + fn (t ) dt + ~ηn (t0 ) dt0
ζ t ζ t
d~rn (t)
= ~vn (t)
dt ≈ ~rn (t) + δt ~vn (t) c1
2
d~vn (t) (δt) ~
m = f~n (t) − ζ ~vn (t) + ~ηn (t) (89) + fn (t) c2 + δ~rnG (94)
dt m
 23

G G
where: δrn,α and δvn,α can then be sampled from the bivariate
   2 Gaussian distribution [13]:
1 − c1 1 1 ζ 1 ζ
c2 ≡ ζ ≈ − δt + δt + ...
m δt
2 6 m 24 m
1 t+δt
Z  
ζ 0
δ~rnG ≡ ~ηn (t0 ) 1 − e− m (t+δt−t ) dt0 (95)
ζ t

We notice that each pair of vectorial components of

δ~rnG and δ~vnG , namely δrn,α
G G
and δvn,α (α = x, y, z), are
Gaussian variables with statistical averages:
G
hδrn,α i = 0 (96)

G
hδvn,α i = 0 (97)

ζ ζ
!
G 2 kB T 3 − 4 e− m δt + e−2 m δt
h(δrn,α ) i = δt 2− ζ
ζ m δt
   
2 ζ kB T 3 ζ
≈ δt3 1− δt + ...
3m m 4 m
(98)
kB T  ζ

G 2
h(δvn,α ) i = 1 − e−2 m δt
m    
2 ζ kB T ζ
≈ δt 1− δt + ... (99)
m m m
kB T  ζ
2
G
hδrn,α G
δvn,α i = 1 − e− m δt
ζ
   
2 ζ kB T ζ
≈ δt 1− δt + ... (100)
m m m

( !2 !2 ! !!)
G G G G

G G
 1 1 δrn,α δvn,α δrn,α δvn,α
ρ δrn,α , δvn,α = × exp − + − 2 cr,v
2π σr σv (1 − c2r,v )1/2 2(1 − c2r,v ) σr σv σr σv
(101)

where: σr2 = h(δrn,α

G 2
) i, σv2 = h(δvn,α
G 2
) i, and cr,v = G 2
mean and variance h(δrn,α ) i = 2D δt.
G G
hδrn,α δvn,α i
σr σv . Details for generating a bivariate distri-
bution are provided in Appendix G of Ref. [13]. We 2. Force-field parametrization
just conclude this section by noticing that, in the formal
ζ
limit m → ∞ (the so called over-damped or high-friction
limit), we can just neglect the time evolution of velocities In the following, I present typical functions used to
(Eq. 92), while positions (Eq. 94) evolve according to the describe monomer-monomer interactions in MD simu-
“simplified” dynamics: lations of coarse-grained polymer models. The intra-
polymer interaction energy was introduced in Ref. [14].
δt ~ The same framework was adopted by me in a series of
~rn (t + δt) ≈ ~rn (t) + fn (t) + δ~rnG
ζ works on the modeling of interphase chromosomes [15–
D ~ 17].
= ~rn (t) + δt fn (t) + δ~rnG (102) Quite generally, monomer-monomer interactions con-
kB T
sist of two terms:
where D = kBζ T is the diffusion coefficient, and δ~rnG can (1) Bonded or permanent interactions describe the phys-
be extracted from the Gaussian distribution with zero ical link between nearest neighbour monomers along the
 24

polymer backbone. Interactions modelling the bending

stiffness of the chain belong also to this category.
(2) Non bonded or non permanent interactions are those
interactions at work only when monomers come into
physical proximity as the consequence of their erratic mo-
tion.
This is formalized into the following energy function
Hint :
Hint = Hbonded + Hnon−bonded (103)
where
N
X −1
Hbonded = [UF EN E (i, i + 1) + ULJ (i, i + 1)]
i=0
N
X −2
+ Ubend (i, i + 1, i + 2) (104) FIG. 11: Illustration of the model energy functions UF EN E ,
i=0 Eq. 106 and ULJ , Eq. 108. The sum UF EN E + ULJ is em-
ployed to describe bonded interactions between connected
and
monomers. It displays a minimum at d ≈ 0.96σ.
N
X −1 N
X
Hnon−bonded = ULJ (i, j) . (105)
i=0 j=i+1 100
N= 50
As usual, N + 1 is the total number of monomers consti- N=100

φ(t) = <R(t)R(0)>/<R(0) >

2
tuting the chain, and i and j run over the indices of the
monomers, numbered consecutively along the chain from
0
one chosen end monomer. 10

By taking the nominal monomer diameter = σ, the 10-1

vector position of the ith monomer, ~ri , the pairwise vec- φ(t) 10
-1

tor distance between monomers i and j, d~i,j = ~rj − ~ri ,

and its norm, di,j , the energy terms in Eqs. 104 and 105 10-2
10-4 10-3 10-2 10-1
are given by the following expressions [14]: t / (<R2> N)
1) The chain connectivity term, UF EN E (i, i + 1) is ex- 10-2
pressed as: 101 102 103
  t [τMD]
 2 
k 2 di,i+1
 − 2 R0 ln 1 − R0 , di,i+1 ≤ R0


UF EN E (i, i+1) = FIG. 12: Time correlation function, φ(t) (see Sec. IV A 4), of

 the end-to-end vector R~ = ~rN − ~r0 for 2 linear chains made
0, di,i+1 > R0

of N = 50 and N = 100 bonds. Inset: The two curves
(106) collapse onto a single master curve after rescaling of time
where R0 = 1.5σ, k = 30.0/σ 2 and the thermal en- “t → hR2ti N ”, as predicted by the Rouse model, see [10].
ergy kB T equals 1.0. The complete bonded interaction Chains trajectories were obtained by using the open source
(which includes the Lennard-Jones term, see Eq. 108) is package LAMMPS [18, 19]
shown in Fig. 11 (blue line).
2) The bending energy has the standard Kratky-Porod
form (discretized worm-like chain), see Eq. 28: Fig. 12 shows an example of relaxation dynamics for
two linear chains made of N = 50 and N = 100 bonds
!
kB T l p d~i,i+1 · d~i+1,i+2
Ubend (i, i + 1, i + 2) = 1− whose time evolution was simulated by using the MD
σ di,i+1 di+1,i+2 implementation provided in LAMMPS [18, 19].
(107)
where lp is the nominal persistence length of the polymer
chain, expressed in units of σ. Appendix A: The Langevin equation
3) The excluded volume interaction between distinct
monomers (including consecutive ones, see Eq. 104) cor-
responds to a purely repulsive Lennard-Jones potential: The Langevin equation for a single point particle of
 mass m and space coordinate r in one dimension (corre-
 4[(σ/di,j )12 − (σ/di,j )6 + 1/4], di,j ≤ σ21/6 sponding generalisation to 3 dimensions being straight-
ULJ (i, j) = forward) is:
0, dij > σ21/6

(108) m r̈(t) = −ζ ṙ(t) + f (t) , (A1)
 25

where the stochastic force f is defined through the time where we have made use of the cross-correlation rela-
correlations hf (t)i = 0 and hf (t)f (t0 )i = A δ(t − t0 ). The tion hf (t)v(0)i = 0. In the long-time limit, v(t) loses
former average specifies that f (t) has no preferential di- memory of the initial condition and hv(t)2 i = 2 ζAm . On
rection, while the latter implies that values of f taken at the other hand, at equilibrium, hv(t)2 i = kBmT (which is
different times are – on average – uncorrelated. For the just the standard equipartition theorem of Statistical Me-
moment, we leave the constant A unspecified: we will chanics) which implies A = 2 ζ kB T . From Eq. A4 and
see later in this section that it can be fixed by simple by similar arguments, it is also easy to deduce the veloc-
physical arguments. ζ
ity cross-correlation function hv(t)v(0)i = hv(0)2 i e− m t :
Eq. A1 is a second order differential equation: it can be m
thus, for time-scales t  ζ (a.k.a., the over-damped or
written in the form of a system of 2 first order differential
equations: high-friction limit), one can just neglect velocity corre-
lations and Eq. A1 simplifies to the single, first-order
ṙ(t) = v(t) (A2) differential equation: 0 = −ζ ṙ(t) + f (t).
ζ 1
v̇(t) = − v(t) + f (t) , (A3)
m m
where v(t) is particle velocity at time t. Eq. A3 has the Particle position at time t, r(t), is obtained by time
following solution: integration of Eq. A4 which leads to:
Z t
1 ζ 0 ζ
v(t) = e− m (t−t ) f (t0 ) dt0 + v(0) e− m t (A4)
m 0

where v(0) is particle velocity at time t = 0. In particu-

lar, the mean square velocity hv(t)2 i reads:
t t t0
Z 
1
Z
ζ
(t−t0 )
ζ
(t−t00 )
Z t Z
2
hv(t) i = e
−
m
0
f (t ) dt ×
0
e
−
m
00
f (t ) dt
00
1 ζ 0 00
m2 0 0 r(t) = e− m (t −t )
f (t00 ) dt00 dt0
2 −2
ζ
t
m 0 0
+hv(0) i e m
m  ζ

2
Z t
−
ζ
(t−t0 )
0 0 −
ζ
t

+ v(0) 1 − e− m t + r(0) , (A6)
+
m 0
e mf (t ) dt × v(0) e m ζ
Z t Z t
1 0 00 −
ζ
(2t−t0 −t00 ) 0 00
= dt dt e m hf (t )f (t )i
m2 0 0
ζ
2 −2 t
+hv(0) i e m
Z t
2 0 −
ζ
(2t−t0 ) 0
+ dt e m hf (t )v(0)i
m 0
Z t Z t
A 0 00 −
ζ
(2t−t0 −t00 ) 0
= 2
dt dt e m δ(t − t )
m 0 0

2 −2
ζ
t where r(0) is particle position at time t = 0. In general,
+hv(0) i e m
A
Z t ζ
one is interested to monitor the time behaviour of how far
0 −2 (t−t0 )
=
m 2
dt e m – on average – the particle has moved from its original
0

2 −2
ζ
t
position. To this purpose, a central role is played by
+hv(0) i e m
  the so called time mean-square displacement δr2 (t) ≡
A ζ ζ
=
−2
1−e m
t 2 −2
+ hv(0) i e m .
t
(A5) h(r(t) − r(0))2 i. By employing Eq. A6, one gets:
2ζ m

t0 t000
*Z + 2 
t Z Z t Z  2
1 ζ
−m (t0 −t00 ) 00 00 0 ζ
−m (t000 −t0000 ) 0000 0000 000 m ζ
2
hδr (t)i = e f (t ) dt dt × e f (t ) dt dt + hv(0) i 2
1 − e− m t
m2 0 0 0 0 ζ
0 000 2 
Z t Z t Z t Z t  2
1 ζ 0 00
+t000 −t0000 ) m ζ
= dt0 dt00 dt000 dt0000 e− m (t −t hf (t00 )f (t0000 )i + hv(0)2 i 1 − e− m t
m2 0 0 0 0 ζ
Z t Z t0 Z t Z t000  2  2
2 ζ kB T 0 00 000 ζ
0000 − m (t0 −t00 +t000 −t0000 ) 00 0000 m ζ
= dt dt dt dt e δ(t − t ) + hv(0) i 2
1 − e− m t
m2 0 0 0 0 ζ
Z t Z t0 Z t  2  2
2 ζ kB T ζ 0 000
−2t00 ) m ζ
= dt0 dt00 dt000 e− m (t +t + hv(0)2 i 1 − e− m t . (A7)
m2 0 0 t00 ζ

By elementary integrations, one thus gets the final ex- pression:

m ζ

hδr2 (t)i = 2D t − 2D 1 − e− m t
ζ
 2 !
m m kB T  ζ
2
2
+ hv(0) i − 2
1 − e− m t .
ζ ζ
(A8)
 26

At large times, the particle follows the standard law The first case is quite general: simply, friction is negli-
of random Brownian motion with δr2 (t) = 2Dt, where gible and the particle behaves ballistically. The second
D = kBζ T is the diffusion coefficient. This clearly holds case holds instead only in the peculiar situation when
in general, in particular regardless of the detailed value of the initial velocity of the particle is exactly = 0. In this
hv(0)2 i. Conversely, the small time regime (t → 0) hides case, the Langevin equation seems to predict an unusual,
some surprise. In fact, by Taylor expansion of Eq. A8 for super-diffusive behaviour. Interestingly, an experimental
t → 0 we get: validation of this behaviour was reported recently [20],
! thus confirming the general validity of the Langevin ap-
 2
2 2 2 2 ζ 2 ζ
proach.
hδr (t)i ≈ hv(0) i t + D − hv(0) i t3
3 m m
2 2 2

 hv(0) i t , if hv(0) i =
 6 0
≈  2 (A9)
 2D ζ

t3
, if hv(0)2
i = 0
3 m

[1] In these notes, average values are indicated by brackets [11] D. Frenkel and B. Smit, Understanding molecular simu-
“h...i”. lation: from algorithms to applications (Academic Press,
[2] M. Rubinstein and R. H. Colby, Polymer Physics (Oxford San Diego, 2002).
University Press, New York, 2003). [12] Notice, that here we consider only classical systems.
[3] M. Doi and S. F. Edwards, The Theory of Polymer Dy- [13] M. Allen and D. Tildesley, Computer simulation of liq-
namics (Oxford University Press, New York, 1986). uids (Oxford University Press, Oxford, 1987).
[4] J. F. Marko and E. D. Siggia, Macromolecules 28, 8759 [14] K. Kremer and G. S. Grest, J. Chem. Phys. 92, 5057
(1995). (1990).
[5] P. J. Flory, Principles of Polymer Chemistry (Cornell [15] A. Rosa and R. Everaers, Plos Comput. Biol. 4, e1000153
University Press, Ithaca (NY), 1953). (2008).
[6] S. M. Bhattacharjee, A. Giacometti, and A. Maritan, J. [16] A. Rosa, N. B. Becker, and R. Everaers, Biophys. J. 98,
Phys.: Condens. Matter 25, 503101 (2013). 2410 (2010).
[7] J.-C. Walter and G. T. Barkema, Physica A 418, 78 [17] M. Di Stefano, A. Rosa, V. Belcastro, D. di Bernardo,
(2015). and C. Micheletti, Plos Comput. Biol. 9, e1003019
[8] A. D. Sokal, arXiv:hep-lat/9405016 (1994). (2013).
[9] P. E. Rouse, J. Chem. Phys. 21, 1272 (1953). [18] S. Plimpton, J. Comp. Phys. 117, 1 (1995).
[10] We remark that, in the continuum limit, the Rouse time [19] Full documentation, including how to down-
can be written in the transparent form τRouse = τ1 ≈ load and install LAMMPS, is provided at:
hR2 (N )i hR2 (N )i b2 [Link]
Dchain
= D mon/N
= Dmon N 2 , highlighting the fact
that, in order to relax, the polymer chain has to explore [20] J. Duplat, S. Kheifets, T. Li, M. G. Raizen, and E. Viller-
a region of its own size. maux, Phys. Rev. E 87, 020105 (2013).