Transcript Rump-Shake Connection

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
[1]
Jake Pushie
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]
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]
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
A=
(
cosdfcosdy-sindfcosdqsindy
sindfcosdy+cosdfcosdqsindy
sindfsindy
-cosdfsindy-sindfcosdqcosdy
-sindfsindy+cosdfcosdqcosdy
sindfcosdy
-cosdfsindq
sindfsindq
cosdq
* 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]
)
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 4space - 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]
Section 8.5.1
* This relationship between the quarternion components and Euler
angles allows us to redefine the rotation matrix as:
A=
(
q02 + q12 - q22 - q32
2(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
)
* The only real “added” overhead here is that this method requires the
generation of 4 random components for the 4D unit vector.
[6]
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]
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]
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]
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 ]
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 ]
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 ]
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.
1. End rotation
2. Kink-jump
3. Crankshaft rotation
[ 13 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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.
dri =
Afi
kBT
+ driG
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 ]
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 ]