Particle methods for Mesoscopic Simulations

1
Particle methods for Mesoscopic Simulations

• C.P. Lowe, University of Amsterdam

• Contents
• What do we try to achieve in a mesoscopic
  simulation?
• How do we simulate realistic dynamics?
• What techniques are at our disposal?
• Dissipative particle dynamics
• The Lowe-Andersen thermostat
• Stochastic rotational dynamics
• Stokesian dynamics
• Langevin dynamics
• What are their strength and weaknesses?

2
Mesoscopic Simulation

• There is no way to do it it still requires the
  scientists input
• To simulate on mesoscopic time and length scales
  you must throw
• out a load of stuff. But you must also keep the
  most important
• characteristics of the original system.
• What should you keep and why. How to you keep it?
• Are you solving a well defined set of equation to
  some degree
• of accuracy?
• How do you relate your simulation to the real
  world?
• (dimensional analysis)

3
Particle Models

• Disadvantages
• Computationally inefficient (?)
• Difficult to completely specify their properties
  (thermodynamic, transport coefficients) a priori.
• Difficult to introduce solid fluid interface
  accurately
• Advantages
• Intuitive, can have basic physical principles
  built in.
• Positions are continuous in space, good for
  interfaces
• Computationally efficient (?)
• Possible to have correct fluctuations (impotant
  on the mesoscopic scale)

4
Case study The dynamics of polymer molecules
Polymers are long molecules consisting of a
large number (up to many millions) of repeating
units. For example polyethylene
Longest time scale in the problem is the time it
takes the polymer to diffuse a distance of its
own size. How long is this? From the diffusion
equation root mean squared displacement D as a
function of time t is
so
where lp is a length of the polymer
5
A Tractable Simulation Model
I Modelling The Polymer Step 1 Simplify the
polymer to a bead-spring model
Simplest modelThe ideal chain The beads do not
interact with each other.
is the vector connecting bead i to bead i1
6
Properties of single polymer chains
Notice that terms involving will, on average
cancel, because they are equally likely to be
positive or negative. So, if the springs between
all monomers are the same we will have where
ltl2gt is the mean square separation between beads.

7
Properties of single polymer chains
Now we have an estimate for the size of the
polymer,
Experimentally b (polyethylene) 5.10-10m
b(DNA)
5.10-8m So for N106 lp(polyethylene)
5.10-7m (1/2m) lp(DNA)
5.10-5m (50m)
Use Stokes-Einstein to estimate D
kBoltzmanns constant TTemperature hshear
viscosity of solvent
kT(room temp.)4.10-14 gcm2/s2 h(water)0.01 g/cm
s
So
what is it for DNA?
8
A Tractable Simulation Model
We still need to simplify the problem
because simulating even this at the atomiclevel
needs ?t 10-9 s and we need to simulate for t
gt 1 s.
Step 2 Simplify the bead-spring model
further to a model with a few beads keeping the
essential (?) feature of the original long polymer
Rg0 , Dp0
Rg Rg0 Dp Dp0
9
Properties of single polymer chains

• The size of the polymer does not depend on the
  kind of springs we have between
• beads, only the mean separation

Can we say more about the distribution of the
end-to-end vector?
The Central limit theorem if we have a
quantity y, which is the sum of a set if random
variabls xi
Then
(average)
(width)
And what is more, for large N the distribution of
y is a Gaussian, regardless of the distribution
of xi
10
Properties of single polymer chains
Now consider the x component of the end to end
vector in terms of the connector vectors
We can apply the central limit theorem to this to
give
where nN-1, with N the number of monomer. So
(now using b2 to denote ltli2gt)
11
Properties of single polymer chains
The same applies for the other three components
so for the full end-to-end vector

• Notice that this does not depend on the form of
  the potential between monomers
• only on the mean squared separation b.
• If we take just a dimer (N2) then

