University of Southern California Dissertations and Theses

Nonlocal Monte Carlo approaches: helium clusters and nucleic acids

(USC Thesis Other)

Nonlocal Monte Carlo approaches: helium clusters and nucleic acids

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content NONLOCAL MONTE CARLO APPROACHES: HELIUM CLUSTERS AND
NUCLEIC ACIDS
by
Nikolay Markovskiy
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulllment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(CHEMISTRY)
December 2008
Copyright 2008 Nikolay Markovskiy
Dedication
to my parents and my sister Irina
ii
Acknowledgements
I would like to thank my scientic advisor Professor Chi H. Mak for his constant
support, patience, and guidance as well as his invaluable advices. While I was con-
ducting my research studies, we had a lot of insightful and inspiring discussions.
Thanks to him I have gained many valuable programming skills, developed a
steady programming approach towards scientic problems, and greatly improved
my ability to think and work independently. He has taught me a great deal of
how to present scientic results in a clear and accessible manner and I thank him
for that.
I am grateful to Professor Stephen Bradforth whom I communicated with
regarding my enrollment to the University of Southern California while I was still
in Russia and who was very supportive answering all my questions and providing
valuable recommendations regarding my future. I would like to thank professor
Hanna Reisler for her guidance during my rst year of study and her willingness
to help while I was adjusting to this new for me eld of physical chemistry. I am
iii
also very grateful to professor Anna Krylov, who was very kind to let me stay in
her group as soon as I had arrived. I had the great opportunity to participate
in her group's research work and get to know her students Lyudmila Slipchenko
and Sergey Levchenko, whom I have learned a lot from.
My special thanks go to professor Howard Taylor for his advice and interest
in my scientic career.
I would like to thank my group members Sergey Zakharov for his help when
I just joined the group, Brian Cheney for his suggestions on my presentations,
and Arun Sharma for all his support on dierent projects and proofreading our
manuscripts during my years of study.
There is number of graduate students whom I would like to thank. I am
very grateful to Russell Sliter from the group of Andrey Vilesov for proof-reading
some parts of this work as well as Chris Nemirow and George Kumi from the
group of Curt Wittig. I had a lot of insightful discussions with Prof. Vilesov's
group members Dmitry Skvortsov and Kirill Kuyanov regarding experimental
side of our research, and I thank them for that. I am also very grateful to Sergey
Malyk, Igor Fedorov, Lukasz Koziol, Anton Zadorozhnyy, Piotr Pieniazek, Askat
Jailaubekov, Roman Rabinovich, Misha Ryazanov and Vadim Mozhayskiy for
their help and many fruitful discussions.
A word of appreciation goes to the department of Chemistry at the University
of Southern California, particulary to Yuki Yabuta, Valery Childress, Heather
iv
Meunier Connor and Michele Dea for all their assistance with paperwork, and
formalities, making my stay at the university enjoyable, scholastic experience.
v
Table of Contents
Dedication ii
Acknowledgements iii
List of Figures viii
Abstract xii
Chapter 1: Introduction 1
1.1 Scope of the Work and Motivation . . . . . . . . . . . . . . . . . 1
1.2 Organization of the Thesis . . . . . . . . . . . . . . . . . . . . . . 5
Chapter 2: Rotational Dynamics of Spherical Top Molecules in
Helium Clusters 7
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Theory and Numerical Implementation . . . . . . . . . . . . . . . 10
2.2.1 Path Integral . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.2 Multilevel Algorithm for Rotations . . . . . . . . . . . . . 16
2.2.3 Details of Implementation . . . . . . . . . . . . . . . . . . 19
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.5 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . 48
Chapter 3: A generalized kinematic formulation of loop closures in
nucleic acids, proteins and other linear polymers 52
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.2 Boundary Constraints in the Loop Closure Problem . . . . . . . . 55
3.3 Solutions to the Generalized Closure Problem . . . . . . . . . . . 63
3.3.1 Part 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.3.2 Part 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.3.3 Generalized Closure Involving More Than Six Degrees of
Freedom . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
vi
3.4 Numerical Implementation of the Generalized Loop Closure Solution 81
3.4.1 An Example . . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.5 Jacobian of the Transformation . . . . . . . . . . . . . . . . . . . 85
3.6 Monte Carlo Step and Preliminary Results . . . . . . . . . . . . . 91
3.7 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . 97
Bibliography 99
Appendices 106
Appendix A: Rotations in Four Dimensions . . . . . . . . . . . . . . 107
Appendix B: Extremum's Search of the Rotational Propagator . . . . 109
Appendix C: Rejection Sampling . . . . . . . . . . . . . . . . . . . . 111
vii
List of Figures
Figure 2.1: Examples of the orientation correlation function for normal
CH
4
, 2CH
4
and 50CH
4
the gas phase, as well as 50CH
4
inside
a shell of twelve
4
He atoms. . . . . . . . . . . . . . . . . . . . . . 26
Figure 2.2: Orientation correlation function of 50CH
4
in a shell of
twelve
4
He atoms (solid line) and the best t to a rigid spherical top
model (dashed line), compared to the non-exchange result (dotted
line) and the gas-phase (dotted dashed line). . . . . . . . . . . . . 27
Figure 2.3: Angular density prole of the
4
He shell in a (a) 50CH
4
-
(He)
12
and (b) 1CH
4
-(He)
12
cluster. H show the positions of the
hydrogen atoms. The
4
He density is strongly templated by the
heavy probe, while it is almost completely isotropic for the light
probe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
Figure 2.4: Eective B constant for (a) 50CH
4
and (b) 1CH
4
in
clusters with dierent number of
4
He atoms, from 1 to 20. For
50CH
4
, the renormalization of the B constant by the
4
He seems
to reach an asymptotic limit at around 6
4
He atoms, with only
minor variations thereafter. For 1CH
4
, the change inB with the
number of
4
He atoms is much more gradual. . . . . . . . . . . . . 31
Figure 2.5: EectiveB constant for SF
6
in clusters with dierent num-
ber of
4
He atoms, from 1 to 8 computed by several methods and
compared with experimental nanodroplet value [Hartmann et al.,
1995], [Harms et al., 1997] B
exp
= 0:034 cm
1
. The renormaliza-
tion of the B constant by the
4
He seems to reach an asymptotic
limit at around 8
4
He atoms.
1
[Kwon et al., 2000] . . . . . . . . 32
Figure 2.6: The super
uid fraction of the helium shell for 50 to 1CH
4
with 12
4
He atoms. Due to the symmetry of the underlying probe,
the super
uid fraction of the shell is close to unity. . . . . . . . . 34
viii
Figure 2.7: The two-particle angular distribution function for (a) 1CH
4
and (b) 50CH
4
with full exchange (solid lines) and no exchange
(dotted lines). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Figure 2.8: Orientation correlation function of 1CH
4
in a shell of
twelve
4
He atoms (solid line)compared to the gas-phase (dotted
line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
Figure 2.9: The percentage change in the eective B constant as a
function of the intrinsic moment of inertia of the probe with full
exchange (circles) and no exchange (squares). . . . . . . . . . . . 39
Figure 2.10: A planar model describing the interaction between the
rotating probe and a rigid shell. The left panel shows a schematic
picture of the the spectrum of the probe rotations. The right panel
shows the low-energy spectrum of the shell and how the density
of states of the shell is dilated when bose symmetry is imposed on
the shell wavefunctions. . . . . . . . . . . . . . . . . . . . . . . . . 42
Figure 2.11: 2D Toy Model. The percentage change in the eective B
constant, obtained from nonlinear t with 2D free rotor orientation
correlation function, as a function of the intrinsic moment of inertia
of the probe with full exchange (circles) and no exchange (squares). 47
Figure 3.1: (a) Drawing showing boundary conditions for the closure
problem. Atoms (i 1), (i) and (i + 1) on the left are xed in
space. The black dashed line indicates that the rest of the chain
to the left of atom (i 1) are also xed. The grey line represents
the movable segment whose conformation is
exible. The atoms
(j 1), (j), (j + 1) and all the atoms to the right represented
by the black dashed line form a rigid segment. (b) When the
conformation of the movable segment is changed, the position and
orientation of the right rigid segment will change, while the left
boundary segment remains xed in space. (c) The minimal number
of torsion angles needed in the movable segment is six. The picture
shows the minimal length if every bond along the movable segment
is torsionable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
ix
Figure 3.2: (a) Similar to Fig. 3.1(c), gure (a) shows details of the
movable segment. But in this more general case, the parts of the
chain inside the movable segment indicated by the thin solid lines
are rigid. (b) To account for the rigid sections inside the movable
segment, we connect atoms (i+1) and (m1), (m+1) and (n1),
and (n + 1) and (j 1) by three ctitious bonds indicated by the
white sticks. (c) An equivalent representation to (b), but using
the three ctitious bonds to represent the constraints introduced
by the rigid sections. . . . . . . . . . . . . . . . . . . . . . . . . . 60
Figure 3.3: Denition of the four RETO coordinates using the ctitious
bond for a section with rigid backbone indicated by the thin solid
lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
Figure 3.4: Denition of the RETO
X
coordinates to represent the
boundary conditions in the closure problem shown in Fig. 3.2(c). 64
Figure 3.5: Drawings showing dierent steps in Part 1 of the solution.
See text for details. . . . . . . . . . . . . . . . . . . . . . . . . . 67
Figure 3.6: Drawing showing the rotation of the C-element in detail.
The bond axis between atoms (2) and (4) becomes the z
C
axis.
Thex
C
-axis is dened so that the atom (1
0
) is on thex
C
-z
C
plane.
The atom (3) is initially placed on the x
C
-z
C
plane dening the
rotation angle
C
= 0. After a rotation of the C-element by
C
,
the nal bond angle (1
0
)-(2)-(3) is . See text for further details. 75
Figure 3.7: Drawing showing two dierent solutions of the generalized
closure problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
Figure 3.8: Drawing showing how the torsion angle
1
is related to

X
1
and
L
. The picture is drawn looking down the (1)-(0) bond
axis. Atom (0) is behind atom (1) (both in grey).
X
1
is dened as
the torsion angle (1)-(0)-(1)-(5), whereas the original
1
is the
(1)-(0)-(1)-(1
0
) angle. . . . . . . . . . . . . . . . . . . . . . . . 80
Figure 3.9: Two solutions of the generalized closure problem for a 4-
nucleotide RNA segment (to reduce clutter, O1P and O2P are not
shown). Atoms (0) through (2) form the rigid L-element, atoms
(2) through (4) the C-element and atoms (4) through (6) the R-
element. The L- and R-element consist of one nucleotide each.
The C-element consists of two nucleotides. . . . . . . . . . . . . . 84
x
Figure 3.10: Solutions of the generalized closure problem for the exam-
ple in Sect. 3.4.1. (a) Solution for
R
from the rst step in Part 1 of
the solution. For every
L
there are two solutions for
R
in general.
Squares denote domain boundaries (see text). (b) Solution for
C
from the second step in Part 1. For each
R
branch in (a), there
are in general two solutions for
C
. The solid and dashed lines in
(b) correspond to the two solutions from the
R
branch associated
with the solid line in (a), and the dotted and dotted dashed lines
in (b) to the
R
branch associated with the dotted dashed line in
(a). (c) Solution for(
L
)

as a function of
L
from the second
step in Part 1 of the solution. Each line (c) corresponds to a
C
branch in (b). The roots are shown as the open circles. The root
corresponding to the native structure is shown as the closed circle. 86
Figure 3.11: Schematic drawing showing the L-, C- and R-elements and
the denitions of the lengths q
i
, bond angles
i
and torsion angles

i
needed to calculate the Jacobian of the RETO
X
$ torsion angle
transformation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
Figure 3.12: RNA example. Drivable backbone torsion angles in a
single residue are in solid. Dashed torsion angles can be driven
when they are not part of the generic closure problem, e.g. when
a RETO element consists of three residues and this residue is in
the middle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
Figure 3.13: Original conguration of tetrahymena ribozyme [Guo et al.,
2004] molecule (a). An example of a new conguration after the
MC move is applied (b). Encircled area is the range where the
boundary conditions are set up and the generic closure meth-
ods is applied. Right and left RETO elements consist of up to
two nucleotides, the central RETO element (stem) consists of 75
nucleotides. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
Figure C.1: Rejection method for generating a random sample x from
a known probability distributionp(x) that is everywhere less than
some other functionf(x). The transformation method is rst used
to generate a random sample x of the distribution f. A second
uniform sample is used to decide whether to accept or reject that
x. If it is rejected, a new sample of f is found; and so on. The
ratio of accepted to rejected points is the ratio of the area under p
to the area betweenp andf. This gure and caption are borrowed
from [Press et al., 1995]. . . . . . . . . . . . . . . . . . . . . . . . 112
xi
Abstract
One of the main challenges in the simulations of many body systems in statistical
physics is the problem of generating sucient coverage of phase space in a rea-
sonable amount of computer time to produce accurate estimates of equilibrium
properties. This work is focused on constructing and applying two dierent non-
local Monte Carlo techniques which incorporate large-scale changes in a number
of degrees of freedom. These methods were then used to study the rotational
dynamics of spherical top molecules in
4
He clusters and the problem of confor-
mational sampling of RNA molecules.
Path integral simulations have been carried out to study the rotations of a
methane inside a single shell of
4
He atoms at 0.3 K to address the question of
whether dopant molecule rotations can be used to probe the quantum statistics
and super
uidity of the shell. We examined the eects of the probe molecule on
the
4
He exchanges and their counter eects on the renormalized rotation constant
of the probe systematically by varying the intrinsic moment of inertia of the
xii
methane. The observed eects show strong dependence on the instrinsic moment
of inertia of the rotating probe, with a heavy probe favoring stronger templating
of the
4
He density and a corresponding suppression of exchanges in the shell, as
well as a large renormalization in the probe's eective rotation constant, while a
light probe shows almost no eect on the shell density or the eective rotation
constant. These results can be rationalized in terms of a rotational smearing
eect and suggest that there is no clearly quantiable relationship between the
super
uid fraction of the shell and the renormalized rotation constant of the
probe for cases where the probe molecule has weak anisotropic interactions with
the
4
He atoms.
The problem of loop closure is important for the conformational sampling
of chain molecules like nucleic acids and proteins. Solutions to the loop closure
problem can be used to generate alternative structures for an internal segment of a
chain when the atomic coordinates of the rest of the chain are constrained at xed
positions. This work describes a new loop closure formulation for linear polymers
based on the kinematic geometry of rigid bodies. This formulation makes use of
a novel representation of the boundary constraints, allowing the loop closure
problem to be reduced from the conventional 6-variable/6-constraint problem to
a 4-variable/4-constraint problem. With this new formulation, loop closure can
then be recast into a one-variable root search, permitting a simple analytical
or numerical solution. The new formulation provides a natural framework for
xiii
closing loop structures of arbitrary sizes, from as small as a ve-member ring to
those containing multiple segments of arbitrary lengths.
xiv
Chapter 1
Introduction
1.1 Scope of the Work and Motivation
The Monte Carlo (MC) method is applied in various elds of science and engi-
neering to determine the properties of a broad variety of systems. In this work we
focus on simulation of two dierent systems one described by quantum and the
other by classical equilibrium statistical mechanics. Unlike Molecular Dynamics
(MD) methods, a MC simulation does not propagate a system in a dynamical
sense, i.e. by numerically solving Newton's equations. Two adjacent MC cong-
urations are not necessarily connected by any dynamical physical process. The
only requirement is that they are representatives of a proper statistical ensem-
ble, or using the probability language, are unbiased samples of the corresponding
statistical distribution. The transition rule between two congurations along the
Markov chain is arbitrary as long as detailed balance is obeyed. This freedom
allows us to employ any stochastic process to sample multidimensional congu-
rational space, as long as it visits the statistically most important states with the
highest contributions according to the to the correct statistical distribution. In a
1
canonical ensemble this ensures that the energy states which have a statistically
signicant population at a given temperature should have the largest contribu-
tions to the Boltzmann sum. The challenge in nding such states is that the
number of possible congurations rises exponentially with the number of degrees
of freedom as well as the number of conformers in case of molecular simulations,
and this leads to a highly "rugged" potential energy prole with multiple local
minima possibly widely separated by large barriers.
It is not very surprising that the natural importance criteria of choosing
(accepting) the next state in the random walk is the relative probability which
depends on the energy dierence between the states, and it favors going downhill
in the energy landscape but also allows going uphill if the energy change is not
very large relative to the temperature [Metropolis et al., 1953]. The simplest
transition rule (move) between congurations is a random change of coordinates.
If the change is large and nonlocal, the importance condition is rarely met and
new congurations are hardly ever accepted. With small random changes, one
risks getting stuck in one local minima during the simulation and obtaining unre-
liable statistics. This is similar in a way to the notorious "critical slowing down"
in discrete case of model magnets, which was initially solved by [Swendsen and
Wang, 1987] via cluster Monte Carlo algorithm, and was generalized later towards
more complicated cases [Dikovsky, 2003, Mak and Zakharov, 2004, and references
therein]. The main idea is to design a transition rule, which independently moves
2
groups of particles (clusters) on multiple length scales as a whole with minimal
eect on the energy criteria. The nature of such nonlocal moves can be very sys-
tem specic and depends on the symmetry, type of interactions, dimensionality,
physics, etc. The a priori knowledge about a system in question and the ability to
solve some part of it rigorously can be of great use, since it does not only facilitate
designing such nonlocal moves, but can also eectively reduce dimensionality of
a problem. Chapter 3 will discuss the implementation of such nonlocal moves in
nucleic acids utilizing the generic loop closure algorithm, which, for given bound-
ary conditions for a selected part of a molecule, yields specic discrete values
of some of its backbone torsion angles independently from the rest of the chain.
Chapter 2 will particulary talk about implementation of nonlocal moves with
respect to rotations of interacting spherical top molecules, making advantage of
the fact that quantum statistical physics of noninteracting rigid bodies has an
analytical solution.
This work is based on two projects. The rst project is devoted to simulation
of spherical top molecules in helium clusters. One of the fundamental goals
of cluster research is the interpretation of the properties of condensed phases
in terms of their constituent atoms and molecules, and to nd out how these
various bulk properties emerge from the nite systems. The property of particular
interest for us is the phenomenon of super
uidity in helium. We have conducted
studies of rotations of spherical top molecules such as CH
4
, CD
4
, and SF
6
in
3
helium clusters which where previously used as molecular probes in spectroscopic
experiments [Hartmann et al., 1996, Nauta and Miller, 2001a, Skvortsov et al.,
2007, Rudolph et al., 2007] to investigate the onset of super
uidity in helium.
The objective of the research is to obtain the rotational spectral features and
investigate the relation of the intrinsic moments of inertia of the probe molecules
to the onset of super
uidity as well as to study the eect of the super
uidic
He shell on their rotations. This led to development of a new Path Integral
Monte Carlo (PIMC) technique which has not been applied to rotations before.
In addition, to the best of our knowledge, it is the only PIMC study explicitly
including rotations of molecules with spherical top symmetry do date.
The second project is dedicated to the development of a generic Monte Carlo
method for more ecient samplings of the conformations of bio-polymers and its
application to nucleic acids specically. With the advancement of supercomput-
ers, molecular dynamics simulations are able to go up to the microsecond time
scale in tracking the motions of small biopolymer systems. But many biologi-
cally interesting phenomena such as folding of bio-polymers or enzyme catalysis
occur over several orders of magnitude greater time scales [Henzler-Wildman and
Kern, 2007, Henzler-Wildman et al., 2007]. As it was mentioned before, MC
methods have the advantage to explore conguration space rapidly by exploring
largescale rearrangements of the system. The simplest example of such is a small
change of a backbone torsion angle of a polymer, which may lead to a very large
4
change in cartesian coordinates of a section of the chain following the rotated
bond. Unfortunately, such rearrangements may also result in many energetically
prohibitive steric clashes [Dodd et al., 1993, and references therein] making such
moves highly inecient. Therefore, it is necessary to develop MC moves which
are able to modify the internal degrees of freedom in a regulated fashion over
selected parts of the chain leaving the rest unchanged. The generic loop clo-
sure algorithm [Mak, 2007, 2008] allows implementation of such MC moves at
arbitrary length scales.
1.2 Organization of the Thesis
Chapter 2 is organized as follows. First, we brie
y review the path integral
method and give a computationally ecient treatment of the path integrals
involved in the rotation of a spherical top molecule. After the theoretical back-
ground, we discuss the importance of a multilevel approach for simulating the
path integrals, its formulation, and its implementation. In Sect. 2.2, we explain
the methodology of extracting rotational constants from imaginary time orienta-
tion correlation function and provide several examples. In Sect. 2.3, starting with
a \heavy" methane with a moment of inertia 50 times larger than normal CH
4
, we
measured the eective rotational constants and the corresponding angular helium
density proles with and without exchange. Then, the moment of inertia of the
5
methane is gradually increased in order to investigate the eects of the probe
on the shell and the counter-eects of the shell on the probe's rotational con-
stant. We also study cluster size eects on methane,"heavy" methane, and SF
6
molecules. The observed eects are discussed and rationalized by a qualitative
two-dimensional model in Sect. 2.4.
In chapter 3, a Monte Carlo technique based on the generic loop closure algo-
rithm is discussed. Section 3.2 provides detailed description of the loop closure
problem and discusses the necessary condition for the problem not to be ill posed.
It also introduces and illustrates the idea of its generic formulation using RETO
coordinates, which can be applied to the loops of arbitrary size. Section 3.3
is devoted to the derivation of the solution to this problem which can serve as
a recipe for its further implementation. Section 3.4 claries some issues in the
numerical implementation and elucidates the method with an RNA example. Sec-
tion 3.5 is devoted to derivation of the jacobian of the transformation from the
torsion angles' space into RETO coordinates' space. And nally, in section 3.6 we
provide an RNA example of Monte Carlo step, which can be applied to polymer
segments of arbitrary size.
6
Chapter 2
Rotational Dynamics of Spherical
Top Molecules in Helium
Clusters
2.1 Introduction
The discovery of super
uidity in liquid
4
He has attracted the interest of many.
Early theoretical work, neutron scattering experiments and more recently, quan-
tum simulations have greatly facilitated the understanding of super
uid behaviors
in the bulk [Ceperley, 1995, Feynman and Hibbs, 1965]. The peculiar properties
of super
uid
4
He is largely captured by the two
uid theory of Tisza and Lan-
dau [Tisza, 1938, Landau, 1941, 1947]. In this picture, the super
uid fraction is
dened by the nonclassical response of the
uid to a slowly rotating macroscopic
probe. Within linear response theory, the structure and the precise nature of the
probe is clearly unimportant.
7
In the area of super
uidity, molecular spectroscopy experiments have been
impossible in bulk helium liquids, because impurity molecules either tend to
aggregate, which leads to phase separation, or they become absorbed onto the
walls of the container [Silvera, 1984]. The recent advent of helium nanodroplet
isolation (HENDI) spectroscopy [Goyal et al., 1992, Toennies and Vilesov,
1998], allowed to spectroscopically probe rotational motion of single molecules
[Grebenev et al., 1998]. A molecular version of the rotating disk experiment
of Andronikashvili [Andronikashvili, 1946] have been carried out using
4
He nan-
odroplets [Toennies and Vilesov, 1998]. In these experiments, the eects of the
helium atoms on a rotating probe embedded inside the droplet was measured
spectroscopically. The earliest experiments were carried out using heavy probe
molecules, like SF
6
and OCS [Toennies and Vilesov, 1998, Grebenev et al., 1998].
The observations of sharp rotational features in helium nanoclusters were ini-
tially explained using ideas borrowed from bulk super
uidity [Kwon et al., 2000].
In general, the eective rotation constants of these heavy probes inside helium
droplets are modied by a factor of about three to ve from their gas-phase
values. These experiments have also been performed more recently employing
lighter probes [Behrens et al., 1998, Nauta and Miller, 2000, 2001b, Slipchenko
and Vilesov, 2005, Skvortsov et al., 2007]. The recent experiments reveal a much
higher degree of rotational coherence in the light probes, with their eective rota-
tional constant nearly identical to their gas-phase values. The experimentally
8
observed renormalization in the rotation constants therefore seems to depend
on the nature of the probe, a feature that is obviously contradictory to linear
response theory. Another question raised by these molecular experiments is how
much of a liquid layer is needed for the observed rotational eects to be seen.
Recent experiments using smaller bosonic clusters [Tang et al., 2002, Tang and
McKellar, 2003, Xu et al., 2003, Tang et al., 2004] suggest that similar eects are
observable with a single shell or even a partial shell of
4
He atoms.
The key question is whether \super
uidity" on the molecular level can be
explained by ideas derived from bulk super
uidity. For a molecular-level descrip-
tion, it is clear that a more complete theory must treat both parts of the quantum
system { the rotating probe and the helium shell { on equal footing instead of
relying on linear response theory applied to the helium alone. If one takes the
rotating probe to be the system whose rotational motions can be directly interro-
gated and the helium shell to be the bath, then we can look for solutions to this
system-bath quantum system. The quantum theory of system-bath interactions
is not completely solved, but since the helium shell is limited in size, one has
the advantage of being able to treat a nite bath without having to resort to an
ad hoc description of quantal dissipative eects. A detailed nite-temperature
solution of the probe+shell system can be obtained from a path integral Monte
Carlo (PIMC) simulation. Using CH
4
-He as a model, we have carefully studied
the eects of the probe's intrinsic moment of inertia and quantum statistics in the
9
4
He separately on the observed renormalization in the rotation constant, while
keeping the interaction between the probe and the shell xed. We found that the
structure of the shell depends sensitively on the moment of inertia of the probe
{ in the heavy-probe limit, the helium shell turns out to be fairly rigid and a
single solvation shell is sucient to observe super
uid-like eects. The results
also suggest that the rotational response of the probe to exchanges in the shell
does not simply follow linear response. In general, the renormalization of the
rotational constant does not directly re
ect the amount of super
uid activities in
the shell.
2.2 Theory and Numerical Implementation
2.2.1 Path Integral
The equilibrium statistics of a quantum system within the canonical ensemble can
be represented by path integrals. In the path integral picture, the density matrix
operator ^ = exp (
^
H) is rst written as a product of n high-temperature
propagators using the Trotter product formula [Trotter, 1959],
^ =
h
e

^
H
i
n
;
10
where = 1=k
B
T , T is the temperature,
^
H is the Hamiltonian, and = =n.
In analogy to real-time path integrals, is often called the \imaginary-time",
and the high-temperature propagators are thus referred to as the \short-time"
propagators.
When the system has both translational and rotational degrees of freedom,
the path-integral form of the density matrix becomes
^ (; R; R
0
) =
Z
^ (; R; R
1
) ^ (; R
n
; R
0
)dR
1
dR
n
; (2.1)
where vectors R
i
consists of rotational,

i
, and translational, r
i
, degrees of free-
dom and(; R
i
; R
i+1
) are short-time propagators of the system. In the absence
of interactions, Eqn.(2.1) is exact for any n. But when interactions are included,
the short-time propagators must be rst approximated by a \short-time approx-
imation" and the path integral is written as the limit
^ (; R; R
0
) = lim
n!1
Z
(; R; R
1
)(; R
n
; R
0
)dR
1
dR
n
: (2.2)
11
In order for Eqn.(2.2) to converge to the original density matrix, the short-time
approximation (; R; R
0
) must be correct to rst order in . One commonly
used approximation is the \primitive" discretization [Feynman and Hibbs, 1965]
(; R
i
; R
i+1
) = ^
rot
(;

i
;

i+1
)^
trans
(; r
i
; r
i+1
)
e

2
[V(R
i
)+V(R
i+1
)]
; (2.3)
whereV is the potential, and ^
rot
and ^
trans
are the exact rotational and transla-
tional propagators of the free system. The translational propagators in Eqn.(2.3)
for the center-of-mass motions can be written in the usual form [Feynman and
Hibbs, 1965]
^ (; r
i
; r
i+1
) =

m
2 h
2

1=2
exp

m
2 h
2
jr
i
r
i+1
j
2

;
where m is the mass. On the other hand, the rotational propagator takes the
form [Muser and Berne, 1996]

rot
(;
;

0
) =
X
j;k;m
e
E
jm
k
h
jjkmihjkmj

0
i (2.4)
whereh
jjkmi are asymmetric top wave functions, which can be expressed in
the basis of eigenvectors of a symmetric top molecule
q
2j+1
8
2
D
j
k;m
(;;) , and
D
j
k;m
(;;) are Wigner's functions,, and the Euler angles [Zare, 1988]. For
12
a spherical top molecule, such as CH
4
, the full rotational propagator in Eqn.(2.4)
reduces to the standard form

rot
(;
;

0
) =
1
8
2
X
j;k;m
(2j + 1)e
Bj(j+1)
D
j
k;m
(;;)D
j
k;m
(
0
;
0
;
0
); (2.5)
The three Euler angles can more conveniently be parameterized as the angles
describing points on a unit sphere embedded in four dimensions, taking the forms
[Goldstein et al., 2002]:
e
0
= cos

2
cos
+
2
e
1
= sin

2
sin

2
e
2
= sin

2
cos

2
e
3
= cos

2
sin
+
2
; (2.6)
such that e
2
0
+e
2
1
+e
2
2
+e
2
3
= 1. The rotational path integrals can therefore be
more conveniently represented as trajectories on the surface of a four-dimensional
hypersphere, spanned by these vectors. As in the case of spherical harmonics,
13
a sum rule exists for the hyperspherical harmonics, leading to the Chebyshev
polynomials of the second kind [Bander and Itzykson, 1966]:
X
k;m
D
j
k;m
(
)D
j
k;m
(

0
) =U
(2j)
(e e
0
)
In terms of the Chebyshev polynomials, the rotation propagator can nally be
written as [Bander and Itzykson, 1966, Norcli, 1972]
^
rot
(;
;

0
) =
X
j
(2j + 1)
8
2
e
Bj(j+1)
U
(2j)
(e e
0
) (2.7)
This form allows for ecient computation of the rotational density matrix during
the path-integral simulations.
For interacting systems with a nonzero potential V , the propagator can be
approximated by the primitive discretization in Eqn.(2.3) above. This primitive
approximation is straightforward to implement in a path integral simulation, but
a large n is often required for convergence. This in turn leads to slow sampling
problems which are often encountered in path integral simulations [Mak and
Zakharov, 2004]. As mentioned above, any approximate propagator correct to
rst order in ==n will generate the exact density matrix in the limitn!1.
Therefore, it is also possible to use alternative approximations which are more
accurate and require fewer discretizations n to converge. This approach was
14
taken by Ceperley et al. in a number of previous path integral simulations of
low-temperature helium [Ceperley, 1995], where a \two-particle density matrix
approximation" was used to keepn small. This approach approximately accounts
for the eect of the pair interaction potential by using the two-particle propagator
as the basis for the Trotter breakup, and the problem of two quantum particles
interacting via a central potential is solved numerically and stored before the
simulation. This approach accelerates the convergence considerably so that a
much small n may be used in the simulation, but it also becomes dicult to
implement for anisotropic systems due to operating memory constraints.
For problems with rotational degrees of freedom, the primitive approach is
preferred in order to simplify the inclusion of rotation motions into the path
integral simulations. But if one chooses to use the primitive approach, more
elaborate Monte Carlo moves must be incorporated in order to overcome the
typical slow sampling problems. A possible approach was given by Mak and
Zakharov [Mak and Zakharov, 2004]. This method uses multilevel moves to
update the paths, allowing nonlocal changes to the paths. This approach has
been implemented for the translational degrees of freedom in a simulation of
bosonic para-hydrogen [Mak et al., 2005]. However, multilevel moves have not
been formulated for rotational paths.
15
2.2.2 Multilevel Algorithm for Rotations
As mentioned above, the rotational paths can be represented by trajectories on
the surface of a four-dimensional unit sphere. In imaginary-time, these paths
are closed, i.e. the beginning and end points of a path must coincide due to the
trace operation inherent in the partition function. The objective of the multilevel
algorithm is to nd ways to sample these closed paths eciently.
The most straightforward way to carry out Monte Carlo sampling is probably
the Metropolis method [Metropolis et al., 1953]. When applied to rotational
path sampling, the Metroplis algorithm would sample a new set of the three
Euler angles
i
,
i
and
i
on time slicei from the Jacobiand
=ddd cos and
accept or reject them based on the acceptance probability
P = min[1;(; R
new
i
; R
new
i+1
)=(; R
old
i
; R
old
i+1
)];
where is the short-time propagator in the primitive discretization from
Eqn.(2.3). In order to obtain converged results, the primitive discretization
requires a largen and hence a small. In this limit, the rotational propagator is
sharply peaked at small displacementsR = R
i+1
R
i+1
causing extremely nar-
row rotational dispersion in a single time slice. This in turn leads to slow sampling
problems when the Metropolis algorithm is used to simulate path integrals.
16
The slow sampling problem may be cured if one is able to generate large
sections of the rotational path spanning several time slices at a time. This is
because longer-time path segments are more delocalized in space, and sampling
them will reduce the correlation between the samples. But to do this, we will need
additional knowledge about the distribution in question. We will now outline a
procedure that can be used to generate long-time rotation paths eciently.
The objective is to sample long-time paths parameterized by the set of discrete
points (R
1
; R
2
; R
N1
) with xed endpoints R
0
, R
N
and imaginary time t =
N from the distribution
Q
N1
i=0
(; R
i
; R
i+1
). We take N = 2
M
, and refer to
M as the number of \levels" in calculation. One starts at the top level m =M,
where the path consists the entire time slice t with the two points R
0
and R
F
.
At the next lower levelm =M1, the path is bisected in time and the midpoint
R
N=2
is sampled from the conditional probability distribution with the ner levels
integrated out:
(R
N=2
;
t
2
jR
0
; R
F
) =
(R
0
; R
N=2
;t=2)(R
N=2
; R
F
;t=2)
(R
0
; R
F
;t)
(2.8)
where (R
0
; R
F
;t) provides the proper normalization. At level m =M 2 with
time slicet=4, we sample two more points, conditional to the end points and the
midpoint R
N=2
from the previous level. In other words, one applies the exact
same procedure as in the previous level (level M 1)to the end points R
0
, R
N=2
17
to sample R
N=4
, and R
N=2
, R
F
to get R
3N=4
. The procedure can then be repeated
recursively down to level m = 0 until the whole path is constructed. The only
requirement is that the exact analytical expression for the conditional distribution
Z
(R
i
; R;=2)(R; R
i+1
;=2)dR =(R
i
; R
i+1
;) (2.9)
is known for arbitrary time scale up to time t. For a free particle or a rigid
free rotor, this is no problem. The method of bisecting the free-particle path
with respect to the translational degrees of freedom is the basis of the bisection
algorithm described by D.M. Ceperley in [Ceperley, 1995] and the multilevel path
sampling algorithm of Mak et al. [Mak and Zakharov, 2004]. With respect to
rotations, the challenge comes from the fact that the conditional probabilities
Eqn.(2.8) do not have the free particles gaussian nature and the corresponding
cumulative probabilities
(;;;

i
;

i+1
;) =(

i
;
;)(
;

i+1
;) (2.10)
need to be integrated numerically.

i
and

i+1
are initial and nal reference
frames and enter the integration as parameters.
18
2.2.3 Details of Implementation
First we outline how to sample from a multivariate distribution by a series of one
dimensional integrations.
In the case of 1D distribution functionp(x)2 [a;b], in order to draw a sample,
one draws a random number R from a uniform distribution between 0 and A
and compares it against P (X), where P (X)
R
X
a
p(x)dx is dened to be the
comulative probability and A is the total area under p(x). Then one solves for
X from equation R = P (X), to obtain a sample from the distribution p(x).
Sampling in 2D and higher dimensions, can be reduced to the aforementioned
scheme. In the simplest multivariate case of 2D p(x;y) with x;y2 [a;b], one
can deneq(y) =
R
b
a
p(x;y)dx and corresponding comulative probabilityQ(Y ) =
R
Y
a
q(y)dy which allows to sample Y variable. Then X can be sampled from
P (X) =
R
X
a
p(xjY )dx. Generalization for higher dimensions can be induced
easily from the above procedure.
Diculty comes from the fact that when analytical expression is not easily
available, numerical integration is usually too computationally expensive for the
method to be of practical use and it happens in our case. Therefore, we have
implemented a rejection sampling scheme, described in appendix C.
In our case, the three dimensional distribution(;;;

i
;

i+1
;) is approx-
imated by its covering piecewise constant function f(;;;

i
;

i+1
;) and;
19
and can be put on an adaptive grid such that f(;;;

i
;

i+1
;) covers the
original function and the points are concentrated around its peaks. Extremum
analysis of(;;;

i
;

i+1
;) is discussed later on. The covering function, the
way it is written, depends on seven parameters: six from

i
and

i+1
and one
more from . This form does not allow to precompute corresponding comulative
probabilities in advance before the simulation due to the operational memory
constraint and is, thus, still not practical.
Physics of the system is not the reference frame specic and the distribution
can be precomputed for the nal positions

0
i+1
only, i.e.

i+1
in the frame of
initial positions

i
for each imaginary time level . This leaves one with only
four external parameters, which is still a lot.
Number of parameters can be reduced even further if one takes into account
the symmetry of the spherical top propagator Eqns. (2.7),(2.10):
(;;; 0;

0
;) =(e
0
e;)(e e
0
;) (2.11)
Where vector e is dened by the four Cayley-Klein parameters Eqn. (2.6). Here,
we dene initial frame, 0, to be e
0
(1; 0; 0; 0) with
0
= 0,
0
= 0,
0
= 0.
Corresponding 4D rotational matrix, R
0
, is discussed in appendix A. Moreover,
after the initial rotation to the origin, we rotate the system in 4D arounde
0
axis
such that the nal coordinate frame, e
0
, lies in e
0
e
2
plane with corresponding
20
rotation matrix R
e
0
e
2
such that e
0
= (cos

0
2
; 0; sin

0
2
; 0) (
0
= 0 or ;
0
= 0;

0
=
0
). Then, the nal coordinate frame can be parameterized by a single
parameter, cos
0
, where
0
2 [0;].
Thus, for each time level , a ne grid in cos
0
is created. For each time
level and cos
0
, extremum analysis of the propagator is performed (appendix
B) and the conditional cumulative probabilities f
c
(;; cos
0
), f
c
(;; cos
0
j)
andf
c
(;; cos
0
j;) are precomputed and stored in advance. Transformations
R
1
0
(

i
;

i+1
) and R
1
e
0
e
2
(

i
;

i+1
) are obtained from the endpoints during the
simulation, just before or after; and are sampled from the approximate rigid
spherical top distribution by the rejection method.
This allows one to generate long time free rotor path between two arbitrary
end points. After the path is constructed, it is nally accepted based on the total
potential in the Metropolis fashion.
2.3 Results
A number of PIMC simulations have been carried out in the past to investigate
the nature of microscopic super
uidity in
4
He droplets with impurities. The
earliest such studies dealt with impurities which are heavy rotors like SF
6
and
OCS. Most of them did not treat the rotations of the probe explicitly [Kwon
et al., 1996, Kwon and Whaley, 1999, Kwon et al., 2000]. In general, it was
21
found that these heavy and highly anisotropic molecules strongly perturb the
helium density around them, and Whaley et al. used the \microscopic two-
uid
model" to describe the super
uid response of the
4
He atoms [Kwon et al., 2000].
In this picture, the renormalization of the rotational constant was thought to be
a re
ection of the \adiabatic following" of the non-super
uid helium atoms to
the rotations of the probe molecule.
Even though the path integral formulation for the rotations of asymmetric top
molecules had been known for a while [Muser and Berne, 1996], the incorporation
of explicit rotations into PIMC simulations of doped helium droplets was not
carried out until very recently because of various computational issues. A number
of such PIMC simulations have appeared recently for linear molecules where all
the rotational degrees of freedom of the probe were taken into account [Paesani
and Whaley, 2004, Zillich et al., 2005, Blinov et al., 2004]. In these simulations,
the renormalized rotational constant of the probe in helium was obtained directly
from the orientational correlation of the probe. It was found that the emergence
of free-rotor-like spectral features for heavy molecules like OCS embedded in
4
He clusters did require exchange, but light molecules like acetylene showed very
similar spectral features to the gas phase with or without exchanges.
Clearly, the question of how the moment of inertia of the rotor aects the
exchange of the
4
He is not completely understood, nor is the counter eect of
how the exchange of the
4
He in turn in
uences the rotation of the molecule.
22
Also unclear is whether the symmetry of the probe makes any dierence. In the
current study, we have used methane (CH
4
) as the probe.
Normal CH
4
is a light spherical top molecule, with a gas-phase rotational
constantB = 5:24 cm
1
[Herzberg, 1991]. Experiments with normal and deuter-
ated CH
4
have been performed [Nauta and Miller, 2001a, Skvortsov et al., 2007,
Rudolph et al., 2007], both showing very weak renormalization inB. The advan-
tage of theory is that we can change the moment of inertia of the CH
4
molecule
while holding everything else xed. This allows us to isolate the eect of the
moment of inertia of the rotor and examine its in
uence on the rest of the prop-
erties of the system. In our simulations, we have measured the super
uid response
and the renormalization of the B constant for a variety of intrinsic moments of
inertia of the probe, from one time to fty times the normal magnitude of natural
CH
4
. In the rest of this paper, we will simply denote these `heavy" CH
4
by their
masses, e.g. 2CH
4
would correspond to deuterated methane. For each mass, we
have also investigated the nature of the bosonic eects by turning the exchanges
among the
4
He atoms on or o.
In each simulation reported here, a single CH
4
molecule was simulated with
a xed number of
4
He atoms around it. The interactions between the CH
4
and
heliums were described by the Buck potential [Buck et al., 1985] and the He-He
interactions by the Aziz potential [Aziz et al., 1987]. All simulations were done at
a constant temperature of 0.3 K with the imaginary time discretization of 1024
23
beads for both translational and rotational degrees of freedom. Approximately
twelve He atoms are sucient to complete the rst solvation shell; consequently,
most of the detailed calculations have been performed for this cluster size. For
twelve He and at 0.3K, the total
4
He super
uid density is close to unity. This
allows us to isolate the eects of just the rst solvation shell from the bulk
4
He.
According to the experiments, dynamical measurables like the ro-vibrational
spectra of the probe embedded in
4
He clusters show signs of quantum statistics.
Extracting such dynamical properties from the simulations would entail either a
real-time path integral simulations or an analytic continuation of imaginary-time
correlation function , both of which are notoriously dicult. In the recent works
of [Paesani et al., 2003] and [Moroni et al., 2003], maximum entropy analysis and
multi-exponential ts were applied to the orientation correlation function in order
to extract the spectral features. Recently, an alternative method was suggested
by several groups to estimate the renormalized B constant using imaginary-time
data alone, based on the orientational correlation function of the probe [Paesani
and Whaley, 2004, Zillich et al., 2005, Blinov et al., 2004]. [Blinov et al., 2004]
suggested using the imaginary time correlation function
h^ n() ^ n(0)i =
1
Z
Trfe

^
H
e

^
H
^ ne

^
H
^ ng; (2.12)
24
where ^ n is a unit vector describing the orientation of the molecule relative to
the laboratory frame and Z =h^ n(0) ^ n(0)i. They tted this to the analytical
expression for a rigid top molecule to extract an eective B constant, allowing
them to estimate the renormalized rotational constant of the OCS molecule based
on the Levenberg-Marquardt algorithm [Press et al., 1995]. We have adopted the
same approach in our work. For our case, the probe and the average density
has overall spherical symmetry. It is therefore logical to t the simulation results
to a rigid spherical top orientation correlation function, for which the following
analytical expression can be derived straightforwardly
h^ n() ^ n(0)i
CH
4
=
1
Z
(
e
2B
+
X
j>0

X
kj1
j
2
k
2
j
e
Bj(j+1)+2Bj
+
X
kj
e
Bj(j+1)

(2j + 1)k
2
j(j + 1)
+
(j + 1)
2
k
2
j + 1
e
2B(j+1)

!)
; (2.13)
where B is the eective rotational constant of a spherical top molecule. An
example of the orientation correlation functions is showed in Fig. 2.1 for CH
4
(B
= 5.24 cm
1
), 2CH
4
, 50CH
4
in the gas phase, as well as a 50CH
4
surrounded
by twelve bosonic
4
He as a function of the imaginary-time =.
To extract the eective B constant for 50CH
4
embedded in twelve
4
He, we
tted its correlation function to that of a free spherical top molecule. The t is
shown in Fig. 2.2. The statistical error of the orientation correlation function is
comparable to the width of the line in the plot. The tted value corresponds to
25
0 0.2 0.4 0.6 0.8 1
τ/β
0
0.2
0.4
0.6
0.8
1
C(τ)
CH
4
2xCH
4
50xCH
4
50xCH
4
in 12
4
He
Figure 2.1: Examples of the orientation correlation function for normal CH
4
,
2CH
4
and 50CH
4
the gas phase, as well as 50CH
4
inside a shell of twelve
4
He atoms.
an eective B constant of 0.091 cm
1
, compared to 0.105 cm
1
for a 50CH
4
molecule in the gas phase.
To assess the eects of quantum statistics, we have carried out the same
simulation turning o all the exchanges among the He atoms. The correlation
function for this non-exchange case is shown in Fig. 2.2 as the dotted line. The
non-exchange correlation function lies above the exchange result, indicating that
exchange facilitates a stronger rotational coherence in the 50CH
4
rotation. The
26
0.2 0.4 0.6 0.8 1
τ/β
0.8
0.85
0.9
0.95
1
C(τ)
50xCH
4
50xCH
4
, gas phase
50xCH
4
, no exchange
fit to effective B
Figure 2.2: Orientation correlation function of 50CH
4
in a shell of twelve
4
He
atoms (solid line) and the best t to a rigid spherical top model (dashed line),
compared to the non-exchange result (dotted line) and the gas-phase (dotted
dashed line).
eective B constant in non-exchange helium is 0.086 cm
1
. Relative to the gas-
phase value, the renormalizedB constant for 50CH
4
is 15.4% and 21.5% smaller
in He, for the exchange and non-exchange cases, respectively. These correspond
to an eective increase in the moment of inertia of the probe from 160.85 amu

A
2
to 185.65 amu

A
2
and 190.35 amu

A
2
for the exchange and non-exchange cases,
respectively. Overall, these results are qualitatively very similar to what was
27
Figure 2.3: Angular density prole of the
4
He shell in a (a) 50CH
4
-(He)
12
and
(b) 1CH
4
-(He)
12
cluster. H show the positions of the hydrogen atoms. The
4
He
density is strongly templated by the heavy probe, while it is almost completely
isotropic for the light probe.
observed for linear rotors previously, showing that bosonic eects promote rota-
tional coherence of a heavy rotor, though the magnitude of theB constant renor-
malization is much smaller in this case (15.4% here versus 177% in OCS [Grebenev
et al., 2000]).
Next, we examine the eects of the probe on the He atoms in the shell.
Figure 2.3(a) shows the angular density of the He atoms in a 50CH
4
-(He)
12
cluster on the surface of the probe molecule, integrated over all distances in
the radial direction. The picture reveals four regions of low densities, with one
of the four hydrogen atoms of the methane positioned at the center of each.
The rest of the helium density is characterized by interconnecting high-density
channels on the surface of the CH
4
. Fig. 2.3(a) shows that there are 10 density
28
maxima inside the channels and 4 weaker local maxima centered on the hydrogen
atoms. The helium atoms sample the potential surface and self-assemble into
a close-packed-like structure, indicating strong templating eects by the probe
molecule underneath. The helium density variation from minimum to maximum
in Fig. 2.3(a) is approximately a factor of two for 50CH
4
-(He)
12
. We also studied
the shell density in the non-exchange case where all exchanges were turned o
and found no discernable dierences in the helium density between the exchange
and non-exchange cases.
Previous studies on SF
6
in
4
He droplets found similar strong templating eects
of the probe on the helium atoms in the rst solvation shell [Kwon et al., 2000].
Based on those results, a microscopic two-
uid model was proposed to rationalize
the super
uid response of the helium atoms. It would therefore be interesting to
examine whether this strong templating eect is re
ected in the super
uid density
of the shell in the case of CH
4
too. For 50CH
4
-(He)
12
, the super
uid fraction
of the shell turns out to be surprisingly large. The simulations show a 961%
super
uid fraction, despite the heavy templating by the probe. Even though there
is a distinct closed-packed-like structure in the helium density, the exchanges can
apparently occur quite eciently along the inconnected high-density channels on
the surface of the CH
4
molecule, encompassing almost the entire shell. Thus, the
ability of the helium atoms to undergo exchange appears to be intimately related
to the symmetry of the underlying probe molecule.
29
If we assume that the non-super
uid fraction of the helium shell follows the
rotation of the probe adiabatically, one can estimate the resulting change in
the moment-of-inertia of the probe. Based on a 4% non-super
uid fraction and
an average separation of approximately 4.05

A between the helium atoms and the
probe, the eective moment-of-inertia of the probe is expected to increase by 13%,
assuming the non-super
uid heliums are uniformly distributed on the shell. On
rst sight, this estimate appears to be in good agreement with the renormalization
in the observedB constant in Fig. 2.2 for the exchange case, which is about 4%..
But if one applies the same logic to the non-exchange case where 100% of the shell
is non-super
uid, then one would have expected a 2500% renormalization in the
B constant. Figure 2.2 shows only a 20% renormalization for the non-exchange
case; therefore, while the renormalized B constant appears to be consistent with
the amount of non-super
uid helium atoms in the exchange case, this agreement
must be taken to be fortuitous and the concept of adiabatic following in the case
of CH
4
does not seem to apply, where the anisotropic interaction between the
methane and the shell is weak.
The next question is how much of a complete helium shell is necessary for
the eects of quantum statistics to show up in the probe rotations. Partial shell
results for the renormalizedB constants are shown in Fig. 2.4(a) for one to twenty
4
He atoms on a 50CH
4
for the full exchange case. The renormalized B constant
seems to have reached an asymptotic limit at around six
4
He atoms. Beyond
30
0 5 10 15 20
Number of
4
He
5.0
5.1
5.2
B, cm
-1
0.09
0.10
B, cm
-1
(b) 1xCH
4
(a) 50xCH
4
Figure 2.4: EectiveB constant for (a) 50CH
4
and (b) 1CH
4
in clusters with
dierent number of
4
He atoms, from 1 to 20. For 50CH
4
, the renormalization
of theB constant by the
4
He seems to reach an asymptotic limit at around 6
4
He
atoms, with only minor variations thereafter. For 1CH
4
, the change in B with
the number of
4
He atoms is much more gradual.
31
Figure 2.5: Eective B constant for SF
6
in clusters with dierent number of
4
He atoms, from 1 to 8 computed by several methods and compared with
experimental nanodroplet value [Hartmann et al., 1995], [Harms et al., 1997]
B
exp
= 0:034 cm
1
. The renormalization of the B constant by the
4
He seems to
reach an asymptotic limit at around 8
4
He atoms.
1
[Kwon et al., 2000]
that, the eective B constants show only minor variations up to 20
4
He atoms.
Clearly, the dominant eects of the
4
He shell on the rotations of the 50CH
4
are
already present with a small number of
4
He atoms, indicating that the observed
eects do not require even a complete solvation shell.
SF
6
molecule has been studied extensively by both experimentalists [Hart-
mann et al., 1995], [Hartmann et al., 1996], [Grebenev et al., 1998], [Harms et al.,
1997] and theoreticians [Kwon et al., 1996], [Kwon and Whaley, 1999], [Kwon
32
et al., 2000]. Path Integral studies were performed by [Kwon et al., 1996] long
time ago. Since it's moment of inertia is very large, in the classical limit,the
rotational propagator was not considered. Later on its rotational constant was
estimated from the super
uid response assuming that
4
He shell adiabatically fol-
lows this heavy rotor [Kwon et al., 2000]. Zero temperature DMC calculations of
SF
6
were previously reported, rotational probe wave functions were introduced,
and the rotational constants were obtained from rotational spectrum [Kwon et al.,
2000]. For the sake of comparison we have performed nite temperature partial
4
He shell calculations taking into account rotational degrees of freedom of SF
6
.
This kind of simulation has not been done before. Rotational constants, as in
the case of methane, are obtained from the orientation correlation function. The
interactions between the SF
6
and the helium are described by the Pack poten-
tial [Pack et al., 1984] and the He-He interactions,as before, by the Aziz potential
[Aziz et al., 1987]. All SF
6
simulations were performed at a constant temperature
of 0.3 K with the imaginary time discretization of 512 beads for both translational
and rotational degrees of freedom. At these sizes of clusters, super
uid fraction
is close to zero, so the results are expected to comply with the previous works.
Fig. 2.5. The curve is expected to saturate at the nanodroplet limit and its shape
resembles one obtained in the case of heavy methane. But, due to its very large
anisotropy, SF
6
with up to 8 He atoms behaves like a rigid top, i.e. a spherical
shell at the equilibrium distance from the molecule, rigidly attached to it. Heavy
33
0 30 60 90 120 150 180
I
o
, amu A
2
0.90
0.92
0.94
0.96
0.98
1.00
1.02
f
S
01 3 6 12 25 50
xCH
4
Figure 2.6: The super
uid fraction of the helium shell for 50 to 1CH
4
with 12
4
He atoms. Due to the symmetry of the underlying probe, the super
uid fraction
of the shell is close to unity.
methane is not even close to be rigid. Due to its much weaker interaction with the
shell, super
uid fraction is about 0.96. Much stronger potential of SF
6
localizes
4
He atoms resulting in near zero super
uid fraction. To what extend this
4
He
shell slipping around the methane is Bose statistics eect still remains unclear.
We now turn to the results for the lighter probes. Figure 2.6 shows how the
super
uid fraction varies with the intrinsic moment of inertia of the probe, from
34
50 to 1CH
4
. Starting at 96% for 50CH
4
, the super
uid fraction increases
quickly to 100% when the intrinsic moment of intertia of CH
4
is decreased. The
reason for this increase in the super
uid fraction for lighter probes is immediately
obvious when we examine the helium atom density on the shell for the lighter
CH
4
, an example of which which is shown in Fig. 2.3(b) for 1CH
4
. In this
case, the helium angular density is almost completely uniform on the surface of
the 1CH
4
, with a minimum-to-maximum variation of only 7%. Clearly, the
light probe provides only very weak templating for the helium density, and for
all practical purposes the helium atoms seem to be distributed uniformly around
the methane. This of course does not imply that there is no structure in the
helium shell itself. Fig. 2.3(b) only shows the one-body density. Any order in
the arrangement of the helium atoms on the shell can only be discerned if we
measure the two-body angular density (similar to the radial distribution function
in bulk liquids). But Fig. 2.3(b) does conrm that there is almost no correlation
between the orientation of the underlying probe molecule and the orientation of
the shell for light CH
4
. This extremely weak templating by the probe enables
rapid exchanges among the
4
He atoms on the shell, leading to a super
uid fraction
close to unity.
While the single-particle density in Fig. 2.3 reveals no noticeable dierence
between the exchange and the non-exchange cases, the two-particle density does.
35
g(θ)
1xCH
4
no exchange
0 π/4 π/2 3/2π
π
θ
g(θ)
50xCH
4
no exchange
(a)
(b)
Figure 2.7: The two-particle angular distribution function for (a) 1CH
4
and (b)
50CH
4
with full exchange (solid lines) and no exchange (dotted lines).
Figure 2.7 (a) shows the two-particle angular density for 1CH
4
for the non-
exchange case (dotted line), compared to the full-exchange results (solid line).
36
There are minor but noticeable dierences between the two, with the peaks in
the non-exchange case being sharper in general. Apparently quantum exchange
induces more delocalization in the
4
He atoms, causing them to tunneling into the
interparticle spacing more eectively. Also clear in the two-particle density is
the solid-like structure in the helium shell, revealed by the distinct peaks in the
angular distribution function. Similar behaviors are seen in Fig. 2.7 (b) for the
50CH
4
.
0 0.2 0.4 0.6 0.8 1
τ/β
0
0.2
0.4
0.6
0.8
C(τ)
1xCH
4
1xCH
4
, gas phase
Figure 2.8: Orientation correlation function of 1CH
4
in a shell of twelve
4
He
atoms (solid line)compared to the gas-phase (dotted line).
37
Figure 2.8 shows the orientation correlation function for the 1CH
4
in twelve
4
He atoms with full exchange, compared to the free CH
4
molecule in the gas phase.
The t to a free rotor orientation correlation function yields a renormalized B
constant of about 5.10 cm
1
. This corresponds to an eective increase in the
moment of inertia of the probe from 3.217 amu

A
2
in the gas phase to 3.30
amu

A
2
in helium, a mere 3% change. This result is unaltered whether exchages
are turned on or o. This suggests that the lack of rotational slowing in the case
of a light CH
4
bears no direct relationship with the super
uidity of the shell {
the molecule has the same eectiveB constant with or without exchanges in the
helium. In addition, both cases show almost identical values to the gas-phase
B constant, indicating that the rotations of the molecule are almost completely
decoupled from the shell.
The renormalizedB constant for CH
4
molecules of dierent intrinsic moments
of inertia are summarized in Fig. 2.9. The renormalized rotational constant of
normal CH
4
is 5.10 cm
1
and for CD
4
is 2.48 cm
1
. These are in good agreement
with the experimental nanodroplet values of 5.0 cm
1
and 2.5 cm
1
[Skvortsov
et al., 2007]. To put these small changes on a more convenient scale to display the
eects of the mass of the probe, we have replotted the percentage increase in the
renormalized B constant relative to their gas-phase moment of inertia's values
in Fig. 2.9. Open circles correspond to results with full exchange and squares
to their non-exchange counterparts. The eects of quantum statistics seem to
38
0 30 60 90 120 150 180
I
o
, amu A
2
0
5
10
15
20
25
(B
o
-B)/B
o
, %
01 3 6 25 50 12
xCH
4
Figure 2.9: The percentage change in the eective B constant as a function of
the intrinsic moment of inertia of the probe with full exchange (circles) and no
exchange (squares).
be signicant only for heavy probes, with only very minor dierences between
the exchange and non-exchange results for a probe mass 25 times that of normal
methane or less.
To summarize the key results, the rotations of a light probe are largely decou-
pled from the helium shell, leading to almost no noticeable dierence between
the measured rotational constant and its intrinsic gas-phase value, regardless of
39
whether the helium atoms are bosonic or not. In contrast, the rotations of a heavy
probes are more sensitive to the quantum statistics in the helium shell, resulting
in a lighter eective moment of inertia when bosonic exchanges are operative.
While this is true, the super
uidity of the shell seems to have no quantiable
relationship with the relative reduction in the B constant, with the super
uid
fraction being close to unity for light or heavy probes. These results suggest
that the major determinant of how theB constant behaves in the case where the
anisotropic interactions between the probe and the shell is weak is the intrinsic
moment of inertia of the probe but not the super
uid fraction of the shell.
2.4 Discussion
The fact that the intrinsic moment of inertia of the probe has such profound
eects on the observed behaviors can be qualitatively understood with a simple
picture. Classically speaking, a light rotor rotates rapidly. From the perspective
of the helium shell, a light molecule would therefore appear to be less anisotropic
than a slowly rotating molecule. This is consistent with the observation that
a heavy rotor produces strong templating of the helium density on the shell,
but a light rotor produces almost none. We will refer to this as the \rotational
smearing" eect.
40
The simple classical rotational smearing analog given above is of course
incomplete, because for very light rotors like normal CH
4
none of the rotational
excited states are signicantly occupied at the typical experimental temperature
of around 0.3 K. A normal CH
4
molecule would therefore be in its ground rota-
tional state most of the time. But since the ground rotational state of CH
4
is
completely symmetric, the shell sees an eectively isotropic molecule and the
helium density shows almost no templating by the probe. While the classical
rotational smearing analog given above is incomplete, its implications are largely
correct in spirit. Furthermore, this rotational smearing eect exists independent
of whether the heliums on the shell are bosonic or not and is therefore unrelated
to quantum statistics.
While the in
uence of the probe on the helium shell can be largely explained
by rotational smearing, the counter eects of the shell on the probe's rotational
constant are more dicult to understand, because they involve not just the
ground but also the excited rotational states of the molecule. Apparently, a
secondary consequence of rotational smearing is a lowering of the dynamical cou-
pling between a light probe and its shell, causing the shell to slip increasingly
relative to the rotational motions of the molecule when the intrinsic moment of
inertia of the probe is lowered. Of course, the origin of the coupling between
the probe and the shell is the anisotropy of the probe-helium potential, whose
magnitude is obviously unchanged regardless of the intrinsic moment of inertia
41
of the probe. Since theB constant shows large renormalization for heavy probes
while almost none for light one, the degree of the shell's slippage apparently
does depend on the intrinsic moment of inertia and there must be a dierence
between light and heavy probes in their eective dynamical coupling to the shell.
In the following, we will try to provide a plasible explanation for this dierence
by appealing to a simple two-dimensional model.
Figure 2.10: A planar model describing the interaction between the rotating probe
and a rigid shell. The left panel shows a schematic picture of the the spectrum
of the probe rotations. The right panel shows the low-energy spectrum of the
shell and how the density of states of the shell is dilated when bose symmetry is
imposed on the shell wavefunctions.
In a planar model rst proposed by Lehmann [Lehmann, 2001], the probe is
taken to be a linear rotor embedded inside a rigid ring ofN
4
He atoms uniformly
distributed around the ring. Figure 2.10 shows a cartoon of the model with six
helium atoms. We assume the only motions for both the probe and the ring are
rotations about the center-of-mass of the system, and the rotations of the probe
42
interact with the ring's via an anisotropic potential, causing hindered motions of
the probe relative to the shell.
The hamiltonian of the model is given byH =H
1
+H
2
+V , whereH
1
=L
2
1
=2I
1
andH
2
=L
2
2
=2I
2
are the hamiltonians of the probe and the ring, respectively, and
they are coupled via an interaction potentialV =V
0
cos(N), where =
1

2
.
Here
1;2
are the rotation angles, I
1;2
are the moments of inertia, and L
1;2
are
angular momentum operators, of the molecule and the ring, respectively. The
rotational constants are related to the moments of inertia by B
1;2
= h
2
=2I
1;2
.
The individual spectra of the probe and the shell are represented pictorially in
Fig. 2.10. The eigenstates of the uncoupled probe and the ring are labeled by
their angular momentum quantum numbers J
1;2
= 0;1;2;::: with eigenener-
gies B
1
jJ
1
j, and B
2
jJ
2
j. These are used to form the basis for the coupled system
fjJ
1
;J
2
ig. SinceV conserves total angular momentum, the interaction only cou-
ples states with the same total angular momentum J
tot
=J
1
+J
2
. The hamilto-
nian in thejJ
1
;J
2
i basis is therefore block diagonal, with each block belonging
to a dierent J
tot
.
We rst consider the consequence of quantum exchange on the spectrum of
the ring. Bosonic symmetry requires a many-particle
4
He wavefunction to be
totally symmetric upon exchange. The only states for which this is true in a
N-helium ring are those where J
2
is an integer multiple of N. This requirement
removes states from the spectrum of the ring as illustrated in Fig. 2.10, leading
43
to a widening of the energy gaps among the allowed states. The resulting dilation
in the density of states leads to a smaller eective coupling between the ring and
the probe. This is the primary reason for the suppression of momentum transfer
between them when bosonic exchanges are in eect.
We now tackle the question of why lighter probes tend to show smaller renor-
malization in their eective B constants. The measured B constant is related
to the energy gap between the lowest J
tot
= 0 state and the lowest J
tot
=1
state of the total system, and these states come from dierent subblocks of the
hamiltonian. We will call the lowest state in the J
tot
= 0 subblock
0
, and

0
is a mixture ofjJ
1
;J
2
i =j0; 0i;j + 6;6i;j 6; +6i, etc.
1
will be used
to denote the lowest state in the J
tot
= +1 subblock, and
1
is a mixture of
j + 1; 0i;j + 7;6i;j 5; +6i, etc.
1
denotes the lowest state in the J
tot
=1
subblock, and it is given by a similar mixture. Since the total hamiltonian has
no preference for the overall direction of the rotation,
1
is necessarily degener-
ate with
1
, The eective B constant of the system is then given by the energy
dierence between
1
(or
1
) and
0
.
For a very light probe with a large B
1
, the energy gap betweenj0; 0i and
j + 6;6i orj 6; +6i is large; therefore,
0
is almost purelyj0; 0i. Similarly,

1
inherits almost all of its character fromj1; 0i when the probe is very light. In
this limit, the eective B constant would of course be just equal to B
1
.
44
When the intrinsic moment of inertia of the probe is increased, the smaller
B
1
causes the higher-energy states in each subblock to mix into
0
and
1
.
0
takes on increasing characters of both thej+6;6i andj6; +6i states, and this
mixing results in a lowering of the energy of
0
. For
1
, a similar thing happens
and the mixing ofj + 7;6i andj 5; +6i intoj1; 0i causes the
1
energy to also
decrease. But the amount of lowering experienced by
0
and
1
are not identical.
Since the energy dierence betweenj 5; +6i andj1; 0i is smaller than between
j+6;6i andj0; 0i (4B
1
+6B
2
versus 6B
1
+6B
2
), the lowering of the energy of
1
is larger than it is for
0
. This produces a smaller energy gap between
0
and
1
compared to the uncoupled system, and as a result, a smaller eectiveB constant.
This largely explains why theB constant is diminished when the probe is coupled
to the ring through the anisotropy and why the observed renormalization in B
depends on the intrinsic moment of inertia of the probe.
Notice that in the picture above, the quantum statistics of the
4
He manifests
itself only through the spectrum of the shell. In general, the symmetry require-
ment for the many-particle bose wave functions removes states from the shell's
spectrum. This dilatation of the density of states in the shell's spectrum then
results in a lowering of the eective coupling between the probe and the shell,
causing the probe to slip increasingly with respect to the shell. Beyond this mod-
ication, no explicit reference to bosonic eects needs to be invoked in order to
explain the observed renormalization in the B constant.
45
As a nal illustration, we computed spectrum of the model in the basis of wave
functions of noninteracting probe-ring system. We set V = 15cm
1
, B
ring
=
0:078cm
1
, and varied moment of inertia of the probe from 1 methane till 25
times moment of inertia of methane ( B
probe
from 7.55 cm
1
till 0.302 cm
1
).
Then, we calculated the orientation correlation function Eqn.(2.12) and obtained
eective rotational constants from the t to the orientation correlation function
of the free 2D rotor:
h^ n() ^ n(0)i
2D
=
1
Z
j=N
b
X
j=N
b
e
Bj
2
()B(j1)
2
; (2.14)
where N
b
is number of basis functions. The renormalized B constants for the
probes of dierent intrinsic moments of inertia are summarized in Fig. 2.11 with
and without of inclusion of the permutational symmetry of the ring . We imme-
diately see a qualitative agreement with Fig. 2.9.
While the planar model seems to be able to explain the qualitative origin of
the counter eects of the shell on the probe's rotations, its implications must
not be taken literally. The reason is that in the planar model, the elementary
excitations of the shell are assumed to be rigid rotations, and this allows us to
rationalize the renormalization of the B constant in terms of the changes in the
spectrum of the shell and its eects on the dynamical coupling with the probe. In
three dimensions, the precise nature of the elementary excitations of the shell is
46
0 30 60 90
I
0
, amu*A
2
0
40
80
120
160
200
240
%(B
0
-B)/B
1 2 3 6 25 12
xCH
4
(2D)
Figure 2.11: 2D Toy Model. The percentage change in the eective B constant,
obtained from nonlinear t with 2D free rotor orientation correlation function,
as a function of the intrinsic moment of inertia of the probe with full exchange
(circles) and no exchange (squares).
not at all clear. We are uncertain whether rigid rotations of the shell will remain
to be the dominant excitations responsible for the low-energy spectrum of the
shell. At present, the precise eects of quantum exchange on the shell's spectrum
in the three-dimensional case is still unclear. We plan to address this question
more fully in future path integral studies.
47
2.5 Conclusions and Future Work
We have presented a novel PIMC technique for symmetric top molecules. The
rotational and translational degrees of freedom are sampled simultaneously, allow-
ing for a better understanding of rotational dynamics in
4
He clusters of probe
molecules with this symmetry. To the best of our knowledge, it is the rst time
that multilevel, nonlocal moves have been applied to the the rotational paths,
which is rather relevant to spherical tops with three rotational degrees of free-
dom. In addition this method can be applied to linear rotors without major
modications. However, with respect to symmetric and asymmetric tops, addi-
tional work is required. The covering distribution from the rejection part of the
method, discussed in section 2.2.3, needs a maxima prole of the original dis-
tribution, which is not clear in more general cases. Memory requirements are
also expected to be considerably higher due to the lower symmetry as well as the
rotational propagator is more expensive to compute. Extension of our approach
would allow us to investigate such probes as NH
3
and H
2
O for which experimental
data is available.
We have applied this technique to methane molecule inside
4
He
12
. The angu-
lar density prole of helium around methane is nearly
at as a consequence of its
very delocalized rotational wave function, whereas \heavy\ methane produced
48
close packed structure of
4
He around the probe molecule. In both cases the den-
sity prole was not sensitive to exchange eects. We have measured rotational
constants of the \articial\ methane molecules in the cluster while increasing
masses of its hydrogens. Results show that taking into account the Bose sym-
metry of
4
He does not produce any visible eect one the eective moment of
inertia of methane, CD
4
, and CT
4
, but does show up in articial C(
50
H)
4
. We
also illustrated this eect using a simple 2D model [Lehmann, 2001]. Though the
model captures some of the main features of the system, it needs considerable
improvement to be realistic.
We have neglected the eect of nuclear symmetry of methane. This molecule
has tetrahedral symmetry, and consequently, the rotational wave function can
be decomposed into states belonging to the identity representation A, the rep-
resentation T , and the representation E, where E can be further decomposed
into the two complex conjugate representationsE
1
andE
2
[Tinkham, 1964]. The
nuclear spin I in the case of CH
4
is I = 2 for A-CH
4
, I = 1 for T -CH
4
, and I
= 0 for E-CH
4
[Huller and Press, 1981]. In order to consider nuclear exchange
eects, the full kernel Eqn.(2.4) has to be symmetrized accordingly. It should
be noted that only the identity and full density representations correspond to
positively valued density matrices. All other representations have negative con-
tributions [Muser and Berne, 1996], which lead to the notorious "sign problem"
in path-integral simulations [Loh et al., 1990]. The methods to treat the sign
49
problem must be then applied [Marx, 1994], [Loh et al., 1990], [Mak and Chan-
dler, 1990], [Mak et al., 1998]. Which spin state is relevant might depend on
the cooling regime and the temperature. Good agreement with experiment sug-
gests that nonsymmetrized propagator is good enough for the properties we have
considered. Smooth character of Fig. 2.9, which shows how eective rotational
constant changes with the mass increase of methane's hydrogens also suggests
that the nuclear statistic eects do not seem to matter.
We have investigated a weakly interacting methane probe while increasing
its intrinsic moment of inertia. The next logical step will be to consider SF
6
molecule. It is a very \heavy\ rotor (approximately moment of inertia of C(
50
H)
4
)
with much higher anisotropy. Articially decreasing its moment of inertia could
provide additional insight to eect of anisotropy on the onset of super
uidity.
We have faced problems with ecient sampling of a SF
6
probe in larger helium
clusters (more than eight helium atoms). As mentioned before, our current sam-
pling scheme is based on the free particle rotational propagator, which is not
a good choice when the interaction with the shell is so strong. We have tried
a Multilevel Monte Carlo scheme for rotations, described by [Ceperley, 1995],
for translational moves. We used primitive density approximation for the level
density which did not seem to improve the convergence. It was pointed out in
[Ceperley, 1995] for the case of translational degrees of freedom, that primitive
approximation is not good enough for this cause. Unfortunately, two particle
50
density matrix approach, developed for central eld interactions, is not feasible
in our case because of a larger amount of degrees of freedom of the anisotropic
interaction resulting in computer memory constraints. Further development of
better density matrix approximations incorporating rotational degrees of freedom
can prove to be useful.
51
Chapter 3
A generalized kinematic
formulation of loop closures in
nucleic acids, proteins and other
linear polymers
3.1 Introduction
The loop closure problem was rst considered by Go and Scheraga in 1970 for
polypeptide chains [Go et al., 1970]. For linear polymers like proteins and nucleic
acids, the bond lengths and bond angles are largely xed and the
exibility of
these molecules are primarily derived from chain torsions. For very long chains,
motions in even a single torsion angle can aect the coordinates of many atoms
simultaneously and may lead to excessive steric overlaps. In order for torsion
angle moves in the conformational sampling of linear polymers to be practical,
multiple torsion angles must be moved at the same time so that the atomic
52
positions of the majority of the chain can remain relatively unperturbed. The loop
closure problem seeks solutions for the possible sets of torsion angles inside a chain
segment that have to be moved simultaneously in order to maintain the atomic
coordinates of the rest of the chain xed. Various formulations and extensions
of the loop closure problem have been reported since Go and Scheraga's work
[Wakana et al., 1984, Dodd et al., 1993, Knapp, 1992, Knapp and Irgensdefregger,
1993, Coutsias et al., 2004]. When applied to Monte Carlo simulations of polymer
conformational sampling, these loop closure solutions give rise to the concerted
rotation or 'conrot' algorithm and its variants [Bruccoleri and Karplus, 1985,
Sartori et al., 1998, Favrin et al., 2001, Bashford and Case, 2000, Ulmschneider
and Jorgensen, 2003, 2004b,a, Dinner, 2000, Deem and Bader, 1996].
In the original formulation of Go and Scheraga, they considered the case
of a protein where the rotatable (; ) angles are separated by xed peptide
bonds in the trans conformation. They concluded that at least six torsion angles
belonging to a tripeptide sequence must be moved simultaneously in order for the
atomic coordinates of the rest of the chain to remain xed. Later, Dinner [Dinner,
2000] extended this formulation to allow arbitrary values for the xed intervening
torsion angles between the rotatable ones. This formulation and extensions of
it also allow for applications to nucleic acids [Dinner, 2000, Ulmschneider and
Jorgensen, 2004a]. Dierent loop closure solutions were also reported recently
53
for proteins by Coutsias et al. [Coutsias et al., 2004] and for single nucleotides
and their ribose rings by Mak [Mak, 2007].
In general, the loop closure problem seeks solutions for the set of torsion
angles inside a segment along a linear chain polymer with the condition that the
atomic coordinates for the rest of the chain, including the two boundary atoms,
remain unchanged. To satisfy the bond length and bond angle constraints at the
two atoms on the boundaries, Go and Scheraga showed that a minimum of six
degrees of freedom are needed in the movable segment [Go et al., 1970]. For the
loop closure problem to be well-posed, at least six torsion angles must therefore
be moved at the same time. The corresponding 6-variable/6-constraint problem
has been solved in dierent ways under various formulations of the problem [Go
et al., 1970, Dodd et al., 1993, Deem and Bader, 1996, Wedemeyer and Scheraga,
1999, Coutsias et al., 2004].
In this work, we consider a generalized formulation of the loop closure problem
and provide a solution for it. This formulation applies to loop structures of
arbitrary sizes, from as small as a ve-member ring to those containing multiple
segments of arbitrary lengths with arbitrary xed intervening torsion angles inside
each segment. This generalized formulation reduces the original 6-variable/6-
constraint problem to a 4-variable/4-constraint problem, and a simple solution is
derived using geometrics of rigid bodies well known in the eld of the kinematics
54
of mechanisms [Hunt, 1978, Duy, 1980]. The generalized formulation allows for
an ecient numerical or analytical solution for loops of arbitrary lengths.
3.2 Boundary Constraints in the Loop Closure
Problem
In the loop closure problem, we consider the chain to be rigid except for chain
torsions, i.e. all bond angles and bond lengths are xed. The problem is to
nd a way to move a set of torsion angles while holding the rest of the chain
xed. Since the motions are in torsion-angle space but the constraints are in
coordinate space, the loop closure problem requires a transformation between
the torsion-angle representation of the chain conformation and its coordinate
space representation. In this section, we will examine the relationship between
the torsion angles inside the movable segment and the boundary constraints.
In Fig. 3.1(a), we illustrate the basic geometry, using a grey dashed line to
represent the movable segment. The movable segment is connected to the rest
of the chain on the left by atom (i + 1) and to the rest of the chain on the right
by atom (j 1). The positions of atoms (i + 1), (i), (i 1) and every atom to
the left indicated by the black dashed line are taken to be xed in coordinate
space. We call this the left boundary segment. On the other hand, we take the
right side of the chain containing atoms (j 1), (j), (j + 1) and beyond to be
55
rigid, i.e. the relative positions of all the atoms are xed, but the entire rigid right
segment is not xed in coordinate space. We call this the right boundary segment.
The precise position of the rigid right segment is of course determine together
by the xed left segment and the conformation of the movable grey segment.
Figure 3.1(b) shows how the position of the rigid right boundary segment may
change if the conformation of the grey segment is moved. Since the left and right
boundary segments are rigid, their relative positions are specied by six degrees
of freedom: three translations and and three rotations. This can be seen easily
as the variables needed to specify the position of the second triangle formed by
atoms (j 1), (j) and (j + 1) relative to the xed position of the rst triangle
formed by atoms (i 1), (i) and (i + 1). These six variables can be taken to be
the translation vector from atom (i) to (j), plus the three Euler angles specifying
the orientation of the second triangle relative to the rst.
Since six variables are required to specify the position and orientation of the
rigid right boundary segment relative to the left, the movable segment connecting
them must have at least six degrees of freedom; otherwise, the closure problem
will be ill-posed. In the simplest case where every torsion angle along the chain
is movable, the situation is illustrated in Fig. 3.1(c). We label the atoms consec-
utively starting with the black atom on the left as (1). The movable segment
in Fig. 3.1(a) is now drawn showing the atoms explicitly. Starting with the xed
positions of atoms (1), (0) and (1) on the left segment, specifying the torsion
56
Figure 3.1: (a) Drawing showing boundary conditions for the closure problem.
Atoms (i 1), (i) and (i + 1) on the left are xed in space. The black dashed line
indicates that the rest of the chain to the left of atom (i 1) are also xed. The
grey line represents the movable segment whose conformation is
exible. The
atoms (j 1), (j), (j + 1) and all the atoms to the right represented by the black
dashed line form a rigid segment. (b) When the conformation of the movable
segment is changed, the position and orientation of the right rigid segment will
change, while the left boundary segment remains xed in space. (c) The minimal
number of torsion angles needed in the movable segment is six. The picture shows
the minimal length if every bond along the movable segment is torsionable.
angle
1
xes the position of atom (2). In the same way, specifying the torsion
angles
i
, i = 2 to 6, denes the positions of all atoms up to (7). Therefore, as
57
we have expected, specifying six torsion angles allows us to uniquely specify the
position and orientation of the rigid right boundary segment relative to the left.
In the more general case, some of the torsion angles along the chain may
take on frozen but arbitrary values. Such is the case for proteins, for example,
where the ! angle in each residue is often close to planar. A similar situation is
encountered in nucleic acids, where the angle in each nucleotide is conned to a
limited range of values due to the ring structure of the ribose. In the most general
case, we can imagine the movable segment in the middle is made up of several
sections, each of which may have any number of frozen internal torsion angles,
but the connections between these sections are torsionable. An example is shown
in Fig 3.1(a). As in Fig. 3.1(a), atoms (i 1), (i), (i + 1) and all atoms to the
left are xed in space, and atoms (j 1), (j), (j + 1) and all atoms to the right
form the rigid boundary segment. But instead of allowing all torsion angles to
be
exible in the movable segment in the middle, we now freeze the bonds along
certain parts of the backbone (indicated by the thin lines in Fig. 3.2(a)) inside this
segment such that the sections from atom (i) to (m), (m) to (n) and (n) to (j)
are individually rigid. As such, we can imagine that atoms (i+1) and (m1) are
connected by a ctitious bond, shown in white in Fig. 3.2(b), and we can forget
about the frozen structure of the backbone in between. Similarly, we connect
atoms (m + 1) with (n 1) and (n + 1) with (j 1) by two other ctitious bonds
as shown in Fig. 3.2(b). Replacing the eects of the constraints being imposed
58
by the frozen backbone structures on the atoms by the three ctitious bonds, we
can now reduce the problem to the one shown in Fig. 3.2(c) where the structural
details of the frozen sections can be ignored. Using a similar atom numbering
scheme as in Fig. 3.2(c), we see that specifying the torsion angle
1
now xes
the positions of both atoms (1
0
) and (2). Specifying
2
xes the position of (3).
Specifying
4
xes the positions of (3
0
) and (4), etc. Consequently, specifying six
torsion angles again allows us to uniquely specify the position and orientation of
the right rigid segment relative to the left.
To further examine the geometry of each of the frozen structures in Fig. 3.2(b),
we show details of one of the sections in Fig. 3.3. The ctitious bond shown in
white connects atoms (3) and (3
0
). Since the bond lengths between atoms (2)
and (3) and atoms (3
0
) and (4) are xed, the geometry of this section can be
uniquely specied by the four variables shown. First, r measures the distance
between atoms (3) and (3
0
). Second, describes the bond angle between the
bond (2)-(3) and the ctitious bond (3)-(3
0
). Third, describes the bond angle
between the ctitious bond (3)-(3
0
) and the (3
0
)-(4) bond. Finally, ! describes
the dihedral angle between the (2)-(3) and (3
0
)-(4) bonds along the ctitious
bond (3)-(3
0
). In the rest of this work, we will refer to these four variables as the
\RETO coordinates" (R for r, E for , T for and O for !) of any rigid section
on the chain.
59
Figure 3.2: (a) Similar to Fig. 3.1(c), gure (a) shows details of the movable
segment. But in this more general case, the parts of the chain inside the movable
segment indicated by the thin solid lines are rigid. (b) To account for the rigid
sections inside the movable segment, we connect atoms (i + 1) and (m 1),
(m + 1) and (n 1), and (n + 1) and (j 1) by three ctitious bonds indicated
by the white sticks. (c) An equivalent representation to (b), but using the three
ctitious bonds to represent the constraints introduced by the rigid sections.
60
Figure 3.3: Denition of the four RETO coordinates using the ctitious bond for
a section with rigid backbone indicated by the thin solid lines.
Interestingly, not only are the RETO coordinates useful for dening the struc-
tures of the frozen sections inside the movable segment in the generalized closure
problem depicted in Fig. 3.2, they are also useful for dening the boundary con-
ditions in the closure problem. At the beginning of this section, we have shown
that the boundary constraints imposed by the xed coordinates of the chain out-
side the movable segment can be dened by six variables specifying the position
and orientation of the right boundary segment relative to the left. In Fig. 3.4, we
show the same generalized closure problem from Fig. 3.2(c), but we have added
a fourth ctitious bond between atoms (1) and (5
0
). Using this additional cti-
tious bond, we can now dene the four boundary RETO coordinates shown in
61
Fig. 3.4 using the atom positions of (0), (1), (5
0
) and (6): r
X
being the (1)-(5
0
)
distance,
X
the (0)-(1)-(5
0
) angle,
X
the (1)-(5
0
)-(6) angle and !
X
the (0)-(1)-
(5
0
)-(6) torsion angle. We call these the RETO
X
coordinates. Since a total of six
variables are needed to specify the boundary conditions, two addition variables
are needed beyond the four RETO
X
coordinates. From Fig. 3.4 it is clear that
these two extra variables must correspond to the rotation of the (0)-(1) bond
together with all the atoms on the left segment around the (0)-(1) bond axis and
the rotation of the (6)-(7) bond together with all the atoms on the right segment
around the (5
0
)-(6) bond axis. These two rotation angles are the torsion angles

1
and
6
. Therefore, the RETO
X
coordinates are directly related to the four
torsion angles
2
to
5
. The fact that the four torsion angles
2
to
5
are uniquely
specied by the RETO
X
coordinates can also be recognized in a dierent way.
As shown in Fig. 3.2(c), the atoms (2), (3), (3
0
) and (4) form a frozen and hence
rigid section. Imagine for the moment that this section is xed in space. Then
specifying the torsion angle
3
determines the position of atom (1
0
). Specifying

4
also determines the position of atom (5). Since the atoms (0), (1), (1
0
) and (2)
also form a rigid section, specifying torsion angle
2
will then x the positions of
both (1) and (0). Similarly, since atoms (4), (5), (5
0
) and (6) form another rigid
section, specifying
5
will x the positions of both (5
0
) and (6). Therefore, with
the four torsion angles
2
to
5
, the positions of atoms (0) to (6) are xed. The
62
four torsion angles
2
to
5
are therefore associated with a unique set of RETO
X
coordinates.
Using the RETO
X
coordinates, we have therefore reduced the original 6-
variable/6-constraint closure problem to a simpler two-part problem. In the rst
part, one must nd the four torsion angles
2
to
5
that would satisfy the bound-
ary constraints imposed by the RETO
X
coordinates leading to a 4-variable/4-
constraint problem. This is followed by a second part in which one must nd the
remaining two torsion angles
1
and
6
that would provide the nal rotations
necessary to bring the left and right boundary segments to exactly match the
xed atomic positions outside the movable segment. As we will see later, the
rst part constitutes the bulk of the calculations in the closure problem, while
the second part turns out to be trivial. Thus, the complexity of the original gen-
eralized closure problem is simplied by also employing the RETO
X
coordinates
to represent the boundary conditions.
3.3 Solutions to the Generalized Closure
Problem
In this section, we will describe the solution to the generalized closure problem in
detail. As we have mentioned above, the solution will be carried out in two parts.
The rst part is to solve for the four torsion angles
2
:::
5
inside the movable
63
Figure 3.4: Denition of the RETO
X
coordinates to represent the boundary
conditions in the closure problem shown in Fig. 3.2(c).
segment so that the geometry of the structure formed by atoms (0), (1), (5
0
) and
(6) has the prescribed RETO
X
coordinates to interface with the rigid boundary
segments to the left and right of the movable one. The second part is to nd the
value of
1
for the rotation of the left segment around the (0)-(1) bond and
6
for the rotation of the right segment around the (5
0
)-(6) bond in order to match
the coordinates of the xed boundary segments.
3.3.1 Part 1
The rst part of the problem is the bulk of the calculation. For its geometry,
we refer to Figs. 3.2(c) and 3.4. To proceed, we again note that the generalized
closure problem assumes three rigid sections inside the movable segment. Specif-
ically, the rigid sections are formed by atoms (0), (1), (1
0
) and (2), which we will
call the left or L-element, atoms (2), (3), (3
0
) and (4), the center or C-element,
64
and atoms (4), (5), (5
0
) and (6), the right or R-element. Each of these three
rigid elements has geometry described by its RETO coordinates: RETO
L
for
the L-element, RETO
C
for the C-element and RETO
R
for the R-element. The
boundary conditions are specied by a fourth set of RETO coordinates RETO
X
,
i.e. those depicted in Fig. 3.4. Finally, the four bond lengths (1
0
)-(2), (2)-(3), (3
0
)-
(4) and (4)-(5) as well as the two bond angles (1
0
)-(2)-(3)

and (3
0
)-(4)-(5)

internal to the movable segment are also assumed known and xed.
We will solve Part 1 of the generalized closure problem using basic kinemat-
ics that are familiar to the studies of mechanisms and mechanical manipulators
[Hunt, 1978, Duy, 1980]. The solution will be carried out in two separate steps
as follow.
First, we disconnect the C-element from the L- and R-elements at atoms (2)
and (4) and use the RETO
X
coordinates to dene the orientations of the L- and
R-elements relative to the ctitious bond (1)-(5
0
) [Chung, 2005]. This geometry is
shown in Fig. 3.5(a). We dene a right-handed \external" (X) coordinate system
such that its y
X
-axis lies along the (1)-(5
0
) bond, the (0)-(1) bond lies on the
y
X
-z
X
plane and the origin is on the atom (1).
The RETO
X
geometry uniquely denes the orientation of the (0)-(1) bond
as well as the position of atom (5
0
) and the orientation of the (5
0
)-(6) bond in
this X coordinate system. Since the L-element is rigid and the (0)-(1) bond must
lie in the direction xed by RETO
X
shown in Fig. 3.5(a), the only degree of
65
freedom the L-element has in this frame is a rotation around the (0)-(1) bond.
We call this rotation angle
L
. We dene a \left" (L) coordinate system for the
L-element, where the z
L
axis lies along the (0)-(1) bond and (1) is the origin.
The rotation angle
L
is dened such that
L
= 0 corresponds to the atom (2)
lying on the y
L
-z
L
plane. The coordinates of all four atoms (0), (1), (1
0
) and
(2) in the L coordinate system are xed by the RETO
L
coordinates. Any point
in the L coordinate system can be transformed into the X coordinate system
by the matrix R
x
[(
X
=2)]R
z
[
L
], where R
x;y;z
[
] is the rotation matrix
for a right-handed rotation about the x, y and z axis, respectively, by angle
.
Similarly, since the R-element is rigid and the orientation of the (5
0
)-(6) bond is
xed by the RETO
X
geometry, the only degree of freedom the R-element has is a
rotation around the (5
0
)-(6) bond, with a rotation angle
R
. We dene a \right"
(R) coordinate system for the R-element, where thez
R
axis lies along the (6)-(5
0
)
bond and (5
0
) is the origin. The rotation angle
R
is dened such that
R
= 0
corresponds to the atom (4) lying on they
R
-z
R
plane. The coordinates of all four
atoms (4), (5), (5
0
) and (6) in the R coordinate system are xed by the RETO
R
coordinates. Any point in the R coordinate system can be transformed into
the X coordinate system by the matrix R
y
[!
X
]R
x
[(
X
=2)]R
z
[
R
] followed
by the translation r
X
^ y
X
, where ^ y
X
is the unit vector along the y
X
-axis in the
X coordinate system.
66
Figure 3.5: Drawings showing dierent steps in Part 1 of the solution. See text
for details.
In Fig. 3.5(b), we show details of the L-element alone. The rotation of the
entire L-element occurs around the (0)-(1) bond, which is dened to be the z
L
-
axis in the L coordinate system. We dene two additional vectors.
~
K
L
is the
projection of the vector from atom (1), the origin of the L coordinate system, to
(2) onto the z
L
-axis.
~
J
L
is the projection of the same vector onto the y
L
-axis.
Similarly in Fig. 3.5(c), we show details of the R-element alone. The rotation of
the entire R-element occurs around the (5
0
)-(6) bond, which is dened to be the
z
R
-axis in the R coordinate system.
~
K
R
is the projection of the vector from atom
(5
0
), the origin of the R coordinate system, to (4) onto the z
R
-axis, and
~
J
R
is its
projection onto the y
R
-axis.
67
After setting up the geometries of the L- and R-element in the X coordinate
frame, we proceed to complete the rst step. In this step, we want to solve for
all pairs of (
L
;
R
) such that the resulting distance between atoms (2) and (4)
comes out correctly. This is necessary because the C-element is connected to the
L- and R-elements by these two atoms and since the C-element is also a rigid
segment the distance between the atoms (2) and (4) must also be xed. Since
there are two degrees of freedom (
L
;
R
) but only one constraint, the solution
will be a curve on the
L
-
R
plane. Analogous to robotic systems, it is easiest to
think of
L
as the \input" variable and resulting
R
as the \output" angle which
is considered to be a function of
L
.
To compute the distance between atoms (2) and (4) upon rotations of the L-
and R-elements, we make use of the vectors dened in Figs. 3.5(b) and (c). We
call the (2)-(4) distance prescribed by the C-element H

. Denoting the vector
between atoms (2) and (4) by
~
H =
~
J
L
+
~
I +
~
J
R
where
~
I =
~
K
L
+r
X
^ y
x
+
~
K
R
,
the square distance between them must satisfy
H
2
= j
~
J
L
+
~
I +
~
J
R
j
2
= (j
~
J
L
j
2
+j
~
Ij
2
+j
~
J
R
j
2
)
+2(
~
J
L

~
I +
~
J
R

~
I
~
J
L

~
J
R
): (3.1)
68
Since only rotations are involved, none of the lengths of the vectors are altered;
therefore, only the second term on the right hand side of Eq.(3.1) will depend on

L
and
R
. Explicitly, in the X coordinate system,
~
J
L
=j
~
J
L
j
0
B
B
B
B
B
B
@
cos
L
sin
X
sin
L
cos
X
sin
L
1
C
C
C
C
C
C
A
; (3.2)
~
J
R
=j
~
J
R
j
0
B
B
B
B
B
B
@
cos!
X
cos
R
sin!
X
cos
X
sin
R
sin
X
sin
R
sin!
X
cos
R
cos!
X
cos
X
sin
R
1
C
C
C
C
C
C
A
; (3.3)
~
I =r
X
0
B
B
B
B
B
B
@
0
1
0
1
C
C
C
C
C
C
A
j
~
K
L
j
0
B
B
B
B
B
B
@
0
cos
X
sin
X
1
C
C
C
C
C
C
A
+ j
~
K
R
j
0
B
B
B
B
B
B
@
sin!
X
sin
X
cos
X
cos!
X
sin
X
1
C
C
C
C
C
C
A
: (3.4)
69
The dot products of these three vectors will contain terms linear in the sine and
cosines of
L
and
R
as well as bilinear terms involving products like cos
L
cos
R
,
cos
L
sin
R
, etc. After rearrangement, Eq.(3.1) can be written in the form
0 = (a cos
L
+b sin
L
+c) cos
R
+(d cos
L
+e sin
L
+f) sin
R
(g cos
L
+h sin
L
+i); (3.5)
where the coecients a through i can be worked out easily using Eqs.(3.2) to
(3.4). Treating
L
as the input angle and
R
as the dependent variable, we arrive
at the following equation for
R
:
0 =u cos
R
+v sin
R
+w; (3.6)
where u = a cos
L
+b sin
L
+c, v = d cos
L
+e sin
L
+f, and w = g cos
L
+
h sin
L
+ i. Using the half-angle formula for cos
R
= (1 t
2
)=(1 + t
2
) and
sin
R
= 2t=(1 +t
2
) where t = tan
1
(
R
=2), Eq.(3.6) can be converted into a
quadratic equation in t:
0 = (wu)t
2
+ 2vt + (w +u) (3.7)
70
whose solutions are easily found.
For the purpose of numerics, it is convenient to know the domains of
L
values
on which real-valued solutions to Eq.(3.7) exist. Such domains are given by the
conditionq
2
v
2
w
2
+u
2
0. q
2
is of course only a function of
L
and can be
written as
q
2
= c
1
cos
2

L
+c
2
sin
2

L
+c
3
cos
L
sin
L
+c
4
cos
L
+c
5
sin
L
+c
6
; (3.8)
where the coecients c
1
to c
6
can be expressed in terms of the coecients a to i
in Eq.(3.5) above. Using the half angle formula for cos
L
= (1s
2
)=(1 +s
2
) and
sin
L
= 2s=(1 +s
2
) wheres = tan
1
(
L
=2), Eq.(3.8) becomes a quartic equation
in s:
0 = s
4
(c
1
c
4
+c
6
) +s
3
(2c
3
+ 2c
5
)
+s
2
(2c
1
+ 4c
2
+ 2c
6
)
+s(2c
3
+ 2c
5
) + (c
1
+c
4
+c
6
): (3.9)
This can be solved numerically by standard quartic equation solvers [Glassner,
1990]. Once the domain boundaries for
L
have been identied by the roots to
Eq.(3.9), whether each is an entry or an exit point can easily be determined by the
71
sign of the derivative ofd(q
2
)=d
L
there. As we will see later, inside the domain(s)
identied by Eq.(3.9), there are additional constraints from the second step of the
solution which may further reduce the feasibility region in
L
. These subdomains
boundaries can be found from Eq.(3.15) which subsumes the condition q
2
0,
making the explicit solution of Eq.(3.9) unnecessary in practice.
After the output angle
R
has been determined as a function of the input

L
, we proceed to the second step. From the rst step, we have determined the
possible pairs of (
L
;
R
) values which permit the C-element to be attached to the
L- and R-elements through atoms (2) and (4). Next we need to determine how
to position the C-element so that it yields the correct values for the bond angles
(1
0
)-(2)-(3) and (3
0
)-(4)-(5). We call these prescribed values

and

,
respectively. Since the C-element is itself rigid and it must be attached through
atom (2) on the left and atom (4) on the right, it only has one degree of freedom
corresponding to a rotation of the entire element around the axis connecting (2)
and (4). We call this rotation angle
C
. For each pair of possible (
L
;
R
) values,
we can therefore calculate the value of
C
that would yield the correct

and use
it to compute the resulting. A solution to the loop closure problem is given by
those triplets (
L
;
C
;
R
) which satisfy all the above constraints and at the same
time yield the correct =

for the (3
0
)-(4)-(5) bond angle.
For each pair (
L
;
R
), we determine
C
as follows. The geometry is shown in
Fig. 3.5(d). After rotation of the L-element by
L
, the positions of the atoms (1
0
)
72
and (2) in the X coordinate system are xed. Similarly, the positions of atoms
(4) and (5) are also xed after the rotation of the R-element by
R
. Using these
positions, we can determine the vector
~
M from atom (2) to (1
0
) and the vector
~
N from atom (4) to (5), with corresponding unit vectors
^
M and
^
N, respectively.
A local \C" coordinate system, shown in Fig. 3.6, is dened for the C-element
with its origin on atom (2) and the z
C
-axis lying along
~
H. We denote the vector
from atom (2) to (3)
~
P and that from atom (4) to (3
0
)
~
Q, with corresponding
unit vectors
^
P and
^
Q, respectively. We dene the rotation angle
C
such that

C
= 0 corresponds to
^
P lying on the same plane as both
^
M and
^
H, which we
will take as thex
C
-z
C
plane for the C coordinate system. With the C coordinate
system dened, the transformation between it and the X coordinate system can
easily be obtained by computing the dot products ^ x
C
^ x
X
, ^ x
C
y
X
, ^ x
C
z
X
, etc.
Figure 6 shows details of the C element geometry within the C coordinate system
and the relative positions of all the atoms involved. The (3)-(2)-(4) angle, which
we call , is determined by the RETO
C
geometry. Because we dene
C
= 0
by placing atom (3) on the x
C
-z
c
plane which also contains the atom (1) (shown
on the right side of Fig. 3.6), the value of the (1
0
)-(2)-(3) angle is at a minimum
when
C
= 0. We will call this angle
min
. When the entire C element is rotated
about thez
C
-axis, the value of this angle will approach the proper value

. This
73
rotation is carried out by the matrixR
z
[
C
] in the C coordinate system. In terms
of
^
M and
^
P , the angle is thus given by
cos =
^
M (R
z
[
C
]
^
P ): (3.10)
The rotation angle
C
that is needed to satisfy the condition =

yields the
desired solution. The positions of atoms (3) and (3
0
) after this rotation are
shown on the right side of Fig. 3.6 in grey. If this rotation is continued beyond
this value of
C
, the (1
0
)-(2)-(3) angle will take on its maximum value
max
when

C
= . Clearly,
max
=
min
+ 2. Eq.(3.10) can easily be solved using again
the half-angle formula for cos
C
and sin
C
, leading to a quadratic equation in
tan
1
(
C
=2). Once
C
is found, we can nally calculate the value of , which is
the resulting (3
0
)-(4)-(5) angle, by cos =
^
N (R
z
[
C
]
^
Q).
As we mentioned previously, one can think of
L
as the input (independent)
variable and both
R
and
C
as output (dependent) variables . The solution to
Eq.(3.7) in the rst step gives at most two real solutions for
R
for every
L
. For
each of these, the solution to Eq.(3.10) for
C
gives another two real solutions
at most. Therefore, for every input value of
L
, there are a maximum of four
real solutions for the angle. The loop closure problem is solved if any one of
these solutions matches the prescribed value

for the (3
0
)-(4)-(5) angle for some
74
Figure 3.6: Drawing showing the rotation of the C-element in detail. The bond
axis between atoms (2) and (4) becomes the z
C
axis. The x
C
-axis is dened so
that the atom (1
0
) is on the x
C
-z
C
plane. The atom (3) is initially placed on the
x
C
-z
C
plane dening the rotation angle
C
= 0. After a rotation of the C-element
by
C
, the nal bond angle (1
0
)-(2)-(3) is . See text for further details.
value(s) of
L
. Fig. 3.7 shows two possible solutions for the movable segment
based on the solutions described above.
Just as in the rst step, it will be useful to be able to determine the domain(s)
of
L
on which the solution for
C
is real. A real solution to Eq.(3.10) exists as
long as
min

max
. Thus, the domain boundaries occur at values of
L
for which
min
=

or cos
max
=

. Referring to Fig. 6 again, we see that
75
Figure 3.7: Drawing showing two dierent solutions of the generalized closure
problem.

min
is given by
^
M
^
H = cos(
min
+ ) whereas
max
=
min
+ 2 is given
by
^
M
^
H = cos(
max
). Notice that the angle is completely determined
by the RETO
C
geometry and does not depend on
C
. The angle between
^
H
and
^
M, however, depends on both
L
and
R
because
~
H is determined by the
positions of atom (2) on the L-element and of atom (4) on the R-element, while
^
M is controlled by
L
alone. Therefore, the domain boundaries can be found by
solving the following equation
^
M
^
H = cos(

): (3.11)
76
To do this, we make use of Eqs.(1)-(3). If
~
H has the prescribed length H

,
Eq.(3.11) can be rewritten as
^
M
~
H =H

cos(

): (3.12)
Therefore, the solution to Eq.(3.11) can be found by simultaneously satisfy-
ing Eq.(3.1) and Eq.(3.12). Again taking
L
as the input variable and
R
as the output, Eq.(3.1) can be rewritten in the form of Eq.(3.6) as shown
above. Eq.(3.11) can be simplied by noticing that since
~
H =
~
J
L
+
~
I +
~
J
R
,
^
M
~
H =
^
M
~
J
L
+
^
M
~
I +
^
M
~
J
R
. The second term on the right hand side of the
last equation depends on
L
alone since
~
I is constant. The third term on the right
is a function of
L
and
R
, containing terms bilinear in the sine and cosine of
L
and
R
, like cos
L
cos
R
, cos
L
sin
R
, etc. The rst term on the right, however,
is independent of either
L
or
R
, since both
^
M and
~
J
L
are transformed by the
same L-element rotation and their dot product is thus invariant. The resulting
structure of Eq.(3.12) is very similar to Eq.(3.6) and it can be represented by
0 =u
0
cos
R
+v
0
sin
R
+w
0
; (3.13)
77
where u
0
, v
0
and w
0
are linear functions of the sine and cosine of
L
. Using the
half-angle formula again for cos
R
and sin
R
, Eq.(3.13) becomes
0 = (w
0
u
0
)t
2
+ 2v
0
t + (w
0
+u
0
) (3.14)
To satisfy both Eqs.(3.6) and (3.12) simultaneously, we must nd the common
root to Eqs.(3.7) and (3.14). This can be accomplished by setting the determinant
of the Sylvester matrix detS = 0, where
S =
0
B
B
B
B
B
B
B
B
B
B
@
S
2
S
1
S
0
0
0 S
2
S
1
S
0
S
0
2
S
0
1
S
0
0
0
0 S
0
2
S
0
1
S
0
0
1
C
C
C
C
C
C
C
C
C
C
A
; (3.15)
and S
2
= (w u), S
1
= 2v, S
0
= (w + u), S
0
2
= (w
0
u
0
), S
0
1
= 2v
0
and
S
0
0
= (w
0
+u
0
). Each element of the S matrix is a linear function of the sine and
cosine of
L
. Making use of the half-angle formula for sin
L
= 2s=(1 +s
2
) and
cos
L
= (1s
2
)=(1+s
2
), detS can nally be rewritten as a 8-th order polynomial
in the variable s. The domain boundaries can now be located by solving for the
real roots of detS = 0 using the Sturm chain [Glassner, 1990]. For each root, we
can further determine whether it is an entry or an exit point by evaluating the
derivative d(
^
M
~
H)=d
L
there.
78
3.3.2 Part 2
In Part 1 of the solution, we have determined the rotations of of the L-, C- and
R-elements needed to close the loop with the two boundary segments xed in
space. In Part 2 of the solution, we will now take these rotations and calculate
the proper backbone torsion angles from them.
Again the six torsion angles
1
to
6
are shown in Fig. 3.2(c). The RETO
X
coordinates are xed. The three rotation angles
L
,
C
and
R
allow us to recon-
struct the positions of atoms (1
0
) through (5) inside the movable segment given
the xed boundary condition. Once these atom positions are determined, one can
easily measure the four torsion angles
2
to
5
internal to the movable segment.
To nally determine
1
and
6
, we note that they are directly related to
L
and
R
by
1
=
L
+C
L
and
6
=
R
+C
R
, where C
L
and C
R
are independent
of
L
or
R
. To demonstrate why this is the case, we show the geometry of the
L-element in Fig. 3.8, looking down the z
L
axis (the (0)-(1) bond axis) in the
L-coordinate system, with all atoms projected onto the x
L
-y
L
plane. The atom
positions of (1
0
) and (2) after the
L
rotation are shown in white. The
L
angle is
dened such that
L
= 0 corresponds to atom (2) on they
L
-z
L
plane. The original
torsion angle
1
is (1)-(0)-(1)-(1
0
). Using the ctitious bond (1)-(5
0
), we dene
a new torsion angle
X
1
as (1)-(0)-(1)-(5
0
). Since the boundary conditions is
xed, the
X
1
torsion angle does not change with
L
. From Fig. 3.8, it is then
79
obvious that
1
and
X
1
are related to each other by
1
=
X
1
+
L
, where is
the (2)-(1)-(1
0
) angle when the atoms are projected onto the x
L
-y
L
plane. Since
atoms (0), (1), (1
0
) and (2) form the rigid L-element, is a constant and it must
not depend on
L
. Thus,
1
and
L
are simply related to each other by a constant
oset. Through the same argument, the
6
torsion angle can also be shown to
be related to
R
in a similar way.
Figure 3.8: Drawing showing how the torsion angle
1
is related to
X
1
and
L
.
The picture is drawn looking down the (1)-(0) bond axis. Atom (0) is behind
atom (1) (both in grey).
X
1
is dened as the torsion angle (1)-(0)-(1)-(5),
whereas the original
1
is the (1)-(0)-(1)-(1
0
) angle.
80
3.3.3 Generalized Closure Involving More Than Six
Degrees of Freedom
In Sect. 3.2, we showed that the boundary conditions in the loop closure problem
correspond to six constraints, and if there are exactly six
exible torsion angles
inside the movable segment, the resulting 6-variable/6-constraint problem have
well-dened discrete solutions. Of course, if there are actually more than six

exible torsion angles inside the movable segment, the problem becomes under-
determined. If there aren
exible torsion angles, then 6 undetermined torsion
angles can take on any value in the range [0; 2]. Once these are xed, the last
six torsion angles can be easily determined using the generalized closure solution
above.
3.4 Numerical Implementation of the
Generalized Loop Closure Solution
The solution to the generalized loop closure problem described above can be
carried out using a simple numerical scheme.
First, the RETO
X
coordinates are used to dene the geometry of the boundary
conditions. To describe the geometries of the L-, C- and R-elements, the RETO
L
,
RETO
C
and RETO
R
coordinates are also needed as inputs. In addition, the bond
81
lengths (1
0
)-(2), (2)-(3), (3
0
)-(4), (4)-(5), as well as the two bond angles (1
0
)-(2)-
(3)

and (3
0
)-(4)-(5)

need to be provided as inputs.
Treating
L
as the independent variable, we seek values of
R
and
C
that
would yield a close loop with the prescribed value of

according to the solution
in Sect. 3.3. To do this, we rst determine the domain(s) on
L
for which a
real solution for
R
exists from the solution to Eq.(3.9). This is easily done
using the quartic formula [Glassner, 1990]. Within each of these domains, we
then determine the subdomain(s) on which a real solution for
C
exists from the
solution of detS = 0, with the Sylvester matrix given by Eq.(3.15). Since each
of the coecients of the Sylvester matrix is a quadratic function of tan
1
(
L
=2),
this corresponds to nding the real roots of an 8-th order polynomial equation,
which can be done numerically by making use of a Sturm chain [Glassner, 1990].
After identifying the domains of
L
on which a real solution to the closure
problem is expected to exist, we use standard root nding algorithm to look for
the roots of the nal closure equation(
L
) =

. We rst divide each
L
domain
into small intervals (typically smaller than 0.002 radians) and bracket the real
root(s) [Press et al., 1992]. The bracketed root(s) are then polished by standard
bisection [Press et al., 1992]. The precision of the bracketing can be improved
by choosing smaller intervals. Particularly dicult for bracketing to identify are
double roots or two very close-by roots which may occur within the same interval.
82
To help improve the accuracy of the root search in these case, we also look for
local extrema within the domains and check for existence of root(s) there.
In some rare cases, a root may actually be on one of the domain boundaries.
In this case, neither bracketing nor extrema search will be able to locate this kind
of roots. To circumvent this, we also added a check in the numerical procedure
to look for possible roots on the domain boundaries.
3.4.1 An Example
The generalized closure solution described here can be used to close loops in any
linear polymers, including proteins and nucleic acids. In this section, we will give
an example for RNAs.
Figure 3.9(a) shows a 4-nucleotide segment of a RNA molecule. We number
the atoms as in Fig. 3.2(c). We dene the L-element to be a rigid section con-
sisting of a single nucleotide from atom (0) to atom (2) and similarly the rigid
R-element from atom (4) to (6) consisting of a single nucleotide. On the other
hand, the rigid C-element is taken to be a two-nucleotide section from atom (2)
to (4). The generalized formulation of the closure problem does not require the
movable torsion angles to be contiguous; therefore, it is possible to dene any
rigid element in the generalized closure problem with an arbitrary number of
atoms. We imagine that the atoms (1) and (1
0
), (3) and (3
0
), (5) and (5
0
) in
the L-, C- and R-elements are connected by the three ctitious bonds shown in
83
Figure 3.9: Two solutions of the generalized closure problem for a 4-nucleotide
RNA segment (to reduce clutter, O1P and O2P are not shown). Atoms (0)
through (2) form the rigid L-element, atoms (2) through (4) the C-element and
atoms (4) through (6) the R-element. The L- and R-element consist of one
nucleotide each. The C-element consists of two nucleotides.
84
Fig. 3.9(a). The generalized loop closure algorithm described above was applied
to this 4-nucleotide sequence to look for closure solutions.
The results are shown in Fig. 3.10. Figure 3.10(a) shows the solution from the
rst step in Part 1 described in Sect. 3.3.1, plotting the output value
R
versus
the input
L
obtained by solving Eq.(3.7). The domain boundaries identied by
Eq.(3.15) are indicated by the open squares. In general, there are two solutions
for
R
for each input value of
L
. Each of these two branches will further gen-
erate two solutions in general for
C
according to the simultaneous solutions of
Eqs.(3.12) and (3.1). These are shown in Fig. 3.10(b). Finally, the dierence
between cos(
L
) and the target value cos

is plotted in Fig. 3.10(c). The
roots found by the numerical algorithm described above are indicated by the
open circles. The root corresponding to the native structure is indicated by the
closed circle. Figure 3.9(a) shows the native structure, while Fig. 3.9(b) shows
the reconstructed structure corresponding to one of the other seven solutions.
3.5 Jacobian of the Transformation
To use the generalized loop closure solution for molecular simulations, the Jaco-
bian of the transformation between the torsion angles internal to the movable seg-
ment and the RETO
X
coordinates for the boundary conditions may be required.
85
-1
-0.5
0
τ
R
/π
-1
0
1
τ
C
/π
0 0.5 1 1.5 2
τ
L
/π
-1
-0.5
0
0.5
1
1.5
cosσ(τ
L
) − cosσ∗
(b)
(a)
(c)
Figure 3.10: Solutions of the generalized closure problem for the example in
Sect. 3.4.1. (a) Solution for
R
from the rst step in Part 1 of the solution.
For every
L
there are two solutions for
R
in general. Squares denote domain
boundaries (see text). (b) Solution for
C
from the second step in Part 1. For
each
R
branch in (a), there are in general two solutions for
C
. The solid and
dashed lines in (b) correspond to the two solutions from the
R
branch associated
with the solid line in (a), and the dotted and dotted dashed lines in (b) to the
R
branch associated with the dotted dashed line in (a). (c) Solution for (
L
)

as a function of
L
from the second step in Part 1 of the solution. Each line (c)
corresponds to a
C
branch in (b). The roots are shown as the open circles. The
root corresponding to the native structure is shown as the closed circle.
86
In this section, we will brie
y describe how this can be done. The development
here parallels that previously given by Dodd, et al. [Dodd et al., 1993].
In Sect. 3.2, we show that the 6-variable/6-constraint problem of general-
ized loop closure can be reduced to a 4-variable/4-constraint problem solving for
(
2
;
3
;
4
;
5
) given RETO
X
, followed by a linear oset in the last two tor-
sion angles
1
and
6
. Therefore, only the Jacobian for the transformation
(
2
;
3
;
4
;
5
)$ (r
X
;
X
;
X
;!
X
) needs to be computed, while the Jacobian for
the rest is unity.
To illustrate how the Jacobian is calculated, we display the geometry schemat-
ically in Fig. 3.11. For clarity, we have employed a dierent atom numbering
scheme in this section compared to the rest of the chapter. Atom (0) is the rst
atom on the L-element. This was also denoted atom (0) in Fig. 3.2(c). But start-
ing with the next atom, we number the rest sequentially up to (9) for the last
atom on the R-element. In the current numbering scheme, the L- and C-elements
are connected by atom (3), which was (2) in Fig. 3.2(c). The C- and R-elements
are connected by atom (6), which was (4) in Fig. 3.2(c). (Notice that atoms
(1) and (7) in Fig. 3.2(c) are not needed for calculating the Jacobian.) For
convenience, we have drawn Fig. 3.11 showing all the relevant bond and torsion
angles schematically; in reality, the atoms are of course not necessarily on the
same plane. We dene the unit vector along the (1)-(0) bond axis to be ^ a and
that along (8)-(9) to be
^
b. We also denote the bond angle at atomj
j
, the bond
87
length between (j) and (j + 1) q
j
and the torsion angle around the (j)-(j + 1)
bond

j+1
. Under this numbering scheme,

3
,

4
,

6
and

7
correspond to the
original torsion angles
2
to
5
, and

2
,

5
and

8
to!
L
,!
C
and!
R
, respectively.
On the other hand, the bond angles have the following correspondences:
1
and

2
to
L
and
L
,
4
and
5
to
C
and
C
,
7
and
8
to
R
and
R
,
3
to and
6
to .
Figure 3.11: Schematic drawing showing the L-, C- and R-elements and the
denitions of the lengths q
i
, bond angles
i
and torsion angles

i
needed to
calculate the Jacobian of the RETO
X
$ torsion angle transformation.
For every atomj, we dene aj-th right-handed coordinate system with origin
at atom j and its x
j
axis along the (j)-(j + 1) bond. We dene the y
j
axis such
that the atom (j1) lies on thex
j
-y
j
plane. The corresponding basis vectors are
^ x
j
, ^ y
j
and ^ z
j
. We denote the vector position of the origin of the j-th coordinate
system by a vector
~
C
j
. Starting with the 1-st coordinate system on the left,
88
the coordinates for atoms (0) and (1) in this local coordinate system are simply
(cos
1
; sin
1
; 0) and (0; 0; 0). The procedure now is to transform the coordinate
system step by step by using the bond lengths, bond angles and torsion angles
along the contour of the chain to arrive at the total transformation needed to
go from coordinate system 1 to 8. It is easy to see that the origin and the basis
vectors of the j + 1-st coordinate system are related to those in j by:
~
C
j+1
=
~
C
j
+q
j
^ x
j
(3.16)
0
B
B
B
B
B
B
@
^ x
j+1
^ y
j+1
^ z
j+1
1
C
C
C
C
C
C
A
= M
j+1

0
B
B
B
B
B
B
@
^ x
j
^ y
j
^ z
j
1
C
C
C
C
C
C
A
; (3.17)
whereM
j+1
=R
x
[

j+1
]R
z
[
j+1
] andR
x
andR
z
are the right-handed rotation
matrices dened in Sect. 3.2. At the end, the vector
~
D
~
C
8

~
C
1
is given by
~
D = (q
1
+q
2
M
2
+q
3
M
2
M
4
+:::
+q
7
M
2
M
3
M
4
M
5
M
6
M
7
) ^ x
1
; (3.18)
and
^
b is related to ^ x
1
by
^
b = (M
2
M
3
M
8
) ^ x
1
.
89
In terms of ^ a,
^
b and
~
D, the RETO
X
coordinates are given by
r
X
= (
~
D
~
D)
1=2
cos
X
= ^ a
~
D=(
~
D
~
D)
1=2
cos
X
=
^
b
~
D=(
~
D
~
D)
1=2
cos!
X
=

^ ar
1
X
cos
X
~
D

^
b +r
1
X
cos
X
~
D

sin
X
sin
X
=

^ a
^
b
(^ a
~
D)(
^
b
~
D)
~
D
~
D
!
=(sin
X
sin
X
): (3.19)
While it is possible to evaluate the Jacobian J =
@(r
X
;
X
;
X
;!
X
)=@(

3
;

4
;

6
;

7
), it is easier to dene the four new vari-
ables p
1
=
~
D
~
D, p
2
= ^ a
~
D, p
3
=
^
b
~
D, and p
4
= ^ a
^
b, and compute
J
0
=@(r
X
;
X
;
X
;!
X
)=@(p
1
;p
2
;p
3
;p
4
) instead. It is trivial to show thatJ andJ
0
are related to each other by a multiplicative constant for all conformations with
the same RETO
X
boundary conditions. To calculate J
0
, we need the partial
derivatives @p
k
=@

l
, k = 1,2,3,4 and l = 3,4,6,and 7. These involve various
derivatives of
~
D and
^
b with respect to

l
. For instance,
@
~
D
@

3
= M
2
dM
3
d

3
(q
7
M
4
M
5
M
6
M
7
+q
6
M
4
M
5
M
6
+ +q
3
) ^ x
1
;
@
^
b
@

3
= M
2
dM
3
d

3
M
4
M
5
M
6
M
7
M
8
x
1
;
90
wheredM
3
=d

3
=dR
x
[

3
]=d

3
R
z
[
3
]. The other partial derivatives can be
obtained in a similar way. We have checked the validity of the analytical formulae
by computing dRETO=d derivatives numerically during a simulation.
3.6 Monte Carlo Step and Preliminary Results
Now we move on to constructing an elementary Monte Carlo step based on the
generic closure algorithm we discussed in the previous sections. Several consider-
ations have to be taken into account. We design a move, which will be reversible
between the initial and nal states, i.e. satisfying detailed balance. All extra
torsional degrees of freedom, which are not used by the closure algorithm can
be "driven", i.e. randomly changed. Obviously, such random changes break the
chain loop, so, the generic closure algorithm needs to be applied to make it up
using available degrees of freedom, and it does not warrantee that the loop can
always be closed again.
Fig. 3.12 provides an RNA illustration of drivable backbone coordinates in
the case when left, central, or right element consist of a single residue only.
Nondrivable angles are also illustrated in Fig. 3.2(c). If a RETO element consists
of several residues one can drive all the backbone angles in intermediate residues.
The maximum change per single driver angle per MC step is denoted by
max
.
As mentioned in section 3.3.3 each driver angle has the potential to change from 0
91
O P
O
-
O
O
OH O
H H
H H
P O O
-
Base
Figure 3.12: RNA example. Drivable backbone torsion angles in a single residue
are in solid. Dashed torsion angles can be driven when they are not part of the
generic closure problem, e.g. when a RETO element consists of three residues
and this residue is in the middle.
to 2. We label variables corresponding to original state as "orig" and destination
state as "dest" and use the following scheme for the Monte Carlo move, similar
to the one proposed by [Dodd et al., 1993]:
I Choose number of residues for RETO
L
, RETO
C
, and RETO
R
element (can
be dierent for each element)
II For each driver anglej, randomly select
j
from within (
max
;
max
),
and set its trial value to
dest
j
=
orig
j
+
j
92
III Obtain RETO
L
, RETO
C
, and RETO
R
parameters corresponding to this
trial move. RETO
X
parameters should be acquired from the original cong-
uration
IV Find all solution sets corresponding to the proposed move. Possible outcomes
are the following:
(a) There is no solution and attempted move is rejected.
(b) There are N
dest
solutions. Based on uniform random distribution, we
pick one ofN
dest
solution sets of torsion angles of the movable segment.
V Original and trial states may have dierent number of solutions. Trial solu-
tion is randomly picked from N
dest
solutions, i.e. with
orig;dest
= 1=N
dest
probability. The move should be reversible, thus, attempted backward prob-
ability should be
dest;orig
= 1=N
orig
[Allen and Tildesley, 1987]
VI The closure problem requires temporary change of variables from torsion
angles to RETO coordinates, and consequently, the Jacobian needs to be
calculated for both original and trial states, we denote it as J(orig) and
J(dest) correspondingly.
VII Calculate the Boltzmann factor for both original and destination states:
exp(V (orig)=k
B
T ) and exp(V (dest)=k
B
T )
93
VIII Accept the destination state with probability
A
orig;dest
= min

1;

dest;orig
exp(V (dest)=k
B
T )J(dest)

orig;dest
exp(V (orig)=k
B
T )J(orig)

It can be shown that the algorithm with this acceptance probability satises
detailed balance [Dodd et al., 1993, Allen and Tildesley, 1987]
We have implemented this method in FORTRAN and tested it in the TINKER
molecular modeling software. We applied the MC move to a part of Tetrahymena
ribozyme group I intron [Guo et al., 2004], consisting of 247 residues. This
molecule is representative of a rare group of introns which are self-splicing. This
method allowed us to move long stem-like part of the molecule on a large scale,
not accessible by traditional algorithms. Fig. 3.13(a,b) shows an example of such
transition. Fig. 3.13(a) corresponds to the original conguration and Fig. 3.13(b)
is a sample from a short MC simulation at room temperature using the Amber94
forceeld and GBSA solvent model.
Though this method provides correct geometries of the three backbone RETO
elements, it does not have any information about the conguration of the side
chains and the rest of the backbone. As result, it still leads to multiple steric
clashes during the simulation and these clashes are partially responsible for the
low acceptance ratios at the room temperature. But since we made the backbone
changes relatively local, this problem happens at considerably smaller scale than
94
(a)
(b)
Figure 3.13: Original conguration of tetrahymena ribozyme [Guo et al., 2004]
molecule (a). An example of a new conguration after the MC move is applied
(b). Encircled area is the range where the boundary conditions are set up and the
generic closure methods is applied. Right and left RETO elements consist of up
to two nucleotides, the central RETO element (stem) consists of 75 nucleotides.
95
a simple random torsion angle move, when one propagates its random change
throughout the whole chain. Moreover, the constraint of xed bond angles and
bond distances is too restrictive. Originally proposed by [Vangunsteren and
Karplus, 1982, Bruccoleri and Karplus, 1985, Sartori et al., 1998] to improve
conformational sampling by neglect of uninteresting high-frequency motions, it
has been noted that applications of such constraints to angles and bond distances
considerably slows down MD dynamics. It was also pointed out that, the use of
traditional force elds results in unwanted artifacts, since the force elds are not
parameterized for rigid models [Sartori et al., 1998]. The reasons for these arti-
facts are that small
uctuations of bond lengths, and especially bond angles, may
considerably reduce rotational barriers. Thus, eective torsional potentials could
be used to compensate for the increased torsional barriers in the presence of con-
straints [Sartori et al., 1998]. These problems are present for all torsion angle
based Monte Carlo techniques and have some alternative solutions [Ulmschnei-
der and Jorgensen, 2003, 2004a]. Currently, in order to overcome aforementioned
diculties and improve acceptance ratio of global moves, the Van der Waals
repulsion terms have to be scaled down by increasing the temperature.
96
3.7 Conclusions and Future Work
We have presented a new MC move which modies internal degrees of freedom
on local parts of polymer chains. The main advantage of this algorithm over
other methods lies in its
exibility and generality. One can formulate the closure
problem for any number of residues per RETO element. It can also be used to
close single residues or sugar rings [Mak, 2007] without major modications in
the code.
One possible application of this method could be to choose the size of a
movable segment on the
y during simulation based on its conformational freedom
with respect to surroundings and the secondary structure, according to a mapping
scheme or a rule created or set in advance before the simulation as we did in the
example in section 3.6. In general, it is not trivial to identify such mapping, so,
random changes of the elements' size have to be introduced.
These moves, which do not generally modify the secondary structure, but
still perform large collective conformational changes can be useful while explor-
ing the congurational space near the equilibrium structure. The fact that the
biological function of bio-polymers, such as proteins and nucleic acids is not
completely dened by the static picture of their folded state has received much
attention recently. They should be rather seen as 'fuzzy' structures that sample
97
a large ensemble of conformations around the average structure as result of ther-
mal
uctuations. Obviously this conformations happen on dierent timescales
since the energy landscape is extremely diverse. For example, protein motions
on a 'slow' timescale dene
uctuations between the states which are separated
by energy barriers of several kT , corresponding to timescales of microseconds or
slower at the room temperature [Henzler-Wildman and Kern, 2007]. Such pro-
tein 'dynamics' is very interesting, because many biological processes - including
enzyme catalysis, signal transduction and protein-protein interaction - occur on
this timescale. In realistic simulations this generic closure algorithm should not
be used alone but rather combined with traditional, small scale MC moves. For
example, they could be side chain rotations and even MD like small changes in
cartesian coordinates of individual/selected atoms. Such local changes should
facilitate the conformational search of the larger scale moves [Ulmschneider and
Jorgensen, 2004b].
Currently, Monte Carlo backbone sampling techniques are limited to usage of
continuum solvent models. Explicit solvent models, popular in MD simulations,
are not feasible here. Each MC step would have to be followed by equilibration
process of individual solvent molecules, which is not practical. Incorporation of
explicit solvent eects into MC simulation might be a challenging task.
98
Bibliography
M.P. Allen and D.J. Tildesley. Computer Simulation of Liquids. Clarendon Press,
Oxford, 1987.
E. L. Andronikashvili. Journal of Physics USSR, 10:201, 1946.
R. A. Aziz, F. R. W. McCourt, and C. C. K. Wong. A new determination
of the ground-state interatomic potential for he-2. Molecular Physics, 61(6):
1487{1511, 1987.
M. Bander and C. Itzykson. Group theory and hydrogen atom .1. Reviews of
Modern Physics, 38(2):330, 1966.
D. Bashford and D. A. Case. Generalized born models of macromolecular solva-
tion eects. Annual Review of Physical Chemistry, 51:129{152, 2000.
M. Behrens, U. Buck, R. Frochtenicht, M. Hartmann, F. Huisken, and
F. Rohmund. Rotationally resolved ir spectroscopy of ammonia trapped in
cold helium clusters. Journal of Chemical Physics, 109(14):5914{5920, 1998.
N. Blinov, X. G. Song, and P. N. Roy. Path integral monte carlo approach for
weakly bound van der waals complexes with rotations: Algorithm and bench-
mark calculations. Journal of Chemical Physics, 120(13):5916{5931, 2004.
R. E. Bruccoleri and M. Karplus. Chain closure with bond angle variations.
Macromolecules, 18(12):2767{2773, 1985.
U. Buck, K. H. Kohl, A. Kohlhase, M. Faubel, and V. Staemmler. Rotationally
inelastic-scattering and potential calculations for he+ch4. Molecular Physics,
55(6):1255{1274, 1985.
D. M. Ceperley. Path-integrals in the theory of condensed helium. Reviews of
Modern Physics, 67(2):279{355, 1995.
W. Y. Chung. Mobility analysis of rssr mechanisms by working volume. Journal
of Mechanical Design, 127(1):156{159, 2005.
99
E. A. Coutsias, C. Seok, M. P. Jacobson, and K. A. Dill. A kinematic view of
loop closure. Journal of Computational Chemistry, 25(4):510{528, 2004.
M. W. Deem and J. S. Bader. A congurational bias monte carlo method for
linear and cyclic peptides. Molecular Physics, 87(6):1245{1260, 1996.
M. Dikovsky. Quantum monte carlo methods for fermionic systems in chemical
physics. PhD Thesis, 2003.
A. R. Dinner. Local deformations of polymers with nonplanar rigid main-chain
internal coordinates. Journal of Computational Chemistry, 21(13):1132{1144,
2000.
L. R. Dodd, T. D. Boone, and D. N. Theodorou. A concerted rotation algorithm
for atomistic monte-carlo simulation of polymer melts and glasses. Molecular
Physics, 78(4):961{996, 1993.
J. Duy. Analysis of mechanisms and robot manipulators. Wiley, New York,
1980.
G. Favrin, A. Irback, and F. Sjunnesson. Monte carlo update for chain molecules:
Biased gaussian steps in torsional space. Journal of Chemical Physics, 114(18):
8154{8158, 2001.
Richard P. Feynman and Albert R. Hibbs. Quantum Mechanics and Path Inte-
grals. McGraw-Hill, New York, 1965.
A.S. Glassner, editor. Graphics gems. Academic press, San Diego, 1990.
N. Go, P. N. Lewis, and H. A. Scheraga. Calculation of conformation of pen-
tapeptide cyclo(glycylglycylglycylprolylprolyl) .2. statistical weights. Macro-
molecules, 3(5):628, 1970.
H. Goldstein, Charles Poole, and John Safko. Classical Mechanics. Addison
Wesley, 3 edition, 2002.
S. Goyal, D. L. Schutt, and G. Scoles. Vibrational spectroscopy of sulfur-
hexa
uoride attached to helium clusters. Physical Review Letters, 69(6):933{
936, 1992.
S. Grebenev, J. P. Toennies, and A. F. Vilesov. Super
uidity within a small
helium-4 cluster: The microscopic andronikashvili experiment. Science, 279
(5359):2083{2086, 1998.
100
S. Grebenev, M. Hartmann, M. Havenith, B. Sartakov, J. P. Toennies, and A. F.
Vilesov. The rotational spectrum of single ocs molecules in liquid he-4 droplets.
Journal of Chemical Physics, 112(10):4485{4495, 2000.
F. Guo, A. R. Gooding, and T. R. Cech. Structure of the tetrahymena ribozyme:
Base triple sandwich and metal ion at the active site. Molecular Cell, 16(3):
351{362, 2004.
J. Harms, M. Hartmann, J. P. Toennies, A. F. Vilesov, and B. Sartakov. Rota-
tional structure of the ir spectra of single sf6 molecules in liquid he-4 and he-3
droplets. Journal of Molecular Spectroscopy, 185(1):204{206, 1997.
M. Hartmann, R. E. Miller, J. P. Toennies, and A. Vilesov. Rotationally resolved
spectroscopy of sf6 in liquid-helium clusters - a molecular probe of cluster
temperature. Physical Review Letters, 75(8):1566{1569, 1995.
M. Hartmann, J. P. Toennies, A. F. Vilesov, and G. Benedek. Spectroscopic evi-
dence for super
uidity in liquid he droplets. Czechoslovak Journal of Physics,
46:2951{2956, 1996.
K. Henzler-Wildman and D. Kern. Dynamic personalities of proteins. Nature,
450(7172):964{972, 2007.
K. A. Henzler-Wildman, M. Lei, V. Thai, S. J. Kerns, M. Karplus, and D. Kern.
A hierarchy of timescales in protein dynamics is linked to enzyme catalysis.
Nature, 450(7171):913{U27, 2007.
Gerhard Herzberg. Molecular Spectra and Molecular Structure, volume 3. Krieger
Publishing Company, Malabar, Florida, 1991.
Alfred Huller and Werner Press. Rotational tunneling in solids: Theory of neutron
scattering. ii. Physical Review B, 24(1):17, 1981.
K. H. Hunt. Kinematic geometry of mechanisms. Oxford, 1978.
E. W. Knapp. Long-time dynamics of a polymer with rigid body monomer units
relating to a protein model - comparison with the rouse model. Journal of
Computational Chemistry, 13(7):793{798, 1992.
E. W. Knapp and A. Irgensdefregger. O-lattice monte-carlo method with con-
straints - long-time dynamics of a protein model without nonbonded interac-
tions. Journal of Computational Chemistry, 14(1):19{29, 1993.
Y. Kwon, P. Huang, M. V. Patel, D. Blume, and K. B. Whaley. Quantum
solvation and molecular rotations in super
uid helium clusters. Journal of
Chemical Physics, 113(16):6469{6501, 2000.
101
Y. K. Kwon and K. B. Whaley. Atomic-scale quantum solvation structure in
super
uid helium-4 clusters. Physical Review Letters, 83(20):4108{4111, 1999.
Y. K. Kwon, D. M. Ceperley, and K. B. Whaley. Path integral monte carlo study
of sf6-doped helium clusters. Journal of Chemical Physics, 104(6):2341{2348,
1996.
L. D. Landau. The theory of super
uidity of he ii. Journal of Physics USSR, 5:
71, 1941.
L. D. Landau. On the theory of super
uidity of helium ii. Journal of Physics
USSR, 11:91, 1947.
K. K. Lehmann. Rotation in liquid he-4: Lessons from a highly simplied model.
Journal of Chemical Physics, 114(10):4643{4648, 2001.
E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino,
and R. L. Sugar. Sign problem in the numerical-simulation of many-electron
systems. Physical Review B, 41(13):9301{9307, 1990.
C. H. Mak. Journal of Computational Chemistry, (in press), 2007.
C. H. Mak and D. Chandler. Solving the sign problem in quantum monte-carlo
dynamics. Physical Review A, 41(10):5709{5712, 1990.
C. H. Mak and S. Zakharov. A multigrid algorithm for sampling imaginary-time
paths in quantum monte carlo simulations, 2004.
C. H. Mak, R. Egger, and H. Weber-Gottschick. Multilevel blocking approach
to the fermion sign problem in path-integral monte carlo simulations. Physical
Review Letters, 81(21):4533{4536, 1998.
C. H. Mak, S. Zakharov, and D. B. Spry. Super
uidity in ch4-doped h-2 nan-
oclusters. Journal of Chemical Physics, 122(10):104301, 2005.
Markovskiy N.D. Wem-Yeuan Chung Mak, C.H. A generalized kinematic for-
mulation of loop closures in nucleic acids, proteins and other linear polymers.
Journal of Computational Chemistry, to pe published, 2008.
D. Marx. Rotational motion of linear-molecules in 3 dimensions - a path-integral
monte-carlo approach. Molecular Simulation, 12(1):33{48, 1994.
Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, and
Augusta H. Teller. Equation of state calculations by fast computing machines.
Journal of Chemical Physics, 21(6):1087, 1953.
102
S. Moroni, A. Sarsa, S. Fantoni, K. E. Schmidt, and S. Baroni. Structure, rota-
tional dynamics, and super
uidity of small ocs-doped he clusters. Physical
Review Letters, 90(14):143401, 2003.
M. H. Muser and B. J. Berne. Path-integral monte carlo scheme for rigid tops:
Application to the quantum rotator phase transition in solid methane. Physical
Review Letters, 77(13):2638{2641, 1996.
K. Nauta and R. E. Miller. Metastable vibrationally excited hf (v=1) in helium
nanodroplets. Journal of Chemical Physics, 113(21):9466{9469, 2000.
K. Nauta and R. E. Miller. Rotational and vibrational dynamics of methane in
helium nanodroplets. Chemical Physics Letters, 350(3-4):225{232, 2001a.
K. Nauta and R. E. Miller. The vibrational and rotational dynamics of acetylene
solvated in super
uid helium nanodroplets. Journal of Chemical Physics, 115
(18):8384{8392, 2001b.
A. Norcli. Relating classical and quantal theories of angular-momentum and
spin. Journal of Physics Part a General, 5(2):231, 1972.
R. T. Pack, E. Piper, G. A. Pfeer, and J. P. Toennies. Multiproperty empirical
anisotropic intermolecular potentials .2. hesf6 and nesf6. Journal of Chemical
Physics, 80(10):4940{4950, 1984.
F. Paesani and K. B. Whaley. Rotational excitations of n2o in small helium
clusters and the role of bose permutation symmetry. Journal of Chemical
Physics, 121(11):5293{5311, 2004.
F. Paesani, A. Viel, F. A. Gianturco, and K. B. Whaley. Transition from molec-
ular complex to quantum solvation in (henocs)-he-4. Physical Review Letters,
90(7):073401, 2003.
W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical
Recipes. Cambidge, New York, 1 edition, 1992.
W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical
Recipes in C: The Art of Scientic Computing. Cambidge Univercity Press, 2
edition, 1995.
S. Rudolph, G. Wollny, K. von Haeften, and M. Havenith. Probing collective exci-
tations in helium nanodroplets: Observation of phonon wings in the infrared
spectrum of methane. Journal of Chemical Physics, 126(12):{, 2007.
103
F. Sartori, B. Melchers, H. Bottcher, and E. W. Knapp. An energy function
for dynamics simulations of polypeptides in torsion angle space. Journal of
Chemical Physics, 108(19):8264{8276, 1998.
I. F. Silvera. Ultimate fate of a gas of atomic-hydrogen in a liquid-helium chamber
- recombination and burial. Physical Review B, 29(7):3899{3904, 1984.
Dmitry S. Skvortsov, Hiromichi Hoshina, Boris G. Sartakov, and Andrey Vilesov.
Rotation of methane and silane molecules in helium droplets. In 62nd Inter-
national Symposium on Molecular Spectroscopy, Ohio State University, 2007.
M. N. Slipchenko and A. F. Vilesov. Spectra of nh3 in he droplets in the 3 mu
m range. Chemical Physics Letters, 412(1-3):176{183, 2005.
R. H. Swendsen and J. S. Wang. Nonuniversal critical-dynamics in monte-carlo
simulations. Physical Review Letters, 58(2):86{88, 1987.
J. Tang and A. R. W. McKellar. Cluster dynamics in the range n=2-20: High
resolution infrared spectra of he-n-co. Journal of Chemical Physics, 119(2):
754{764, 2003.
J. Tang, Y. J. Xu, A. R. W. McKellar, and W. Jager. Quantum solvation of
carbonyl sulde with helium atoms. Science, 297(5589):2030{2033, 2002.
J. Tang, A. R. W. McKellar, E. Mezzacapo, and S. Moroni. Bridging the gap
between small clusters and nanodroplets: Spectroscopic study and computer
simulation of carbon dioxide solvated with helium atoms. Physical Review
Letters, 92(14):145503, 2004.
Michael Tinkham. Group Theory and Quantum Mechanics. McGraw-Hill, New
York, 1964.
L. Tisza. Natire, 141:913, 1938.
J. P. Toennies and A. F. Vilesov. Spectroscopy of atoms and molecules in liquid
helium. Annual Review of Physical Chemistry, 49:1{41, 1998.
H.F. Trotter. Proceedings of American Mathematical Society, 10, 1959.
J. P. Ulmschneider and W. L. Jorgensen. Monte carlo backbone sampling for
polypeptides with variable bond angles and dihedral angles using concerted
rotations and a gaussian bias. Journal of Chemical Physics, 118(9):4261{4271,
2003.
104
J. P. Ulmschneider and W. L. Jorgensen. Monte carlo backbone sampling for
nucleic acids using concerted rotations including variable bond angles. Journal
of Physical Chemistry B, 108(43):16883{16892, 2004a.
J. P. Ulmschneider and W. L. Jorgensen. Polypeptide folding using monte carlo
sampling, concerted rotation, and continuum solvation. Journal of the Ameri-
can Chemical Society, 126(6):1849{1857, 2004b.
W. F. Vangunsteren and M. Karplus. Eect of constraints on the dynamics of
macromolecules. Macromolecules, 15(6):1528{1544, 1982.
H. Wakana, H. Wako, and N. Saito. Monte-carlo study on local and small-
amplitude conformational
uctuation in hen egg-white lysozyme. International
Journal of Peptide and Protein Research, 23(3):315{323, 1984.
W. J. Wedemeyer and H. A. Scheraga. Exact analytical loop closure in proteins
using polynomial equations. Journal of Computational Chemistry, 20(8):819{
844, 1999.
Y. J. Xu, W. Jager, J. Tang, and A. R. W. McKellar. Spectroscopic studies of
quantum solvation in he-4(n)-n2o clusters. Physical Review Letters, 91(16):
163401, 2003.
R. N. Zare. Angular Momentum. Wiley, New York, 1988.
R. E. Zillich, F. Paesani, Y. Kwon, and K. B. Whaley. Path integral methods
for rotating molecules in super
uids. Journal of Chemical Physics, 123(11):12,
2005.
105
Appendices
106
Appendix A
Rotations in Four Dimensions
Rotations in four space dier from the rotations in three space, because one
axis is not enough to dene the rotation. More general approach to rotations in
any space is to dene them as a series of rotations in corresponding coordinate
planes. In two space it is just one plane, XY , in three space it can be dened
by a sequence of rotations in three coordinate planes: XY , YZ, ZX. In
four dimensions one has XY , YZ, ZX, XW , YW and ZW planes. Thus,
any rotation in 4D can be thought as a product of six rotational matrixes
RR
xy
R
yz
R
zx
R
xw
R
yw
R
zw
:
R
xy
=
2
6
6
6
6
6
6
6
6
6
6
6
6
4
cos (
xy
) sin (
xy
) 0 0
sin (
xy
) cos (
xy
) 0 0
0 0 1 0
0 0 0 1
3
7
7
7
7
7
7
7
7
7
7
7
7
5
R
yz
=
2
6
6
6
6
6
6
6
6
6
6
6
6
4
1 0 0 0
0 cos (
yz
) sin (
yz
) 0
0 sin (
yz
) cos (
yz
) 0
0 0 0 1
3
7
7
7
7
7
7
7
7
7
7
7
7
5
107
R
zx
=
2
6
6
6
6
6
6
6
6
6
6
6
6
4
cos (
zx
) 0 sin (
zx
) 0
0 1 0 0
sin (
zx
) 0 cos (
zx
) 0
0 0 0 1
3
7
7
7
7
7
7
7
7
7
7
7
7
5
R
xw
=
2
6
6
6
6
6
6
6
6
6
6
6
6
4
cos (
xw
) 0 0 sin (
xw
)
0 1 0 0
0 0 1 0
sin (
xw
) 0 0 cos (
xw
)
3
7
7
7
7
7
7
7
7
7
7
7
7
5
R
yw
=
2
6
6
6
6
6
6
6
6
6
6
6
6
4
1 0 0 0
0 cos (
yw
) 0 sin (
yw
)
0 0 1 0
0 sin (
yw
) 0 cos (
yw
)
3
7
7
7
7
7
7
7
7
7
7
7
7
5
R
zw
=
2
6
6
6
6
6
6
6
6
6
6
6
6
4
1 0 0 0
0 1 0 0
0 0 cos (
zw
) sin (
zw
)
0 0 sin (
zw
) cos (
zw
)
3
7
7
7
7
7
7
7
7
7
7
7
7
5
In this matrix notation, a four component vector is dened as (Eqn. 2.6):
e =
0
B
B
B
B
B
B
B
B
B
B
@
e
x
e
y
e
z
e
w
1
C
C
C
C
C
C
C
C
C
C
A

0
B
B
B
B
B
B
B
B
B
B
@
e
1
e
2
e
3
e
0
1
C
C
C
C
C
C
C
C
C
C
A
Transformations R
1
0
(

i
;

i+1
) and R
1
e
0
e
2
(

i
;

i+1
) (see section 2.2.3) are
obtained in terms of these six rotation matrices during the simulation (routines
sample.f ang4d and rot4d of the "Cluster" program).
108
Appendix B
Extremum's Search of the
Rotational Propagator
This appendix implies the completion of reading of the section 2.2.3.
Here we assume, that before our analysis, initial position of the propagator
is rotated towards (1; 0; 0; 0) and the nal position is in e
1
e
2
plane, as it was
discussed in section 2.2.3. Then, for the simplicity sake,as it was done in the
program, we redene the Euler angles:
+
2
,

2
and

2
Hawing it
in mind the scalar products of the propagator Eqn. 2.11 can be written as
e
0
e = cos cos
e e
0
= sin
0
sin cos cos
0
cos cos (B.1)
Here, one has two parameters:
0
and ("+" if
0
= 0 and "" if
0
= ). In
order to nd the extrema, one needs to dierentiate Eqn. (2.11). Having done
that, in the case of and , the extremum points are set to be at 0;

2
;;
3
2
and
2. As soon as and are put on the grid, numerical one dimensional maximum
109
search in is performed. This procedure does not insure 100% of the maxima are
recovered, but it does not seem to aect noticeably the free particle distribution,
which can be compared against analytically, and the properties been measured
so far.
110
Appendix C
Rejection Sampling
The rejection method of sampling from probability distribution p(x) can be
described as follows. First, the distribution in question is covered with the more
computationally feasible comparison function f(x) which lies always above p(x)
Fig. C.1. Then one picks a random number between 0 and A, where A is the
total area under f(x), and use it to get corresponding x. Then one picks a uni-
form sample between 0 and f(x). This would generate uniformly distributed set
of points under f(x). Next, one compares f(x) with p(x). If f(x) > p(x) one
rejects the point, otherwise accepts it, creating a uniform set of points underp(x).
This should work for any number of dimensions. The closer the ratio of volumes
under f(x) and p(x) is to unity the fewer rejections the method will produce.
111
Figure C.1: Rejection method for generating a random sample x from a known
probability distribution p(x) that is everywhere less than some other function
f(x). The transformation method is rst used to generate a random sample x of
the distributionf. A second uniform sample is used to decide whether to accept
or reject that x. If it is rejected, a new sample of f is found; and so on. The
ratio of accepted to rejected points is the ratio of the area under p to the area
betweenp andf. This gure and caption are borrowed from (Press et al., 1995).
112

Abstract (if available)

Abstract One of the main challenges in the simulations of many body systems in statistical physics is the problem of generating sufficient coverage of phase space in a reasonable amount of computer time to produce accurate estimates of equilibrium properties. This work is focused on constructing and applying two different nonlocal Monte Carlo techniques which incorporate large-scale changes in a number of degrees of freedom. These methods were then used to study the rotational dynamics of spherical top molecules in 4He clusters and the problem of conformational sampling of RNA molecules.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Spectroscopy of hydrogen and water clusters in helium droplets

PDF

Rovibrational spectroscopy in helium droplets

PDF

Location, ionization and shielding of impurities in helium nanodroplets

PDF

Infrared and Raman spectrosopy of molecules and molecular aggregates in helium droplets

PDF

Atomic and molecular clusters at ultra-low temperature

PDF

Vorticity in superfluid helium nanodroplets

PDF

Imaging superfluid helium droplets

PDF

Evaporative attachment of slow electrons to free sodium clusters

PDF

Infrared spectroscopy of carbocations in helium droplets

PDF

A technique for assessing water cluster sizes by pickup momentum transfer and a new source with loading system for lithium clusters

PDF

Dissociation energy and dynamics of HCl and aromatic-water clusters

PDF

Temperature-dependent photoionization and electron pairing in metal nanoclusters

Asset Metadata

Creator Markovskiy, Nikolay (author)

Core Title Nonlocal Monte Carlo approaches: helium clusters and nucleic acids

Contributor Electronically uploaded by the author (provenance)

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Chemistry (Chemical Physics)

Publication Date 12/02/2008

Defense Date 10/17/2008

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag CH4,helium clusters,methane,Monte Carlo,nucleic acids,OAI-PMH Harvest,rotational propagator,sampling,SF6,spherical top

Language English

Advisor Mak, Chi (committee chair), Kresin, Vitaly V. (committee member), Vilesov, Andrey (committee member)

Creator Email markovsk@usc.edu,markovskiy@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m1857

Unique identifier UC1188630

Identifier etd-Markovskiy-2501 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-130554 (legacy record id),usctheses-m1857 (legacy record id)

Legacy Identifier etd-Markovskiy-2501.pdf

Dmrecord 130554

Document Type Dissertation

Rights Markovskiy, Nikolay

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu