Accelerated Molecular Dynamics through stochastic iterations to strengthen yield of path hopping over upper states (SISYPHUS)

5.7 MB·English

by  Pratyush Tiwary

Accelerated Molecular Dynamics through stochastic iterations to strengthen yield of
path hopping over upper states (SISYPHUS)
Pratyush Tiwary and Axel van de Walle
School of Engineering, Brown University, Providence, Rhode Island 02912, USA
(Dated: January 3, 2013)
We present a new method, called SISYPHUS (Stochastic Iterations to Strengthen Yield of Path
Hopping over Upper States), for extending accessible time-scales in atomistic simulations. The
method proceeds by separating phase space into basins, and transition regions between the basins
based on a general collective variable (CV) criterion. The transition regions are treated via tra-
ditional molecular dynamics (MD) while Monte Carlo (MC) methods are used to (i) estimate the
expectedtimespentineachbasinand(ii)thermalizethesystembetweentwoMDepisodes. Inpar-
3 ticular, an efficient adiabatic switching based scheme is used to estimate the time spent inside the
1 basins. The method offers various advantages over existing approaches in terms of (i) providing an
0
accuraterealtimescale,(ii)avoidingrelianceonharmonictransitionstatetheoryand(iii)avoiding
2
the need to enumerate all possible transition events. Applications of SISYPHUS to low temper-

n ature vacancy diffusion in BCC Ta and adatom island ripening in FCC Al are presented. A new
a CVappropriateforsuchcondensedphases,especiallyfortransitionsinvolvingcollectivemotionsof
J several atoms, is also introduced.

