Efficient exploration of configuration space and identification of metastable structures in condensed phase systems are challenging from both computational and algorithmic perspectives. In this regard, schemes that utilize a set of pre-defined order parameters to sample the relevant parts of the configuration space [L. Maragliano and E. Vanden-Eijnden, Chem. Phys. Lett. **426**, 168 (2006); J. B. Abrams and M. E. Tuckerman, J. Phys. Chem. B **112**, 15742 (2008)] have proved useful. Here, we demonstrate how these order-parameter aided temperature accelerated sampling schemes can be used within the Born-Oppenheimer and the Car-Parrinello frameworks of *ab initio* molecular dynamics to efficiently and systematically explore free energy surfaces, and search for metastable states and reaction pathways. We have used these methods to identify the metastable structures and reaction pathways in SiO_{2} and Ti. In addition, we have used the string method [W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B **66**, 052301 (2002); L. Maragliano *et al.*, J. Chem. Phys. **125**, 024106 (2006)] within the density functional theory to study the melting pathways in the high pressure cotunnite phase of SiO_{2} and the hexagonal closed packed to face centered cubic phase transition in Ti.

## I. INTRODUCTION

In different material and chemical systems, processes such as electronic excitations, defect formation, chemical reactions, and phase transitions often occur on time scales that are orders of magnitude separated from each other. In addition, the relaxation time scales associated with metastable states are often much smaller than the waiting times associated with transition between different metastable states.^{1} This disparity, between the relaxation and waiting times, stems from the stiffness of the normal modes: the softer modes (which are generally associated with transitions between metastable states) decay much slower than the relatively stiffer modes (which are generally associated with relaxations in a metastable basin). As a consequence, transitions between metastable states are often statistically rare and are generally difficult to resolve experimentally or simulate using traditional techniques. As such, a considerable amount of effort has been invested over several decades to develop computationally efficient schemes for scanning the configuration space.^{2–12}

The commonly used approaches of exploring the potential energy landscape using genetic algorithms, random searches, or predictions of phase boundaries based on quasi-harmonic free energy calculations often fail to properly account for the anharmonic contributions of the vibrational free energy and configurational entropy that affect the overall stability of a phase.^{5,7,10} To systematically incorporate such contributions, efficient free energy exploration is necessary and is generally a challenging proposition. Recently, a set of free energy sampling approaches, such as the adiabatic free energy dynamics (AFED)^{13} and its variants like d-AFED,^{12,14} crystal-AFED,^{11,12} temperature accelerated molecular dynamics (TAMD),^{12,15–19} and gentlest ascent dynamics in the space of collective variables (CVs),^{20,21} to name a few, were proposed for the discovery and thermodynamic ranking of stable and metastable states in a system.^{22} The central theme of these methods is to define a set of collective variables (CVs) that can capture the interesting transition events taking place. Ideally, we would like these CVs to correspond to the slow decaying soft normal modes of the system, but the soft normal modes are often not known *a priori*. Hence, a suitable strategy is to select CVs such that the soft normal modes are a small subset of the space spanned by these CVs. Once identified, the CVs can be used to systematically drive the system by imposing masses to the CVs that are much higher than the atomic masses, and at the same time assigning the CVs a temperature that is orders of magnitude higher than the physical temperature. To date these CV assisted exploration techniques have been successfully used to study complex transitions with empirical potentials.^{12,21,23,24}

Classical potentials can be computationally inexpensive and are routinely used to study the mechanistic details of complex transition events, but they often fail to capture the effects related to changes in chemistry or electronic structure and their transferability to different conditions is questionable. Thus, the objective of this paper is to demonstrate that these CV assisted free energy exploration ideas (i.e., d-AFED, crystal AFED) can be seamlessly coupled to *ab initio* molecular dynamics (AIMD) simulations. The popular approaches for AIMD simulations are (a) the Born-Oppenheimer molecular dynamics (BOMD)^{25} (wherein the eigenvalue problem of the electronic Hamiltonian is solved explicitly to obtain the interatomic interaction potential *on the fly* and then the ions are subsequently evolved), and (b) Car-Parrinello molecular dynamics (CPMD)^{26} (wherein the electronic wavefunctions and the ions are evolved simultaneously); the methods differ in the way the time scale separation between the electronic and the ionic degrees of freedom is treated. The goal of this paper is to propose suitable schemes to combine the existing order parameter aided sampling techniques that utilize the near-adiabatic response of the atomic degrees of freedom to the CVs, within the BOMD and CPMD approaches of AIMD simulations.

This paper is organized in the following way: in Section II, we illustrate how the order parameter aided sampling can be performed within the framework of BOMD, and Section III illustrates how this sampling can be performed using CPMD. In Section IV, we present a scheme to reconstruct the free energy surface (FES) using information extracted from the sampling trajectories. The applicability of the proposed schemes is illustrated in Sections V and VI using two examples—predicting the equilibrium phase boundaries in SiO_{2} phase diagram under extreme conditions and evaluating the metastability of the face centered cubic (FCC) phase of Ti. A specific case of exploring transition paths using the string method in AIMD is discussed in Appendix A. Another specific case of AIMD sampling when the collective variable is the $h-matrix$, i.e., the matrix formed by the lattice vectors of the simulation cell, is discussed in Appendix B. In Appendices C and D, we have discussed the convergence of the equilibrium distribution of the microscopic variables and the CVs obtained from this order parameter aided sampling scheme.

## II. TEMPERATURE ACCELERATED SAMPLING USING BORN-OPPENHEIMER MOLECULAR DYNAMICS

Let us consider a system with *N* atoms and *N _{e}* electrons. The positions of the nuclei are denoted by

**R**= {

**R**

_{1},

**R**

_{2}, …,

**R**

_{N}}, where

**R**

_{i}∈ ℝ

^{3}is the position of the nucleus

*i*. Following Kohn and Sham, the energy of this system within the Born-Oppenheimer approximation is

^{27}

