数理解析研究所講究録
第 1576 巻 2008 年 84-104

84

Dimension Reduction and Identification of Dynamic
Barriers in Structural Transitions of Clusters *Tomohiro Yanao
Fukui Institute for Fundamental Chemistry, Kyoto University, Kyoto 606-8103, Japane-mail: yanao@fukui.kyoto-u.ac.jp

Abstract

This report highlights the author’s recent study with W. S. Koon (Caltrh),
J. E. Marsden (Caltei), and I. G. Kevrekidis (Princeton University). Molecular
reactioo as well as functional motions of biopolymers are typically large-amplitude
collective motioo that involve alarge number of $degr\infty$ of heedom in acoherent
manner. It has been asignificant challenge to understand the mrhtism for sui
collective motions. By taking up asix-atom cluster as an illustrative example, wedevelope amethodology to understand $\infty 1lective$ motions of molecules with many
rad\"u, whii serve
degrees of&eedom in terms of the reduced dynamioe ofas predominant collective variables. We highlight the competition between two
radii: One
different kinds of forcae that aesentially mediate the dynamic8 ofthe ordinary potential foroe that originates&om the potential enery function of
the system, and the other is an internal centrifugal force that $ori\dot{g}nataebom$ theintrinsic non-Euclidean nature of molecular internal space. While the potential
force generally works to keep the mass distribution of the system compact andsymmetric, the internal centrifugal force works to inflate and elongate it. We showthat the internal centrifugal force induces asignificant dynamic barrier for reaction,
whii can often overshadow the $ori\dot{g}nal$ potential energy barrier. This hnd of$g\infty metry$ -related dynamical effect8 should shed new light on the conventional picture$\mathfrak{B}^{r}ration$$\mathfrak{H}^{r}ration$$i\epsilon$

of molecular reactioo.

1

Introduction

Large-amplitude collective motions play an essential role in chemical reactions as wellas in the functional motions of biomolecules. .Understanding the mechanism for suchcollective motions has been a significant challenge in current molecular science. Sincesuch collective motions involve a large number of degrees of freedom in a coherent manner,
reducing the dimensionality by using a small number of appropriate collective variables$\overline{*Based}$on the paper, T. Yanao, W. S. Koon, J. E. Marsden, and I. G. Kevrekidis, J. Chem. Phys.
126, 124102 (2007).

85

is crucially important. Indeed, dimensionality reduction has been along-standing issue
in nonequilibrium statistical mechanics [1-3], molecular dynamics [4-6], Monte Carlo
methods [7], and in reaction path theories [8-10].
The first issue of this report is to find out appropriate collective variables (or equivreaction coordinates) that play predominant roles rather ubiquitously in awide
class of large-amplitude motions of molecular systems. We specifically highlight the threegyration mdii of amolecule as such collective variablae in this study. $Physi_{Ca}g_{y}$ thegyration radii are the measure of mass distribution of asystem along the three principal axes. We investigate the three-dimensional dynanics of gyration radii of ageneralnumerically at the level of their classical equations-atom molecule both talyticallythree rotationalof motion. To this end, we first reduce the thrae translational
grees of freedom $bom$ the $3n$ total degrees of freedom configuration space correctly using
the bamework of standard reduction theory for symmetric systems foUowing $[11, 12]$ andthe gauge tfeory [13-19]. In order to carry out further reduction in dimensions, beyond
symmetry reduction, we employ the principal-axis hyperspherical formalism, which was
early paper by Eckart [20] and reformulated in the hyperspherical context
lnitiated in
Kuppermrn [23-25].
by Chapuisat et al. $[21, 22]$
formalism, the instantaneous principal axes ofhypersphericalprincipal-axis
In thethe system are used as abody frame. Then the $(3n-6)$ internal degrees of freedomthethe $(3n-9)$ hypertglae.
are parametrized by the three yration radii
gyration radii characterize the mass distribution of the system, hyperangles are associatedwith cyclic and democratic deformation of the system called kinematic rotation, whii
Aquilanti et al. [29-31]. It was shown in
was scrutinized by Littlejohn et al. [26-28]
Refs. [20-22] that the use of $(3n-9)$ quasivelocities instead of the hypertgles simplifiae
the expression of kinetic energy of -body systems remarkably. By taking the advrtageof this concise expression for kinetic energy, we here derive classical equations of motion,
physically appealing.
which are also concisethethree $\Re^{r}ration$ radii are regarded as slow collective variablae
present
study,In thethat essentially dominate large-amplitude motions of the system, while the hypertgularvariablae are regardd as fast “bath” modes. This distinction of variablae is based mainly
on time scale separation and kinetic energy partitioning as wif be shown in this report.
Based on this claesification of the variables, we average out the hypertgular variablae thatappear in the equations of motion for the yration radii to obtain aset of approximately
closed equations of motion for the yration radii.
The averaged equations of motion for the $ae^{r}ration$ radii shed light on the meianism of
motions
adynamic force in larcompetition between the potential force
of polyatomic molecul\’e. The dynamic force hae its origin in the kinematic coupling ofgyration radii with the hyperangular modes (kinematic rotations) via the non-Euclidet
intemal centrifirgal force
metric of the internal space. This dynamic force is essentiaUy
the potential force
speaking,
Generallyarising from the internal motions themselves.
compact. On
works to keep the internal maes distribution of amolecule symmetricelongatingtd
the other htd, the internal centrifugal force has the remarkable effect ofinflating the internal maes distribution of the molecule. In other words, the molecularvibrations induoe aspontteous tendency of deformation. In this way, the internal centrifugal force ct be the critical driving force for amolecule to move from one potential$alen\dot{t}ly$$\bm{t}d$$n$$\bm{t}d$$d\triangleright$$\bm{t}$$\bm{t}d$$Wl\dot{u}1e$$\bm{t}d$$\bm{t}d$$n$$\bm{t}d$$g\triangleright amplitude$$\bm{t}d$$\bm{t}$$\bm{t}d$86well to another. We show that this dynamic force induces asignificant dynamic barrier
often overshadow the original potential energy barrier. The effects
for reaction, which
four-atom reactions have been studiedof the internal centrifugal force in three-atom
ageneralized
present
report,wemethod of these results to larger
in Refs. [32-34]. In this-atom moleculae involving syvtematic averaging based on our recent publication [35].
As an iUustrative example, we take up the structural trtsition dynamics of asix-atom-like cluster [36-38], with constrt enery and withMorse cluster, which represents
zero total tgular momentum. The cluster possessae two kinds of geometrically distinct
isomers, one of which hae aregular octahedron structure, whii is highly symmetric,
the other hae aelongated structure. The octahedron structure lies at the bottomof amuch daeper potential $weU$ . Interestingly, this cluster shows aclear switching ofstructural preference depending on its total energy. In the low ener $gy$ range (solid-like
phase), the cluster spends more time in the potential well of the octahedron structure,
while in the high ener $gy$ rt $ge$ (liquid-like phase), the cluster spends more time in thecompeting potential well. This kind of two-state switching is ageneral feature of the firstproteins $[42, 43]$ . We
order phaee trtsitions as is observed in larger clusters [39-41]
will provide anovel expltation for such switching in terms of the competition between
the internal centrifugal force.
the potential force
This report is orgtized foUows. In Sec. 2, after introducing the model system, thesix-atom Morse cluster, we pr\’eent how the structural trtsition dynamioe of the cluster
is coarsely characterized by the yration radii. In Sec. 3, we introduce the principal-ris
radii
hyperspherical coordinatae,
derive general equations of motion for theof isolated -atom $system8$ . By averaging these equations of motion for the gyration rad\"u,the competition between akinematic force and the potential force is characterized in Sec.
4. Meianism of switiing of structural preference of the cluster is explained in terms of
adynanic barrier. The report concludae in Sec. 5with some remarks regarding futurestudies.
$C\bm{t}$$\bm{t}d$$n$$\bm{t}Ar_{6}$$\bm{t}d$$\bm{t}d$$\bm{t}d$$\epsilon s$$\bm{t}d$$\mathfrak{B}^{r}ration$$n$2Collective Variables in Structural Transition Dynamics of Clusters