2

] Achieving usefully long timescales (seconds or longer) coarse grained dynamics typically provided by Metady-
ci inatomisticsimulationsofmaterialsisaproblemofgreat namics methods1,15. The method also provides an accu-
s interest, and the search for a practical and general solu- raterealtimescalethatdoesnotdeterioratewithsystem
-
l tion has generated intense activity in the field over last size and that does not rely on harmonicity assumptions.
r several decades1–9 (see10–12 for excellent reviews of these There is no need to construct a priori or in situ a cat-
t
m and other efforts). The problem arises because, as the alog of possible transition mechanisms. The method is
. system moves from one energy basin to another through specially suited for exploiting massive parallel comput-
t
a infrequent rare events, it stays trapped in some energy ing.
m basinforextendedperiodsoftime. Alongwiththesmall Theproposedalgorithmgeneralizesourpreviouswork3
- time steps (on the order of femtoseconds) needed for the along multiple dimensions. First, we use a general CV
d totalenergystayingconserved, thisseverelyrestrictsthe χ (instead of the system’s potential energy) to discrimi-
n
time scales accessible in MD simulations and also leads nate between basins and transition regions and propose
o
c to limited phase-space exploration. a novel type of CV suitable for this purpose in con-
[ Although many methods exist to increase the rate of densed phases. (Recently proposed dimensionality re-

rare events and efficiently explore the rugged free energy duction algorithms16,17 that discover CV automatically
2
landscape, a frequent limitation is the inability to effi- could be used as well.) Second, we introduce an adi-
v
9 ciently obtain an accurate estimate of the “real” time abatic switching scheme to efficiently calculate the real
4 scale of the simulation in general physical systems. Ex- timespentinsidewells. Finally,weuseamorerobustcri-
6 isting schemes to achieve this either terion to determine when the system has been trapped
6 in a basin for sufficiently long time to have equilibrated
. 1. require cataloguing all possible transitions paths
2 therein.
out of given basin4–6, which can be computation-
1 Let the system be characterized by position x =
2 ally prohibitive in low-symmetry systems, or
(r ,...,r ) and velocity v = (v ,...,v ). The CV χ
1 1 N 1 N
is a function of x (we give an example of such a function
: 2. rely on harmonic approximations to the system’s
v energy surface13, or laterinthisletter). Forauser-specifiedcut-offvalueχcut,
i
X we define the basins, or wells, W in this χ-space as a set
r 3. involve computing averages that converge slowly, of connected states for which χ<χcut. In the wells, the
a especially for large system sizes, because they in- method does not follow the system’s exact trajectory in
volve exponentials of the system’s total energy.2,14 phase space, but instead provides the expected amount
of time spent in each well. In contrast, the χ ≥ χ
cut
In this letter, we propose a new mixed Monte Carlo- regionofthe phasespacecontains theinterestingbutin-
Molecular Dynamics method that uses a collective vari- frequently occurring events whose dynamics is fully de-
able (CV) χ solely to discriminate between basins and scribed. The choice of χ thus specifies the level of de-
cut
transition regions, thus placing very weak requirements tails one wishes to retain in the simulation. Large values
onthechoiceofCV,lessstringentthanwhatrequiredin of χ may cause wells to merge and limit the ability to
cut
metadynamicsmethods31. Themethodstillprovidesthe resolvetheprecisedynamicsofsomeevents. Themethod
system’s dynamics via conventional atomic coordinates is still formally correct but the definition of the wells W
and thus provides more detailed information than the maynothaveananobviousphysicalmeaning. Neverthe-
2
less, the method is very robust to changes in χ that T (Kelvins)
do not change the topology of the well structurec(uwte will 10−101250.0 1000.0 833.3 714.3 625.0 555.6 500.0
Brute−force MD
provide numerical evidence of this fact). SISYPHUS
When the CV χ of the system is above the thresh- c) 10−15
otrldueχHcuatm,itlthoenisaynstwemithecvoonlvsetasnvt-iapreMssDureac(coorrdvionlugmteo)iotsr 2m/se
D ( 10−20
constant-temperature (or energy) ensemble as entailed
by the simulation. When χ < χ , the algorithm con-
cut
tinues performing MD until either (i) χ≥χcut (in which 10−25 0.8 1 1.2 1.4 1.6 1.8 2
case the system is considered to have exited the well inverse T (Kelvins−1) x 10−3
and standard MD continues as described above) or (ii)
a time equal to the system’s decorrelation time τ has FIG.1: DiffusionconstantforvacancydiffusioninTaatvar-
c
elapsed. In that latter case, the system is considered to ious temperatures as through brute-force MD (circles) and
havebeentrappedinthewellforlongenoughthatithas SISYPHUS (stars). Errorbars (roughly same as marker size)
are also provided as obtained over 16 independent runs.
reached a local thermodynamic equilibrium. This crite-
rion is similar to the one used in other methods8,18 to
define a transition event. At this time, we launch a MC
bias inside the well changes the fraction of time spent at
simulation (called MCa) whose aim is to generate a new
the boundary but not the ratio of the times spent at any
random starting position at the well’s boundary to ini-
twopointsoftheboundary. Theparametermdetermines
tializethenextMDepisode. MCaisrunforlongenough
howsharptheboundaryofthewellis(avaluearound0.5
that the system loses memory of how it entered the well
is found to perform well in practice). V is kept around
and visits the boundary of the well a few times. MD is 0
thestandarddeviationinpotentialenergyofthesystem.
restarted with the position x where the system last vis-
As we demonstrate numerically later, the algorithm is
itedthewell’sboundaryandwithvelocitiesvdrawnfrom
very robust with respect to choice of parameters in Eq.
atruncatedMaxwell-Boltzmanndistributionconditional
(2).
on v·∇χ(x) > 0 (i.e. we only consider velocities in the
Havingdescribedawaytoacceleratetheexplorationof
half-space pointing outwards of the well).
various wells, we now turn to the question of calculating
In parallel to the first MC (MCa) we perform a sec-
(via a Monte Carlo labelled MCb) the expected time the
ond MC run (called MCb) that calculates the expected
system would have spent inside well W if there was no
time t the system would have spent inside the well.
W
acceleration of the dynamics. This time, denoted t ,
ThisseparationoftwoMCrunsmakesouralgorithmex- W
can be calculated as the reciprocal of the flux of states
tremelyparallelizableandespeciallyamenabletobeused
exiting the well:
on loosely coupled cluster of computers. We can launch
asmanyMCbrunsaswehaveprocessorsavailable. These
v
runs do not need to communicate with each other, and t = lim((cid:104) 1(x∈S )(cid:105))−1 (3)
W w→0 w w
because of the system’s ergodicity, we can make a quick
estimate of the quantity tW by averaging over these in- where the average (cid:104)···(cid:105) is taken over x drawn from
dependentruns. Thisparallelizationisevensimplerthan the well W with a probability density proportional to
for the Parallel Replica method18. e−V(x)/(kBT). kB is Boltzmann’s constant, T is the tem-
Before we describe how we calculate the time the sys- perature and 1(A) equals 1 if the event A is true and 0
tem would have spent in whichever well it visited, we otherwise. v denotes the mean projection of a Maxwell-
describe specifically how MCa is implemented. We de- Boltzmann-distributed velocity along the unit vector u
fine the boundary of the well W with thickness w: parallel to ∇χ(x), conditional on v.∇χ(x) < 0. When
(cid:112)
all atoms have the same mass m, v = k T/2πm (a
S ={x=(r ,...,r ):|χ(x)−χ |<w} (1) B
W 1 N cut general expression can be found in Ref. 3).
Trying to visit S with an unbiased potential would be AstraightforwardimplementationofEq. (3)willagain
W
of no avail, since we will rarely visit states in S . We suffer from a rare event problem. Even an importance
W
thus do MC with a biased potential V∗(x) defined as sampling scheme, as suggested in Ref. 3 is not very ef-
follows: ficient. Consider what happens if we used the biased
potential as defined in Eq. (2) to calculate t :
(cid:40) W
∞ χ(x)≥χ
cut
V∗(x) = V(x)+ V (cid:16)χcut−χ(cid:17)m χ(x)<χ (2) (cid:104)e−β(V(x)−V∗(x))(cid:105)∗
0 χcut cut tW =wli→m0(cid:104)wve−β(V(x)−V∗(x))1(x∈Sw)(cid:105)∗ (4)
PerEq. (2),MCanevervisitstheoutsideofthewellW
and biases states with a penalty function that increases where (cid:104)···(cid:105)∗ denote expectations taken under a density
with their depth inside the well. Note that the bias is proportionaltoe−βV∗(x),inwhichβis1/(k T). Thisap-
B
zero at the well boundary, which is important to obtain proach is exact in the limit of an ensemble average, but
the correct sampling distribution of the boundary2. The thereisafundamentaltrade-offthatlimitsitsusefulness:
3
10−14 D106 10−16
M
o
2D (m/sec)1100−−1186 −up relative t110024 2D (m/sec)1100−−1187
eed 10−19
p
10−020. 12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 s1000. 12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 1 .8 2 2.2 2.4 2.6 2.8 3 3.2
χ − <χ> (units of 0K lattice parameter) χ − <χ> (units of 0K lattice parameter) V (eV)
cut cut 0
FIG.2: (Top)Insensitivityofdynamicstochoiceofχ (relativetoaverageχatthattemperature)forvacancydiffusioninTa
cut
across temperatures. (cid:5), ◦,(cid:63) and ∗ denote 500,600,700 and 800K temperatures respectively. (Centre) Corresponding speed-ups
relative to physical time achieved in brute-force MD run in the same wall-clock time. (Bottom) Insensitivity of dynamics to
choice of V . Errorbars over 16 independent runs for each data point.
0
A large bias V∗(x)−V(x) leads to a more rapid conver- only when x ∈ S , and whenever it is nonzero, the
w
gence of the denominator (due to an increased sampling difference Vˆ (x,0)−Vˆ (x,1) is very small(see Eqs.(1-2)).
rate of the boundary) but a slower convergence of the Since this average is calculated with the maximally
numerator (due to an increase in V∗(x)−V(x)). biased potential Vˆ (x,1), the boundary x∈S is visited
w
To avoid this problem, we now propose a technique frequently, and thus the first term in Eq.(5) can be
that bears some resemblance to adiabatic switching evaluated very quickly. The second part in Eq.(5) is R,
methods19,20, in which the system is continuously, adi- where the average (cid:104)∂Vˆ(x,α)(cid:105) = (cid:104)V∗(x)−V(x)(cid:105) does
abatically switched from V (x) (the true potential) to ∂α α α
not contain any exponentials that could cause a slow
V∗(x) (identical to the potential used in MCa). Let
convergence. In SM, we prove rigorously that for this
Vˆ (x,α) smoothly interpolate between Vˆ (x,0) ≡ V (x) switching scheme and for the choice of biasing potential
andVˆ (x,1)≡V∗(x). Thenwecanexpresstheensemble in Eq.(2) we can use a non-uniform grid to evaluate R
averageinEq.(3)asbelow(workingintermsofratet−1): which can be made finer as α → 0 but kept coarse for
W
larger α, leading to further computational efficiency.
(cid:82) v¯1(x∈S )e−βVˆ(x,0)dx For solid-state systems where bond-breaking is the
t−1 = lim w w
W w→0 (cid:82) e−βVˆ(x,0)dx dominant mechanism of interest, we take χ to be the
bond distortion function (BDF), defined below for a N-
(cid:28) (cid:29)
v1(x∈S )
≡ lim w e−β(Vˆ(x,0)−Vˆ(x,1)) R (5) particle system:
w→0 w
1  1/p
wheredxdenotesadifferentialvolumein3-Ndimensional χ(x)=a0 (cid:88) |rijr−eqriejq|p (9)
configuration space for N particles, the integration being ∀i,j rij=|ri−rj| ij
performedoverentireconfigurationspacewithinthewell
wherep>1anda isthe0Kelvinlatticeparameter. For
W and the expected value (cid:104)···(cid:105) in Eq.(5) is defined by 0
α eachbond,req istheequilibriumbondlengththatcanbe
ij
(cid:82) (···)e−βVˆ(x,α)dx obtained by a few conjugate gradient steps each time a
(cid:104)···(cid:105) = (6) transitionisdetected. Inthelimitthatp→∞,theBDF,
α (cid:82) e−βVˆ(x,α)dx which is a p-norm over fractional bond distortions, ap-
proachesthemaximumnorm,i.e. thestrain(timesa )in
0
Below we define the term R in Eq.(5) and re-express it
themaximallystrainedbond. Forp→∞wethusrecover
in a computationally tractable form (see Supplemental
the so-called bond-boost function21. We pick p around
Materials (SM) for a more detailed derivation):
8-12 for the systems studied in this letter. Not taking
(cid:82) e−βVˆ(x,1)dx (cid:32) (cid:90) 1(cid:42)∂Vˆ (x,α)(cid:43) (cid:33) p=∞allowsustotreatonasimilarfootingcaseswhere
(a)onebondisdistortedbyalargeamount,or(b)several
R = =exp −β dα
(cid:82) e−βVˆ(x,0)dx ∂α bonds are collectively distorted by a significant amount
0
α
which is however less than the amount in (a). As soon
(7)
as either (a) or (b) happens, the BDF detects it through
We pick a linear switching scheme for Vˆ(x,α), i.e. an a spike and thus we can then switch back to doing com-
pletely unbiased MD. This way we can treat transition
interpolation scheme between Vˆ(x,0) and Vˆ(x,1):
mechanisms involving small but concerted and collective
motion of several atoms (see Fig. 3 g-h)) - for example
Vˆ(x,α)=(1−α)V(x)+αV∗(x) (8)
in glasses, shear transformation zones22 involve several
atoms moving displacing together by a small amount.
We now make a few observations regarding
We now describe applications of the algorithm that
Eq.(5). It involves 2 parts. The first is
(cid:68) (cid:69) validate it and demonstrate its insensitivity to choice of
limw→0 v1(xw∈Sw)e−β(Vˆ(x,0)−Vˆ(x,1)) . This is nonzero parameters. The first example is vacancy assisted lattice
1
4
(a)1adatom (b)1adatom (c)2adatoms (d)2adatoms (e)Ad&sub- (f)Ad& (g)3adatoms (h)3adatoms
(before) (after) (before) (after) strate(before) substrate(after) (before) (after)
FIG. 3: Mechanisms seen by SISYPHUS for adatom movement on Al(001) surface. Atoms colored per z-coordinate.
Blue/orange = substrate atom, green/red=adatom, red/orange = atom with maximum movement. Solid lines are periodic
boundaries.
diffusion at low temperatures in BCC Tantalum. Lat-
tice diffusion at low temperatures is a problem impor-
tant in a spectrum of sciences from Materials Science to
Geology23,24, but is beyond the time scales one can ac-
cess in current MD simulations, requiring times longer
than milliseconds3,25. We describe the parameters for
(a)t=0ns(0eV) (b)t=1µs(-1.6eV) (c)t=2ms(-1.8eV)
the MD part26 of this simulation in SM. In Fig.1 we
demonstrate how the diffusion constant through brute-
force MD and SISYPHUS for vacancy assisted diffusion
lieonthesameArrheniusplot, givinganArrhenius-type
activation energy of 0.9(+/-0.1) eV (in rough agreement
with Harmonic Transition State Theory(HTST) calcula-
tion of 1.1 eV). Fig.2(top) demonstrates the lack of sen- (d)t=4ms(-2.6eV) (e)t=9ms(-3.7eV) (f)t=15ms(-5.0eV)
sitivity of the dynamics to what χ and V values we
cut 0
pick. As expected the speed-up is higher for higher χcut FIG. 4: Room-temperature adatom ripening on Al (001)
(Fig.2(middle)). At higher temperature we find slightly surface as a function of time. Each structure was quenched
highersensitivitytoχ (stillwithinorderofmagnitude to find its corresponding local minima, and the energies af-
cut
accuracy). This is because for a low χ , the t values ter quenching (relative to (a)) are provided. (a) shows the
cut W
becomesclosertotimeτ thesystemtakestoequilibrate startinggeometry. After1µs,wehavetwodisconnectedclus-
c
in a well. Insensitivity to V values can be seen from ters (b-c). Corresponding brute-force MD runs were found
0
trapped in similar configuration in the fraction of microsec-
Fig.2(bottom). A smaller V leads to slower sampling in
0 ond time they could achieve. At around 4ms these two clus-
Eq.4, and MCa also takes longer to converge. With too
ters join (d). At 9ms (e) there is further joining and we have
high a V however, the periphery of the well W can be
0 effectively one long chain of adatoms. At 15ms (f) the chain
too steep and one might again face sampling issues since
has coarsened into one entity across simulation cells. Color
the system can be trapped in some regions of the well scheme same as in Fig. 3.
boundary. In SM, we provide a back-of-the-envelope es-
timateofhowtopickanoptimumV foragivensystem.
0
For our second application, we studied the room- used to get rough estimate of the time-scale for adatom
temperature dynamics of Al adatoms on Al (001) thin islandripeningthatwecancompareSISYPHUSwith. In
film (625 atoms, roughly 3 times larger than the Ta va- SM we provide details of the simulation parameters for
cancy diffusion example). We picked this problem be- this system.
causefirstly,fromatechnologicalperspective,itisofim- Fig. 3 and the movies in SM illustrate the most com-
mense importance for fabrication processes in nanoscale mon mechanisms seen through SISYPHUS. We have the
devices involving growth of thin films from deposited single adatom hop (a-b), the concerted two adatom hop
adatoms27,28. This problem is very interesting from a (c-d), and the concerted event involving an adatom and
theoretical perspective too, given that it is an inherently a substrate atom as the latter moves to the adatom
non-equilibrium phenomenon dictated by the interplay layer(e-f). Occasionallyweseemorecomplicatedmecha-
between kinetics and thermodynamics. Being able to nismslikethe3-atomhop(g-h),andeventswithcreation
modelandcontrolthegrowthandpropertiesofsuchfilms of a vacancy in the top substrate layer (see SM). The
is hugely desirable - the time-scales needed are however first three mechanisms are the most common and are in
far beyond MD. Accurate 0 Kelvin saddle point search fact the lowest energy transitions found using the dimer
methods5,29 have shown the existence of a large number method for saddle point search. In Fig. 4 and movie in
oftransitionmechanismswithlowandsimilaractivation SM we show the typical evolution of the island ripening
energies(smallerthan0.4eV).Suchlowlyingbarrierscan at room temperature over several milliseconds of physi-
behardtodealwithinmostacceleratedMDmethods10. cal time (exact time can vary from run to run well but
On-the-flyKMC5 calculationshavebeenusedpreviously stillwithinanorderofmagnitude). Shownalongwithare
5
correspondingenergiesobtainedbyquenchingeachstruc- of atoms. The method works well irrespective of system
turetoitslocalminima,illustratingfurthertheeffective- size, and can be applied in the general setting of any
ness of the algorithm in escaping and exploring various collective variable. We also introduced a new CV appro-
local energy minima. The overall time can be compared priate for solid state systems especially for transitions
with the analogous on-the-fly KMC work5 for same sys- involving collective motions of several atoms.
tem and interatomic potential30 where 20 adatoms (half
as many as current work) formed one compact cluster We would like to thank Dr. Arthur Voter and Prof.
around 1ms. As a verification that our proposed BDF is Graeme Henkelman for helpful discussions and stimulat-
effective to decomposing phase space into disconnected ingthetrajectoryofthisresearch,Dr. ChristopherWein-
wells, we show, in a movie accompanying the SM, a su- bergerfortheHTSTcalculationforTavacancydiffusion,
perpositionofsnapshotsofthesystemduringaMCarun, andProf. JuliaGreerforhospitalityduringPT’svisitto
illustrating that the system does not jump from one well Caltech. ThisresearchwassupportedbytheUSNational
toanotherwithinoneMCarun(sinceitrejectsallmoves Science Foundation through XSEDE computational re-
with χ≥χ ). sources provided by NCSA under grant DMR050013N
cut
In conclusion, we have shown SISYPHUS to be an ex- and DMR120056, and NSF Condensed Matter and Ma-
tremely parallelizable and robust set of algorithms that terials Theory program DMR0907669, and the Office of
help achieve fraction of second timescale for thousands Naval Research through grant N00014-11-1-0886.
1 A. Laio and M. Parrinello, Proceedings of the National 17 G.A.Tribello,M.Ceriotti,andM.Parrinello,Proceedings
Academy of Sciences 99, 12562 (2002). of the National Academy of Sciences 109, 5196 (2012).
2 A. F. Voter, Phys. Rev. Lett. 78, 3908 (1997). 18 A. F. Voter, Physical Review B 57, 13985 (1998).
3 P. Tiwary and A. van de Walle, Phys. Rev. B 84, 100301 19 D. Frenkel and B. Smit, Understanding molecular simu-
(2011). lation : from algorithms to applications (Academic Press,
4 Y. Fan, A. Kushima, S. Yip, and B. Yildiz, Phys. Rev. 2002), 2nd ed., ISBN 0122673514.
Lett. 106, 125501 (2011). 20 C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
5 G. Henkelman and H. Jonsson, The Journal of Chemical 21 K. A. Fichthorn, R. A. Miron, Y. Wang, and Y. Tiwary,
Physics 115, 9657 (2001). Journal of Physics: Condensed Matter 21, 084212 (2009).
6 G. T. Barkema and N. Mousseau, Phys. Rev. Lett. 77, 22 M.L.FalkandJ.S.Langer,Phys.Rev.E57,7192(1998).
4358 (1996). 23 S. Mrowec and S. Marcinkiewicz, Defects and diffusion in
7 K.Lindorff-Larsen,S.Piana,R.O.Dror,andD.E.Shaw, solids: an introduction (Elsevier, 1980).
Science 334, 517 (2011). 24 E. B. Watson and E. F. Baxter, Earth and Planetary Sci-
8 C.-Y.Lu,D.E.Makarov,andG.Henkelman,TheJournal ence Letters 253, 307 (2007).
of Chemical Physics 133, 201101 (pages 4) (2010). 25 M. I. Mendelev and Y. Mishin, Phys. Rev. B 80, 144111
9 X.-M. Bai, A. F. Voter, R. G. Hoagland, M. Nastasi, and (2009).
B. P. Uberuaga, Science 327, 1631 (2010). 26 G. Ackland and R. Thetford, Philosophical Magazine A
10 A.Voter,F.Montalenti,andT.Germann,AnnualReview 56, 15 (1987).
of Materials Research 32, 321 (2002). 27 Z. Zhang and M. G. Lagally, Science 276, 377 (1997).
11 A.LaioandF.L.Gervasio,ReportsonProgressinPhysics 28 K.E.Becker,M.H.Mignogna,andK.A.Fichthorn,Phys.
71, 126601 (2008). Rev. Lett. 102, 046101 (2009).
12 D. Perez, B. Uberuaga, Y. Shim, J. Amar, and A. Voter, 29 G. Henkelman and H. Jonsson, The Journal of Chemical
AnnualReportsinComputationalChemistry5,79(2009). Physics 111, 7010 (1999).
13 M. Sørensen and A. Voter, The Journal of Chemical 30 A. Voter and S. Chen, in Mater. Res. Soc. Proc (1987),
Physics 112, 9599 (2000). vol. 82, pp. 175–180.
14 C. Jarzynski, Physical Review E 73, 046105 (2006). 31 In Metadynamics, the CV should distinguish between ini-
15 G.Bussi,A.Laio,andM.Parrinello,Phys.Rev.Lett.96, tialbasin,finalbasinandtransitionregions.OurCVhow-
090601 (2006). ever is allowed (but not required) to take the same value
16 G.A.Tribello,M.Ceriotti,andM.Parrinello,Proceedings in two basins.
of the National Academy of Sciences 107, 17509 (2010).