where $\psi ={\psi 1x,\psi 2x,\u2026}$ are the single particle orthonormal orbitals, $nex=\u2211i\psi ix2$, *m _{e}* is the mass of an electron, the functional

*U*contains the electron-electron and electron-ion interactions, part of the electronic kinetic energy excluding the single particle electronic kinetic energy, and

*U*contains the ion-ion interactions. The Euler-Lagrange equation corresponding to Eq. (1) is

_{R} where *ϵ _{i}* are the eigenvalues and $\psi ix$ are the eigenfunctions. In BOMD, the eigenvalue problem for the Euler-Lagrange equation is solved at each step of evolution of the nuclei and the interatomic interactions for fixed positions of nuclei are obtained

*on the fly*. Correspondingly, the evolution equations for the nuclei can be obtained from the following Lagrangian:

where *m _{j}*,

*j*= 1, 2, …,

*N*are the masses of the nuclei.

At finite temperatures, the electronic density is modified to $nex=\u2211ifi\psi ix2$, where *f _{i}* is the occupancy of the eigenlevel

*i*specified by the Fermi-Dirac equation. In addition, following the finite temperature formulation of the density functional theory by Mermin,

^{28,29}the energy functional in the above equation is modified to

where *μ* is the chemical potential of the electrons. After we obtain the free energy $FKSR$, the nuclei in the BOMD approach are evolved (within the isothermal-isochoric ensemble) using the following equations of motion:

Here, the “bath” refers to the coupling to some heat bath maintained at *T _{R}*, e.g., Nosé-Hoover chains,

^{30}generalized Gaussian moment thermostating,

^{31}or a Langevin bath.

^{32–34}In Eq. (5), the force on the nuclei is obtained from the electronic free energy $FKS$, but at lower temperatures the occupancies

*f*behave like step functions and $FKS\u2192EKS$. The details of error incurred in this semi-classical evolution are elucidated in Ref. 35.

_{i}For efficient sampling of the configuration space, let us consider a set of CVs $\theta \u0304(R)=\theta 1(R),\theta 2(R),\theta 3(R),\u2026,\theta M(R)$ and a set of coarse-grained variables (CGVs) $z=z1,z2,z3,\u2026,zM$. We assume that this specific set of CVs can capture the important transition events taking place in the system and that the CVs are smooth functions of the microscopic variable **x**, i.e., the first and second order derivatives of the CVs are well defined. The free energy surface (FES) of the system within the isothermal-isochoric ensemble is^{15,21}

where *β* = 1/*k*_{B}*T _{R}*,

*T*is the physical temperature,

_{R}*Z*= ∫

_{ℝ3N}exp[ −

*β*ℱ

_{KS}]d

**R**is the configurational partition function, and

*δ*

^{(M)}(⋅) is the

*M*-dimensional Dirac-delta function. In the spirit of TAMD and d-AFED methods, a natural way to efficiently explore the free energy surface is to consider the modified potential energy $U\kappa BO$,

Using this potential, we obtain the equations of motion that govern the dynamics of **R**, **z** under an isothermal-isochoric ensemble,

Here, *T _{z}* ≫

*T*is the temperature at which the heat bath that is coupled to the CGVs is maintained and

_{R}*m*

_{zj}≫

*m*,

_{i}*i*= 1, 2, …,

*N*and

*j*= 1, 2, …,

*M*are the masses of the CGVs. Assigning masses to the CGVs that are sufficiently larger than those of the nuclei means that the CGVs evolve much slower than the nuclei. Consequently, the faster variables (i.e., the nuclei) have sufficient time to equilibrate. This means that the CGVs only feel the ensemble averaged effect of the

**R**variables. It has been shown by Rosso

*et al.*

^{13}and separately by Maragliano and Vanden-Eijnden

^{15}that as a consequence of this time scale separation, the CGVs evolve on the FES. Since

*T*≫

_{z}*T*, the CGVs can overcome all barriers and systematically explore the FES within a computationally accessible time frame.

_{R}The masses of the CGVs that govern the dynamics of **z** in Eq. (8) are fictitious parameters. To illustrate the effect of the choice of these masses on the dynamics of CGVs, let *ϵ*^{2} = *m*_{zj}/*m _{j}* and

*τ*=

*t*/

*ϵ*. Now expressing the masses of the CGVs in terms of physically relevant quantities such as the masses of the atoms (the effective mass associated with a CV is defined in Ref. 36), the evolution of

*z*in Eq. (8) is modified to

_{j} When the masses of the CGVs are very large, say, for example, *m*_{zj} → ∞ then *ϵ* → ∞, and under these conditions the **z** variables do not evolve while the **x** variables evolve under a fixed set of constraints imposed by the CGVs. On the other hand, when the CGVs are finite and *ϵ*^{2} ≫ 1, the evolution of the CGVs is much slower than the atomic variables. Thus, depending on the choice of masses of the CGVs, the space of collective variables can be efficiently explored using two different paradigms: in the first approach (protocol A), the coupled dynamics of **R**, **z** can be separated such that the CGVs are held fixed and the nuclei are evolved to obtain a sufficiently long trajectory. The CGVs are then evolved using the average of the forces on the structures obtained by evolving the nuclei, i.e.,

Thus, the time scale separation is systematically imposed so that the slow variables *feel* an ensemble averaged effect of the fast variables. By repeating this process, we can sample the free energy surface. Methods like the string method in the space of collective variables^{23} and gentlest ascent dynamics in the space of collective variables^{20,21} are based on this philosophy and can identify metastable states and saddle points on the FES *on the fly* without actually calculating the FES. These techniques fall under the category of heterogeneous multi-scale modeling (HMM) schemes wherein the information that is required to evolve the coarse-grained macroscopic variables is obtained on the fly by evolving the underlying microscopic degrees of freedom (DOFs).^{37,38} In Appendix A, we have illustrated how the string method in collective variables can be used within the framework of BOMD.