Where A is the normalization constant and U the
potential between the beads. If we set
the distribution of the end-to-end vector of the
dimer is the same as the polymer
12
Effective potentials for polymer chains
The potential
between two beads in a dimer model of the polymer
is an effective potential. This is
an artificial potential that reproduces some
desired fearture of the original system (in this
case the end to end vector)
In general if there are monomers in the
original polymer and in the model, the
monomer/ monomer stretching potential
Reproduces the correct end to end
distribution and the correct bead/bead
distribution.
13
Effective potentials for polymer chains
In terms of the structure factor, S(k)
S(k)
Increasing
k
Note that S(k) is the response function for
density fluctuations (it determines how the
density in the system change with a an
external petrurbation). As such, the more beads
we have in the model the shorter wavelength
peturbations the model will resolve.
14
Summary

• To model the long time behaviour of polymers we
  need a very simple model.
• We necessarily throw away a lot of information
  about the true long polymer
• In this case we introduce effective potentials
  between beads that reproduces the true structure
  factor at long wavelengths.
• The more beads we put in the model, the shorter
  the wavelengths we resolve (we can put
  information back)
• This is almost a best case scenario for
  mesoscopic modelling,
• but what about the dynamics?

15
A model for the Dynamics
To simulate the dynamics we need to model the
effect of the solvent we can do this by 1)
Just including the thermal effects of the solvent
(i.e the fluctutations that jiggle the polymer
around) 2) Including the thermal and fluid-like
like behaviour of the solvent. This will include
the hyrodynamic interactions between the
monomers.
16
Approach 1
Solve a Langevin equation for the beads
Force on particle i
-gv is the friction force, here the friction
coefficient is related to the monomer diffusion
coeffiient by DgkT
is a random force with the property
is the sum of all other forces (in our case, the
bead force)
Many ways to solve this equation but often
pointless
Forbert HA, Chin SA. Phys Rev E 63, 016703 (2001)
Takes in into account particle inertia but not
solvent inertia which is inconsistent
17
Approach 1
Use an Andersen thermostat A method that
satisfies detailed balance (equilibium properties
correct)
Integrate the equations of motion with a normal
velocity Verlet algorithm
Then with a probability GDt (G is a bath
collision probability) set
Where qi is a Gaussian random number with zero
mean and unit variance. (i.e. take a new velocity
component from the correct Maxwellian)
Gives a velocity autocorrelation function
C(t)ltv(0)v(t)gt
Identical to the Langevin equation with g/mG
18
Approach 1-2
Use Brownian/Stokesian dynamics
Integrates over the inertial time in the Langevin
equation and solve the corresponding Smoluchowski
equation (a generalized diffusion equation). As
such, only particle positions enter.
are random displacements that satisfy
and is the mobility tensor
D.L. Ermak and J.A. McCammon, J. Chem. Phys. 969,
1352 (1978) Computer simulations of liquids,
M.P. Allen and D.J. Tildesley, (O.U. Press, 1987)
19
Approach 1-2
If the mobility tensor is approximated by
The algorithm is very simple. This corresponds to
neglecting hydrodynamic interactions
(HI) Including HI requires the pair terms. A
simple approximation based on the Oseen tensor
(the flow generated by a point force) is.
For a more accurate descrition it is much more
difficult but doable, see the work of Brady and
co-workers.
A. J. Banchio and J. F. Brady J. Chem. Phys.
118, 10323 (2003)
20
Hydrodynamics of polymer diffusion
b
a
b is the kuhn length a is the hydrodynamic
radius Note that this is not a real physical
radius. It is related to the diffusion
coefficient a bead would have were it not part of
the chain, Dmon It is still, however, a
parameter we must input into our model
21
Hydrodynamics of polymer diffusion
Exact result for any N (Kirkwood-Reissman)
hydrodynamic
bead
f(N) is a known function related to the structure
factor
For any long chain (n ? 8)
The second term becomes irrelevant, as does the
ratio a/b so
22
Hydrodynamics of polymer diffusion
Note that If we neglect hydrodynamic
interactions we have which is completely wrong
(so a method without hydrodynamics is useless)
The large N result really requires large N
so that f(N), a quantity we have no control over,
reaches its asymptotic scaling. Although the
ratio a/b is irrelevant for long chains it is not
for short model chains. Notably if a/bltlt1 there
will effectively be no hydrodynamic interaction
between beads for small N (see almost all DPD
polymer simulations). What value should we
choose for it?
23
Dynamic scaling
Choosing the Kuhn length b For a value a/b ¼
the lowest order finite N correction to f(N)
cancels the first (non-hydrodynamic term) and
long polymer scalingholds to a good approximation
even for small N.
Pure renormalization also gives 1/vN scaling
for small N to a good approximation.
Since the hydrodynamic radius is fixed by the
properties of the solvent this then fixes the
value we use for b
24
A Hydrodynamic model for the Dynamics
An explicit solvent model The solvent is
modelled explicitly as an ideal gas coupled to a
Lowe-Andersen thermostat More later but this
thermostat - Satisfies detailed balance -
Gallilean invariant (looks the same for a fluid
in uniform motion viewed in a co-moving frame of
reference as for a stationary fluid) -
Conservation of momentum - Isotropic (no
preferred direction) fluctuations
fluctuating hydrodynamics
Correct equilibrium properties
Ingredients for correct hydrodynamics
25
A Tractable Model for the Dynamics
II Modelling The Solvent Using an ideal gas
coupled to a Lowe-Andersen thermostat, in
parctice
(1) For all particles identify neighbours within
a distance rc (using cell and neighbour lists)
(2) Decide with some probability if a pair will
undergo a bath collision
(3) If yes, take a new relative velocity from a
Maxwellian, and give the particles the new
velocity such that momentum is conserved
(4) Advect particles
26
A Tractable Model for the Dynamics
Is it reasonable to use a gas to model the
solvent? Advantage the density distribution is
always uniform Disadavantage the dynamics of a
gas are not generally the same as the dynamics
of a liquid. Notably the Schmidt number Sc,
Here D is the diffusion coefficent and n the
kinematic viscosity
v
F
t
(nt)1/2
Gas Sc1
27
A Tractable Model for the Dynamics
Is it reasonable to use a gas to model the
solvent? Advantage the density distribution is
always uniform Disadavantage the dynamics of a
gas are not generally the same as the dynamics
of a liquid. Notably the Schmidt number Sc,
Here D is the diffusion coefficent and n the
kinematic viscosity
t
(nt)1/2
Liquid Scgtgt1
28
A Tractable Model for the Dynamics
But the thermostat transfers momentum so
contributes an additional viscosity. Roughly

