Simulation of polymer systems

Bull. Mater. Sci., Vol. 7, No. 1, March 1985,pp. 35-41. © Printed in india. Simulation o f p o l y m e r s y s t e m s S L NARASIMHAN, P S GOYAL and B A DASANNACHARYA Nuclear PhysicsDivision.BhabhaAtomicResearchCentre,Trombay,Bombay400085,India Abstract. The present paper describes an algorithm which can generate, even on a small computer,arbitrarily long polymerchains, making sure that the configurationsgenerateddo not suffer from boundary effects.This has been achieved by employing the concept of a window,whichis an analogueof virtual memoryscheme.The algorithmhas been tested for the case of dilute polymer solution. Keywords. Polymerchains; self avoiding random walks; Monte-Carlo methods. I, Introduction The configurational properties of polymer chains in polymer solutions or in polymer melts can be studied experimentally using light scattering or neutron scattering techniques. While the light scattering technique is employed for studying dilute polymer solutions, small angle neutron scattering (SANS) technique is employed for studying dilute as well as concentrated polymer solutions, and also polymer melts (Higgins and Stein 1978). The configurational properties of linear polymer chains have also been studied by computer-simulating self avoiding random walks (sAw) on various lattices (Mckenzie 1976) using either the exact enumeration method or the Monte Carlo method. Configurational properties of short saws consisting of about 24 steps or less have been studied by enumerating all possible configurations of the saw of a given length. But it is not possible to list all the configurations of a long saw which consists of, say, hundreds of steps. Therefore one has to be content with generating as large a set of configurations as one can, using Monte Carlo simulation technique. One of the attempts at generating long saws using Monte Carlo technique is due to Wall (1954). His method consists in generating saws on a finite lattice whose size is determined by the available main memory of the computer used. One usually starts generating a walk from the centre of the lattice. Steps are generated by choosing, at random, one of all the nearest neighbours of the currently occupied site and remembering that backward step is not allowed. Should the site chosen be already occupied, the attempt is discarded and a fresh attempt to generate the chain is made from the centre of the lattice. The effect of this phenomenon, called "attrition", becomes more and more severe as one tries to generate longer and longer walks. It has been reported that out of 140000 attempts made at generating a 121-step walk, only one was successful (Wall 1954), But the successful attempts lead to configurations which are equally probable. T o reduce the effect o f attrition, Wall has proposed an enrichment scheme (Wall et a11963) with the help of which walks consisting of about 800 steps have been generated. Rosenbluth and Rosenbluth (1955) independently proposed a method which reduces 35 36 S L Narasimhan, P S Goyal and B A Dasannacharya the effect of attrition to a large extent but at the cost of biasing the walks. In this method steps are generated by choosing, at random, one of the unoccupied nearest neighbours--not one of all the nearest neighbours as in Wall's m e t h o d - - o f the currently occupied site. This way, generation o f the walk proceeds until a site which has no unoccupied nearest neighbours is encountered. In such a situation, generation of the walk is terminated, and a fresh attempt is made. As this method generates walk with some bias the configurations generated are weighted so as to make them equally probable. In either of these methods, attrition is not zero. Moreover, for a given size of the lattice, walks which consist of more than a certain number of steps will suffer from the effects of the boundary (Wall et al 1976). To give an example, suppose we have a 16-bit microcomputer which provides a program space of about 20 K words to the user, and we want to generate SAWS on a cubic lattice. Assume that, in this program environment, we are constrained to generate saw s within a cubic array, LAX, of dimensions 15 x 15 x 15. Clearly, even in principle, we cannot generate as long a SAWas we want to. Also, walks which consist of more than, say 30 or 40 steps, will have finite probability o f hitting the boundary of LAX,and therefore will have to be coiled back into the interior of LAT. In the present paper we report an algorithm which treats LAr as merely a 3-dimensional window through which we can look into the 3d-real space. Starting from the centre of the window, it generates a segment of the walk, stores the coordinates of the sites visited, and uses the same window after proper initialization, to generate another segment of the ,same length. This process is repeated till a walk of desired length has been generated. If and when it encounters a walk which enters a blocked state and cannot proceed further, it enforces a back-tracking algorithm so that the walk gets out of its blocked state and proceeds further. This way it can generate as long a walk as is desired even in a restricted program environment, and every attempt made at generating such a walk is successful. The walks generated are weighted using a procedure due to Rosenbluth and Rosenbluth (1955). The concept of a window is discussed in § 2. Section 3 describes the algorithm and § 4 gives the results we have so far obtained. A discussion of results and the capabilities of the algorithm are given in § 5. 2. Concept of a window Let LATbe a finite cubic array, having dimensions L x L x L, so that any location can be accessed by specifying the array indices I, J, K, where I, J, K = l . . . . . L. Imagine that EATis kept in 3d-real space. Let Xo, Yo, Zo be the space coordinates in 3d-real space of the centre of EAT.If b is the given step length, then the site (I, J, K) of EATwill map in to a point with coordinates (X, Y,Z) in 3d-real space according to the transformation: X -- [ I - (L + 1)/2] b + Xo Y = [ J - (L + 1)/2]b + Yo (1) Z = [K - (L + 1)/2]b + Z o LA'rcan be moved in space in such a way that its centre coincides with the point (X, Y,Z). Then any site (1', J', K') of LAT will map into a point with coordinates (X', Y', Z') according to the transformation. Simulatwn of polymer systems }'~ 2 1 5 37 5 4 2 3 ZoL. (x, .... i , / 4 II/ ~ X 5 I / / [ / ', ,~" ...../" 11 I ., ..-- 'I /+ i " I I t~ .~' I I I I I I I I I ,,, ,,, ,I r , k___L___L___L___J it" Figure i. ', A ' I , . . . . 4I - - - - tj - ~ - - +I - - - a t t I::> X Illustration of a two dimensional 5 x 5, moving window. X' = [I' - (L + l)/2]b + X Y' = [ J ' - (L + 1)/2]b + Y (2) Z' = [K' - (L + l)/2]b + Z We will call (1) the "window transformation" and LATthe "window". Figure 1 illustrates a 2d-window, for L = 5. 3. Algorithm Suppose we want to generate a 4095-step SAW and are constrained to work with a minimal data base consisting of the following three arrays: (a) EAT (]5, 15, 15)--the window; (b) COORD(5 ] 2)--an array in which to store the space coordinates of the sites generated and (c) NRAND (512)--array in which to store the number of unoccupied neighbours each site has when the walk proceeds from there. Assume that our machine configuration includes a floppy subsystem. (i) Choose some point (Xo, Yo, Zo) in space to start the walk from; start with a clean window, i.e., initialize array EAT with zeros. Assume that in the window site with coordinates (I -- 8, J = 8, K -- 8), the centre of the window, maps on to (Xo, Yo, Zo)Let (I,J,K) be the window indices of the current site in LaW. Occupy tile site (Xo, Y0,Zo) by entering data 1 in location (I,J,K) of LAT. Store the coordinates (X0, Yo, Zo) in array COORD.Let p be the number of unoccupied neighbours of(I, J, K). Store p in array NRAND. Choose one of the unoccupied neighbours of (I,J,K) at random and occupy it. Using the window-transformation (1) calculate the space coordinates of the point into which the chosen site maps into. Store these coordinates in cOOeD and so on. Repeat this process until a walk segment consisting of 6 steps is generated. We restrict ourselves to generating a walk segment consisting of 6 steps in one window because a segment consisting of more than 6 steps has a non-zero S L Narasimhan, P S Goyal and B A Dasannacharya 38 r" . .......... . . . .i--. . . . . . . . . . . . . . . . -I I I E-,i l . F-~i i Figure 2. S c h e m a t i c r e p r e s e n t a t i o n o f c h a i n g e n e r a t i o n p r o c e s s in t w o d i m e n s i o n s . T h r e e s e g m e n t s o f the c h a i n g e n e r a t e d u s i n g w i n d o w s I, II a n d III a r e s h o w n b y d a s h e d , solid a n d d a s h e d lines respectively. IIr I T -J 1_. . . . . . . . . . . . . . . . . probability of hitting the boundary of LAXand therefore would have to be artificially coiled into LAT.Thus we eliminate boundary effects. (ii) Define the next window, which is first wiped clean, such that its centre maps into the last point of the previous segment generated, as shown in figure 2 for the case of two dimensional walk. Initialize the window with respect to the previous segment(s) generated. This is done by mapping the coordinates stored in COORD back into the window-indices and checking if they are within the bounds of LAX.If they are, then those locations of LAXare initialized with data 1. Otherwise, they are ignored. One more segment is generated in this window. (iii) Process (ii) is repeated till the arrays COORDand SRANDare full. Once they are full, they are dumped into a floppy file, and are initialized for fresh use. (iv) Process (iii) is repeated till a walk of desired length has been generated. (v) If in the process of generating this walk, we encounter a situation where the walk cannot proceed further because it has reached a site which has no unoccupied neighbours, then we back-track on the walk, keeping track of which window we are in, to a point from where the walk could proceed along a different path. This is done by scanning the array NRANOfrom the current location towards the top till we encounter a location which contains p >/2. We determine the window in which this site was generated. Proceed the walk from that site making sure that we avoid those locations of this window which mapped onto points which got the walk blocked. It may be noted that window initialization and some time back-tracking also may require that we retreive the coordinates and values from the floppy files back into the arrays COORDand NRAND. Our algorithm can take care of these data movements back and forth between the buffers COOROand NRANDand their corresponding files on floppy diskette. Once a walk is generated, it is weighted according to the scheme of Rosenbluth and Rosenbluth (1955) in the following way: Let the walk consist of N steps. Let p(i) be the number of options available for generating the ith link. (2 <~ i <~ N; 0 ~ p(i) < 5 for a cubic walk). These values are retreived from the array NRAND.Then the weight associated with the walk is given by: W = I-I p(i) i=2 (3) 5 The algorithm also calculates R 2, the square of end-to-end distance of the walk. To get an average value ( R 2 ) of R 2, we generate a sample consisting of large number, say n of such walks. Simulatton of polymer systems 4. 39 Results The above algorithm has been used to calculate* the mean square end-to-end distance ( R 2 ) as a function of the number, N, of monomers in the chain for a cubic lattice**. This calculation essentially gives information about the dependence of ( R 2 ) on N in a dilute polymer solution. Calculations have been done for N = 8, 16, 32, 48, 64, 80, 128. The window used in these calculations had a size of 31 x 31 × 3I. For a given position of the window, only a small segment (16 monomers) of the walk was generated to ensure that the boundaries o f the window are not encountered in the process of generating the walks. The mean value ( R 2 ) was evaluated using a large number, n, of walks by either giving equal weightage to all the walks or by weighting them as per the procedure outlined in earlier section. To get an accurate value of ( R 2 ), it is desirable to have as large a sample as possible. Restriction to sample size comes from the computation time. The optimum value of n was arrived at by calculating ( R 2 ) as a function of n. Figure 3 shows ( R 2 ) as function n for a 128 monomers chain for the case when all the chains were given equal weightage. It may be noted that though the value of ( R 2 ) fluctuates for small values of n, it stabilizes to within 2% for n = 750. For smaller walks (N < 128), the accuracy on ( R 2 ) is expected to be better than 2 ~/o for n = 750. In view of this, when all the walks had same weightage we used sets of 750 walks to compute the mean value of R 2. In the case of weighted chains, we have used larger samples as it was seen that some walks acquired very low weightage in the process of weighting. The sample sizes used in various calculations are given in table 1. Also given in the same table are the calculated values of ( R 2 ) for different values of N both for weighted and unweighted situations. The dependence of ( R 2 ) on N is shown in figure 4, where ( R 2 ) is plotted against N p V 8 c s • • • • e B • B • c Choin length N = 128 Lattice : Cubic No weighting a) N - - - - 0 [ I _.A__ _ _ _ 1000 [ I 20oo Somple size n Figure 3. Variation of mean square end-end distance, ( R 2 ), as a function of the sample size n. * Calculations performed on the super-mini computer, PRIME 450. ** Dependence of (R2) on N is known (Domb 1963) to be independent of the lattice used. 40 S L Narasimhan, P S Goyal and B A Dasannacharya Table 1. Calculated values of mean square end-to-end distance (R 2 ) for a polymer chain for different values of N. Unweighted Weighted Chain length (N) Sample size (n) ( R2) 8 16 32 48 64 80 96 112 128 750 750 750 750 750 750 750 750 24.14 51.74 80.68 104.08 131.91 16t.64 188.45 215.23 Sample size (n) ( R2) 9999 8000 8000 5000 5000 6000 10.93 27.26 66.97 101.91 129.73 167.83 I000 333.10 / 6,0 @ UnIweighted samples - 5-0 A - / / 1.11 N 1'09 / / j _ 0 Weighted samples _ _ N1.192 4.0 / 0C V c 3.0 2.0 10 L I 2.0 3:0 .... I 4.0 I 5.0 In N Figure 4. Variation of mean square end to end distance,(R2), for a single polymer chain as a function of number of monomers N. Filled circles correspond to unweighted average and open circles correspond to weighted average as discussed in the text. on a log-log scale. It m a y be m e n t i o n e d that for real p o l y m e r systems, to get an a p p r o p r i a t e value o f ( R 2 ) using M o n t e C a r l o m e t h o d , one has to m a k e sure that the walks used in c o m p u t i n g ( R 2 ) are generated with equal probability. In view o f this, the open circles represent a m o r e relevant situation. 5. Discussion It is well established b o t h experimentally a n d theoretically (De G e n n e s 1979) that the mean square end to end distance ( R 2 ) for a linear p o l y m e r chain having N m o n o m e r s is given by an expression o f the type: (R 2 > = AN ~ Simulation of polymer systems 41 where the exponent v depends on the environment in which polymer chain resides. For example, the value of v is 1.2 for dilute polymer solutions (De Gennes 1979). Our calculations, where we have generated an isolated polymer chain in an infinitely large vessel, essentially corresponds to dilute polymer solution. It is seen that calculated (open circles in figure 4) values of ( R 2 ) lie on the dashed line ( R 2 ) = N 1.192. That is, our algorithm gives v = 1' 192 in agreement with earlier simulation studies (Rosenbluth and Rosenbluth 1955) as well as light and neutron scattering experiments on dilute polymer solutions (Flory 1953; Cotton et al 1974). We also note that the solid points, corresponding to unweighted samples, give v = 1"09 again in agreement with earlier simulation work (Rosenbluth and Rosenbluth 1955). In short, we find that the above algorithm for examining the chain configurations in polymer systems using Monte Carlo method gives results similar to those obtained using existing algorithms. Using the present algorithm, however, it is possible to examine very long polymer chains (N >/1000) even on a small computer without encountering boundary effects, which is not possible with the other existing algorithms. Moreover, as this algorithm keeps track of all the occupied neighbours of each monomer of the chain, unlike earlier algorithm, it can be easily used for examining concentrated polymer solutions and polymer melts. In particular, using the above algorithm we are examining the configurations of a long (N monomer) polymer chain in a melt of small (P monomer) chains (Joanny et al 1981). References Cotton J P, Decker D, Benolt H, Farnoux B, Higgins J S, Jannik G, O ~ r R, Picot C and des Coiseaux J 1974 J. Macromol. 7 863 De Gennes P G I979 Scaling concepts in polymer physics (Ithaca, New York: Cornell University press) Domb C 1963 J. Chem. Phys. 38 2957 Flory P J 1953 Principles of polymer chemistry (Ithaca, New York: Cornell University press) Higgins J S and Stein R S 1978 J. Appl. Crystallogr. 11 346 Joanny J H, Grant P, Turkevich L A and Pincus P 1981 J. Phys. 42 1045 Mckenzie D S 1976 Phys. Rep. C27 35 Rosenbluth M N and Rosenbluth A W 1955 J. Chem. Phys. 23 356 Wall F T 1954 J. Chem. Phys. 22 1036 Wall F T, Windwer S and Guns P J 1963 Methods in computational physics (New York: Academic Press) vol. 1 pp. 217 Wall F T, Mandel F and Chin J C 1976 J. Chem. Phys. 63 4592