The second approach, which we will refer to as protocol B, is based on evolving both {**R**, **z**} simultaneously in a coupled manner as shown in Eq. (8), i.e., the ratios of masses of the CGVs to that of the atoms are finite quantities. In this case, the scale separation is not directly evident but is incorporated by the difference in the masses of these variables. Methods like the TAMD and d-AFED are based on protocol B.^{11–15} One particular variant of AFED that is suitable to study phase transitions is crystal-AFED^{11} in which the $h-matrix$ (the matrix formed by the lattice vectors of the simulation cell) is used as a collective variable. As shown in Appendix B, an implementation of crystal-AFED in AIMD involves a simple modification of any existing NPT routine (i.e., equations of motion in which the number of atoms *N*, temperature, and external pressure are kept constant). In Appendix C, we have analyzed the convergence of the equilibrium densities of the CGVs in sampling using protocols A and B convergence.

## III. TEMPERATURE ACCELERATED SAMPLING USING CAR-PARRINELLO MOLECULAR DYNAMICS

The order parameter aided sampling technique can also be implemented within the Car-Parrinello approach.^{26,39} As proposed by Car and Parrinello, one way to write the evolution equations for the orbitals and the ionic positions is to consider the following Lagrangian:^{26}

where $m\u0304e$ is a fictitious electron mass parameter and Λ_{ik} are the Lagrange multipliers that ensure the orthogonality constraints on the wavefunctions. By selecting a value of $m\u0304e$ that is much smaller than *m _{i}*, an adiabatic time scale separation between the two systems can be achieved. The equations of motion that govern the dynamics of the electrons and ions are

^{26}

In Eq. (12), the Born-Oppenheimer approximation is not directly imposed but the separation of time scales between the electrons and the nuclei is incorporated by the difference in the masses $m\u0304e$ and *m _{i}*. In the spirit of TAMD and d-AFED methods, we consider the modified potential energy $U\kappa CP$,

Thus, the Lagrangian in Eq. (11) is modified to

Using this Lagrangian, the equations that describe the evolution of the wavefunctions, ionic positions, and the CGVs are

where *T _{R}* is the physical temperature,

*T*is the temperature of the CGVs, and the “bath” refers to the thermostat that regulates the temperatures in the system.

_{z}The order parameter aided sampling can be implemented in CPMD in a similar fashion as the BOMD protocols A and B described previously in Section II. In the CPMD protocol A approach, the CGVs are held fixed while the atomic and electronic DOFs are allowed to evolve simultaneously in order to obtain the mean forces on the CGVs, while in the protocol B approach, all of the DOF are evolved together in a coupled manner using Eq. (15). Additional details on the application of the CPMD approach in the context of the string method and the use of the simulation cell matrix as the CGVs are provided in Appendices A and B, respectively. In addition, a general discussion of the use of the CPMD protocols A and B for a simple model system is presented in Appendix D.

## IV. FREE ENERGY SURFACE RECONSTRUCTION

A free energy surface provides a pictorial representation of the metastable, stable and unstable states in a system at physically relevant conditions from which the probable reactive pathways connecting these states can be calculated. The data from sampling trajectories can be used to reconstruct the free energy surface by accumulating the CGVs in a multidimensional histogram. Since the temperature at which the CGVs evolve is generally a few orders of magnitude higher than the physical temperature of the system, a FES with reasonable accuracy from histogram binning can only be obtained when the equilibrium density is converged to a very high accuracy. Consequently, the computational cost of FES construction from histogram binning increases significantly with an increase in the number of CGVs.

For this reason, we propose a modification to the single sweep^{40} (SS) method that treats the FES reconstruction as an inverse problem. The SS method involves a variational reconstruction of the FES and, in principle, is not limited by system dimensions. The method follows two simple steps: (a) first, a sufficiently long TAMD or d-AFED trajectory is obtained and a set of centers $C={z1,z2,\u2026,zK}$, in the space of *m* collective variables, is chosen along this trajectory such that the Euclidean distance *d _{ij}* between any two centers

*i*,

*j*is more than a cutoff value, and then (b) the free energy is expressed as a superposition of Gaussian radial basis functions (RBFs) placed at the centers

**z**

_{j}, i.e.,

where | ⋅ | is the *L*_{2} norm in ℝ^{m}, 𝒜 : Ω → ℝ, Ω ∈ ℝ^{m} is a bounded open set, *ϕ*_{σ}(*r*) = *e*^{−r/2σ2}, and *σ* is the width of the Gaussian. Optimal values of the coefficients *a _{i}* are obtained by optimizing the cost function $EA$ which is the standard mean squared error,

^{40}

where **f**_{j}, the mean force on the CGV at **z**_{j}, is the raw data used to fit the FES, *K* is the number of centers **z**_{j} at which **f**_{j} is known, and the gradient of the free energy with respect to **z** is

The cost function $E$ is convex with respect to **a**, and for a given *σ*, the optimal coefficients are obtained from $a\sigma *=argminaE\sigma ,a$. The optimal width of the Gaussian RBFs is obtained from the following optimization:

Substituting *σ*^{*}, **a**^{*} back into Eq. (16), we recover the FES. Being an inverse problem, the FES reconstruction is generally ill-posed because the input data are not sufficient to uniquely determine the coefficients. A simple way forward is to restrict the class of admissible solutions and to provide a method for ranking the possible solutions. In this regard, constraints such as smoothness of the solution have been introduced to act as a stabilizing term in the regularization theory pioneered by Tikhonov and others.^{41} Tikhonov regularization involves selecting a solution set that minimizes |**a**|^{2}, i.e., the *L*_{2} norm of the RBF expansion coefficients. Thus, the cost function in (17) is modified to

which, upon substituting the functional form of $EA$ from Eq. (16), becomes

We have used regularization along with multi-fold cross-validation to fine tune the model predictions.^{42,43} Let us assume that *V _{i}*,

*T*are the validation and model training data subsets, respectively, obtained from the partitioning of the raw data into

_{i}*p*-subsets for the multi-fold cross-validation. The subsets

*V*,

_{i}*T*satisfy the condition that $Vi\u222aTi=C$,

_{i}*i*∈ [1,

*p*]. Let us further assume that for fixed values of the parameters