• Note that this is independent of temperature
  whereas DkT
• So if we fix G and reduce kT, D goes down while n
  stays roughly constant
• This way we can make the Schmidt number of our
  gas as high as we want
• More technically, we have a parameter that we
  must set

and
so we must set Lltlt1
29
A Tractable Model for the Dynamics
What about the density? If the density is low
There will be few thermostat collisions and a
very small value of L will be necessary to
ensure Scgtgt1
If the density is high

• easier to get Scgtgt1 but
• rc becomes large and we only expect hydrodynamic
  behaviour on lengths great than rc
• You have to evaluate a lot of interactions, which
  is computationally demanding
• Compromise

30
A Tractable Model for the Dynamics
Locating the interactions
Construct a neighbour list for each particle of
all particles within a distance rcd. To do so
use a cell list for cubic cells of side rcd.
Update cell list when maximum relative
displacement exceeds d.
Construct a Cell list for cubic cells of side rc.
Search the 27 cells within which all interactions
must lie.
A neighbour list reduces the number of
interactions evaluated by a factor 6 if, with
small d, it only needs updating infrequently.
Otherwise its an overhead Lltlt1 pays to use a
neighbour list Lgt1 does not pay to use a
neighbour list. Check!
31
A Tractable Model for the Dynamics

• Isnt and ideal gas too compressible?
• Compressibility is a matter of length scales.
• Dynamically it is the time it takes sound to
  propagate
