Exploring the free energy surface using ab initio molecular dynamics

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.¹ 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)¹³ 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.²² 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)²⁵ (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)²⁶ (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₂ 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.

Let us consider a system with N atoms and N_e electrons. The positions of the nuclei are denoted by R = {R₁, R₂, …, R_N}, where R_i ∈ ℝ³ is the position of the nucleus i. Following Kohn and Sham, the energy of this system within the Born-Oppenheimer approximation is²⁷ 

where $ψ={ψ1x,ψ2x,…}$ are the single particle orthonormal orbitals, $nex=∑iψ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_R contains the ion-ion interactions. The Euler-Lagrange equation corresponding to Eq. (1) is

where ϵ_i are the eigenvalues and $ψ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=∑ifiψ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,³⁰ generalized Gaussian moment thermostating,³¹ 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_i behave like step functions and $FKS→EKS$⁠. The details of error incurred in this semi-classical evolution are elucidated in Ref. 35.

For efficient sampling of the configuration space, let us consider a set of CVs $θ̄(R)=θ1(R),θ2(R),θ3(R),…,θM(R)$ and a set of coarse-grained variables (CGVs) $z=z1,z2,z3,…,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_BT_R, T_R is the physical temperature, Z = ∫_ℝ^3Nexp[ − βℱ_KS]dR 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κ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_R is the temperature at which the heat bath that is coupled to the CGVs is maintained and 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.¹³ and separately by Maragliano and Vanden-Eijnden¹⁵ that as a consequence of this time scale separation, the CGVs evolve on the FES. Since T_z ≫ T_R, the CGVs can overcome all barriers and systematically explore the FES within a computationally accessible time frame.

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 ϵ² = 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_j in Eq. (8) is modified to

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 ϵ² ≫ 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²³ 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¹¹ 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.

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:²⁶ 

where $m̄e$ 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̄e$ 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²⁶ 

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̄e$ and m_i. In the spirit of TAMD and d-AFED methods, we consider the modified potential energy $Uκ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_z is the temperature of the CGVs, and the “bath” refers to the thermostat that regulates the temperatures in the system.

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.

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⁴⁰ (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,…,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₂ 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,⁴⁰ 

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σ*=argminaEσ,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.⁴¹ Tikhonov regularization involves selecting a solution set that minimizes |a|², i.e., the L₂ 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_i are the validation and model training data subsets, respectively, obtained from the partitioning of the raw data into p-subsets for the multi-fold cross-validation. The subsets V_i, T_i satisfy the condition that $Vi∪Ti=C$⁠, 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_i, are then used to validate these predictions and calculate the validation error $Eraσλ,z∈Vi,σ,λ$⁠. This process of training and validation is repeated 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 $Ω̄σ,λ$⁠. Finally, the cross-validation estimate of the RBF coefficients is given by

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

At ambient conditions, SiO₂ exists as stishovite and with an increase in pressure, it transforms to a CaCl₂-type structure followed by a transformation to an α-PbO₂ type structure which persists to approximately 200 GPa in pressure.⁴⁴ Above this pressure, the pyrite structure becomes the stable phase and at pressures above 600 GPa, the cotunnite phase is stable.⁴⁴ Phase transitions in SiO₂ at ultra-high pressure conditions are relevant to our understanding of the structure, formation, and evolution of giant planets.⁴⁵ Here, we report our theoretical predictions of the melting line of high pressure cotunnite phase of SiO₂ 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]3s²3p²) and O ([He]2s²2p⁴). 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₂ units with supercell axes along (100), (010), (001).

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.⁴⁸ For an atom i, let us define

Here, $Nbi$ is the number of nearest neighbors of particle i, l is an integer parameter, $m=−l,−l−1,…,l−1,l$⁠, R_ij is the vector from particle i to j, and $YlmRij$ are the spherical harmonics. The orientation order parameter is⁴⁸ 

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

We have explored the high temperature-pressure configuration space of SiO₂ 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₆ as CVs. The parameter W in crystal-AFED was set to W = 10⁴W₀ ≫ W₀, where $W0=NkBT/γ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⁻¹ and that of the $h-matrix$ was set to γ_h = 2 ps⁻¹.

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 $γh−1$⁠, $γSi−1$⁠, $γO−1$⁠. In particular, $γh−1>γSi−1$⁠, $γO−1$⁠, 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 $γh−1$⁠. 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⁶ K and the transition from solid to liquid happens in about 2 × 10⁴ steps (which is much greater than $γh−1$⁠).