2.1

Model cluster and residence probabilities

We first introduce our prototypal model for the study of molecular conformational transitions; an atomic cluster composed of six identical atoms that mutually interact throughcluster. The dimensionless Hamiltonian
the pairwise Morse potential, This is called anof the system is given by$M_{6}$$\frac{\mathcal{H}}{\epsilon}=\frac{1}{2}\sum_{:\approx 1}^{6}(\dot{r}_{\dot{\alpha}}\cdot\dot{r}_{\epsilon i})+\sum_{i\triangleleft}[e^{-2(d_{ij}-d_{0})}-2e^{-(\phi_{j}-d_{0})}]$,(1)$(i=1, \cdots , 6)$ is the three-dimensional position vector of the atom .
The
where
subscript represents the quantity with respect to the space-fixed frame. (This rule isrepresents the time derivative. An thealso applied to other quantities.) The dot over$i$$r_{\dot{n}}$$s$$r_{\epsilon i}$

87

(a)

$t1.8$$Saddle\Phi$$\wedge\omega-\{\{.\S$$7 \epsilon\frac{a}{\S}$$12$$\sim\S_{\Leftrightarrow}12.2$

12.4$\gamma^{/’\backslash }\backslash _{\backslash _{\backslash }\Phi}-CTBP$$\propto T/$ノ$\theta$$1.S- 1$

0.5

$0$

0.$

1

1.5

Energy (e)Arc-length

Figure 1: (a) Isomerization scheme of the six-atom Morse cluster. The cluster has two ge
ometrically distinct isomers, OCT and CTBP. The potential energy curve along the steepest
descent path connecting the saddle point and the two potential
ima is shown. (b) Energy
dependence of the residence probabilities of thecluster in the two isomers. Open squares
represent the residence probability for the OCT isomer, while the open triangles represent theresidence probability for the CTBP isomer.
$m\dot{i}$$M_{6}$

masses of the atoms are set to unity, $m_{1}=1(i=1, \cdots 6)$ . The paxameter represents$\epsilon$

the depth of the Morse potential
is the inter-particle disttce between atom td
atom . The parameter , whii corresponds to the equilibrium disttce of the $p\dot{u}rwise$
Morse potential, is set to $d_{0}=6.0$ . This provides apotential topography similar to thatof the Lennard-Jones potential [37]. Therefore, this system can be regarded as amodelof thecluster. Since the main interest of this study is in the internal dynanics ofpolyatomic systems, the total $\bm{t}1^{1ar}$ momentum of the system is assumed to be zerothroughout the report. In what follows, our numerical results are presented in absolute$\bm{t}dd_{ij}$

$j$

$i$

$d_{0}$$Ar_{6}$

units.
The isomerization (structural transition) scheme of thecluster is shown in Fig.
l(a). This cluster has two gmmetrically distinct isomers: One is the regular octahedron(OCT) and the other is the capped trigonal bipyramid (CTBP) [36-38]. The potential
energy minimum of the OCT isomer is $V=-12.49\epsilon,$that of the CTBP isomer is$V=-12.13\epsilon$ . Thaee two isomers are connected through asaddle point whose potential
ener$y$ is $V=-11.83\epsilon$ . In Fig. l(a), the potential energy topography along the steepaet$de8cent$ path is shown. The horizontal axis is the arc-len $gth$ of the path in the space of
yration radii, which will be introduced later. Note that the OCT isomer has ahighly
symmetric spherical structure with adeep potential well, while the CTBP isomer has
elongated (coUapsed) structure with ashallow potential $weU$ .
In this study, we employ amicroctonical, constant ener $y$ simulation. When the
total enery of the cluster is igher than that of the saddle point, the isomerization
reaction between the two isomers is energetically possible. Occurrence of the structural
isomerization reaction is detected by the queniing method [36], which solves the firstorder equations,, at each instrt until the system arrives at one of
the minima on the potential energy surface. The parameter is arbitrary. With thismethod, every point along aclassical trajectory is attributed to one of the two isomers$M_{6}$$\bm{t}d$$\bm{t}$$dr_{\epsilon i}/d\tau=-\partial V/\partial r_{\epsilon i}$$\tau$

88

except for the saddle points.
Fig. l(b) shows the energy dependence of the percentages of the total residenoe times
in the two isomer wefs obtained through along-time simulation. These are essentiallythe probabllIties of finding the system in the respective configurations. In the low energy
range, thecluster spends more time in the OCT isomer tht in the CTBP isomer.
This is understtdable from the view point of the potential ener $y$ topography since
the potential well of the OCT isomer is much deeper tht that of CTBP as shown in
Fig. 1(a). However, as energy increases, the structural preferenoe of the cluster switches
dramaticafy,
the system prefers to be in the CTBP isomer than in the OCT. In
this high energy range, the potential energy curve in Fig. l(a) is obviously not sufficientfor explaining residence probabilities. In this report, we will present the mechanism of
the switching of structural preference observed in Fig. l(b) in terms of the effect of a
dynamical and geometrical force that often counteracts with the potential force.
$M_{6}$$\bm{t}d$2.2

Coarse dynamics of gyration radii

Here, we coarsely characterize the structural transition dynamics of thecluster using
the three gyration radii, our collective variables. For generality, we consider an n-atom
molecule whose constituent atoms have masses $m_{i}(i=1, \cdots n)$ . After eliminating theoverall translational degrees of freedom of the molecule via the mass-weighted Jacobi
vectors,$M_{6}$$\rho_{\iota i}=\sqrt{\frac{m_{1+1}\sum_{k--1}^{1}m_{k}}{\sum_{k=1}^{i+1}m_{k}}}(\frac{\sum_{k--1}^{i}m_{k}r_{l}k}{\sum_{k=1}^{i}m_{k}}-r_{\iota(i+1))}(i=1, \cdots, n-1)$,(2)

the singular-value decomposition theorem [44] can be applied to decompose the $3x(n-1)-$
dimensional matrix
into the product of the three matrices$W_{s}\equiv(\rho_{\epsilon 1}\cdots\rho_{\epsilon(n-1)})$$W_{*}=RN\cup^{T}$,(3)

where $R=(e_{1}e_{2}e_{3})$ , and
$N$

$=$

$\cup$$=$($a_{0}o^{1}$$a_{2}00$$a_{3}00|_{0}^{0}0\ldots 000)\equiv($$(u_{1}u_{2}u_{3}|u_{4}\cdots u_{n-1})\equiv(P|Q)$

A $|0)$ ,

(4)

.

(5)

The symbol $T$ on the matrix in Eq. (3) indicates the matrix trtspose. The matrix, are orthogonal td
is a3 $x3$ orthogonal matrix $(SO(3))$ , whose column vectors,
normalized. The matrix is a $3\cross(n-1)$ diagonal matrix, which is split up
the leftthe right 3 $x(n-4)$ zero matrix $0$ for later use. The diagonal
3 $x3$ diagonal matrix A, are the gyration radii, which are the key qurtities
elements (singular values),of this study. Although there are mry alternative definitions for the Jacobi vectors,
gyration radii are independent of the choice of Jacobi vectors. They are non-negative tdordered according to $a_{1}\geq a_{2}\geq a_{3}\geq 0$ . The matrix is $\bm{t}(n-1)x(n-1)$ orthogonal$\cdots u_{n-1}$ , are orthogonal
matrix $(SO(n-1))$ , whose column vectors,
normalized.
$U$

$R$

$e_{1},$ $e_{2},$ $e_{3}$$N$$\bm{i}to$$\bm{t}d$$a_{1},a_{2},$ $a_{3}$$U$$u_{1},$$\bm{t}d$89

Figure 2; (a) Typical time evolution of the three gyration radii
, and
of
the
cluster at total energy $E=-11.0\epsilon$ . The cluster quenches to the OCT isomer $bomt=0$
to $t=468$ and $bomt=884$ to $t=$ 1200, while the cluster is quenched to the CTBP isomer
ffom $t=468$ to $t=884$ . $(b)$ The same trajectory as in (a) projected onto the three-dimensional
space of gyration radii. The two trapping regions correspond to the two isomers, OCT (solidline) and CTBP (dashed line).
$a_{1},$$a_{2}$$a_{3}(a_{1}\geq a_{2}\geq a_{3})$$M_{6}$The matrix is split up into the left $(n-1)\cross 3$ matrix
the right $(n-1)x(n-4)$matrix for later convenience.
The physical meaning of the decomposition Eq. (3) is as follows. It is shown thatthe matrix coincides with the principal axis hame, which specifiae the $instant\epsilon n\bm{m}us$orientation of the system, while and determine the internal structure (sizesymmetry) of the system. Gyration radii,, are the measure of the mass-weighted
lengh of the system along the respective principal axes. They essentially characterize the mass distribution, and are related to the principal moments of inertia,,
and $M_{3}(M_{3}\geq M_{2}\geq M_{1}\geq 0)$ , by,,.
The values of the yration radii of thecluster at the potential minimum of the
OCT isomer are $(a_{1}, a_{2}, a_{3})=(5.97,5.97,5.97)$ , while those of the CTBP isomer are$(a_{1}, a_{2}, a_{3})=(8.19,5.89,4.27)$ . The coincidence of the three yration radii at the OCT
structure Indicates that it has isotropic mass distribution (spherical top). The CTBP
structure has an aeymmetric maes distribution.
Fig. 2(a) shows atypical time evolution of the threerffiii of thecluster at
total ener $yE=-11.0\epsilon$ . In the time period of Fig. 2(a), the system first transited $bom$
the OCT isomer to the CTBP isomer. After awhile it returned back to the OCT isomer.
The isomerization reaction is clearly marked by the change of the three
radii.
In eai isomer, the three yration radii take values close to those at the corresponding
potential minimum. Fig. 2(b) shows the same trajectory as in Fig. 2(a) projected onto thethree-dimensional space of $ae^{r}ration$ radii, which captures the chtges of mass distributionof the syst $em$ . In Fig. 2(b), the two regions corresponding to the OCT isomer
the
CTBP isomer are distinguished. Both of the two regions are small when the total ener $y$
is low. As the total energy increasae, the two regions become larger. The essential fact
is that certain regions of the three-dimensional space of yration radii
certain isomerwells of the cluster have good one-to-one correspondence. Therefore, the trtsition of the$U$$P_{\bm{t}}d$

$Q$

$R$

$N$

$U$

$\bm{t}d$$a_{1},$ $a_{2},$ $a_{3}$$M_{1},$ $M_{2}$$M_{1}=a_{2}^{2}+a_{3}^{2}$$M_{2}=a_{3}^{2}+a_{1}^{2}$$M_{3}=a_{1}^{2}+a_{2}^{2}$$M_{6}$$M_{6}$$\mathfrak{H}^{r}ration$$\mathfrak{B}^{r}ration$$\bm{t}d$$\bm{t}d$90

trajectory in the three-dimensional space of gyratim radii from one region to the otheris indicative of astructural chtge of the full system. This is the most vital condition
for the three $ae^{r}ration$ radii to be good collective variables. In general, gyration radii cancharacterize structural transitions of amolecule in this way as long as the system chtges
its mass distribution significantly through the trtsitions.
The number of total lnternal degrees of freedom of thecluster is twelve. Therefore
nine other internal degrees of beedom are “hidden” behind the dynamics in Fig. 2(b).
If these
degrees of freedom ct be regarded as “bath” modes, one ct think ofaclosed dynanical system in terms of only the three gyration radii. If this is the case,one will be able to extract essential information about the structural transition of thefull system from the low-dimensional dynamics of gyration radii. Thaee are the main
issues of this report in the foUowing sections, where we will finally extracteffective
barrier structure that separatae the different isomer regions as in Fig. 2(b). Thaee barrier
structures will in turn explain the true dynamical stability of the respective isomers.
$M_{6}$$h\ddagger dden$$\bm{t}$3

Equations of Motion for the Collective Variables

3.1

Kinetic energy of an n-atom molecule in the hypersphericalcoordinates

The goal of this section is to derive classical equations of motion for the $a^{r}ration$ radii, ourcollective variables, using the pricipal-axis hyperspherical coordinatae. For this purpose,
we first present the general expression for the total kinetic energy of-atom molecule,such as Eq. (27)
Eq. (29), in this subsection.
Let amatrix $R\in SO(3)$ be abody (body-fixed) frame. Since thevectors withrespect to the body hame $\rho_{i}(i=1, \cdots n-1)$ are related to the Jwobi vectors with, the time-derivative
respect to theframe $\rho_{\epsilon i}(i=1, \cdots n-1)$ byof
ct be expraesed as.
(6)
Similarly, the total angular momentum with respect to the spacefixed frame
and thatwith respect to the body frame $L$ are related by $L,$ $=RL$ . Since$\bm{t}n$$\bm{t}d$$J\infty obi$

$spac\triangleright fixed$$\rho_{\dot{n}}=R\rho_{i}$$\rho_{\iota i}$$\dot{\rho}_{\ell\dot{*}}=\dot{R}\rho_{i}+R\dot{\rho}_{i}$$L_{f}$

$L$

$L$

.

$= \sum_{1=1}^{n-1}\rho_{\ell 1}x\dot{\rho}_{\epsilon i}$,(7)can be expressed as.

$L= R^{T}\sum_{:-1}^{n-1}\rho_{\dot{\alpha}}x\dot{\rho}_{\ell 1}=M\Omega+\sum_{i=1}^{n-1}\rho_{i}x\dot{\rho}_{1}$(8)The matrix in Eq. (8) is the moment of inertia tensor with respect to the body frame,
whose components are defined by$M$$M_{\alpha} \rho=\sum_{1=1}^{n-1}[(\rho_{i}\cdot\rho_{i})\delta_{\alpha\beta}-\rho_{i\alpha}\rho_{i\beta}]$

,

(9)

91

is theand $\beta(\alpha, \beta=1,2,3)$ represent the axes of the body frame andwhere
Kronecker delta. The vectorin Eq. (8) is the angular velocity withrespect to the body frame whose components are defined through the anti-symmetric$\delta_{a\beta}$$\alpha$$\Omega\equiv(\Omega_{1}, \Omega_{2}, \Omega_{3})^{T}$

matrix,$R^{T}\dot{R}\equiv(\begin{array}{lll}0 -\Omega_{3} \Omega_{2}\Omega_{3} 0 -\Omega_{l}-\Omega_{2} \Omega_{l} 0\end{array})$

.

(10)

The kinetic energy of the system in a space fixed frame is given by$K= \frac{1}{2}\sum_{i=1}^{n}m_{i}$(箇 .

$\dot{r}_{\epsilon i}$)$= \frac{1}{2}\sum_{i=1}^{n-1}\dot{\rho}_{\dot{n}}$.

$\dot{\rho}_{\iota*}\cdot$

.

(11)

After applying Eq. (6) to Eq. (11), it can be expressed in terms of the body-frame
quantities as$K= \frac{1}{2}\Omega^{T}M\Omega+\Omega^{T}\sum_{i=1}^{n-1}\rho_{i}x\dot{\rho}_{i}+\frac{1}{2}\sum_{:<1}^{n-1}\dot{\rho}_{i}\cdot\dot{\rho}_{i}$

.

(12)

Now we reformulate the kinetic energy of Eq. (12) in terms of the principal-axis
hyperspherical coordinates. In what follows, we identify the body ffame in the above
argument with the principal-axis bame introduced in Eq. (3). Correspondingly, theangular velocity vector , the moment of inertia tensor , and the Jacobi vectors
in Eq. (12) are all regarded as the ones with respect to the $principal-\alpha is$ frame.
The second $term$ of Eq. (12) can be reformulated as follows. Based on Eq. (3), theJacobi vectors
with respect to the principal-axis frame are expressed collectively ina matrix fom as
$R$

$R$

$M$

$\Omega$$\{\rho:\}$$\{\rho_{i}\}$$(\rho_{1}\cdots\rho_{n-1})=N\cup^{T}\equiv W$

.

(13)

wher the 3 $x(n-1)$ matrix $W$ is a body-frame counterpart of$W$ ,the components ofare given by$e$$W_{e}$

of Eq. (3). By using$\sum_{i=1}^{n-1}\rho_{i}x\dot{\rho}_{1}\equiv(\xi_{1},\xi_{2}, \xi_{3})^{T}$$(\begin{array}{lll}0 -\xi_{3} \xi_{2}\xi_{3} 0 -\xi_{l}-\xi_{2} \xi_{l} 0\end{array})=\dot{w}w^{\tau_{-W}\tau_{=-2AP^{T}\dot{P}A}}$

,

(14)

where A and are the sub-matrices defined in Eq. (4) and in Eq. (5). After introducing
quasivelocity components$P$$(i,j=1,2,3, i\neq j)$ ,: $u_{j}=-u$: .
the anti-symmetric matrix
can be expressed as$\dot{u}$(15)$\dot{u}_{j}\equiv\omega_{1j}$$P^{T}\dot{P}$$P^{T}\dot{P}=(\begin{array}{lll}0 -\omega_{l2} \omega_{3l}\omega_{l2} 0 -w_{23}-w_{31} \omega_{23} 0\end{array})$

.

(16)

By substituting Eq. (16) into Eq. (14) and comparing the leftmost part with the rightmost
part of Eq. (14), we obtain$\xi_{1}=-2a_{2}a_{3}w_{23},$$\xi_{2}=-2a_{3}a_{1}w_{31},$$\xi_{3}=-2a_{1}a_{2}\omega_{12}$

.

(17)

92

These formulas allow us to express the second term on the right hand side of Eq. (12)finally as$\Omega^{T}\sum_{:=1}^{n-1}\rho_{i}x\dot{\rho}_{i}$(18)$-2a_{2}a_{3}\Omega_{1}\omega_{23}-2a_{3}a_{1}\Omega_{2}\omega_{31}-2a_{1}a_{2}\Omega_{3}\omega_{12}$$=$$\Omega^{T}Bw$

$=$

,

(19)

where the following two matrices are introduced for notational compactness:
一

$B=($

$-2a_{3}a_{1}00$$2a_{2}a_{3}00$$-2a_{1}a_{2}00$),$\omega=(\begin{array}{l}\omega_{23}\omega_{31}w_{l2}\end{array})$

.

(20)

Similarly, the third term in the right hand side of Eq. (12) can also be reformulatedusing the matrix $W$ of Eq. (13) as$\frac{1}{2}\sum_{:=1}^{n-1}\dot{\rho}_{1}\cdot\dot{\rho}_{1}$$=$$\frac{1}{2}b(\dot{W}\dot{W}^{T})$$=$$\frac{1}{2}R(\dot{A}\dot{A}^{T})+\frac{1}{2}$

(21)

‘hr

$( A\dot{P}^{T}PP^{T}\dot{P}A)+\frac{1}{2}$Tt$(A\dot{P}^{T}QQ^{T}\dot{P}A)$

,

(22)

, and
are the matrices
represents the trace of a matrix and where
where
components,
quasivelocity
defined in Eq. (4) and in Eq. (5). After introducing the$A,$$R(\cdot)$$\dot{u}_{i}$. $u_{j}=-u_{i}$ .

the $(n-4)x3$ matrix$Q^{T}\dot{P}$$\dot{u}_{j}\equiv\gamma_{ij}$

$P$

$Q$

$(i=1,2,3, j=4, \cdots n-1)$ ,(23)

in Eq. (22) can be expressed as$Q^{T}\dot{P}=(\begin{array}{lll}\gamma_{l4} \gamma u \gamma_{34}\vdots \vdots \vdots\gamma_{1(n-1)} \gamma_{2(n-l)} \gamma_{3(n-l)}\end{array})\equiv(\gamma_{1}\gamma_{2}\gamma_{3})$

,

(24)

where we have introduced the $(n-4)$ -dimensional column vectors $\gamma_{i}(i=1,2,3)$ . By
applying Eq. (16) and Eq. (24) to Eq. (22), we obtain the final expression for the third
term of Eq. (12) as$\frac{1}{2}\sum_{i=1}^{n-1}i$

$i$

$=$

$\frac{1}{2}(\dot{a}_{1}^{2}+\dot{a}_{2}^{2}+\dot{a}_{3}^{2})+\frac{1}{2}(M_{1}\omega_{23}^{2}+M_{2}w_{31}^{2}+M_{3}w_{12}^{2})$$+ \frac{1}{2}a_{1}^{2}\sum_{k=4}^{n-1}\gamma_{1k}^{2}+\frac{1}{2}a_{2}^{2}\sum_{k=4}^{n-1}\gamma_{2k}^{2}+\frac{1}{2}a_{3}^{2}\sum_{k=4}^{n-1}\gamma_{3k}^{2}$$\frac{1}{2}\dot{a}\cdot\dot{a}+\frac{1}{2}\omega^{T}Mw+\frac{1}{2}\sum_{i=1}^{3}\gamma_{i}^{T}a_{1}^{2}\gamma_{i}$where $a\equiv(a_{1},a_{2}, a_{3})^{T}$ , andprincipal-axis frame.

$M$

,

,

(25)
(26)

is the moment of inertia tensor with respect to the93

Substituting Eq. (19) and Eq. (26) into Eq. (12), we obtain the desired general
expression for the total kinetic energy of an n-atom molecule as$K= \frac{1}{2}\Omega^{T}M\Omega+\Omega^{T}B\omega+\frac{1}{2}\dot{a}\cdot\dot{a}+\frac{1}{2}w^{T}Mw+\frac{1}{2}\sum_{:=1}^{3}\gamma_{i}^{T}a_{i}^{2}\gamma_{i}$

,

(27)

which is essentially the same as the one given by Eckart [20] and by Chapuisat et al.
[21]. Introduction of the quasi-velocities,
and , has been crucial in making this acompact expression; in it, the $(3n-6)$ internal degrees of freedom are represented by
the three gyration radii and the $(3n-9)$ quasivelocity components. Thes quasivelocities are essentially the “angular velocities” of kinematic (or democratic) rotations, which
are continuous shape changes associated with the permutations (relabellings) among theconstituent atoms [23-28]. They should be distinguished from the angular velocity of theprincipal-axis frame .
$w$$\gamma_{i}$$e$$\Omega$

3.2

Equations of motion for gyration radii under the conditionsof vanishing total angular momentum

We next investigate the equations of motion for the three yration radii,, and ,
under the conditions of zero total angular momentum. Since the total angular momentum
with respect to the body frame is given by$a_{1},$$L=M\Omega+Bw$ ,$a_{2}$$a_{3}$(28)

after using Eq. (8) and Eq. (19), the total kinetic energy can be expressed as follows
after applying Eq. (28) to Eq. (27),
,

(29)

).

(30)

$K= \frac{1}{2}L^{T}M^{-1}L+\frac{1}{2}\dot{a}\cdot\dot{a}+\frac{1}{2}w^{T}\tilde{M}w+\frac{1}{2}\sum_{1=1}^{3}\gamma_{i}^{T}a_{1}^{2}\gamma_{i}$where we have introduced a diagonal matrix,$\frac{(a_{2}^{2}-a_{\theta}^{2})^{2}}{(a_{0}^{2}2^{+a\S)}0}$$\overline{M}\equiv M-BM^{-1}B=($$\frac{(a_{\theta}^{2}-a_{1}^{2})^{2}0}{(a_{\theta}^{2}+a_{1}^{2}),0}$$\frac{(a_{1}^{2}-a_{2}^{2)^{2}}0}{(a_{1}^{2}+a_{2}^{2})}0$

The expression of Eq. (29) is important from the gauge-theoretical point of view [19]:Each term of Eq. (29) is independent of the choice of body frame, while the kinetic energy
expression of Eq. (27) is specific to the principal-axis choice of body frame. By applying
the zero angular momentum condition, $L=0$ , to Eq. (29), we obtain the purely internal

94

kinetic energy,
$K$

,

(31)

$=$

$\frac{1}{2}\dot{a}\cdot\dot{a}+\frac{1}{2}w^{T}\overline{M}w+\frac{1}{2}\sum_{i=1}^{3}\gamma_{l}^{T}a_{i}^{2}\gamma_{i}$$=$$\frac{1}{2}(\dot{a}_{1}^{2}+\dot{a}_{2}^{2}+\dot{a}_{3}^{2})+\frac{(a_{1}^{2}-a_{2}^{2})^{2}}{2(a_{1}^{2}+a_{2}^{2})}w_{12}^{2}+\frac{(a_{2}^{2}-a_{S}^{2})^{2}}{2(a_{2}^{2}+a_{3}^{2})}\omega_{23}^{2}+\frac{(a_{3}^{2}-a_{1}^{2})^{2}}{2(a_{3}^{2}+a_{1}^{2})}w_{31}^{2}$$+ \frac{1}{2}a_{1}^{2}\sum_{k=4}^{n-1}\gamma_{1k}^{2}+\frac{1}{2}a_{2}^{2}\sum_{k=4}^{n-1}\gamma_{2k}^{2}+\frac{1}{2}a_{3}^{2}\sum_{k=4}^{n-1}\gamma_{3k}^{2}$. (32)

metric components fortd
this point, it is important to note that the-dimensional space of yration radii
are all equal to unity. This means that theas Fig. 2(b) is Euclidean even though the $(3n-6)$ -dimensional full internal space
is non-Euclidean. This $wiU$ be of great advantage in the $foUowing$ discussions in reducing-dimensional space of gyration radii.
the full dynamics of amolecule to the
By using the intemal kinetic energy Eq. (32), the Lagrtgit for $the*atom$ moleculewith vtishing total angular momentum is given by $\mathcal{L}=K-V$ . Recall that for our studyof isolated moleculae the potential energy is afunction only of the internal degrees ofrad\"u can be written as
heedom. The Euler-Lagrange equations for the$ma\epsilon s$$At2$$\dot{a}_{1}^{2},\dot{a}_{2)}^{2}$$thr\infty$$a_{3}$$su\bm{i}$$thr\infty$$V$$\mathfrak{B}^{r}ration$$\ddot{a}_{1}=\frac{a_{1}(a_{1}^{2}+3a_{2}^{2})(a_{1}^{2}-a_{2}^{2})}{(a_{1}^{2}+a_{2}^{2})^{2}}w_{12}^{2}+\frac{a_{1}(a_{1}^{2}+3a_{3}^{2})(a_{1}^{2}-a_{3}^{2})}{(a_{1}^{2}+a_{3}^{2})^{2}}w_{13}^{2}+a_{1}\sum_{k=4}^{n- 1}\gamma_{1k}^{2}-\frac{\partial V}{\partial a_{1}},(33)$$\ddot{a}_{2}=\frac{a_{2}(a_{2}^{2}+3a_{1}^{2})(a_{2}^{2}-a_{1}^{2})}{(a_{2}^{2}+a_{1}^{2})^{2}}w_{21}^{2}+\frac{a_{2}(a_{2}^{2}+3a_{3}^{2})(a_{2}^{2}-a_{3}^{2})}{(a_{2}^{2}+a_{3}^{2})^{2}}w$襲$+a_{2} \sum_{k=4}^{\mathfrak{n}-1}\gamma_{2k}^{2}-\frac{\partial V}{\partial a_{2}},(34)$$\ddot{a}_{3}=\frac{a_{3}(a_{l}^{2}+3a_{2}^{2})(a_{3}^{2}-a_{2}^{2})}{(a_{3}^{2}+a_{2}^{2})^{2}}w_{32}^{2}+\frac{a_{3}(a_{3}^{2}+3a_{1}^{2})(a_{3}^{2}-a_{1}^{2})}{(a_{3}^{2}+a_{1}^{2})^{2}}w_{31}^{2}+a_{3}\sum_{k=4}^{n-1}\gamma_{3k}^{2}-\frac{\partial V}{\partial a_{3}}.(35)$The left htd sides of these equations are the components of acceleration in the thraedimensional space of yration radii. The right htd sides are physically the forces inthis space. While the fourth terms on the right htd sides of Eqs. (33)-(35)represent
the force that originates from the potential energy function, the first three terms on theright htd sides of these equations represent adynamic (kinematic) fome that originates$bom$ the dynanical coupling of gyration radii with the kinematic-rotation $modes$ . All.
of these dynamic force terms are quadratic in the quasivelocity components,$bom$$es8entially$
therepreaent the intemal
Therefore, these terms
force arisingkinematic rotations. It has been shown that the internal centrifugal terms proportionalarise in the three- and four-atom dynamics [32]. But the internal centrifugal forc
toaxe intrinsic to systems with more tht four-atom systems interms proportional to
the thr $ee$-dimensional physical space. As the number of atoms increases, the number ofincreases.
the terms proportional to
While the force arising from the potential function is dependent on the system, theinternal centrifugal force terms in Eqs. (33)-(35)are common to the dynamics of general-atom moleculae. Therefore it is quite interaeting to explore the essential properties$w_{ij}\bm{t}d\gamma_{ik}$$cent\dot{n}fi\nu gal$$e$$\omega_{ij}^{2}$$\gamma_{1k}^{2}$$\gamma_{ik}^{2}$

$n$

95

of this kinematic force. The internal centrifugal force terms proportional tohave
an effect of breAing the symmetry of mass distributIon of the -atom molecule. This
property can be understood by noting the sign of the corresponding terms in Eqs. (33)(35): Since $a_{1}\geq a_{2}\geq a_{3}\geq 0$ holds by definition, the firstthe second terms in Eq.
(33) are always positive or zero. This means that these terms have apersist$ent$ effect ofenlarging the largoet yration radius . On the other hand, by similar arguments, thefirst and the second terms in Eq. (35) are always negative or zero. Thaee terms wif
therefore persistently diminish the smallest gyration radius . Finafy, the first termof Eq. (34) is negative or zero and the second term of Eq. (34) is positive or zero.
Therefore, the force ct act both posltively
negatively for . SInoe the yration
radii are the meaeure of the mass-weighted size of the system along the principal axae,
these properties of the internal centrifugal foroe indicate that amolecule is always
to elongate in the most massive dioection and to collapse in the leut massive dioection.
The terms proportional toin Eqs. (33)-(35)are always positive. Therefore, thaee havean effect of inflating the molecule in all three principal axis directions. It is expectd thatthe indating effect becomae greater ae the number of atoms increase since the number of
the terms proportional toincreasae.
The internal centrifugal force discussed above have similaritiae to those of the normal
centrifugal force. However the origins of the internal centrifugal force
the normal centrifugal force are different. The normal centrifugal force is associated with the non-zero
total $an\infty ar$ momentum,is thereby absent in the system with zero total angularmomentum. On the other htd, the internal centrifugal force of this study is associatedwith the kinematic rotations
ct arise even in the dynamics with zero angular momentum. That is, the internal centrifugal force is essentially adynamical force inducedby molecular vibrations themselves. In the next section, we scrutinize the competition
between this kind offorce
force.
the ordinary$\omega_{ij}^{2}$$n$$\bm{t}d$$a_{1}$$a_{3}$$\bm{t}d$$a_{2}$$for\mathfrak{r}ed$$\gamma_{1j}^{2}$$\gamma_{1j}^{2}$$\bm{t}d$$\bm{t}d$$\bm{t}d$$dyna\iota\dot{m}cal$$potent\ddagger al$$\bm{t}d$4Dynamics of Collective Variables: Competition between Dynamic Force and Potential Force

4.1

Distinction between collective variables and “bath” modes

In this section, we scrutinize the low dimensional dynanics of yration radii on the basisof the equations of motion presented in the previous section and numerical experiments.
Competition between the internal centrifugal force
the potential force is highlighted.
We first
the distinction between coUective variablae and “bath” modae. Fig. 3$(b)$ the squarae ofshows the typical time evolution of (a) the yration radii
the quasivelocities
at total energy $E=-10.0\epsilon$ . During
the time interval of Fig. 3, the cluster had undergone isomerization from the OCT isomer
to the CTBP isomer once. In Fig. 3, we, andthat the three yration radii,, itge slowly
smoothly, in contrast to the internal quasivelocity components,, whii osciUate rapidly and sharply. This result supports the tsatz that theradii ct be considered as slow variablae while other internal modes behave like“bath” modes. The rapid oscillations of the quasivelocities induce the rapid osciUations$\bm{t}d$$clari\mathfrak{h}^{\gamma}$$a_{1},$$\omega_{ij}^{2},$$a_{3}$$\bm{t}d\gamma_{1}^{2}$$\mathfrak{B}^{r}ration$$\bm{t}da_{3},$$(c)\gamma_{i}^{2}\equiv\gamma_{u}^{2}+\gamma_{\dot{\iota}5}^{2}(i=1,2,3)$$s\infty$$\bm{t}d$$a_{2},$$a_{1},$$a_{2}$$w_{1j}^{2}$

96

Time

Time

Time

$(b)$ the squares of, and
Figure 3: Typical time evolution of (a) the three gyration radii
at total energy $E=-10.0e$ . The system quenches
the quasivelocities,
to OCT from $t=0$ to $t=21.5$ , and to CTBP from $t=21.5$ to $t=50$ .
$a_{1},$$w_{12}^{2},\omega_{23}^{2},w_{31}^{2},$$(c)\gamma_{1}^{2},\gamma_{2}^{2},$$a_{2}$$a_{3},$$\gamma_{3}^{2}$of the internal centrifugal force in the dynamics of $ae^{r}ration$ radii. This $re8ult$ justifies
our strategy of averaging the rapid oscillations of the forces in Eqs. (33)-(35)to deduce
effective force fields that the gyration radii “feel” as will be shown in the next subsection.
We next praeent tother dynamical evidence that gyration radii play apredominrtrole over other hypertgular variablae in the large-amplitude motions of the cluster. Opensquaroe in Fig. 4(a) represent average kinetic energies in the 12 internal modes of thecluster when the cluster is in the OCT isomer. Here, each term in Eq. (32) represents
the kinetic energy in each $mode$ . Open circlae represent the similar average kinetic enery
distribution when the system is in the CTBP isomer. The total internal ener$y$ of thecluster is $E=-10.0\epsilon$ . It is evident that the kinetic energy is partitioned equaUy among
CTBP isomers. That is, kinetic ener $gy$ is
the 12 internal modes in both the OCT
equipaxtitioned among all the internal modes when the cluster is in each of the two isomer.
On the other hand, Fig. 4(b) shows the kinetic enery distribution among internal modes
averaged over oeactive tmjectories. As opposed to Fig. 4(a), we see that kinetic ener $gy$
modes have markedly (about $1.5\sim 2$ times)
is not equipartitioned, but themore kinetic energy tht other internal modes. Note that the kinetic ener $y$ distributlon
is $almo8t$ the same for the reaction from OCT to CTBP and the reverse reaction bom
CTBP to OCT. This is due to time reversal symmetry of the system. The result of Fig.
must acquire more hnetic ener $y$4(b) indicates that the two gyration radii
than other internal degrees of freedom for both the forward and ba&ward reaction. This
clearly shows that the gyration radii play apredomlnrt role over other internal varlablaein the large-amplitude motion.
$M_{6}$$\bm{t}d$$a_{1}\bm{t}da_{3}$$a_{1}\bm{t}da_{3}$4.2

Characterization of the potential force and the internal centrifugal force via averaging

Based on the results in the previous subsection, we regard the three yration radii as
predominant collective variables (reaction coordinates), and all other hyperangular modes
are regarded as “bath” modes in this study. For the understanding of the mechanismof structural transition of the cluster, the next important step is to characterize thecompetition between the internal centrifugal force and the potential foroe in the dynamics

97

Intemal modes

Intemal modes

Figure 4: (a) Open squares represent the average kinetic energies distributed into the 12 internal
modes of thecluster when the system is in the OCT isomer. Open circles represent the
average kinetic energies distributed among the 12 internal modes when the system is in the
CTBP conformation. (b) Average kinetic energies partitioned by the 12 internal modes overreactive trajectories. Open squares correspond to the reaction &om OCT to CTBP, and open
circles correspond to the reaction from CTBP to OCT.
$M_{6}$of yration radii. We here deduoe
averaged field of the internal centrifugal foroe td
that of the potential force numerically along aone-dimensional reaction path introducedin the space of $ae^{r}ration$ radii. The reaction path is determined by directly averaging the$traj\propto tories$ in the space of gyration
radii as in Fig. 2(b) in the following way:First, thetrajectories are averaged in terms of thecomponents with
fixed to different
values. We then obtaln representative (averaged) points for respective
values. Thaee
representative points are connected successively to obtain asingle path.
The internal centrifugal foroe
the potential force can now be iaracterized along
the averaged $re$ action path. Since theve two forces oscillate rapidly ae compared withthe time scale of variation of yration radii as aresult of the rapid oscillations of thehypertgular variablae (see Fig. 3and Ref. [35]), it is reasonable to average out thedependence of these forces in the dynanics of gyration radii. In Fig. 5, shown with arrows
are (a) the averaged field of the internal centrifugal force (the sum of the first three termson the right hand sides of Eqs. (33) $-(35)),$ $\bm{t}d(b)$ that of the potential force (the fourthterms on the right hand sides of Eqs. (33) $-(35))$ along the reaction path at total $ener_{\Psi}$$E=-10.0\epsilon$ . In the figures, the averaged
reaction path is shown by the thin broken curve.
The averaged field of internal centrifugal foroe is roughly direct $ed$ from the region of the
OCT isomer to that of the CTBP isomer. ffirthermore, the internal centrifugal force Is
stronger in the OCT region tht in the CTBP region. This indicates that the internal
centrifugal force strongly “pushes” the system from the OCT region to the CTBP region.
On the other hand, the averaged potential foroe field along the reaction path in Fig. 5(b)exhibits aitge in direction, roughly pointing, locally, to both points corresponding to
the potential minima of the OCT
the CTBP isomer. In other words, the potential
force is directed towards the minima of the potential wefs in the reduced space of yration
radii as is naturally expected. Thus, the averaged otentIal force works to keep the internalmaes distribution of the cluster close to that of each potential minimum structure, while
the averaged internal centrifugal force has apersistent tendency to bring the system mass$\bm{t}$$a_{2}\bm{t}da_{3}$$a_{1}$$a_{1}$$\bm{t}d$$tim\triangleright$$\bm{t}d$

$p$

98

Figure 5: Averaged field of the internal centrifugal foroe (Panel $(a)$ ) and that of the potential
foroe (Panel $(b)$ ) along the reaction path at $E=-10.0\epsilon$ . The thin broken lines in the two figures

constitute the averaged reaction path at $E=-10.0\epsilon$ . The arrows represent the directionality
and the magnitude of respective forces at each point along the reaction path. The thick solid
curve is the projection of the steepest descent path onto the space for comparison.

distribution from that of the OCT isomer to that of the CTBP isomer.

4.3

Characterization of a dynamic barrier for reaction

The results in Fig. 5(a) indicate that the internal centrifugal force ct be adynamic
forforce for the $re$action from OCT to CTBP, while it ct be adynamic
extra amount of “work”
the reaction from CTBP to OCT. This, in turn, indicates thatradii in addition toyration
chtge
its
to
vibratingcluster
may be needed in order for the$quanti6^{r}$thaee two works
the work necessary to overcome the potential forces. Here, we
path.
This procedure
along
the
reactionforcestwooftheby considering the line integrals$ent$
switiing of the favored mass
gives apossible rationalization for the enery-depend
l(b).
Fig.
in
cluster presented
distribution of the
In Fig. 6, the upper three ptels show the probability distributions along the reaction
paths for different energies, (a) $E=-11.2\epsilon,$ $(b)E=-10.0\epsilon$ , td(c) $E=-8.5\epsilon$ . The
horizontal $ris$ of eai of these flgures is the arc-length along each reaction path in thespace of $ae^{r}ration$ radil, whii serves ae the $re$action coordinate. The left part (the shortarc-length part) of ach $fi_{1}re$ corresponds to the region of the OCT isomer, while theright part corresponds to the region of the CTBP isomer. Since the reaction path becomae
from (a) to (c). In
longer as the total enery increaees, the width of these ptels
distribution,
theleft one of which
each panel, two pe&s are observed in the probability
corresponds to the OCT isomer while the right one corresponds to the CTBP isomer. At
low total energy (Fig. $6(a)$ ), the left $pe$ ak is higher tht the right peak. This indicatesthat the system prefers to be in the mass distribution of the OCT isomer tht that of
CTBP. On the other hand, as the total enery increaees (Fig. 6(b) td $(c)$ ), the left
$peak$ is lowered and the right peak becomae higher. In other words, the preferable mass
distribution of the cluster switches from that of the OCT isomer to that of the CTBP
$ban\dot{\tau}er$$d_{7}\dot{v}\dot{m}ng$$\bm{t}$$M_{6}$$e$$inCre\mathfrak{B}eS$99

represent the probability distribution along the reaction
coordinates introduced at total energy (a) $E=-11.2\epsilon,$ $(b)E=-10.0\epsilon$ , and (c) $E=-8.5\epsilon$(dotted
The lower three panels represent the $\infty rraeponding$ reconstructed potential energyline), and the superposed potential $V_{pot+cent}$ ( $8olid$ line) along the reaction coordinates.

Figure 6: The upper three$pane1_{8}$.

$V_{pot}$

isomer. In this way, the switching of the structural preference observed in Fig. l(b) ischaracterized in terms of the reaction paths. This kind of switching in the probability
distribution is analogous to $first-order$-like phase tramsitions observed in clusters [39-41]
and polymers $[42, 43]$ .
Now, the mechanism of the switching to the favored isomer can be $rationd\dot{u}$ ed in termsof the “work” necessary to bring the system against the averaged force fields obtained in
Figs. 5(a) and (b). The line integral of the negative of the averaged potential force fieldalong the reaction path gives a reconstructed potential energy$(\langle f_{pot,1}\rangle, \langle f_{pot,2}\rangle, \langle f_{pot,3}\rangle)$curve$V_{pot}$as$V_{pot}=- \int_{path}\sum_{\dot{*}=1}^{3}\langle f_{pot,i}\rangle da_{i}$

.

(36)

We carried out this line integral numerically starting at one of the end points of theis the smallest. Since the absolute value of the potential
reaction path, where theenergy is arbitrary, we set
to zero at the starting point of the line integral. The
potentialis essentially the reduced potential energy curve that the gyration radii feelis shown for thealong the reaction coordinate. In the lower three panels of Fig. 6,, with dotted, andthree representative total energy values, $E=-11.2\epsilon,$on thedependingsignificantly
curves. The essential topography ofdoes not change
always has two wells. The left well corresponds to the OCT
total energy. That is,$a_{1}$$V_{pot}$$V_{pot}$$V_{pot}$$-10.0\epsilon$$V_{pot}$$V_{pot}$$-8.5\epsilon$

100

isomer while the right one corresponds to the CTBP isomer. The well for the OCT isomeris always much deeper than that for the CTBP isomer. The potential
resembles theoriginal potential energy curve along the steepest descent path shown in Fig. l(a). The
is obviously not sufficient for explaining the energy-dependent change of thecurve ofprobability distribution along the reaction paths shown in the upper panels of Fig. 6.
We next superpose the effect of the internal centrifugal force over that of the potential. Similarly to Eq. (36), a superposed line integral of the averaged potential
curve) is introduced
force and the averaged internal centrifugal foroe$V_{pot}$$V_{pot}$$V_{pot}$$(\langle f_{cont,1}\rangle, \langle f_{cent,2}),$$\langle f_{cent,3}\rangle$as$V_{pot+ccnt}=- \int_{path}\sum_{:=1}^{3}(\langle f_{pot,i}\rangle+\langle f_{cent,i}\rangle)da:$