• a distance l, ts, relative to other time-scales
  in the problem.
• Notably here, the viscous time tH
• Speed of sound in an ideal gas Cs

So the sonic time is
For mesoscopic l, in reality tsltlttH, but note
that decreasing L to increase Sc also increases
ts/tH. That is, the model solvent is becoming too
compressible Again a trade off is required
32
A Tractable Model for the Dynamics
Bead-Solvent interactions Thermostat
interactions between the beads and the solvent
are the same as the solvent-solvent
interactions. There are no bead-bead
interactions. We are solving a well defined
model in which monomers only act only via the
fluid The model fluid provides the fluctuations
and the time-dependent hydrodynamic interactions
between monomers. The bead hydrodynamic radius
is strictly related to the solvent diffusion
coefficient (solvent and monomers being
equivalent)
33
Time Scales
So how do we do in the end?
time it takes momentum to diffuse l time it
takes sound to travel l time it takes a polymer
to diffuse l
Reality ts lt tH ltlt tD Model (N 2) ts tH lt
tD Gets better with increasing N
34
Relative to the Alternatives
I DPD Very similar but harder to integrate
the equations of motion. Same parametric
considerations and callibration required II
Stochastic Rotational Dynamics Not shown to
work in the correct parameter regime. Same
parametric considerations and callibration
required
35
Relative to the Alternatives
III Lattice-Boltzmann Better control of the
parameters. No rigorous thermodynamics
(fluctuation dissipation). Bead must be mapped in
a continuous way onto the lattice IV
Stokesian Dynamics No callibration required.
Only parametric consideration is the ratio a/b.
Does not need an explicit solvent. Much faster
for this problem. But neglects the
time-dependence of the HI (fluid inertia). Poor
scaling with number of beads/polymers. External
geometries are difficult.
36
Does it work?
N 32 (?)
N 16 (?)
One simulation is a 32 bead chain with only half
the beads shown The other is a sixteen bead
simulation In both cases the time display time is
in units of the polymer diffusion time
37
Does it work?
b 4a requires b solvent particle separation
so maybe not, for small N
Hydrodynamic contribution to the diffusion
coefficient for model chains with varying bead
number N
38
Centre of mass motion
Convergence excellent. Not exponential decay
(time dependence effect).
39
Surprise, its algebraic
40
Solves a more relevant problem viscosity
Time dependent polymer contribution to the
viscosity For polyethylene tp 0.1 s
41
Summary

• This problem is doable.
• We can calculate the true long-time behaviour of
  model
• ideal polymers from simulations of short model
  polymers
• Despite this problems relative simplicity it
  requires a lot of work
• We used a lot of scientific input along the way
• The choice of every parameter has to be carefully
  justified
• For this problem, an explicit solvent may not be
  the
• best method
• Dimensional analysis is vital in relating the
  results to real world
• quantities.

42
Interacting chains
Really we want to simulate interacting chains. We
start with one excluded volume
chain. Question How do we renormalize the
static properties (interactions between blobs
of polymer)?
43
Interacting chains
Problem I Flory excluded volume parameter ?
(effective monomer volume).
? lattice volume x (1/2 ? )
What is ? off-lattice (i.e. in reality)? This is
solved.
44
Interacting chains

• Problem II
• What do we need to reproduce with an effective
  monomer/monomer potential?- Ideal chain size
  (easy)
• Same degree of expansion independent of N (hard)

45
Interacting chains
Is this problem already solved ? R.M. Jendrejack
et al., J. Chem Phys 116, 7752 (2002)
46
Interacting chains
Plot one way
Plot another way
47
Interacting chains
Alternative Florys result So keep
constant. (B2 second virial
coefficient)
48
Interacting chains
This works for depressing large N.
49
Stochastic Rotational dynamics
A. Malevanets and R. Kapral, J. Chem. Phys 110,
8605 (1999).
Advect
collide
random grid shift recovers Gallilean invariance
50
Stochastic Rotational dynamics
Collide particles in same cell
basically rotates the relative velocity vector
where the box centre of mass velocity is
with Ncell the number of particles in a given
cell. R is the matrix for a rotation about a
random axis

