Real and Biased MC Simulations of Molecules and Polymers

1. Real and Biased MC Simulations of Molecules and Polymers CHEM 699.08 2nd Ed.: 8.5, 8.6, 8.7 1st Ed.: 7.5, 7.6, 7.7 20.09.01 Jake Pushie [ 1 ]

2. Section 8.5.1 Within the bounds of the MC simulation the orientations and positions in space of the components in the system are varied during each MC step. * The “components” themselves may be atoms, or molecules. * easy not as easy Translation: described by the center of mass. * Rotation: rotate about one of the Cartesian axes (x,y,z) by an angle. * The angle of rotation, dw, is chosen based on criteria such that dw <= dwmax * [ 2 ]

3. Section 8.5.1 The process of rotation is described by applying trigonometic relationships. If the initial orientation is described by (xi, yj, zk), then the new orientation, (x’i, y’j, z’k), generated following rotation by dw about the x-axis is calculated as: * ( ) ( ) ( ) x’ y’ z’ 1 0 0 0 cosdw sindw 0 -sindw cosdw x y z = The Euler angles (f, q, y) are often invoked to describe orientations of a molecule. * [ 3 ]

4. Section 8.5.1 If all f, q, and y angles are changed by a small amount, say df, dq, and dy then the molecule is moved according to the matrix equation: * vnew = A vold ( ) cosdfcosdy-sindfcosdqsindy sindfcosdy+cosdfcosdqsindy sindfsindy -cosdfsindy-sindfcosdqcosdy -sindfsindy+cosdfcosdqcosdy sindfcosdy sindfsindq-cosdfsindq cosdq A = A note on implementation: Sampling from q directly results in an uneven distribution of data points as q becomes larger. To correct for this (or more appropriately to “correctly” do this) one must sample from cosq. * [ 4 ]

5. Section 8.5.1 Unfortunately, the use of Euler angles results in a computationally expensive calculation of the rotation matrix, as there are a total of six trigonometric functions to evaluate on each iteration (sine and cosine for each of f, q, and y angles). * An alternative approach is to treat the problem as a vector in 4-space - given the condition that the components of the vector sum to 1: q02 + q12 + q22 + q32 = 1 * The quarternion components are related to the Euler angles by: * q0 = cos½qcos½(q + y) q1 = sin½qcos½(q + y) q2 = sin½qsin½(q + y) q3 = cos½qsin½(q + y) [ 5 ]

6. Section 8.5.1 This relationship between the quarternion components and Euler angles allows us to redefine the rotation matrix as: * ( ) q02 + q12 - q22 - q322(q1q2 + q0q3) 2(q1q3 - q0q2) 2(q1q2 - q0q3) q02 - q12 + q22 - q32 2(q2q3 + q0q1) - 2(q1q3 + q0q2) 2(q2q3 - q0q1) q02 - q12 - q22 + q32 A = The only real “added” overhead here is that this method requires the generation of 4 random components for the 4D unit vector. * [ 6 ]

7. Section 8.5.2 MC simulations of molecules are difficult, in that they require long computation times unless some internal degrees of freedom are frozen, or unless the system is very small. * One method used for MC simulations is to translate and rotate the molecule, while making random changes to the Cartesian coordinates of individual atoms. * Very small displacements are necessary: even small changes to atomic coordinates can give rise to large energies and a high frequency of rejection during the simulation. * [ 7 ]

8. Section 8.5.2 Another method is to freeze the bond lengths and/or angles. * Drawbacks: Distribution of internal conformations may be affected. Small bond rotations give rise to large displacements for long molecular chains. [ 8 ]

9. Section 8.6 Polymers may be any macromolecule which is constructed from smaller subunits. Proteins fall into this category as well. * Modeling polymers requires that large systems are needed to accurately describe a polymer’s behaviour. This implies that the duration of the simulation must also be sufficiently long in order to explore the large conformational space. * For such large systems empirical models are used to reduce calculation time -- although present increases in computational resources have allowed QM methods to be applied to large simulations. * [ 9 ]

10. Section 8.6.1 The choice of how the system will be represented is governed by the type of problem at hand. * In the lattice model representation each effective monomer occupies one of the vertices. Upon each successive iteration one of the monomers may occupy an adjacent, unoccupied, vertex. * [ 10 ]

11. Section 8.6.1 In more sophisticated models each vertex may represent several chemical bonds, atoms between the vertices are assumed to be rigid. * One of the most common applications of the lattice systems is to perform a “random walk” where the polymer chain grows until it contains the desired number of monomers. * = subunit Polymer length = 8 [ 11 ]

12. Section 8.6.1 Simple properties which can be determined from the random walk: * Polymer size: the mean square end-to-end distance <Rn2> Radius of gyration: r.m.s. distance from the centre of mass <s2> = <Rn2>/6 Of course once the polymer is ‘grown’ into the lattice it is useful to explore alternate configurations of the system. * [ 12 ]