.

(37)

The solid curves in the lower three ptels of Fig. 6are the superposed ener $gyV_{pot+\infty nt}$ .
, this energy curve $V_{pot+cent}$ shows two welk corresponding to
Similarly to the curve ofthe two isomers. However, the topography of $V_{pot+oent}$ chtges dramatically depending on,
the total energy of the system. At low total energy, $E=-11.2\epsilon$ , the left $weU$ ofwhich corresponds to the OCT isomer, is deeper than the right wef that corresponds to
the CTBP isomer. As the total energy increases, the left $weU$ gets shallower while the rightwell gets deeper. This is because the effect of the internal centrifugal foroe becomae moresignificrt as the total ener$y$ of the dynamics increases, while the potential curve
does not itge significrtly. Note that the internal centrifugal force has the directionality
to “push” the $sy$stem from the OCT region to the CTBP region as confirmed in Fig.
5(a). In this way, the OCT isomer becomae less stable dynamically while the CTBP
total energy increases. This dramatic iange of theisomer becomae more preferred
explains the chtge of the probability distribution insuperposed energy curve
the upper ptels of Fig. 6fairly well. To summarize, the predomintce of the OCT
isomer at low energy is mainly supported by the predominance of the potential force,
while the predomintce of the CTBP isomer at high energy is mainly supported by the
predominance of the internal centrifugal force.
$V_{pot}$$V_{pot+c\bm{m}t}$$V_{pot}$$a\epsilon$$V_{pot+c\varpi t}$