• Advantages
• Trendy
• Computationally simple
• Conserves mometum
• Conserves energy
• Disadvantages
• Does not conserve angular momentum
• Introduces boxes
• Isotropy?
• Gallilean invariance jammed in by grid shift
• Conserves energy (need a thermostat for
  non-equilibrium simulations)

51
Stochastic Rotational dynamics

• Equation of state Ideal gas
• Parametrically exactly the same as all other
  ideal gas models
• must fix
• number of particles per cell (cf r)
• degree of rotation per collision (cf G)
• number of cells traversed before velocity is
  decorrelated (cf L)
• Transport coefficients the same as the
  dissipative ideal gas,
• theoretical results accurate in the range of
  parameters where
• Sc is small. For realistic parameters, must
  callibrate.
• For an analysis see
• J.T. Padding abd A.A. Louis, Phys. Rev. Lett. 93,
  2201601 (2004)

52
Dissipative Particle Dynamics
First introduced by Koelman and Hoogerbrugge as
an off-lattice lattice gas method with discrete
propagation and collision step. P.J.
Hoogerbrugge and J.M.V.A. Koelman, Europhys.
Lett. 19, 155 (1992) J.M.V.A. Koelman and P.J.
Hoogerbrugge and , Europhys. Lett. 21, 363
(1993) This formulation had no well defined
equilibrium state (i.e. corresponded to no known
statistical ensemble). This didnt stop them and
others using it though. The formulation usually
used now is due to Espanol and Warren. P.
Espanol and P.B. Warren, Europhys. Lett. 30, 191,
(1995). Particles move according to Newtons
equations of motion


