Efficient discovery of multiple minimum action pathways using Gaussian process
EEfficient discovery of multiple minimum action pathways using Gaussian process
JaeHwan Shim, ∗ Juyong Lee, † and Jaejun Yu ‡ Department of Physics and Astronomy, Seoul National University1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea Department of Chemistry, Division of Chemistry and Biochemistry, Kangwon National University1, Gangwondaehak-gil, Chuncheon-si, Gangwon-do 24341, Republic of Korea
We present a new efficient transition pathway search method based on the least action principleand the Gaussian process regression method. Most pathway search methods developed so far rely onstring representations, which approximate a transition pathway by a series of slowly varying replicasof a system. Since those methods require a large number of replica images, they are computationallyexpensive in general. Our approach employs the Gaussian process regression method, which takesthe Bayesian inference on the shape of a given potential energy surface with a few observed dataand Gaussian-shaped kernel functions. Based on the inferred potential, we find multiple low-actionpathways by carrying out the action optimization based on the Action-CSA (Conformational spaceannealing). Here we demonstrate a drastic elevation of computing efficiency about five orders ofmagnitude for the system with the M¨uller-Brown potential. Further, for the sake of demonstratingits real-world capabilities, we apply our method to ab initio calculations on alanine dipeptide. Theimproved efficiency of GPAO makes it possible to identify multiple transition pathways of alaninedipeptide and calculate their transition probabilities with ab initio accuracy. We are confidentthat our GPAO method is a powerful approach to investigate the mechanisms of complex chemicalreactions
I. INTRODUCTION
Finding multiple transition pathways of phase transi-tions, chemical reactions, and conformational transitionsremains challenging in physics and chemistry. In general,it takes a considerably long time to simulate an evolutionover an energy barrier between reactants and products.Indeed, for some cases with large energy barriers, it ispractically not feasible to carry out a proper simulationeven on state-of-the-art computers. To overcome sucha limitation, researchers have suggested numerous pathsampling approaches [1–6]. Among various attempts, themethods based on the fixed-boundary-value formulationare commonly used. Unlike the initial-value formula-tion, which awaits a long time to observe rare events,the fixed-boundary-value formulation presumes a virtualtrajectory and evaluates the pathways’ quality based onthe height of an energy barrier or its action value.Among many fixed boundary value approaches, sev-eral methods based on the least action principles wereproposed and considered a promising way to find sad-dle points on a potential energy surface (PES) [5, 7].Passerone and Parrinello [5] suggested a “so-called”action-derived molecular dynamics (ADMD) method,which looks for low-action pathways via the local opti-mization of a modified classical action with an energyconservation restraint term. Later, Lee and coworkers [8]developed a kinetically controlled ADMD method by in-troducing a restraining term to satisfy the equipartitiontheorem. These classical-action-based methods have an ∗ [email protected] † [email protected] ‡ [email protected] inherent limitation that the classical principle of least ac-tion is an extremal, not a minimum principle. In otherwords, the classical action can be either minimized ormaximized, which makes the computational outcome am-biguous.Another class of pathway sampling approaches uses theOnsager-Machlup (OM) action [9] to obtain the mostprobable transition pathway of a diffusive system de-scribed by the Langevin equation [7, 10–13]. For a diffu-sive system, the probability p ( x B | x A ; t ) for a state x A to become another state x B within a given time interval t , can be described as follows [9, 11, 14]: (cid:90) x B x A D x ( s ) e − i ¯ hS OM ( x ( s )) (1)where S OM is the OM action and D is the notation forintegrating all the possible pathways connecting two endstates. The probability of a transition from the state x A to x B can be calculated by summing up the probabil-ity of all pathways (cid:80) s e − S OM ( x ( s )) /k B T [9, 10, 15]. Themost dominant pathway may effectively be determinedby the probability of that state transition. However, formost cases, a state transition usually incorporates multi-ple pathways, such as multiple protein folding pathways.Lee and coworkers utilized an efficient global optimiza-tion method to identify the most dominant transitionpathway corresponding to the global minimum of OMaction and sample multiple less prevalent pathways [13].In their work, the locally optimized transition pathwayswere sampled by using the modified classical action toavoid the Hessian calculation of a potential function. Leeand Brooks tested the efficiency of the direct minimiza-tion of OM action by using the system’s Hessian infor-mation. They calculated the gradient of OM action ana-lytically and directly minimized OM action [7]. The di- a r X i v : . [ phy s i c s . c o m p - ph ] F e b rect optimization of the Onsager-Machlup action (DOA)method guarantees to locate the local minima of OM ac-tion. However, it did not preserve the total energy of thesystem.A major bottleneck of the existing pathway searchmethods is that they require a significant amount ofcomputational resources. Conventional pathway searchmethods optimize the transition pathway by calculatingthe energy and gradient of every configuration image.From this information, they obtain the optimal coor-dinates, which minimize a given action. However, thecalculation of the energies and forces for multiple im-ages increases the computational burden quite seriously.A small enough step size is necessary to generate con-tinuous and smooth trajectories. However, the use of asmaller step size obviously increases computational costs.In other words, there is a trade-off relationship betweenthe accuracy and computational efficiency of trajectories.For instance, to achieve the continuity of a trajectory,Passerone and Parrinello [5] used 300 images to find tra-jectories on a simple 2D M¨uller-Brown (MB) potential.As a way to get around the computational cost is-sue, J´onsson and coworkers [16, 17] suggested the one-image evaluation (OIE) method incorporating the Gaus-sian process (GP) [18, 19] with a nudged elastic band(NEB) [4] to reduce the number of energy and force cal-culations. The OIE method finds the minimum energypathway (MEP) by building a Gaussian PES that ap-proximates the actual PES iteratively. The followingsteps construct the approximate Gaussian PES. First,a posterior distribution is used to evaluate the most un-certain image along the pathway. Second, the most un-certain image is evaluated with actual energy calcula-tion and then added the value to a database to improvethe Gaussian PES’s accuracy. Third, on the newly con-structed Gaussian PES, the MEP is updated to fit thenew Gaussian PES using NEB. After the relaxation ofNEB, the whole process is repeated until the level of un-certainty decreases below a certain threshold. With theMB potential, the number of function calls is effectivelyreduced by two orders of magnitude [17].Another major limitation of existing pathway searchmethods, such as NEB, is that they are local optimizationmethods whose final results heavily depend on the initialguesses on pathways. Generally, conventional pathwaysearch methods start from the linear interpolation be-tween end states and perform local optimization of theinitial pathway. These approaches can find the correctminimum energy or action pathways only when a givenPES is simple and smooth. However, many complex sys-tems have highly rugged PES, which does not guaranteeto locate the most dominant pathway near the linear in-terpolation between two states. To avoid this problem,an extensive search on a pathway space is essential.In this work, to overcome these limitations, we presenta novel transition-pathway-sampling method by combin-ing the state-of-the-art machine learning techniques andthe advantages of both Passerone’s and Lee’s methods. We also carry out direct OM action optimization witha total-energy-conservation restraint without using theoriginal system’s Hessian. One of the core algorithmsin our method, Gaussian process action optimization(GPAO), shares a similar principle with the OIE method.To avoid calculations of all the intermediate images, wecalculate only the most uncertain images as measuredby the variance via GP and then add the calculation re-sult to the database. The potential energies, forces, andHessian of the intermediate images are estimated usingthe Gaussian PES. Within the Gaussian PES, the MEPis directly optimized so as to find its minimum action.Also, GPAO samples multiple sub-optimal pathways, lo-cal minima of action, because it uses the identical sam-pling method with Action-CSA [13].In this work, we also implement a routine that au-tomatically finds suitable target energy throughout theiteration and combines the GP with the conventionalADMD or DOA methods. Target energy is inherentlyhard to know prior to an actual pathway sampling be-cause proper target energy is closely related to the heightof the energy barrier of the MEP. Here we utilize a mech-anism that approximates the target energy to the maxi-mum potential energy of the temporal MEP throughoutthe iteration. As a result, the computation efficiencyof our transition-pathway-sampling method shows a dra-matic enhancement. The conventional ADMD requires,even for the simple MB example, given with suitable tar-get energy, 1,200,000 calls of force calculations to find theMEP.On the other hand, GPAO with classical action, with-out knowledge of the target energy, reduces the numberof 1,200,000 actual function calls down to 12. Encour-aged by the significant computational efficiency gain byGPAO in five orders of magnitude, we take a step for-ward to undertake a computationally challenging exam-ple: finding multiple transition paths of alanine dipeptidethrough ab initio calculations. Action-CSA’s capability,which preserves the diversity of distinct pathways whileoptimizing the OM action globally and efficiently, resultsin multiple physically accessible trajectories between thetwo end-states. To the best of our knowledge, our work isthe first report, which identifies multiple conformationaltransition pathways of alanine dipeptide with an ab initio potential. II. THEORY AND METHODA. Gaussian Process
Significant computational improvements of GPAO overthe conventional ADMD or DOA methods are attributedto the ability to quantify the variance of specific im-ages in a pathway by using the Bayesian inference. Letus suppose a pathway consisting of N images and M evaluated data points. In our GP model, the proba-bility distribution of a given string of images, X ∗ = (cid:2) x (1) , x (2) , · · · , x ( n ) , · · · , x ( N ) (cid:3) , whose energies and forcesare represented as Y ∗ = [ ∗ , ∇ V ∗ ], is given by a multivari-ate Gaussian distribution: p ( Y ∗ ; X ∗ , X , Y ) = N ( X ∗ | µ ∗ , Σ ∗ ) (2)where X = (cid:2) x (1) , x (2) , · · · , x ( m ) , · · · , x ( M ) (cid:3) are the eval-uated data points and Y = [ V, ∇ V ] correspond to theircorresponding potential energies and forces. µ ∗ and Σ ∗ are the mean and the variance of the Gaussian distribu-tion, respectively: µ ∗ = m ∗ + K T ∗ ( K + σ n I ) − ( Y − m ) (3) Σ ∗ = K ∗∗ + K T ∗ ( K + σ n I ) − K ∗ (4)where m , ( K ) mm (cid:48) = k (cid:16) x ( m ) , x ( m (cid:48) ) (cid:17) are the mean andthe kernel matrix of the multivariate Gaussian distribu-tion p ( Y ) = N ( m ( X ) , K ( X , X )). m ∗ and ( K ∗∗ ) nn (cid:48) = k (cid:16) x ( n ) , x ( n (cid:48) ) (cid:17) are the mean and the kernel matrix of N ( m ( X ∗ ) , K ( X ∗ , X ∗ )). ( K ∗ ) mn = k (cid:0) x ( m ) , x ( n ) (cid:1) is thekernel matrix connecting coordinates between the evalu-ated data points and the images in a pathway. σ n is anoise parameter and I is the identity matrix. A detaileddescription of GP is provided in Supplemental Materi-als [20]. B. Classical and Onsager-Machlup action withtotal energy conservation
For a pathway consisting of N images and A atoms foreach image, the discretized representation of the classical action with energy conservation restraint, Θ cls [5], is,Θ cls ( { X ∗ } , E t ) = S cls + µ E N − (cid:88) n =0 (cid:16) E ( n ) − E t (cid:17) (5)where µ E is the coefficient of the total energy conserva-tion restraint term, E t is the target total energy of thesystem. The classical action is defined as, S cls = N − (cid:88) n =0 dt A − (cid:88) a =0 m a (cid:16) x ( n ) a − x ( n +1) a (cid:17) dt − V (cid:16) x ( n ) (cid:17) (6)Total energy at the time step n is defined as follows: E ( n ) = 12 A − (cid:88) a =0 m a (cid:16) x ( n ) a − x ( n +1) a (cid:17) dt + V (cid:16) x ( n ) (cid:17) (7)To minimize the action analytically, the derivatives ofaction are necessary. The derivatives of classical actionand the restraint terms with respect to the coordinatesof the n th image are, ∂S cls ∂ x ( n ) = dt A − (cid:88) a =0 (cid:18) − ∂∂ x ( n ) a V (cid:16) x ( n ) (cid:17) + m a a ( n ) a (cid:19) (8) ∂∂ x ( n ) (cid:32) µ E N − (cid:88) n =0 (cid:16) E ( n ) − E t (cid:17) (cid:33) = 2 µ E dt A − (cid:88) a =0 (cid:20)(cid:18) ∂∂ x ( m ) a V (cid:16) x ( n ) (cid:17) + 1 dt p ( m ) a (cid:19) (cid:18) − E t + V (cid:16) x ( n ) (cid:17) + 12 m a v ( n ) Ta · v ( n ) a (cid:19) (9)+ 1 dt p ( n − a (cid:18) − E t + V (cid:16) x ( n − (cid:17) + 12 m a v ( n − Ta · v ( n − a (cid:19)(cid:21) where v ( n ) a = ( x ( n +1) a − x ( n ) a ) /dt, a ( n ) a = (2 x ( n ) a − x ( n − a − x ( n +1) a ) /dt and p ( n ) a = m a v ( n ) a . Our parameter set for theclassical action on the MB potential consists of dt = 0 . m a = 1, µ E = 0 . E t = − .
45 and N = 300.If a system propagates according to the Langevin dynamics, the dynamics of the diffusive system can be describedby the OM action [9]. The discretized formula of the OM action and its derivatives are given by, S OM = ∆ V N (cid:88) n =0 (cid:20) dt γ (cid:18)(cid:12)(cid:12)(cid:12) ∇ V (cid:16) x ( n +1) (cid:17)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) ∇ V (cid:16) x ( n ) (cid:17)(cid:12)(cid:12)(cid:12) (cid:19) − (cid:16) ∇ V (cid:16) x ( n +1) (cid:17) − ∇ V (cid:16) x ( n ) (cid:17)(cid:17) · v ( n ) + γdt v ( n ) T · v ( n ) (cid:105) (10) ∂S OM ∂ x ( n ) = γdt a ( n ) + ∂ ∇ V (cid:0) x ( n ) (cid:1) ∂ x ( n ) (cid:18) dt γ ∇ V (cid:16) x ( n ) (cid:17) − dt a ( n ) (cid:19) + 14 ∇ V (cid:16) x ( n − (cid:17) − ∇ V (cid:16) x ( n ) (cid:17) + 14 ∇ V (cid:16) x ( n +1) (cid:17) (11)where ∂ ∇ V (cid:0) x ( n ) (cid:1) /∂ x ( n ) is the Hessian of the potential,∆ V is the energy difference between the initial and the fi-nal step, and γ is a damping constant. We used dt = 0 . N = 300 and γ = 1 (unit time) − for the MB potential.The method we propose in this paper imposes the totalenergy conservation restraint on the OM action, that is,Θ OM ( { X ∗ } , E t ) = S OM + µ E N − (cid:88) n =0 (cid:16) E ( n ) − E t (cid:17) (12)We used dt = 0 . N = 300, γ = 1 (unit time) − , µ E = 0 . E t = − .
45 for the MB potential example.The BFGS algorithm in the scipy Python library is usedto minimize the action and the minimization is conductedon the sine components of the actual coordinates. Thatis, we searched ˜x ( k ) : ˜x ( k ) = 2 N − (cid:88) n =0 x ( n ) sin (cid:18) π (2 k + 1)( n + 1)2 N (cid:19) +( − k x ( N − (13)that minimizes a given action. C. Action minimization using Gaussian processbased on full atom coordinates
Our GPAO method handles two tasks simultaneously.First, it minimizes the uncertainty of a Gaussian PESnear the true MEP in an iterative way. Second, themethod determines the total energy ( E t ) of the MEPautomatically. In this study, we assumed that the E t of a system is close to the maximum potential en-ergy along with the MEP, V (max) . To gradually ap-proximate the E t to V (max) , we approximated E t as (cid:16) µ (max) ∗ − E t − . (cid:17) /
2, where µ (max) ∗ is the estimatedmaximum potential energy of MEP on the Gaussian PES.How the GPAO method minimizes the MEP’s uncer-tainty and approximates the E t to − .
45 on the MB po-tential is illustrated in Fig. 1. In Fig. 1(a), plotted arethe Gaussian PES and µ , which are constructed from theevaluated data points marked as × . In the middle plots of Fig. 1(b), the uncertainty maps, Σ, are plotted withthe trajectory indicating the most uncertain point at eachiteration, Σ (max) . Potential energy and force data of themost uncertain point are added to the database, and theGaussian PES near the trajectory is refined close to thetrue PES. The bottom graph of Fig. 1(c), shows the po-tential energy V , total energy H , kinetic energy E k , andtarget energy E t as a function of distance. The detailedprocedure of GPAO is given as follows:1. Find the most uncertain configuration along thepathway. If the maximum uncertainty is lower thancertain tolerance, find the configuration having amaximum potential energy2. Calculate the energy and force of the configurationobtained at step 1 and add them to the database3. Tune the hyperparameters using the GP regressionalgorithm.4. Adopting a new hyperparameter set and data, con-struct a new Gaussian PES5. Set E t ← (cid:16) µ (max) ∗ − E t − . (cid:17) /
26. With new target energy and a Gaussian PES, op-timize the given action7. Check the convergence of uncertainty, action con-vergence, and the maximum energy difference toreal data. If any of the three conditions do notmeet the requirement, reiterate from step 1. Oth-erwise, it is considered converged.Fig. 1(a) and (b) show how Gaussian PES is graduallyapproximated to the true PES (Fig. 3(a)) as data ac-cumulates throughout the iteration. At iteration 1, thepathway is optimized with the three data points: the ini-tial, the final, and a random intermediate states betweentwo. Due to the lack of data points, potentials estimatedalong the trajectory are still uncertain. Fig. 1(c) illus-trates that the uncertain range (light blue intervals) ofthe estimation at iteration 1 is much wider than the one G a u ss i a n P E S (max) : 0.101 S OM : 3.512 E t : 1.47Iteration : 1 (max) : 0.200 S OM : 3.045 E t : 0.81Iteration : 2 (max) : 0.368 S OM : 2.556 E t : 0.53Iteration : 3 (max) : 0.406 S OM : 1.684 E t : 0.43Iteration : 9 U n ce r t a i n t y (max) : 0.536 l : 1.77 × 10 : 1.00 × 10 : 1.00 × 10 (max) : 0.542 l : 1.98 × 10 : 1.00 × 10 : 1.00 × 10 (max) : 0.444 l : 2.11 × 10 : 1.00 × 10 : 1.00 × 10 (max) : 0.049 l : 9.90 × 10 : 6.77 × 10 : 7.41 × 10 Distance T o t a l & P o t e n ti a l Distance
Distance
Distance
VHE t E k K i n e ti ce n e r gy (a)(b)(c) FIG. 1. Demonstrative images of GPAO process. The expectation value of the Gaussian PES map, uncertainty map, andenergy graph along the pathway are plotted from top to bottom. From left to right, map, and graph at iteration 1, 2, 3, and9. The red line at the Gaussian PES indicates trajectory, the black cross position where data acquired and the white arrowsare force data at that coordinates. On the uncertainty map, we marked additional red crosses corresponding to the latestposition where data being acquired in addition to the white crosses marked at the identical locations on all Gaussian PESs. Inthe energy graph, the blue line shows potential energy, sky blue means 95% confidence range of potential energy, the red dotindicates kinetic energy, the orange dashed line is the total energy and the green diamond line expresses target energy. at iteration 9. Furthermore, Fig. 1(b) shows that maxi-mum uncertainty at iteration 1 ( Σ (max) = 5 . × − ) issmaller than that at iteration 9 ( Σ (max) = 9 . × − ).At iteration 2, including the most uncertain data pointof the previous iteration (marked with a red cross inFig. 1(b)), the pathway is optimized with four datapoints. Most uncertain points in the pathway are markedwith red crosses in Fig. 1(b) at iteration 2, which areadded to the next iteration database. At iteration 9, thepathway is optimized with 11 data points and Σ (max) < . E t to V (max) iteratively. At iteration 1, E t was set to − . µ (max) = − .
101 from thetrajectory optimized with three data points. Using thisinformation, at iteration 2, E t was raised to − .
81, whichis a half of E t − µ (max) of the previous iteration. GPAOresulted in µ (max) of − .
20 with four data points. Simi- larly, at iteration 3, E t was set to half of E t − µ (max) ofthe previous iteration, − . S OM , E t , µ (max) , Σ (max) = 1 . × (cid:112) max diag ( Σ ∗ ) and the hy-perparameters of the kernels l , σ En , σ fn of each iterationare displayed on Fig. 1(a) and Fig. 1(b). σ En and σ fn arenoise parameters corresponding to the energy and forcesof data points, respectively. The parameter set for the ac-tion was assigned to the same as the conventional ADMDwith Θ OM calculation.We used the squared exponential kernel as the Gaus-sian kernel for the MB model: k (cid:16) x ( i ) , x ( j ) (cid:17) = σ f exp (cid:18) − l (cid:16) x ( i ) − x ( j ) (cid:17) T (cid:16) x ( i ) − x ( j ) (cid:17)(cid:19) (14)where l and σ f are isotropic real hyperparameters. Weused the multivariate Gaussian distribution with zeromean for our GP regression model, m = m ∗ = 0. D. Action optimization using the Gaussian processand reduced coordinates
150 100 50 0 50 100 15015010050050100150
Alanine Dipeptide PES (eV)(a)(b) C eq C ax φ (deg) ψ φ ψ C α C β NC − CN +1 FIG. 2. (a) The PES of alanine dipeptide as a functionof reduced coordinates, φ and ψ angles, V ( ψ, φ ). (b) The C eq and C ax conformations of alanine dipeptide are drawn.The brown balls correspond to carbon, the light blue onesare nitrogen, the red ones are oxygen, and the ivory ones arehydrogen. Alanine dipeptide is a small molecule consisting of 22atoms, and it has been widely used as a test system forcomputational methods in biophysics [13]. As displayed in Fig. 2, alanine dipeptide has two rotatable bonds asso-ciated with two dihedral angles, the φ and ψ angles. The φ angle is a dihedral angle between C − N C α C , and the ψ angle is defined as a dihedral angle between N C α CN +1 .These two dihedral angles almost exclusively determinealanine dipeptide’s conformation because the bond an-gles and bond lengths are much more rigid than dihedralangles. Alanine dipeptide has two stable conformations,C7 eq and C7 ax .The PES’s shape shows that multiple conformationaltransition pathways from one conformation to the otherare possible (Fig. 2). Lee and coworkers discovered allpossible conformational transition pathways through theaction optimization using the conformation space anneal-ing (Action-CSA) method [7, 13].The potential energies and forces of alanine dipep-tide are calculated using the ab initio based density-functional-theory package VASP [21, 22]. The con-strained MD method of
VASP is used to optimize theconformations of C7 ax and C7 eq with fixed φ and ψ an-gles. With fixed both φ and ψ dihedral angles, MD sim-ulations are performed with the Anderson thermostat.With random initial velocities, the heat bath tempera-ture decreased from 80 K to 0 K for 500 steps. Theprobability that a system interacts with the heat bathwas given to 5%. We use the generalized gradient ap-proximation (GGA) exchanged functional and a van derWaals scheme considering three-body interactions. Theenergy cutoff is set to 520 eV, and a convergence ratefor the SCF scheme is 10 − eV. The system is slowlyannealed to a stable structure with the fixed φ and ψ angles.Since the structures of alanine dipeptide are al-most exclusively defined by φ and ψ angles, its po-tential can be expressed as the function of themas follows: V (cid:0)(cid:8) x ( n ) (cid:9)(cid:1) = V (cid:0) φ ( n ) , ψ ( n ) (cid:1) . This re-duction on coordinates x ( n ) = (cid:2) φ ( n ) , ψ ( n ) (cid:3) , v ( n ) = (cid:2) φ ( n +1) − φ ( n ) , ψ ( n +1) − ψ ( n ) (cid:3) /dt is extended on the for-mulation of the modified OM action for the alaninedipeptide system. Our Θ OM for the alanine dipeptidesystem is,Θ OM = S OM + µ E N − (cid:88) n =0 (cid:16) E ( n ) − E t (cid:17) = ∆ V N (cid:88) n =0 (cid:34) dt γ (cid:32) I ( n ) φ ( φ, ψ ) + 1 I ( n ) ψ ( φ, ψ ) (cid:33) (cid:18)(cid:12)(cid:12)(cid:12) ∇ V (cid:16)(cid:110) x ( n +1) (cid:111)(cid:17)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) ∇ V (cid:16)(cid:110) x ( n ) (cid:111)(cid:17)(cid:12)(cid:12)(cid:12) (cid:19) − (cid:16) ∇ V (cid:16)(cid:110) x ( n +1) (cid:111)(cid:17) − ∇ V (cid:16)(cid:110) x ( n ) (cid:111)(cid:17)(cid:17) · v ( n ) + γdt (cid:32) I ( n ) φ ( φ, ψ ) (cid:18) φ ( n +1) − φ ( n ) dt (cid:19) + I ( n ) ψ ( φ, ψ ) (cid:18) ψ ( n +1) − ψ ( n ) dt (cid:19) (cid:33)(cid:35) + µ E N − (cid:88) n =0 (cid:32) I ( n ) φ ( φ, ψ ) (cid:18) φ ( n +1) − φ ( n ) dt (cid:19) + 12 I ( n ) ψ ( φ, ψ ) (cid:18) ψ ( n +1) − ψ ( n ) dt (cid:19) + V (cid:16)(cid:110) x ( n ) (cid:111)(cid:17) − E t (cid:33) (15)where ∇ V (cid:0)(cid:8) x ( n ) (cid:9)(cid:1) = (cid:2) ∂ φ V (cid:0)(cid:8) x ( n ) (cid:9)(cid:1) , ∂ ψ V (cid:0)(cid:8) x ( n ) (cid:9)(cid:1)(cid:3) indicates the potential difference with respect to φ and ψ , I ( n ) φ ( φ, ψ ) and I ( n ) ψ ( φ, ψ ) indicate the moments of inertiaassociated with φ and ψ angles of n th image, respectively.In case we obtain ∇ V (cid:0)(cid:8) x ( n ) (cid:9)(cid:1) with full atomic coor-dinates, we map them onto the reduced coordinates byprojecting each force vector to a corresponding rotationaxis. The average values of I φ and I ψ are 112.96 u˚A and 75.93 u˚A , respectively. The derivatives of the mo-ments of inertia, ∂ φ I ψ , ∂ ψ I ψ , ∂ φ I φ , and ∂ ψ I φ are calcu-lated with VASP using the data we used for generatingthe PES map. Except for high energy regions, where φ ≈ ◦ and ψ ≈ ◦ , the changes in the moment ofinertia are negligible. About 99% of the data show that | ∂ φ I ψ | is less than 2.211, which is less than 2%, and 95%of the data show | ∂ φ I ψ | < . | ∂ ψ I ψ | , | ∂ φ I φ | and | ∂ ψ I φ | are less than 2.702u˚A , 0.572 u˚A , and 1.258 u˚A , respectively. About 95%of the data are less than 0.948 u˚A , 0.388 u˚A , and 0.540u˚A . We assume that the changes in the moments of in-ertia, ∂ φ I ψ , ∂ ψ I ψ , ∂ φ I φ , and ∂ ψ I φ are negligible becausea trajectory rarely visits high energy regions where themoment of inertial may change rapidly. We use the GPto estimate the moment of inertia with the same databaseto construct the Gaussian PES. The action for the relax-ation was Θ OM with parameters γ = 0.001 (˚A (cid:112) u / eV) − , dt = 0.1, 0.2, 0.3 (˚A (cid:112) u / eV) and µ E = 0 . E. Efficient global action optimization usingAction-CSA
We combined GPAO with the Action-CSA method [13]to search multiple transition pathways because stochas-tic reactions occur through multiple pathways in general.In this work, the Action-CSA is set to start with ran-domly generated 40 initial pathways, followed by localoptimizations using GPAO. Afterward, we proceed with the following steps:1. Pick two pathways from the candidates randomly.2. Generate a new pathway from two pathways usinga crossover operator.3. Apply random mutations/perturbations to the newpathway with a certain probability4. Locally minimize the new pathway.5. Pick the most similar pathway to the new path-way from the candidates. If the distance betweenthe new pathway and a similar pathway is lowerthan a cutoff distance, compare the OM action be-tween them. Otherwise, compare the highest OMpathway in the candidates.6. If a new one has lower action than its counterpart,replace it with the new one.7. Reduce the cutoff distance. If the cutoff distance islower than a certain threshold, stop the iteration.Otherwise, go back to step 1.With two randomly selected pathways, we choose arandom image along a pathway. With the selected posi-tion, we slice both trajectories at the selected image andmerge them. After a new pathway is generated, a ran-dom mutation/perturbation is performed with a proba-bility of 30%. For the pathway’s mutation, the pathwayis modified by adding another random number at thesine components of a pathway. That is, the ˜x ( k ) in theEq. (13) is altered. After the local minimization of thepathway using the GPAO method, we search the nearestneighbor by measuring the Fr´echet distance between thenew pathway and the current candidates. The nearestneighbor is considered similar , readily accessible to eachother if the Fr´echet distance is shorter than the currentcutoff distance d cut .The initial cutoff distance is set to half of the averagedistance between the initial pathways. After an iteration,we reduce d cut by 2% and repeat the whole procedurefrom step 1 until d cut becomes below 0.05. III. RESULT AND DISCUSSIONA. Action minimization on the M¨uller-Brownpotential (a) Muller Brown clsOM S OM time step (unitless) T o t a l E n e r gy ( un itl e ss ) Gap : 3.85 Gap : 0.05 ( OM )Gap : 0.09 ( cls )(b) FIG. 3. A comparison of the minimum action pathways ob-tained with classical action with energy conservation restraint( S cls , blue dashed line), OM action ( S OM , orange dash-dottedline), and OM action with energy conservation constraint(Θ OM , green solid line). (a) the minimum action pathwaysbetween two stable potential energy minima on the M¨uller-Brown potential and a zoomed image near the saddle point.(b) energy profile. To compare the characteristics of pathways obtainedwith different actions, we carry out action optimizationcalculations on the MB potential, defined as follows: V ( x, y ) = (cid:88) µ =1 A µ e a µ ( x − x µ ) + b µ ( x − x µ )( y − y µ ) + c µ ( y − y µ ) , (16)where x µ = (1 , , − . , − y µ = (0 , . , . , a µ = ( − , − , − . , . b µ = (0 , , , . c µ =( − , − , − . , . A µ = ( − , − , − . , . C OMΘ , the pathway with theminimum modified OM action Θ OM , successfully con-serves the system’s total energy throughout the pathway (Fig. 3). The magnitude of total energy fluctuations ofthe pathway was reduced by 98.7% compared to thatof C OM S , the pathway optimized using S OM . The OMvalue of C OMΘ increased only about 19%, and the Fr´echetdistance to C OM S is less than 0.067 (Table I). This re-sult demonstrates that the total energy conservation termdoes not change the shape of the pathway significantlywhile keeping the total energy of the system constant.Fig. 3(b) shows the total energy profile of three trajec-tories and their energy gaps, the differences between themaximum and the minimum total energies of pathways.The energy gaps of C clsΘ , C OM S and C OMΘ are 0.09, 3.85 and0.05, respectively. The large spikes of the total energyoccurred near the region where the trajectory climbs upthe energy barriers of the PES between the timesteps 0.8 ∼ S OM equation (Eq. (10)). The firstterm in Eq. (10) prefers to minimize atom’s forces, andthe third term keeps the velocities between images slow.The second term in Eq. (10) favors following a concavepotential energy surface along its ongoing direction. Tobe more specific, if a pathway has to climb up a PES dur-ing a transition, it prefers to climb up the surface fast.This tendency is clearly manifested in Fig. 3(a). Rightbefore C OM S reaches the saddle point with the maximumpotential energy of the PES, the images are largely sep-arated because of large velocities of images.In contrast, the images of C OMΘ near the saddle pointare almost continuous and uniformly distributed. Thisuniform distance between images is because the total en-ergy spike is mainly due to the second term in Eq. (10).Our results demonstrate that unphysical increases ofatomic velocities to lower the second term can be sup-pressed effectively by the energy conservation restraintwithout changing the trajectory direction.Our method avoids the occurrence of unphysical en-ergy spikes observed with C OM S , but still maintains S OM low. In Table I, S OM values of C clsΘ , C OM S and C OMΘ are14.49, 1.387, and 1.652, respectively. The potential ener-gies at the saddle point, V (max) , of C clsΘ , C OM S , and C OMΘ are − . − . − . C OM S gives only little effect on findinglow V (max) and S OM . Furthermore, Fig. 3(a) shows howtwo pathways of C OM S and C OMΘ overlap each other. TheFr´echet distance between two trajectories is only 0.067.Also, 95% of the intermediate images of C OMΘ are locatedwithin a distance of 0.013 from C OM S . B. Efficiency of GPAO
The most significant advantage of the GPAO methodsis a dramatic reduction of the number of force calls by four-to-five orders of magnitude than the conventionalaction optimization method while preserving their accu-racy (Table I). The numbers of force calls to find C clsΘ , Number of force calls V (max) S OM Distance to C OM S Distance to non-GP C clsΘ × − .
400 14.49 0.250 - C OM S × − .
406 1.387 - - C OMΘ × − .
406 1.652 0.067 -GP - C clsΘ − .
401 11.08 0.250 0.006GP - C OM S − .
405 1.187 0.020 0.082GP - C OMΘ − .
406 1.684 0.083 0.042 C NEB × − .
407 - - -GP- C NEB − .
412 - - 0.0766TABLE I. A comparison of performance between the GP-assisted algorithms and the conventional action optimization algo-rithms was conducted on the M¨uller-Brown potential. The number of force call is the number of actual function evaluationsduring the optimization process, V (max) is the maximum potential energy of a pathway, S OM is the OM action value of a path-way. Distance to C OM S is a Fr´echet distance to C OM S and distance to non-GP is a Fr´echet distance from the pathway optimizedvia non-GP-assisted method S OM (eV) µ (max) (eV) T total A B C D E F G C OM S , and C OMΘ using the conventional NEB methodswere 1,220,400, 586,200, and 993,900, while those of ourGPAO methods required only 12, 14, and 12, respec-tively. In other words, our GPAO methods decreasedcomputational costs to find transition pathways by fiveorders of magnitude.The significant reduction of the computational cost bythe GPAO method allows us to find transition pathwaysof a molecule with a quantum-mechanical (QM) poten-tial. As discussed in the next section, a single force calcu-lation of alanine dipeptide with
VASP requires around1741.4 CPU seconds. Using our method, the optimaltransition pathway of alanine dipeptide within a singleday. However, it will take five years to complete thesame calculation by the conventional action optimizationmethod and the QM potential.The significant computational gain is attainable be-cause the computation cost of GPAO is independent ofthe number of images and atoms. The number of po-tential energy evaluations of GPAO depends only on theconvergence of a Gaussian PES. In contrast to GPAO, thecalculation cost of the conventional action optimization method, such as ADMD and DOA, is strictly propor-tional to the number of images in a pathway. To ensurethe continuity of a pathway, we followed the same numberof step sizes, 300 images for a pathway, with the previousstudy [5] on the MB potential. Thus, for a single actionevaluation, the energy and forces of 300 images must becalculated.On the contrary, the number of force calls required forthe GPAO method depends on the number of data points,which should be large enough to accurately approximatethe PES. Overall, the GPAO required similar numbersof force calls regardless of the kinds of action (Table I).In addition, the results of the GP-assisted NEB methods,GP- C NEB , showed similar numbers of force calls to repro-duce the result of the conventional NEB method. Indeed,it turns out that the number of data points required tofind an accurate MEP on the MB potential is also in arange between 10 and 20. These results imply that thetype of local relaxation method does not significantly af-fect the efficiency of the GPAO methods.Our GPAO methods lead to the transition pathwaysremarkably close to those obtained by the conventionalaction optimization method. Overall, the discrepancybetween the pathways obtained with the two classes ofthe methods is less than 0.08% of the pathways’ totaldistances (Table I). The distance between C clsΘ and GP- C clsΘ was 0.006. Similarly, the distances between C OM S and C OMΘ and the corresponding pathways obtained bythe GPAO method were 0.082 and 0.042, respectively. Ingeneral, the results of all five GP-assisted methods show adeviation of less than a Fr´echet distance of 0.1 comparedto the conventional action optimization result. Even fur-ther, the GP- C clsΘ results deviate from the conventionalresult less than a Fr´echet distance of 0.01. These resultssupport that approximation of PES by GP is accurateenough to locate the transition pathways of low actionon a given PES.GPAO gives us a way to improve the Hessian’s accu-racy by accumulating data points, which allows us toutilize calculation tools that are hard to acquire Hes-0 ( d e g ) AB CD T total = 3 psAlanine Dipeptide PES (eV) ( d e g ) AB CD EFG T total = 6 ps
100 0 100 (deg) ( d e g ) A B CD EFG T total = 9 ps (a)(b)(c) FIG. 4. The distinctive pathways of each time scale gained viaGP-CSA. Pathways gained with T total = 30 ˚A (cid:112) u/eV ≈ T total = 60 ˚A (cid:112) u/eV ≈ T total =90 ˚A (cid:112) u/eV ≈ sian. The conventional way of obtaining the Hessian of aDFT potential is using a finite difference method. Thatis, to acquire a Hessian at an arbitrary point φ ( n ) , ψ ( n ) ,we need to compute ∇ V ( φ ( n ) , ψ ( n ) ) , ∇ V ( φ ( n ) + ε, ψ ( n ) )and ∇ V ( φ ( n ) , ψ ( n ) + ε ). These two additional force callstriples the computation cost to evaluate Θ OM , which pre-vents the use of finite difference methods from being ap-plied to relatively large systems. Thus, such a fast evalu-ation of a potential enables us to use Hessian informationfor practical problems. C. Multiple transition pathways of alaninedipeptide
The PES of alanine dipeptide is periodic because twodihedral angles, φ , and ψ , are chosen as the reaction coor-dinates. We use the following periodic kernel to describethe PES with GP: k (cid:16) x ( m ) , x ( n ) (cid:17) = σ f exp − l D (cid:88) d =1 sin (cid:16) x ( m ) d − x ( n ) d (cid:17) (17)where x = ( φ, ψ ) is a two-dimensional reaction coordi-nate, l, σ f are isotropic hyperparameters.We sampled multiple transition pathways of alaninedipeptide using GP-CSA calculations with three differ-ent total transition times, T total = N dt = 3, 6, and 9ps. Among the pathways obtained with GP-CSA, thedistinctive pathways having low action values are plot-ted in Fig. 4. Path A in Fig. 4 has the lowest S OM valuefor all transition times, suggesting that Path A is themost dominant pathway between C7 eq and C7 ax confor-mations. This result is consistent with previous studiesconducted with molecular mechanics potentials [13, 23].The maximum potential energy of a pathway estimatedwith GP, µ (max) is near -128.48 eV (Table. II). For allthree different T total , our method consistently located thecorrect saddle points.The S OM value of Path D with a T total of 3 ps is 11.282eV, which is high enough to make the path improbable.However, the action value decreases drastically to 5.551eV and 4.887 eV when T total becomes 6 ps and 9 ps. As T total increases, Path D becomes longer, which allows thepathway to detour the high potential region and to seeka lower valley of the PES (see Path D in Fig. 4(b), (c)).This tendency is also observed in the µ (max) values inTable II. µ (max) of Path D at 3 ps is 0.2 eV higher thanthe other transition times.It is noticeable that the relatively longer pathways, i.e,Path E , F , and G , were not observed in simulations with T total = 3 ps. As T total becomes longer, GPAO can ac-cess additional pathways, which require longer T total tooccur. S OM values of Path E obtained with T total = 6 psand 9 ps were 8.547 eV and 6.745 eV, and µ (max) valueswere − .
237 eV and − .
316 eV, respectively. Themaximum potential energy at the saddle point of Path1 E of T total = 6 ps was 0.1 eV higher than Path E with T total = 9 ps. This difference arises because Path E with T total = 6 ps did not pass the saddle point due to a shorttransition time. However, Path E with T total = 9 ps cor-rectly passed the saddle point (see E in the Fig. 4(c)).With T total = 6 ps, we find two additional pathways,Path F and Path G , similar to the most dominant path-way, Path A . These two pathways are longer pathwayswith higher S OM values, 4.56 eV and 5.293 eV, that ofPath A . These results clearly indicate that our GPAOmethod can accurately find multiple transition pathwaysof the conformational change of alanine dipeptide with aQM potential, which have not been reported previously. IV. CONCLUSION
In this study, we demonstrate that the use of the GP al-gorithm dramatically enhances the efficiency of searchingminimum action pathways between potential energy min-ima by the accurate approximation of PES with much-reduced numbers of energy evaluations. Our results im-ply that pathway search using the modified OM actionwith a total energy conservation restraint is a prominentapproach. Compared to the sole optimization of OM ac-tion, large fluctuations of kinetic energies are effectivelyremoved, without affecting the OM action of a pathway.OM action gradient requires the Hessian of a potentialfunction, which is computationally expensive to obtainin general. Combining GP with the action optimizationmethod gives not only superior computational efficiencybut also the gradients of OM action without actual Hes-sian calculations. This advantage is particularly impor-tant for potential energies whose Hessians are challeng-ing to obtain. For example, the Hessian of DFT-basedpotentials is often incorrect. However, with GPAO, theaccumulation of gradient information at multiple datapoints leads to accurate estimation of the Hessian of aPES from a Gaussian PES.Since the Hessian of a Gaussian PES is analyticallygiven, we can achieve a way to obtain the OM action ofDFT-based potentials and enhance computational effi-ciency by five orders of magnitude. Additionally, we alsoshowed that our GPAO methods find suitable externalparameters, such as target energy, automatically.This study also presents a new multiple-pathway-sampling method by combining the Action-CSAmethod [13] with the GPAO methods. Due to the GPAOmethod’s efficiency, multiple conformational transition pathways between the two stable conformations of ala-nine dipeptide can be successfully obtained using a DFT-based potential. This is the first work that sampled mul-tiple low action pathways of alanine dipeptide with aDFT potential to the best of our knowledge.Undertaking higher-dimensional problems in an effi-cient way is still a big challenge in machine learning fieldsand others. This paper proves that our GPAO methodproperly works for small-dimensional problems. How-ever, the method does not guarantee its efficiency onproblems with a large degrees of freedom yet. Fortu-nately, many methods have been suggested to solve prob-lems on a high dimensional space using various types ofdescriptors [24–27]. Finding reaction pathways of thesystems with many degrees of freedom can be applied tovarious problems such as finding drug binding pathways,battery, and protein folding. Thus, we believe that theGPAO methods will open up new possibilities for suchproblems.
ACKNOWLEDGMENTS
We like to thank to In-Ho Lee for the help in the ini-tial stage of this project. JL acknowledges the supportof the National Research Foundation of Korea (NRF)grant funded by the Korean government (MSIT) (No.2018R1C1B600543513). JY acknowledges the support bythe National Research Foundation of Korea (NRF)(No.2020R1F1A1066548). Additional financial support fromSamsung Electronics is also acknowledged. This workwas also supported by the National Supercomputing Cen-ter with supercomputing resources and technical sup-ports (KSC-2020-CRE-0026).
V. DATA AVAILABILITY
The source code of GPAO can be found in theGitHub repository (https://github.com/schinavro/taps).GPAO consists of multiple tools in the open-source pack-age
TAPS . TAPS incorporates tools needed for data-driven pathway search methods such as a database oratomic calculator. For atomic calculations, we use theatomic simulation environment (
ASE ) [28] package,which can easily link with atomic calculators such as
VASP [21, 22],
Quantumespresso [29],
OpenMX [30],and SIESTA [31]. [1] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Optimiza-tion by simulated annealing, Science , 671 (1983).[2] B. H. Good, Y.-A. de Montjoye, and A. Clauset, Perfor-mance of modularity maximization in practical contexts,Physical Review E , 046106 (2010). [3] R. Guimer´a and L. A. Nunes Amaral, Functional cartog-raphy of complex metabolic networks, Nature , 895(2005).[4] H. J´onsson, G. Mills, and K. W. Jacobsen, Nudged elas-tic band method for finding minimum energy paths oftransitions, in Classical and Quantum Dynamics in Con- densed Phase Simulations (World Scientific, Singapore,1998) pp. 385–404.[5] D. Passerone and M. Parrinello, Action-Derived Molecu-lar Dynamics in the Study of Rare Events, Physical Re-view Letters , 108302 (2001).[6] L. R. Pratt, A statistical method for identifying transi-tion states in high dimensional problems, The Journal ofChemical Physics , 5045 (1986), publisher: AmericanInstitute of Physics.[7] J. Lee and B. R. Brooks, Direct global optimizationof Onsager-Machlup action using Action-CSA, ChemicalPhysics , 110768 (2020).[8] I.-H. Lee, J. Lee, and S. Lee, Kinetic energy control inaction-derived molecular dynamics simulations, PhysicalReview B , 064303 (2003).[9] S. MacHlup and L. Onsager, Fluctuations and irre-versible process. II. Systems with kinetic energy, PhysicalReview , 1512 (1953).[10] A. Bach and D. D¨urr, Entropy density in function spaceand the Onsager-Machlup function, Physics Letters A ,244 (1978).[11] P. Faccioli, M. Sega, F. Pederiva, and H. Orland, Dom-inant pathways in protein folding, Physical Review Let-ters , 108101 (2006).[12] P. Faccioli, Characterization of protein folding by dom-inant reaction pathways, The journal of physical chem-istry. B , 13756 (2008).[13] J. Lee, I.-H. Lee, I. Joung, J. Lee, and B. R. Brooks,Finding multiple reaction pathways via global optimiza-tion of action, Nature Communications , 1 (2017).[14] L. Onsager and S. MacHlup, Fluctuations and irre-versible processes, Physical Review , 1505 (1953).[15] H. S. Wio, Path integrals for stochastic processes: Anintroduction (World Scientific, Singapore, 2013).[16] O.-P. Koistinen, F. B. Dagbjartsd´ottir, V. ´Asgeirsson,A. Vehtari, and H. J´onsson, Nudged elastic band calcu-lations accelerated with Gaussian process regression, TheJournal of Chemical Physics , 152720 (2017).[17] J. A. Garrido Torres, P. C. Jennings, M. H. Hansen,J. R. Boes, and T. Bligaard, Low-Scaling Algorithm forNudged Elastic Band Calculations Using a SurrogateMachine Learning Model, Physical Review Letters ,156001 (2019).[18] J. M. Bernardo, J. O. Berger, A. P. Dawid, and A. F. M.Smith, Bayesian model averaging and model searchstrategies, in
Bayesian Statistics , Vol. 6 (Oxford Univer-sity Press, 1999) pp. 157–185.[19] A. O’Hagan, Curve Fitting and Optimal Design for Pre-diction, Journal of the Royal Statistical Society. Series B(Methodological) , 1 (1978).[20] See Supplemental Material at [url will be inserted bypublisher] for detailed explanation on Gaussian processregression. [21] G. Kresse and J. Furthm¨uller, Efficiency of ab-initio totalenergy calculations for metals and semiconductors usinga plane-wave basis set, Computational Materials Science , 15 (1996).[22] G. Kresse and J. Furthm¨uller, Efficient iterative schemesfor ab initio total-energy calculations using a plane-wavebasis set, Physical Review B , 11169 (1996).[23] R. Elber and D. Shalloway, Temperature dependent reac-tion coordinates, The Journal of Chemical Physics ,5539 (2000).[24] A. P. Bart´ok, R. Kondor, and G. Cs´anyi, On representingchemical environments, Physical Review B , 184115(2013).[25] J. Behler, Atom-centered symmetry functions for con-structing high-dimensional neural network potentials,The Journal of Chemical Physics , 074106 (2011).[26] J. Behler and M. Parrinello, Generalized Neural-NetworkRepresentation of High-Dimensional Potential-EnergySurfaces, Physical Review Letters , 146401 (2007).[27] E. Kocer, J. K. Mason, and H. Erturk, Continuous andoptimally complete description of chemical environmentsusing Spherical Bessel descriptors, AIP Advances ,015021 (2020).[28] A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E.Castelli, R. Christensen, Marcin Du(cid:32)lak, J. Friis, M. N.Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C.Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L.Kolsbjerg, J. Kubal, Kristen Kaasbjerg, S. Lysgaard,J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka,Andrew Peterson, C. Rostgaard, J. Schiøtz, O. Sch¨utt,M. Strange, K. S. Thygesen, Tejs Vegge, L. Vilhelmsen,M. Walter, Z. Zeng, and K. W. Jacobsen, The atomicsimulation environment—a Python library for workingwith atoms, Journal of Physics: Condensed Matter ,273002 (2017).[29] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ-cioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris,G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis,A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello,L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentz-covitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of mate-rials, Journal of Physics: Condensed Matter , 395502(2009).[30] T. Ozaki and H. Kino, Efficient projector expansion forthe ab initio LCAO method, Physical Review B ,10.1103/PhysRevB.72.045121 (2005).[31] J. M. Soler, E. Artacho, J. D. Gale, A. Garc´ıa, J. Jun-quera, P. Ordej´on, and D. S´anchez-Portal, The SIESTAmethod forab initioorder-Nmaterials simulation, Journalof Physics: Condensed Matter14