*λ*,

*σ*,

**a**

_{σλ}solves the following minimization:

These coefficients correspond to predictions based on one training data subset, namely, *T _{i}*. The remaining data, i.e.,

*V*, are then used to validate these predictions and calculate the validation error $Era\sigma \lambda ,z\u2208Vi,\sigma ,\lambda $. This process of training and validation is repeated

_{i}*p*-times, so that each data set is used at least once for training. Thus, the cross-validation estimate of the total prediction error is the sum of the validation error for all the validation sets, i.e.,

Now, the optimum values of the regularization parameter *λ*^{*} and the optimum width *σ*^{*} of the Gaussian RBFs are obtained by minimizing $\Omega \u0304\sigma ,\lambda $. Finally, the cross-validation estimate of the RBF coefficients is given by

Substituting *σ*^{*}, **a**^{*} back into Eq. (16), we recover the FES.

## V. APPLICATION: DETERMINATION OF THE MELTING POINT OF COTUNNITE

At ambient conditions, SiO_{2} exists as stishovite and with an increase in pressure, it transforms to a CaCl_{2}-type structure followed by a transformation to an *α*-PbO_{2} type structure which persists to approximately 200 GPa in pressure.^{44} Above this pressure, the pyrite structure becomes the stable phase and at pressures above 600 GPa, the cotunnite phase is stable.^{44} Phase transitions in SiO_{2} at ultra-high pressure conditions are relevant to our understanding of the structure, formation, and evolution of giant planets.^{45} Here, we report our theoretical predictions of the melting line of high pressure cotunnite phase of SiO_{2} and discuss the relevance of the results in the context of the core structure of large planets. We have used sampling techniques to obtain an initial guess path in the CV space joining the solid and the liquid basins. Subsequently, this guess path was evolved using the string method in the CVG space to obtain the melting trajectories from which the Gibbs free energy differences between the solid and liquid phases were obtained.

Our *ab initio* calculations were performed using plane wave basis density functional theory (DFT) implementation in Vienna *Ab initio* Simulation Package (VASP).^{29,46,47} We have used the PBE exchange correlation functional and the projector augmented wave (PAW) method pseudo-potentials for the Si ([Ne]3*s*^{2}3*p*^{2}) and O ([He]2*s*^{2}2*p*^{4}). The wave functions were expanded to a plane wave cutoff of 640 eV and Brillouin zone integrations were performed at the Γ-point. We have used a supercell with 48 SiO_{2} units with supercell axes along (100), (010), (001).

### A. Selection of order parameters

While analyzing the systems that exhibit a solid-solid or a solid-liquid transition, it is often helpful to distinguish between particles belonging to different phases using the bond orientation order parameters proposed by Steinhardt *et al.*^{48} For an atom *i*, let us define

Here, $Nbi$ is the number of nearest neighbors of particle *i*, *l* is an integer parameter, $m=\u2212l,\u2212l\u22121,\u2026,l\u22121,l$, **R**_{ij} is the vector from particle *i* to *j*, and $YlmRij$ are the spherical harmonics. The orientation order parameter is^{48}

An ideal simple cubic lattice has (*Q*_{4}, *Q*_{6}) = (0.764, 0.354), an ideal body-centered cubic (BCC) lattice has (*Q*_{4}, *Q*_{6}) = (0.036, 0.511), an ideal face centered cubic lattice has (*Q*_{4}, *Q*_{6}) = (0.191, 0.574), and an ideal hexagonal closed packed (HCP) lattice has (*Q*_{4}, *Q*_{6}) = (0.097, 0.485).^{49} In practice, the two dimensional parameter space defined by the *Q*_{4}–*Q*_{6} plane is useful for mapping out the different phases in a crystal.^{50} Consequently, we have used the Steinhardt order parameters *Q*_{6}, *Q*_{4} that can quantify the orientational order, along with the volume (*V*) of the sample as the collective variables.^{24}

### B. Selection of sampling parameters

We have explored the high temperature-pressure configuration space of SiO_{2} using *ab initio* crystal-AFED simulations (described in Appendix B) at 500, 600, 700, and 800 GPa. In addition, for sake of comparison, we also used d-AFED within BOMD to sample the free energy landscape using the volume and *Q*_{6} as CVs. The parameter *W* in crystal-AFED was set to *W* = 10^{4}*W*_{0} ≫ *W*_{0}, where $W0=NkBT/\gamma h2$ is the usual value of the weight-like parameter in the Martyna-Tobias-Klein equations of motion. The time step of integration was set to 0.40 fs. We have used Langevin thermostats for the ionic DOF and the $h-matrix$. The friction coefficients for the ions were set to *γ*_{Si}, *γ _{O}* = 15 ps

^{−1}and that of the $h-matrix$ was set to

*γ*

_{h}= 2 ps

^{−1}.

An issue that is central to the sampling schemes described in Sections II and III and *ab initio* crystal-AFED, d-AFED is selecting the temperature of evolution of the CGVs. Within a NPT simulation, the relevant time scales are $\gamma h\u22121$, $\gamma Si\u22121$, $\gamma O\u22121$. In particular, $\gamma h\u22121>\gamma Si\u22121$, $\gamma O\u22121$, which is a manifestation of the fact that the relaxations in the $h-matrix$ are slower than the relaxation time scales associated with the internal degrees of freedom. The CGVs should be maintained at a temperature (*T*_{h}) such that the simulation time needed to observe a transition event is sufficiently longer than $\gamma h\u22121$. This ensures that the sampling trajectories are not polluted by any temperature related relaxations in the $h-matrix$. In our simulations, the $h-matrix$ was maintained at 4 × 10^{6} K and the transition from solid to liquid happens in about 2 × 10^{4} steps (which is much greater than $\gamma h\u22121$).