53
Dissipative Particle Dynamics
So what are the forces? They are three fold and
are each pairwise additive
The conservative force
where
Is a repulsion parameter Is an interaction
cut-off range parameter
54
Dissipative Particle Dynamics
What is the Conservative force? Simple a
repulsive potential with the form
U(r)
aijrc
It is soft in that, compared to molecular
dynamics it does not diverge to infinity at any
point (there is no hard core repulsion.
rc
The dissipative force
Component of relative velocity along line of
centres
55
Dissipative Particle Dynamics

• What is the Dissipative force?
• A friction force that dissipates relative
  momentum
• (hence kinetic energy)
• A friction force that transports momentum
• between particles

?
wd
rc
The random force
56
Fluctuation Dissipation
1
To have the correct canonical distribution
function (constant NVT) the dissipative (cools
the system) and random (heats the system) forces
are related
wd
rc
For historical (convenient?) reasons wd is given
the same form as the conservative force
The weight functions are related
As are the amplitudes
57
DPD as Soft Particles and a Thermostat
Without the random and dissipative force, this
would simply be molecular dynamics with a soft
repulsive potential. With the dissipative and
random forces the system has a canonical
distribution, so they act as a thermostat. These
two parts of the method are quite separate but
the thermostat has a number of nice
features. Local Conserves Momentum Gallilean
Invariant
58
Integrating the equations of motion

• How to solve the DPD equations of motion is
  itself something of an issue.
• The nice property of molecular dynamics type
  algorithms (e.g. satisfying
• detailed balance) are lost because of the
  velocity dependent dissipative force.
• This is particularly true in the parametrically
  correct (Llt1) regime
• Why is this important?
• Any of these algorithms are okay if the time-step
  is small enough
• The longer a time-step you can use, the less
  computational time your
• simulations need
• How long a time step can I use?
• Beware to check more than that the temperature is
  correct
• The radial distribution function is a more
  sensitive test. The temperature
• can be okay while other equilibrium properties
  are severely inaccurate.
• L-J.Chen, Z-Y Lu, H-J ian, Z-Li, and C-C Sun, J.
  Chem. Phys. 122, 104907 (2005)

59
Integrating the equations of motion
Euler-type algorithm
P. Espanol and P.B. Warren, Europhys. Lett. 30,
191, (1995).
And note that, because we are solving a
stochastic differential equation
(Applies for all the following except the LA
thermostat)
60
Integrating the equations of motion
Modified velocity Verlet algorithm
R.D. Groot and P.B. Warren, J. Chem. Phys. 107,
4423, (1997).
Here l is an adjustable parameter in the range
0-1

• Still widely used
• Actually equivalent to the Euler-like scheme

61
Integrating the equations of motion
Self-consistent algorithm I Pagonabarraga,
M.H.J. Hagen and D. Frenkel, Europhys. Lett.
42, 377, (1998).

• Updating of velocities is performed iteratively
• Satisfies detailed balance (longer time-steps
  possible)
• Computationally more demanding

62
Changing the equations of motion?
Lowe-Andersen thermostat (LAT) C.P.Lowe,
Europhys. Lett. 47, 145, (1999).
Bath collision

• Here G is a bath collision frequency (plays a
  similar role to g/m in DPD)
• Bath collisions are processed for all pairs with
  rijltrc
• The current value of the velocity is always used
  in the bath collision (hence
• the lack of an explicit time on the R.H.S.)
• The quantity x is a random number uniformly
  distributed in the range 0-1
• The quantity mij is the reduced mass for
  particles i and j, mijmi mj/(mimj)

63
DPD versus LAT

• The LAT is an alternative to the thermostat part
  of the DPD algorithm
• (the combined random and dissipative forces)
• Similarities as a thermostat
• Conserves linear and angular momentum
• Gallilean Invariant
• Local
• Advantage satisfies detailed balance by
  construction, allowing longer time-steps without
  the use of an iterative scheme for solving the
  equations of motion.
• Disadvantage? It does not use weight functions
  wd and wr (or alternatively
• you could say it uses a hat shaped weight
  functions)
• But, no-one has ever shown these are useful or
  what form they should best take. The form
  wr(1-rij/rc) is only used for convenience (work
  for someone?)
• They could be introduced using a distance
  dependent collision probability
• In the limit of small time-steps LAT and DPD are
  actually equivalent!
• E.A.J.F. Peters, Europhys. Lett. 66, 311 (2004).
• Word of warning in the LAT, bath collisions must
  be processed in a random order

64
Which method should I use?
1) It depends on the conservative force
(interaction potential). The time step must
always be small enough such that the conservative
equations of motion adequately conserve total
energy. To check this, run the simulation
without the thermostat and check total
energy. 2) If this limits the time-step the
methods that satisfy detailed balance lose their
advantage. 3) If not, use the self-consistent
or LAT methods. Never Euler or modified Verlet. 4
) There are some much better methods that still
do not strictly satisfy detailed balance (based
on more sophisticated Langevin-type algorithms). W
.K. den Otter and J.H.R. Clarke, Europhys. Lett.
53, 426 (2001). T. Shardlowe, SIAM J. Sci.
Comput. (USA) 24, 1267 (2003). 5) For a review
see P. Nikunen, M. Karttunen and I. Vattulainen,
Comp. Phys. Comm. 153, 407 (2003).
65
DPD Summary

• The dissipative and random forces combine to act
  as a thermostat
• (Fullfilling the same function as Nose-Hoover or
  Andersen thermostats
• in MD)
• As a thermostat it has a number of advantages
  over some commonly
• used MD thermostats
• The conservative force corresponds to a simple
  soft repulsive harmonic
• potential between particles, but in principle it
  could be anything
• (The DPD thermostat can also be used in MD
• T. Soddemann, B. Dunweg and K. Kremer, Phys. Rev.
  E68, 046702 (2003) )
• The equations of motion are awkward to integrate
  accurately with
• large time-steps. Chose your algorithm and test
  it with care.
• The Lowe-Andersen thermostat has the same
  features as the DPD