5

Summary and Outlook

illustrative
By taking up the structural transition dynamics of asix-atom cluster
example, we have presented ageneral methodolo $gy$ to elucidate the mechanism for largeamplitude $coUective$ motions of molecules with mry degrees of freedom in terms of the
reduced dynamics of the three molecular yration radii. Based on the hamework ofgeometric mechanics, we have first separated the three rotational degrees of freedom td-atom system. Then we have applied thethe $(3n-6)$ internal degrees of freedom of
principal-axis hyperspher.ical coordinatae to further decompose the total $(3n-6)$ internal
degrees of freedom into the three gyration radii and the $(3n-9)hyPer\bm{t}1^{1ar}$ modes.
coUective variablae, while the remaining
The three gyration radii are regarded as slow
hyperangular variablae are regarded ae “bath” modes. This classification of variablae intocollective“bath” variablae has been supported by the numerical observations on timescale separation and kinetic energy partitioning. The time scale separation has justified$\mathfrak{B}\bm{t}$$\bm{t}n$$\bm{t}d$$\bm{t}d$101

ur strategy to average out the hyperangular variables from the equations of motion for
the three gyration radii, our collective variables.
We have scrutinized the averaged equations of motion for the gyration radii to reveal
the mechanism of competition between akinematic force
the potentlal force in the
motions of the cluster. The kinematic force is identified asinternal
centrifugal force that arisae from the dynamical coupling of yration radii with the hyperangular modes via the intrinsic non-Euclidet metric of the internal space. This internal
centrifugal force has the remarkable effect of elongating
inflating the
distribution of the syst $em.$ Thaee effects often counteract the effect of the potential force, whii
generally serves to keep the mass distribution of the system symmetriccompact. Via
the competition between these two forces, the most preferable conformation of the system
is determined.
While the potential function is sp cific to the system of interest, the kinetic enery
expressions such as Eq. (29) td Eq. (32) are common to general -atom moleculae.
Therefore the internal centrifugal force terms in Eqs. (33)-(35) also arise ubiquitously
in the dynamics of general -atom moleculae. Rom this aeon, the elongationtheinflation effects of the internal centrifugal force should be of general importance in awide
class of molecular reactions, where the syst $ems\ddagger gnific\bm{r}tly$ chtges its mass distribution.
For example, because of the directionality of this force, its roles should be importrt indissociation and recombination reactions.
It is an interesting next step to elucidate thespace structure that mediates the
low-dimensional dynamics of gyration radii. Recently, the dynanical systems approai
hae provided amathematical basis for the study of non-statistical reaction procaes\’e ofmoleculae [45-49]. Since the hi$gh$-dimensional nature of phase space has often prevented
from applying the dynamical systems approach to molecular reactions, our strategy tofocus on the low-dimensional dynamics of gyration radii should be very useful. The
space structure of the reduced dynamics of molecular yration radii will be scrutinized in
our future publications. We also plt to extend the present methodolo $y$ to the systems inseveral different environments: Large-amplitude motions of moleculae with non-zero total
angular momentum are an interesting issue for the next step, in which the kinematicforces such as the (normal) centrifugal force and the Coriolis force also come into play inaddition to the internal oentrifugal force of this study. It is ako an interaeting issue toextend the present approach to the systems in thermal environments such as biomolecules.
Since the internal centrifugal force has ageneral tendency to become strong in high enery(or temperature) range and has an effect of inflating
elongating the system, this force
is naturafy expected as adriving force for the unfolding of proteins.
Finally, an important consequenoe of the present study is that molecular vibrations
themselves ct induce adynamical force such as the internal centrifugal force, whii ct
be acritical driving force for large-amplitude motions.of the system. It would be interesting to note that this kind of dynanic forces are quite analogous to the dynamic forcesthat stabilize the inverted pendulum whose suspension point is violently vibrated [50] as
well as the particles in electromagnetic traps [51]. Obviously, the importance of this kindof dynamical effect hae not been fully appreciated in the current reaction-rate theories. In
the current standard reaction-rate theories, molecular $re$actions are usuaUy regard $ed$ as
an event in which asyst $em$ surmounts abarrier of (Born-Oppenheimer) potential ener $y$$o$$\bm{t}d$$l\arg\triangleright amplitude$$\bm{t}$$\bm{t}d$$ma\epsilon s$$\bm{t}d$