13. Section 8.6.1 The Verdier-Stockmayer algorithm (1962) generates new configurations of the system based on combinations of three simple manipulations. * = Polymer subunit 1. 2. 3. End rotation 1. Kink-jump 2. Crankshaft rotation 3. [ 13 ]

14. Section 8.6.1 Another relatively simple algorithm to implement is the “slithering snake” model. * One end of the polymer chain is randomly chosen as the “head”, and it is advanced to one of its adjacent vertices. Each previous subunits is then advanced to that of its successor’s position. * [ 14 ]

15. Section 8.6.2 The simplest of the continuous polymer methods uses a string of connected beads to represent the polymer chain. The beads are the effective monomers, while to links are the effective bonds. * The links between the effective monomers may be fixed, or may be allowed to fluctuate using a harmonic potential. * With the linked bead type of model the system is free to sample from a continuum of possible configurations; the pivot algorithm is one way this is accomplished. * Unfortunately this type of model treats rotation around the individual bonds in such a way that any torsion angle from 0 to 360o is equally likely. * [ 15 ]

16. Section 8.6.2 To increase the complexity (and hopefully the accuracy) of the model, one must move to an atomistic model of the polymer chain: which is the closest model in reality to the actual polymer. * One of the drawbacks to the atomistic representation is the necessity for a realistic initial state for the entire system. This must be “close” to the state from which properties will be derived. * [ 16 ]

17. Section 8.7 All things being equal a “real” MC simulation should not have any preference for exploring only one region of phase space over another. * A “biased” MC simulation, on the other hand, can explore a specific region of phase space that may be of importance to the operator. * Example: solute-solvent interactions, where solvent molecules further away from the solute should be treated differently from those in the primary (or secondary, … etc.) solvation shell. * Preferential sampling can be employed, whereby the solvent molecules in the vicinity of the solute move more frequently than those further away. * [ 17 ]

18. Section 8.7 The implementation in such a case can be done by defining a cutoff region around the solute molecule, wherein any solvent molecule is moved more often than those outside. * The probability that one of the molecules will be moved is determined by a probability parameter p. * randomly choose molecule if molecule is inside cutoff distance, move it otherwise if pi > the random number [0,1] a trial move is attempted for molecule i otherwise don’t move molecule i So the closer p is to zero, the more often the closer molecules are moved than the further away molecules * [ 18 ]

19. Section 8.7 In applying preferential sampling one must be sure that a procedure is followed in deciding whether a move will be accepted or rejected. * We must apply the principle of microscopic reversibility. * Why does this come into play? * pinner = <always move> pouter = 0.2 [ 19 ]

20. Section 8.7 In applying preferential sampling one must be sure that a procedure is followed in deciding whether a move will be accepted or rejected. * We must apply the principle of microscopic reversibility. * Why does this come into play? * pinner = <always move> pouter = 0.2 [ 19 ]

21. Section 8.7 In applying preferential sampling one must be sure that a procedure is followed in deciding whether a move will be accepted or rejected. * We must apply the principle of microscopic reversibility. * Why does this come into play? * pinner = <always move> pouter = 0.2 [ 19 ]

22. Section 8.7 In applying preferential sampling one must be sure that a procedure is followed in deciding whether a move will be accepted or rejected. * We must apply the principle of microscopic reversibility. * Why does this come into play? * pinner = <always move> pouter = 0.2 This must be taken into account when devising acceptance criteria. * [ 19 ]

23. Section 8.7 Another type of biased MC simulation is the force-bias Monte Carlo method. [Pangali et al. 1978; Rao et al. 1979] * Here, the forces are calculated for the randomly selected molecule, then the molecule moves not in a wholly random matter, but in the direction of the forces on it. * The exact direction the molecule moves is determined by a probability distribution whose peak is in the direction of of the calculated force. * The smart Monte Carlo method [Rossky et al. 1978] also requires that forces be calculated, but the displacement of the atom or molecule has more of an element of “randomness” to it * [ 20 ]

24. Section 8.7 There are two components to the smart Monte Carlo method. The first is the force on atom or molecule, the second component is a random vector drG. * Afi dri = + driG kBT fi = force on the atom/molecule A = “a parameter” (fudge factor?) driG = random displacement This method, unlike the force-biased MC method does not restrict the magnitude of displacement that an atom or molecule may undergo. * [ 21 ]

25. Section 8.7 In practice, the force-bias and smart MC methods perform very similarly, and there is little to choose between them. * The necessity to calculate the forces makes these last two methods much more elaborate and computationally expensive. * The “perks” of these methods are that they both lead to much higher acceptance rates for each iteration, and also allow for larger displacements per iteration. * ~ fin ~ [ 22 ]