• thermostat but is computationally more efficient
  as it allows longer
• time-steps.

66
Why this form for the conservative force?
In principle the conservative force can be
anything you like, what are the reasons for this
choice? Some common statements
It is the effective interaction between blobs of
fluid No it isnt, at least not unless you are
very careful about what you mean by effective.
A soft potential that allows longer
time-steps Maybe, but relative to what?
Factually it is not a Lennard-Jones (or
molecular-like) potential. It is the simplest
soft potential with a force that vanishes at some
distance rc
As with any soft potential it has a simple
equation of state in the fluid regime and at high
densities.
67
The equation of state of a DPD Fluid
For a single component fluid with pairwise
additive spherically symmetric interparticle
potentails the pressure P in terms of the
radial distribution function g(r) is
where r is the density. For a soft potential
with range rc at high densities, rgtgt3/(4prc3)
g(r)1 so
g(r) real fluid
g(r) DPD fluid
Where a is a constant. For DPD a.101aijrc4

• Note though that
• If r is too high or kT too low the DPD fluid will
  freeze
• making the method useless.
• And a/kT is not the the true second Virial
  coefficient
• so this does not hold at high densities

EoS
68
Mapping a DPD Fluid to a real fluid
R.D. Groot and P.B. Warren, J. Chem. Phys. 107,
4423, (1997). Match the dimensionless
compressibility k for a DPD fluid to that of real
fluid
For a (high density) DPD fluid, from the
equation of state
For water k-116 so in DPD aij75kT/rrc4 Once the
density is fixed, this fixes the repulsion
parameter. You can use a similar procedure to
map the dimensionless compressibility of other
fluids.
69
Whats right and whats wrong
By setting the dimensionless compressibility
correctly we will get the correct thermodynamic
driving forces FThfor small pressure
gradients (the chemical potential gradient is
also correct)
Technically, we reproduce the structure factor at
long wavelengths correctly. But, other things
are completely wrong, eg the compressibility
factor P/rkT And this assumes on DPD particle is
one water molecule. If it represents n water
molecules the r(real)r(model) so aij must be
naij(n1). That is the repulsion parameter is
scaled with n and if ngt1 the fluid freezes. R.D.
Groot and K.L. Rabone, Biophys. J. 81, 725 (2001).
70
DPD for systems with a given equation of state
I. Paganabarraga and D. Frenkel, J. Chem. Phys
155, 5015 (2001)
The basic idea is to input an equation of
state. To do so a local density is defined
where r is a weight function that vanishes for
rijgtrc The conservative force is the the
derivative of the free-energy (as a function of
r) w.r.t. the particle positions
Where is the excess free energy per partilce
as calculated from the EoS Is the density really
the density? Is it a free energy or a potential
energy?
71
DPD for systems with a given equation of state
Eg a van der Waals fluid
where A and B are parameters (related to the
critical properties of the fluid). Simplest EoS
that gives a gas liquid transition.
72
DPD for immiscible fluids
If the DPD partices are not the all the same
species (i.e. the repulsion parameter/and or
interaction range between different type
particles differs) the equation of state will
still be given to a good approximation by
the mean field approximation at high
densities. Mapped on the the c-parameters of
Flory Huggins theory by Groot and Warren (not
dwelled on here because this equation of state is
poor) Here we look at an immiscible DPD mixture
used to model a droplet. Parameters give
73
DPD for immiscible fluids
Corresponds to oil/1-propanol mixture at room
temperature. Droplet size 173 mm
vs Experiment
74
DPD for immiscible fluids
Coalescence
Minimum distance between droplet and substrate
75
Case Study The dynamics of biofilaments
M. C. Lagomarsino, I. Pagoabarraga and C.P.
Lowe Phys. Rev. Lett. 94, 148104 (2004). M. C.
Lagomarsino, F.Capuani and C.P. Lowe, J. Theor.
Biol. 224, 205 (2003)
Example, tying a knot in Actin
76
Effective potentials?
This was derived for continuum solids and has no
right to work for the biofilaments we are
considering now. It turns out that an energy of
this form is a good approximation for nano-scale
filaments, but not the relation of the bending
modulus to Youngs modulus. It is a
better approximation to assume these filaments
cannot be stretched (Y is infinite) The reason
is, resistance to bending has more to do with
resistance to bending of bonds between the
elements that compression/expansion.
77
Modelling the dynamics of biofilaments