The bond orientational order parameters *Q*_{4}, *Q*_{6} were calculated only for the Si sub-lattice. The values of *r*_{min}, *r*_{max}, obtained from the radial distribution function of Si, were set to 3.62 Å and 3.87 Å, respectively. To account for the electronic entropy, we have used the Fermi smearing with the electronic temperature set to the physical temperature. For d-AFED simulations in AIMD, the coupling constants *κ* for *Q*_{6} and *V* were set to 3000 eV and 2.50 eV/Å^{6}, respectively, and the temperature of evolution of the CVs was set to 4 × 10^{6} K. These parameters satisfy the conditions obtained in Appendices C and D.

Fig. 1 shows the distribution of structures along crystal-AFED and d-AFED sampling trajectories projected on the *Q*_{6}-*V* CV space at 600 GPa and 8600 K. The initial structure for these sampling simulations is a sample with 48 SiO_{2} molecules in the cotunnite phase that has been thermalized at 600 GPa and 8600 K for 4 ps using the isothermal-isochoric ensemble.

### C. Details of string method simulations in the space of collective variables

The string method simulations (described in Appendix A) were performed at 500, 600, 700, and 800 GPa using *Q*_{6} and *V* as the collective variables. The initial guess paths for the string method simulations were obtained from *ab initio* crystal-AFED sampling trajectories. The string joining the solid and liquid basins consists of 18 discretized configurations including the end points. Restrained simulations, at each discretized configuration along the string, were performed for 15 000 steps using a time step of 0.40 fs to calculate the forces on the CGVs. The coupling constants *κ* for the CVs *Q*_{6} and *V* were set to 3000 eV and 2.50 eV/Å^{6}, respectively. These coupling constants were selected such that perturbations in a CGV lead to a change in the respective CV. Physically, this means that the coupling constants should be strong enough to initiate a change in the microscopic system so that the CV follows the respective CGV closely.

At each of iteration of the string, after calculating the mean forces acting on the CGVs, the volume CV was made dimensionless by dividing it with the equilibrium volume of the cotunnite phase. Next, the position of the string was updated in the space of the CGVs using a dimensionless time step of 0.005. For simplicity, we have set the correlation matrix *M* to be an identity matrix. The Gibbs free energy $G\alpha $ along the discretized string is calculated using the trapezoidal rule,

where *α _{k}* ∈ [0, 1] is the scalar parameter used to discretize the string, and $g\alpha k$ is the mean force on the CGV $z\alpha k$. We have taken the normalized arc length of the discretized string as the scalar parameter

*α*.

Fig. 2(a) shows the Gibbs free energy profiles at 600 GPa at 8500, 8700, and 8950 K. At 8950 and 8700 K, the liquid is the stable phase while at 8500 K, the cotunnite phase is the stable phase. Fig. 2(b) shows the Gibbs free energy difference between the liquid and solid basins Δ*G*_{ls} as a function of temperature. The thermodynamic definition of the melting point is the temperature at which the solid and liquid phases have same Gibbs free energy, so at 600 GPa we find that Δ*G*_{ls} = 0 at ∼8543 K. Generally, melting of a solid is accompanied by a large change in volume and density. However, for the melting of cotunnite at 600 GPa, the volume changes from ∼613 Å^{3} to ∼626 Å^{3} in the liquid phase, a marginal volume increase of about 2% relative to the solid phase.

### D. Discussion of systematic and statistical errors

Sampling on a free energy surface is significantly different from sampling on a potential energy surface because in the latter case, the calculated quantities are generally deterministic while in the former scenario the relevant quantities are calculated from ensemble averages. In addition, exploring the FES within the framework of DFT involves two separate adiabatic time scale separations and the effect of approximations made in calculating the free energy surfaces has to be accounted for. The uncertainties in the free energy calculations stem from systematic error, due to certain approximations made, and statistical error that comes from finite length of simulations used to evaluate the derivatives of the free energy.^{23} The statistical error can often be controlled by obtaining longer MD trajectories while the systematic error associated with the choice of the coupling constant that connects the CVs to the CGVs can be controlled by using larger values of the coupling constants.^{21}

The error in the Gibbs free energy along the melting trajectories (see Fig. 2), obtained from the string method, is calculated by accumulating the error in estimation of the mean forces on the CGVs,

where var $g\alpha k$ is the squared error in the estimation of the gradient of the free energy at $z\alpha k$. The uncertainty in the mean force calculated from a molecular dynamics trajectory of finite length is estimated using the protocol proposed by Maragliano *et al.* in Ref. 23. The convergence of the free energy barrier to melting and that of the mean force on a CGV at the cotunnite basin at a representative temperature and pressure of 8500 K and 600 GPa, respectively, are shown in Fig. 3.

The melting temperature at each pressure obtained from our analysis is used to obtain the equilibrium phase diagram of SiO_{2} in Fig. 4. Our predicted melting line is steeper than the extrapolations used in Ref. 44. The melting line in Fig. 4 suggests that at pressures corresponding to the core of planets like Uranus, Neptune, and Saturn, SiO_{2} is likely to be in the solid phase. Furthermore, the proposed pyrite-cotunnite-liquid triple point seems to be shifted to higher pressures. Our analysis of the high pressure melting line in SiO_{2} is only the first step in better understanding of the structure of giant planets. The melting mechanisms of SiO_{2} are not clearly understood: questions like whether the silicon or oxygen sub-lattices melt before the other,^{51,52} or whether both the sub-lattices melt together remain unanswered. We plan to analyze these issues in our future work.

## VI. APPLICATION: STABILITY OF THE FACE CENTERED CUBIC PHASE IN BULK Ti