$e$

$n$

$n$

$re$

$\bm{t}d$$ph_{\mathfrak{B}}e$$ph_{\mathfrak{B}}e$$\bm{t}d$102

along a certain reaction coordinate. However, the result of the present study and anotherrecent one [52] suggest that dynamic forces, such as the internal centrifugal force, modify

the original potential energy barrier significantly and induces another dynamical barrier
for reactions. Thus, in order to improve the current reaction-rate theories, it would bevery important to take into consideration not only the static potential energy barriers
but also such dynamic barriers. This kind of dynamic barriers should definitely shed newlight on the conventional picture of molecular reactions.

Acknowledgments
This report is based on the joint work with Wang S. Koon (Caltech), Jerrold E. Marsden(Caltech), and Ioannis G. Kevrekidis (Princeton University). The author thanks for theirvery stimulating discussions and thoughtful suggestions in developing this research.

References[1] R. Zwanzig, Phys. Rev. 124, 983 (1961).
[2] H. Mori, Prog. Theor. Phys. 33, 423 (1965).
[3] R. Zwanzig, Nonequilibrium Statis tical Mechanics (Oxford University Press, New
York, 2001).

[4] A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. 9912562 (2002).
[5] R. Elber, A. C\’ardenas, A. Ghosh, and H. A. Stern, Adv. Chem. Phys. 126, 93 (2003).