• Motivation
• 1) To test how well this assumption works
• 2) For practical use in non-biological systems a
  way must be found
• to at least point the filaments in a
  particular direction (although they
• have a dipole moment it is too small for them
  to align in an electric
• field)
• 3) We have to do it numerically because the
  equations are too hard

78
Numerical Model
Fb - bending force (from the bending energy for
a filament with stiffness k that we described
earlier) Ft - Tension force (satisfies
constraint of no relative displacement along the
line of the links) Ff - Fluid force (from the
model discussed earlier, with F the sum of all
non hydrodynamic forces) Fx - External
force Solve equations of motion using a Langevin
Equation!!
79
Dynamics we must account for the very sticky
fluid around the filament
At its simplest, resistive force theory
80

• Taking the simple model for the energy and
  applying it to a sperms
• tail (a circular arrangement of nine mircotubules
  with dyneins supplying
• the motions by sliding between filaments)
• Assuming that the filaments do not stretch
• Using this model for the dynamics

Good predictions for the swimming speed of simple
spermatozoa
81
Why might this not give a complete picture and
what about orientation?
A simple model, a chain of rigidly
connected point particles with a friction
coefficient g
82
Why might this not give a complete picture?
A simple model, a chain of rigidly connected
point particles with a friction coefficient g
subject to an external force F
Ff -g (v-vf)
Vf
v
83
The Oseen tensor gives the solution to the fluid
flow equations (on a small scale) for a point
force acting on a fluid. This gives the velocity
of the fluid due to the force on another bead as
These equations are linear so solutions just add
Stokesian dynamics without the the fluctuations
84
Approximate the solution as an integral. For a
uniform perpendicular force.

• s the distance along a rod of unit length
• b is the bead separation

85
Approximate the solution as an integral. For a
uniform perpendicular force.

• s the distance along a rod of unit length
• b is the bead separation

If the velocity is uniform the friction is higher
at the end than in the middle
86
What happens with uniform force acting downwards?
Sed B FL2/k ratio of bending to hydrodynamic
forces If the filament is long enough, the
bending modulus small enough or the force high
enough, the filament bends significantly.
87
B 300
88
B 3000
89
B 15000
90
B 1, filament aligned at 450
91
Why?
Aligned more parallel, lower friction force
F
A component of the force perpendicular to the
force bends it and moves it left.
Aligned more perpendicular, higher friction force
So a torque acts on the fibre to rotate it
towards the prependicular
92
How long does it take to reorientate?

• From this we can
• work out what conditions are necessary
• in the real world to see the effect
• work out when the approximation of
• neglecting diffusion and dipole orientation
• is sensible.

93
Is this practically relevant?

• For a microtubule the bending modulus is known
• and we estimate, B 1 requires F1 pN for a 10
  micron
• microtubule. This is reasonable on the micrometer
  scale.

• For sedimentation (external force is gravity) ,
  no.
• Gravity is not strong enough. Youd need a
  ultracentrifuge

• Microtubules are barely charged and the charge is
  known, we estimate
• an electric field of 100 V/m for B 200 (L30
  microns). So it
• should be doable.

94
But this only orientate the filament in a plane
(perpendicular to the force direction) What if
we apply a force in a direction that rotates?
Circularly polarised electric field
Electric field as a function of time
95
Dimensional Analysis 30m Microtubule in
water, Field 100 V/m Frequecy 1 Hz Movie
timereal time
96
Dimensional Analysis 30m Microtubule in
water, Field 100 V/m Frequecy 1 Hz Movie
timereal time