Due to its excellent strength-to-weight ratio, formability, and corrosion resistance, titanium and its alloys have garnered much attention in a variety of industries where advanced mechanical properties are essential, notably in the transportation industry and in medical devices.^{53} Titanium under room temperature and pressure conditions exists in the *α*-phase which is a hexagonal closed packed (HCP) structure and transforms to a body centered cubic (BCC) structure at temperatures above 1150 K.^{53} Additionally, under pressure *α*-Ti transforms to a hexagonal phase. In the past few years, experimental studies have reported the existence of a face centered cubic (FCC) phase in thin films and nanoparticles of Ti obtained from mechanical attrition.^{54–56} Though some experimental studies have attributed the existence of the FCC phase to the presence of impurities like hydrogen that can lead to the formation of titanium hydrides, an analysis of stability or metastability of this phase is lacking. From the stability point of view, it is possible that there exists a FCC basin in bulk Ti, even though HCP is the stable phase at low temperatures. To understand the metastability of the FCC phase diagram in bulk Ti, we have used the FES exploration scheme to first explore the configuration space for probable metastable structures and then to obtain the FES itself. In addition, we have used the string method to gain insights into the phase transition pathways, mechanisms, and the free energy barrier between the HCP and FCC phases.

Our free energy surface explorations were again performed using the Vienna *Ab initio* Simulation Package (VASP).^{29,46,47} We have used the PBE exchange correlation functional and projector augmented wave (PAW) method pseudo-potentials for the Ti ([Ar]3*d*^{2}4*s*^{2}) along with a plane-wave cutoff of 400 eV. The simulation cell contained 108 Ti atoms with 12 basal planes stacked along the [0001] which is commensurate with both ABABAB… and ABCABC stacking sequences in HCP and FCC phases, respectively. The free energy surface exploration was performed using the $h-matrix$ as the CV by maintaining the atoms at 500 K and 1 atm and setting the temperature of the $h-matrix$ to 2 × 10^{7} K. The time step of integration was 1 fs and the friction coefficients were set to *γ _{Ti}* = 15 ps

^{−1},

*γ*

_{h}= 2 ps

^{−1}. We had initiated two sampling trajectories: one from the FCC structure and another from the HCP structure. Both these initial structures were equilibrated at 500 K and 1 atm for 10 ps within the isothermal-isobaric ensemble. We found that the majority of the sampling trajectories seem to cluster around the HCP basin while very few trajectories cluster around the FCC basin, suggesting that either FCC is unstable or only a shallow local minimum.

Next, we have used structures from the FES sampling trajectories to obtain a 2-dimensional grid in the *Q*_{6}-*Q*_{4} collective variable space. The grid spacings along the *Q*_{6} and *Q*_{4} axes were 0.04. At each point on the grid, the mean forces on the CGVs were calculated from a 12 ps long restrained NPT trajectory at 500 K and 1 atm pressure with coupling constants set to 3000 eV for *Q*_{4} and 2000 eV for *Q*_{6}. Fig. 5(a) shows the reconstructed Gibbs free energy surface at 500 K for bulk Ti. In agreement with the phase diagram of Ti, we see that the stable phase in bulk Ti is the HCP phase and FCC is clearly unstable.

Even though the FES in Fig. 5(a) shows that there is no metastable state corresponding to FCC phase, using a molecular dynamics simulations we find that the structure does not spontaneously transform to the HCP phase at 500 K and 1 atm even in 15 ps while at temperatures above 750 K, a FCC structure spontaneously transforms to a HCP structure in less than 4 ps. To further verify the stability of the FCC phase, we resorted to string method calculations for which the initial path was obtained from a sampling trajectory using the **h** as the CV. In agreement with the reconstructed FES, our results at 500 K and 1000 K (Fig. 5(b)) indicate that the FCC phase is unstable in bulk Ti.

Since our calculations show that the FCC phase is unstable in bulk Ti, an important question is what makes the FCC phase stable in thin films or nanostructures. While it is plausible that surface effects play a dominant role, the presence of local stress fields due to dislocations, impurities, or solute atoms can also significantly affect the phase stability in Ti. We plan to study some of these issues in our future work.

## VII. CONCLUSIONS

In this paper, we have outlined how the CV assisted free energy surface exploration schemes can be implemented within the Car-Parrinello or the Born-Oppenheimer framework and have demonstrated the use of the d-AFED, crystal-AFED sampling schemes by using internal order parameters (*Q*_{6}) and extensive quantities (*V*, $h-matrix$) as collective variables to explore relevant parts of the configuration space in SiO_{2} and Ti. The stability of the high pressure cotunnite phase of SiO_{2} was predicted based on zero temperature quasi-harmonic free energy calculations and our results from free energy sampling at finite temperatures clearly support this conclusion. In addition, using finite temperature AIMD calculations, we have also calculated the phase boundaries in this regime of the (*P*-*T*) phase diagram that had not been studied before. Similarly, in case of Ti we have systematically studied the metastability of the FCC phase in bulk Ti. Our FES and string method results show that the FCC phase is not metastable in bulk Ti.

## Acknowledgments

The authors are grateful to Liang Qi and Qian Yu for motivating us to study the stability of FCC phase in bulk Ti and to Sebastian Hamel for motivating us to study the cotunnite phase of SiO_{2}. A.S. wishes to thank Professor Weinan E, Professor Mark Tuckerman, Professor Giovanni Ciccotti, Professor Eric Vanden-Eijnden, and Professor Sara Bonella for their comments. In addition, we want to thank Alfredo Correa and Xavier Andrade-Valencia for their help in improving the manuscript. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344 and was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.

### APPENDIX A: MINIMUM FREE ENERGY PATH USING THE STRING METHOD IN THE COLLECTIVE VARIABLE SPACE

The string method is an efficient way to calculate the maximum likelihood path associated with a rare transition event.^{21,57} On the free energy surface, the idea of searching for minimum free energy paths (MFEPs) can be implemented within the framework of heterogeneous multi-scale modeling (HMM),^{37,38} wherein one starts with a macro-scale solver (the string method in this case) and the values of the quantities, such as the forces on the CGVs, required by the macro-solver are then computed on the fly using a (constrained or restrained) micro-scale solver.^{23} The explicit form of macro-solver (i.e., the string method) that we have used is slightly different from the original version of the string method.^{58,61}

To obtain a MFEP in the space of CGVs, we consider a curve parametrized by $\gamma ={z\alpha :\alpha \u2208[0,1]}$. The evolution of the curve using the string method is given by