[6] G. Hummer and I. G. Kevrekidis, J. Chem. Phys. 118, 10762 (2003).

[7] D. I. Kopelevich, A. Z. Panagiotopoulos, and I. G. Kevreidis, J. Chem. Phys. 122,
044907 (2005); 122, 044908 (2005).
[8] W. H. Miller, N. C. Handy, and J. E. Adams, J. Chem. Phys. 72, 99 (1980).
[9] P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, Annu. Rev. Phys. Chem.
53, 291 (2002).
[10] W.

$E$, W. Ren, and E. Vanden-Eijnden, Chem. Phys. Lett. 413, 242 (2005).

[11] J. E. Marsden, T. S. Ratiu, and J. Scheurle, J. Math. Phys. 41, 3379 (2000).

[12] J. E. Marsden and T. S. Ratiu, Introduction to Mechanics and Symmetry (Springer,
New York, 1999).
[13] M. Kummer, Indiana Univ. Math. J., 30, 281 (1981).
[14] J. E. Marsden, R. Montgomery, and T. Ratiu, Contemp. Math. 28, 101 (1984).

103[15] A. Guichardet, Ann. Inst. H. Poincar\’e 40, 329 (1984).

[16] T. Iwai, Ann. Inst. H. Poincar\’e 47, 199 (1987).
[17] T. Iwai, J. Math. Phys. 28, 964 (1987).
[18] A. Tachibana and T. Iwai, Phys. Rev. A 33, 2262 (1986).