The bond orientational order parameters Q₄, Q₆ 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₆ and V were set to 3000 eV and 2.50 eV/Å⁶, respectively, and the temperature of evolution of the CVs was set to 4 × 10⁶ 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₆-V CV space at 600 GPa and 8600 K. The initial structure for these sampling simulations is a sample with 48 SiO₂ molecules in the cotunnite phase that has been thermalized at 600 GPa and 8600 K for 4 ps using the isothermal-isochoric ensemble.

The string method simulations (described in  Appendix A) were performed at 500, 600, 700, and 800 GPa using Q₆ 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₆ and V were set to 3000 eV and 2.50 eV/Å⁶, 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α$ 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αk$ is the mean force on the CGV $zα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 ų to ∼626 ų in the liquid phase, a marginal volume increase of about 2% relative to the solid phase.

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.²³ 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.²¹ 

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αk$ is the squared error in the estimation of the gradient of the free energy at $zα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₂ 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₂ 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₂ is only the first step in better understanding of the structure of giant planets. The melting mechanisms of SiO₂ 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.

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.⁵³ 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.⁵³ 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]3d²4s²) 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⁷ K. The time step of integration was 1 fs and the friction coefficients were set to γ_Ti = 15 ps⁻¹, γ_h = 2 ps⁻¹. 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₆-Q₄ collective variable space. The grid spacings along the Q₆ and Q₄ 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₄ and 2000 eV for Q₆. 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.

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₆) and extensive quantities (V, $h-matrix$⁠) as collective variables to explore relevant parts of the configuration space in SiO₂ and Ti. The stability of the high pressure cotunnite phase of SiO₂ 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.

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₂. 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.

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.²³ 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 $γ={zα:α∈[0,1]}$⁠. The evolution of the curve using the string method is given by

where $zλαn$ is the CGV along the discretized string at α_n, $Pμνn$ is a projection matrix, $Mνλn$ is the correlation matrix that accounts for the projection of the MFEP in the configuration space to the CGV space, and $−∂Azα,t/∂zλ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κ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_R 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$⁠) 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.

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,¹¹ 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:⁵⁹ 

where V is the volume of the system, h₀ 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⁵⁹ 

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

where I is the identity matrix, N_f = 3N is the number of degrees of freedom in the system, W = (N_f/3 + 1) kT_Rτ² is a mass-like parameter, τ determines the time scale of the barostat motion, and $P(int)=det(h)−1∑i=1Npi⊗pi/mi+Ri⊗Fi$ is the internal pressure tensor. The evolution of p_g drives the fluctuations of the internal pressure tensor estimator such that 〈P^(int)〉 = PI. The conserved quantity for the dynamics in Eq. (B3) is⁵⁹ 

In the limit W → ∞, from Eq. (B3a) we see that $ḣ→0$⁠. 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τ², 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.¹² 

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

where

is the force on the ion i.

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=θx$ and a modified potential energy surface $Uκx,y=Vx+κ y−θ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 ≪ β_x is the inverse of the temperature at which y is maintained, γ_x, γ_y are the friction coefficients, η_x, η_y are δ-correlated stationary Gaussian processes with zero mean, and m_y ≫ m_x are the masses. Our goal is to explore the following free energy in 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=βx∇yAy2−ΔyAy$⁠.

The equilibrium distribution, $ρ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 $ϵ1=γxmx/γymy$ is a small parameter, $Lxρx,y=∇x∇xUκρx,y+βx−1Δxρx,y$⁠, and $Lyρ=∇y∇yUκρx,y+βy−1Δyρx,y$⁠. We limit ourselves to solutions of the following class where the equilibrium density is expressed as an expansion involving different powers of ϵ₁:¹⁵ 

Substituting this ansatz into Eq. (C5), we obtain the following set of partial differential equations:¹⁵ 

Thus, the first order correction to $ρ0x,y$ is the solution of $Lxρ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₀ 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 ρ₁ 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 ϕ₁ 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 ρ₁ is independent of ϵ₁ and assuming that the coefficients of the higher order terms are finite and bounded, then as m_xγ_x/m_yγ_y → 0, $Aκ,ϵ1→Aκ$⁠.

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_κdx. The evolution of y is given by

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

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_i are the momentum and mass of the variable i = {x, y, z}. The potentials $Vx,y$ and $Wy,z=κ z−θ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¹⁵ and is similar to the Liouville operator based analysis presented by Marx et al. for the ab initio centroid molecular dynamics scheme.⁶⁰ 

We seek to explore the free energy surface $Aκ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