where $z\lambda \alpha n$ is the CGV along the discretized string at *α _{n}*, $P\mu \nu n$ is a projection matrix, $M\nu \lambda n$ is the correlation matrix that accounts for the projection of the MFEP in the configuration space to the CGV space, and $\u2212\u2202Az\alpha ,t/\u2202z\lambda n$ is the mean force on the collective variable. Thus, the implementation of string method in CGV space requires calculation of the forces on CGVs, i.e., the gradient of the free energy for which we use restrained MD simulations as the micro-scale solver using the potential $U\kappa CP$ (in CPMD approach) defined in Eq. (13). The equations of motion for the variables {

*ψ*,

**R**} for fixed

**z**are given by

where *m _{i}* is the mass of nuclei

*i*and

*T*is the physical temperature. The gradient of the free energy is obtained by the average of the instantaneous forces on the CGVs in a sufficiently long trajectory obtained by evolving the above equations. As discussed in Section II, another approach is to solve the eigenvalue problem for the electronic part of the Kohn-Sham Hamiltonian at each step of evolution of the ionic positions and obtain the interatomic interaction potential ($FKS$)

_{R}*on the fly*that is then used to evolve the ions by keeping the CGVs fixed,

The free energy profile along the string parametrized by the scalar parameter *α* ∈ [0, 1], such that the end points correspond to *α* = 0, 1, is given by the following line integral:

In the limit of large {*κ _{i}*}, the error in the above estimate of the free energy is determined by the physical temperature of the ensemble and the discretization of the string.

### APPENDIX B: *AB INITIO* SAMPLING USING THE SIMULATION CELL MATRIX AS THE COLLECTIVE VARIABLE

Structural phase transitions often involve changes in the atomic symmetry and local density. Consequently, such transformations can often be captured by the supercell edge lengths, angles, or the elements of the full supercell matrix, as the collective variables. However, structural changes that are controlled by internal order parameters lie on a hypersurface that is orthonormal to these CVs and, hence, are not captured by using the full cell matrix as CV. Nevertheless, for complex systems, using the cell matrix as CV can provide valuable information about the probable metastable states and transition pathways. The crystal-AFED technique introduced by Yu and Tuckerman,^{11} which is an isobaric variant of d-AFED designed to sample phase transitions, is in line with this philosophy.

To understand how crystal-AFED works and how it can be extended to sample metastable states within the framework of DFT, we start with an ensemble that allows anisotropic cell fluctuations and consider the following partition function:^{59}

where *V* is the volume of the system, **h**_{0} captures the different possible shapes of a supercell at constant volume, $Q2$ is the partition function for a fixed supercell size and shape, $Q1$ is the partition function for a fixed volume of the supercell, and *P*_{ext} is the external pressure. Let **h** be the matrix whose columns are the three vectors describing the shape of the supercell, then the partition function in Eq. (B1) becomes^{59}

If the *N* atoms in the system have momenta **p** ∈ ℝ^{3N} and masses *m*_{1}, …, *m _{N}*, then the Martyna-Tobias-Klein equations of motion within the framework of BOMD that generate the ensemble in Eq. (B2) are

^{59}

where **I** is the identity matrix, *N*_{f} = 3*N* is the number of degrees of freedom in the system, *W* = (*N*_{f}/3 + 1) *kT _{R}τ*

^{2}is a mass-like parameter,

*τ*determines the time scale of the barostat motion, and $P(int)=det(h)\u22121\u2211i=1Npi\u2297pi/mi+Ri\u2297Fi$ is the internal pressure tensor. The evolution of

**p**drives the fluctuations of the internal pressure tensor estimator such that 〈

_{g}**P**

^{(int)}〉 =

*P*

**I**. The conserved quantity for the dynamics in Eq. (B3) is

^{59}

In the limit *W* → ∞, from Eq. (B3a) we see that $h\u0307\u21920$. This means that the effect of **p**_{g} becomes negligible and $Pextdeth$ merely becomes a constant term in Eq. (B4). Under these conditions, the dynamics in Eqs. (B3) reduces to

which samples the partition function for a fixed **h**, i.e., the partition function $Q2$. Now, if instead of taking *W* to infinity, we assign a very large value to *W*, i.e., *W* ≫ (*N*_{f}/3 + 1) *kT _{R}τ*

^{2}, then

**h**evolves very slowly while the ionic DOF relaxes very fast and following a similar logic as described in Sections II and III, the effective dynamics of

**h**explores the free energy surface. All that remains to be done is to assign a very high temperature to the

**p**

_{g}variables so that the FES exploration is tractable and can be performed efficiently in a finite amount of computer time. Recently, it was shown that the internal order parameters can also be incorporated into this sampling scheme.

^{12}

This temperature accelerated sampling scheme can be seamlessly integrated within the Car-Parrinello framework as

where

is the force on the ion *i*.

### APPENDIX C: ANALYSIS OF SAMPLING USING PROTOCOLS A AND B WITHIN BOMD

Here we illustrate that the equilibrium density obtained from sampling within protocol B converges to the true solution as the ratio of the masses of the fast and slow variables, i.e., *m _{i}*/

*m*

_{zj}in Eq. (8), goes to zero. Let us consider a model system defined by a potential $Vx$. We define an order parameter $y=\theta x$ and a modified potential energy surface $U\kappa x,y=Vx+\kappa \u2009y\u2212\theta x2/2$. In the case of BOMD, the potential $Vx$ is replaced by $FKSx$. According to the protocol B, the coupled equations of motion of the variables

**x**and

*y*within the Langevin settings are

where *β _{x}* is the inverse of the physical temperature,

*β*≪

_{y}*β*is the inverse of the temperature at which

_{x}*y*is maintained,

*γ*,

_{x}*γ*are the friction coefficients,

_{y}*η*,

_{x}*η*are

_{y}*δ*-correlated stationary Gaussian processes with zero mean, and

*m*≫

_{y}*m*are the masses. Our goal is to explore the following free energy in

_{x}*y*within isothermal-isochoric ensemble:

In the above, by approximating the Dirac-delta function by a Gaussian distribution, we obtain