[19] R. G. Littlejohn and M. Reinsch, Rev. Mod. Phys, 69 , 213 (1997).

[20] C. Eckart, Phys. Rev. 46, 383 (1934).
[21] X. Chapuisat and A. Nauts, Phys. Rev. A 44, 1328 (1991).
[22] X. Chapuisat, Phys. Rev. A 45, 4277 (1992).

[23] A. Kuppermann, Adv. Mol. Vib. Col. Dyn. $2B$ , 117 (1993).

[24] A. Kuppermann, J. Phys. Chem. 100, 2621 (1996).
[25] A. Kuppermann, J. Phys. Chem. 101, 6368 (1997).

[26] R. G. Littlejohn and M. Reinsch, Phys. Rev. A 52, 2035 (1995).
[27] R. G. Littlejohn, K. A. Mitchell, V. A. Aquilanti, and S. Cavalli, Phys. Rev. A 58,
3705 (1998).

[28] R. G. Littlejohn, K. A. Mitchell, M. Reinsch, V. Aquilanti, and S. Cavalli, Phys.
Rev. A 58, 3718 (1998).
[29] V. Aquilanti, A. Lombardi, and E. Yurtsever, Phys. Chem. Chem. Phys. 4, 5040
(2002).
[30] V. Aquilanti, A. Lombardi, M. B. Sevryuk, and E. Yurtsever, Phys. Rev. Lett. 93,
113402 (2004).

