Kurt Kremer
(4) Addition and subtraction of these two equations yields
3 Polymer simulations
3.2 How to generate conformations
Having demonstrated that the methods which follow the realistic dynamics are not very useful for single isolated chains, we now try to generate the conformations by a purely stochastic process.
3.2.1 Simple sampling techniques
First let us consider SAWs on a lattice. To fulfil the excluded volume requirement each lattice site can only be occupied once; but otherwise each conformation of an N-step walk has the same probability. If we fix the first step, then each new step is taken with probability l / ( q - l), where q is the coordination number of the lattice, and we account for the fact that backward steps are ruled out.
This most naive way of generating conformations with equal probability is called simple sampling. Each time an attempted new bond hits a site which is already occupied, one has to start again at the beginning. Otherwise, different SAW conformations will receive different probabilities, generating an (attractive) effective interaction among the monomers. Each sampled conformation is therefore taken randomly out of the q ( q - l)N-l possible random walk paths, which do not include direct back-folding, whereas the total number of SAWs on a lattice is given by [3]
Z ( N ) = ciJq$N7-1 N >> 1 (14) with qeff
<
q - 1 and CO a number of order unity. The critical exponent y is dimension- but not lattice-dependent (y zz 7/6 for d = 3, and 4/3 for d = 2). Typical numbers for qeff are 2.6385 (square, p = 4), 2.879 (diamond, q = 4) and 4.6835 (simple cubic, q = 6) [3].Thus the success rate of the sampling process, A ( N ) = Ao(qeE/ ( q - l ) ) N NT-', decays exponentially with N . A typical value for A is e.g. A(100) = 0.03 on the diamond lattice.
This illustrates that simple sampling is only useful to get a rough estimate for chain properties like qeff and y.
The first improvement was suggested as early as 1955 by Rosenbluth and Rosen- bluth [26]. Their idea of biased sampling is to look ahead for at least one step in order to overcome the attrition. More modern approaches look several steps ahead or implement this idea within in a dynamical scheme [8, lo]. Here I explain two alternatives, which can mainly be used for isolated chains.
3.2.2 Dimerisation
The dimerisation approach (Figure 3) [27, 281 allows a simultanous determination of both the exponents 7 and U. The idea is to begin by generating many short SAWs of
Computer simulations in soft matter science 155
I I I
t i
J rl
Figure 3. Illustration of the dimerisation method.
length No. These walks can then be combined randomly. For every step the success rate
&in for a binary assembly of walks of length NI, and NZ = N - N I is Pbin(N1, N2) = Z ( N ) / ( Z ( N 1 ) Z ( N 2 ) ) . Using Z c(qNiVy-' we get Pbin(IV1, N z ) = N7-'/(N1N2)7-'; ob- viously &in is optimal for N 1 = Nz = N / 2 . The original procedure, where short walks were generated and then stored, is computationally not very efficient. It is more direct to generate each sub-walk by simple sampling, as follows [3, 291. After a S24W of N1 monomers is generated we continue. If the next N I steps violate the excluded volume condition within the second piece, we start at N1 again. Only an overlap with the first iV1 monomers makes it necessary to start at the very beginning. This defines a very efficient hierarchical procedure.
Batoulis and Kremer generalised this approach to study good-solvent properties of star polymers [29]. In particular, the dependence of y on the number
f
of arms of the star is interesting and very difficult to estimate analytically. Using this method we were able to give precise results for~ ( f )
and also size statistics such as the mean radius of gyration, and the mean hydrodynamic radius (RH) = N - z ~ i & j ( l / ~ z j ) . For R H , very strong corrections to scaling are observed, which explain deviations of experiments such as light scattering from the asymptotic power law behaviour. Figure 4 shows the extrapolation of the hydrodynamic radii as a function of N - ' J 2 , which was thought to be the correction to the scaling for this quantity. The correct value is 1 - v [30].3.2.3 Pivot algorithm
For this algorithm [31, 321, as illustrated in Figure 5 , a point on the chain is chosen randomly and one part of the chain is rotated at random. As with dimerisation, this can easily be done for both lattice and off-lattice systems. Given a chain of length N and a
156 Kurt Kremer
0.27
1
0.25
0 0.02 0.04 0.06 0.08 0.10 0.12
Figure 4. The extrapolated hydrodynamic radii as a function of N1i2 129,’.
t-‘
Figure 5. Illustration of the pivot algorithm.
pivot point at position Nl, then the acceptance rate p is given by the probability that the new conformation has no overlaps. On a lattice this needs a simple U ( N ) check, while in continuous space this is more difficult. On lattices, especially for d = 2, one has to take care that the choice of moves does not contain hidden conservation rules. The approach is then ergodic as shown by Madras and Sokal [31] who claim an No.2 power law for the
‘relaxation time’ of the mean-square end-to-end distance. This method was used to obtain very precise estimates of the exponent v, v = 0.7496 & 0.0007 ( d = 2: exact value 3/4) and v = 0.592 & 0.003 ( d = 3), which is about as accurate as the results coming from dimerisation. For the single isolated SAW there is probably no better way to generate very quickly many conformations that are globally different.
However there is no ‘free lunch’. As soon as the concentration is increased, either due to other chains or monomer-monomer attraction along the same chain, the acceptance rate is dramatically reduced. A more subtle aspect is that this method relaxes large length scales very fast, while short ones need a longer time. For example, nearest-neighbour bond correlations need up to O ( N 2 ) moves to relax completely. For many problems it is essential to cover the equilibration of the short as well as the long distances. This can
Computer simulations in soft matter science 157
be overcome in two ways (depending on the model). One is a hybrid method of pivot moves and molecular dynamics simulations (for bead spring models) [33]. For fixed bond lengths it is sufficient instead to generalise the pivot approach so that not only tails of the chain, but also internal pieces of arbitrary length, are rotated. The choice of lengths is then adjusted to the density, so that the acceptance rate stays above 10%. This version can be applied to ring polymers as well.
3.2.4 Many-chain methods: generalised reptation algorithms
In semidilute solutions or melts the generalised pivot algorithm will only produce accept- able success rates for moves involving a small number of monomers. Eventually it joins the class of dynamic algorithms, which follow realistically the very slow physical dynamics of the system. This is certainly not what we are looking for in the present context.
The generalised reptation algorithm is an effective way to simulate dense solutions, and also for single chains, in the collapsed regime for T
<
0, the Theta temperature (see Khokhlov, this volume). The method is explained here for a single lattice chain. The original idea of the ‘slithering snake’ goes back to Kron in 1965 and Wall and Mandel 1975 [34]. One randomly takes an end monomer and tries to add it at the other chain- end with a random orientation. If the new chain fulfils the SAW condition, the move is accepted, otherwise rejected. In this way the chain moves (forward and backward) like a snake along its own contour. (The algorithm resembles, but should not be confused with, the physically realistic reptation dynamics of polymer melts: see McLeish, this volume, and below.) For this method detailed balance is fulfilled trivially, but on a lattice, only a subspace of all conformations is reached. (This subspace excludes conformations where both chain ends are completely surrounded by other monomers.) In continuous space this problem is absent. The extension to many chains as well as interacting chains is trivial. One can also introduce a grand canonical version by allowing the chain length to fluctuate. The use of pointers allows for very fast and efficient codes.How fast is the relaxation of a chain conformation in CPU time? To estimate this, let us look at the position of the middle monomer of the chain. At first this does not move. A measure of the relaxation time Tgpu is the time until this space point is not on the chain for the first time. For random walks this first passage time follows the same power law as a one-dimensional diffusion along the chains. Since each move only requires a constant N-independent number of operations we get Tgpu c(N Z for ideal chains. (To simulate this is a good check on the program.) This method is O ( N ) faster than algorithms with
‘realistic dynamics’ such as the Rouse model (see below), whose real relaxation time is TN N N 2 giving N 3 in CPU time. tage for simulations, is not completely understood. The explanation is probably in the difference between diffusion times and first passage times [35]. In addition, it should be kept in mind that because of the SAW condition, SAWs are correlated objects!
158 Kurt Kremer