Here, we have neglected an additive constant that depends on the coupling constant and the temperature. The interchange between the limit *κ* → ∞ and the spatial integral has been discussed in Ref. 23. The *κ* dependent free energy is related to the true free energy in Eq. (C2) by

where $Qy=\beta x\u2207yAy2\u2212\Delta yAy$.

The equilibrium distribution, $\rho x,y$, of the fast and the slow variables can be obtained from the solution of the Fokker-Planck equation for the stochastic processes whose dynamics is governed by the overdamped limit of Eq. (C1),

where $\u03f51=\gamma xmx/\gamma ymy$ is a small parameter, $Lx\rho x,y=\u2207x\u2207xU\kappa \rho x,y+\beta x\u22121\Delta x\rho x,y$, and $Ly\rho =\u2207y\u2207yU\kappa \rho x,y+\beta y\u22121\Delta y\rho x,y$. We limit ourselves to solutions of the following class where the equilibrium density is expressed as an expansion involving different powers of *ϵ*_{1}:^{15}

Substituting this ansatz into Eq. (C5), we obtain the following set of partial differential equations:^{15}

*ρ*

_{0}in the second differential equation, we obtain $Wx,y=Ly\rho 0x,y$, where

Thus, the first order correction to $\rho 0x,y$ is the solution of $Lx\rho 1x,y+Wx,y=0$. An analytical solution is probably not available for a multidimensional potential energy landscape, but for a single dimensional potential, the solution takes the following form:

where *x*_{0} is some reference point. In principle, the higher order correction can be obtained by solving the other partial differential equations in Eq. (C7), but they are difficult to solve analytically for a general high dimensional system. An important point to note is that the form of *ρ*_{1} illustrates the complex coupling between the fast and the slow variables and that the higher order corrections depend on the non-local structure of the potential energy landscape.

If we assume that a function *ϕ*_{1} exists for a general multidimensional system (that exhibits time scale separation between the different DOFs), then the above discussion shows that the equilibrium distribution obtained by evolving the coupled dynamics of **x**, *y* in Eq. (C1) is given by

Thus, the free energy in the space of CGV *y* is given by

Since *ρ*_{1} is independent of *ϵ*_{1} and assuming that the coefficients of the higher order terms are finite and bounded, then as *m _{x}γ_{x}*/

*m*→ 0, $A\kappa ,\u03f51\u2192A\kappa $.

_{y}γ_{y}Let us now analyze the dynamics of the CGV *y* in protocol A. Within protocol A, the information required to evolve *y* is obtained by first evolving **x** by keeping *y* fixed,

The equilibrium density generated from such an evolution is *ρ* = *e*^{−βxUκ}/*Z*_{o}, *Z*_{o} = ∫*e*^{−βxUκ}*d***x**. The evolution of *y* is given by

which samples free energy surface $A\kappa $ in Eq. (C3). Thus, both protocols A and B can be used to explore the free energy surface.

### APPENDIX D: ANALYSIS OF SAMPLING USING PROTOCOLS A AND B WITHIN CPMD

To properly understand the sampling capabilities of Eq. (15), let us first consider the Hamiltonian of a simple model system that exhibits two levels of adiabatic time scale separations

where *p _{i}*,

*m*are the momentum and mass of the variable

_{i}*i*= {

*x*,

*y*,

*z*}. The potentials $Vx,y$ and $Wy,z=\kappa \u2009z\u2212\theta y2/2$ define the interactions between the variables. This form of interaction is significant for a system of interacting electrons, ions, and coarse-grained variables: in the potential

*U*

_{κ}, the coupling between the CGVs and the ionic variables (thorough the CVs) does not depend on the electronic DOF while the interaction potential between the electrons and the ions is independent of the CGVs. The analysis presented below is an extension of the approach used by Maragliano and Vanden-Eijnden

^{15}and is similar to the Liouville operator based analysis presented by Marx

*et al.*for the

*ab initio*centroid molecular dynamics scheme.

^{60}

We seek to explore the free energy surface $A\kappa z$ of the slowest variable *z*,

According to the protocol B, the coupled equations of motion of the variables *x*, *y*, and *z* within the Langevin settings are

where *β* is the inverse temperature, *γ _{x}*,

*γ*,

_{y}*γ*are the friction coefficients, and

_{z}*η*,

_{x}*η*,

_{y}*η*are

_{z}*δ*-correlated stationary Gaussian processes. Since

*m*≫

_{z}*m*≫

_{y}*m*,

_{x}*x*evolves faster than

*y*,

*z*. The equilibrium distribution can be obtained from the Fokker-Planck equation which in the limit of over-damped Langevin dynamics takes the following form:

where *ϵ*_{2} = *γ _{y}m_{y}*/

*γ*,

_{z}m_{z}*ϵ*

_{1}=

*γ*/

_{x}m_{x}*γ*and as a consequence of the difference in masses, we have

_{y}m_{y}*ϵ*

_{1},

*ϵ*

_{2}≪ 1. As in case of Eq. (C5), the operators $Lx$, $Ly$, and $Lz$ correspond to the dynamics of

*x*,

*y*, and

*z*, respectively. We limit ourselves to solutions of the type $\rho =\rho 0+\u03f51\rho 1+\u03f512\rho 2+\cdots $. Substituting this ansatz into Eq. (D4), we obtain the following set of partial differential equations:

If we take $\rho 0x,y,z=e\u2212\beta Vx,y+Wy,z$, then the first order correction $\rho 1x,y,z$ is obtained from $Lx\rho 1+W0x,y,z=0$, where, $W0x,y,z=Ly\rho 0x,y,z+\u03f52Lz\rho 0x,y,z$. Similar to the discussion in Section II, the solution is given by

Thus, the equilibrium distribution sampled in the overdamped limit of the coupled dynamics in Eq. (D3) is given by $\rho x,y,z=\rho 0x,y,z+\u03f51\rho 1x,y,z+O\u03f512$. Consequently, protocol B provides an effective strategy to explore the FES in a computationally accessible time frame.