[31] V. Aquilanti, A. Lombardi, and M. B. Sevryuk, J. Chem. Phys. 121, 5579 (2004).

[32] T. Yanao and K. Takatsuka, Phys. Rev. A 68, 032714 (2003).
[33] T. Yanao and K. Takatsuka, J. Chem. Phys. 120, 8924 (2004).
[34] T. Yanao and K. Takatsuka, Adv. Chem. Phys. 130 $B,$ $87$ (2005).

[35] T. Yanao, W. S. Koon, J. E. Marsden, and I. G. Kevrekidis, J. Chem. Phys. 126,
124102 (2007).
[36] R. S. Berry, Chem. Rev. 93, 2379 (1993).

[37] P. A. Braier, R. S. Berry, and D. J. Wales, J. Chem. Phys. 93, 8745 (1990).

[38] T. Komatsuzaki and R. S. Berry, J. Chem. Phys. 110, 9160 (1999); 115, 4105 (2001).

104[39] J. Jellinek, T. L. Beck, and R. S. Berry, J. Chem. Phys. 84, 2783 (1986).
[40] H. L. Davis, J. Jellnek, and R. S. Berry, J. Chem. Phys. 86, 6456 (1987).

[41] P. Labastie and R. L. Whetten, Phys. Rev. Lett. 65, 1567 (1990).
[42] Y. Zhou, C. K. Hall, and M. Karplus, Phys. Rev. Lett. 77, 2822 (1996).
[43] M. Takano, T. Takahashi, and K. Nagashima, Phys. Rev. Lett. 80, 5691 (1998).

[44] G. Strang, Linear Algebra and its Applications (Academic Press, New York, 1976).

[45] N. De Leon, M. A. Mehta, and R. Q. Topper, J. Chem. Phys. 94, 8310 (1991).
[46] T. Uzer, C. Jaff\’e, J. Palaci\’an, P. Yanguas, and S. Wiggins, Nonlinearity 15, 957(2002).

[47] H. Waalkens, A. Burbanks, and S. Wiggins, J. Chem. Phys. 121, 6207 (2004).

[48] F. Gabern, W. S. Koon, J. E. Marsden, and S. D. Ross, Physica $D211$ , 391 (2005).
[49] C. B. Li, Y. Matsunaga, M. Toda, and T. Komatsuzaki, J. Chem. Phys. 123, 184301(2005).

[50] M. Levi, Physica $D132$ , 150 (1999).
[51] W. Paul, Rev. Mod. Phys. 62, 531 (1990).
[52] T. Yanao, W. S. Koon, and J. E. Marsden, Phys. Rev. A 73, 052704 (2006).