Achieving an Order of Magnitude Speed-up in Hybrid Functional and Plane Wave based Ab Initio Molecular Dynamics: Applications to Proton Transfer Reactions in Enzymes and in Solution
aa r X i v : . [ phy s i c s . c h e m - ph ] J a n Achieving an Order of Magnitude Speed-up in Hybrid Functional and PlaneWave based
Ab Initio
Molecular Dynamics: Applications to Proton TransferReactions in Enzymes and in Solution
Sagarmoy Mandal,
1, 2
Vaishali Thakkur, and Nisanth N. Nair a) Department of Chemistry, Indian Institute of Technology Kanpur (IITK), Kanpur - 208016,India Interdisciplinary Center for Molecular Materials and Computer Chemistry Center,Friedrich-Alexander-Universit¨at Erlangen-N¨urnberg (FAU), N¨agelsbachstr. 25, 91052 Erlangen,Germany (Dated: 11 January 2021)
Ab initio molecular dynamics (AIMD) with hybrid density functionals and plane wave basisis computationally expensive due to the high computational cost of exact exchange energyevaluation. Recently, we proposed a strategy to combine adaptively compressed exchange(ACE) operator formulation and multiple time step (MTS) integration scheme to reduce thecomputational cost significantly [
J. Chem. Phys. , 151102 (2019)]. However, it wasfound that the construction of the ACE operator, which has to be done at least once in everyMD time step, is computationally expensive. In the present work, systematic improvementsare introduced to further speed-up by employing localized orbitals for the construction of theACE operator. By this, we could achieve a computational speed-up of an order of magnitudefor a periodic system containing 32-water molecules. Benchmark calculations were carried outto show the accuracy and efficiency of the method in predicting the structural and dynamicalproperties of bulk water. To demonstrate the applicability, computationally intensive freeenergy computations at the level of hybrid density functional theory were performed toinvestigate (a) methyl formate hydrolysis reaction in neutral aqueous medium and (b) protontransfer reaction within the active site residues of class-C β -lactamase enzyme. I. INTRODUCTION
Ab initio molecular dynamics (AIMD) techniquesare extensively used to investigate structural and dy-namical properties of a wide variety of molecular sys-tems, and to predict mechanism and free energetics ofphysiochemical processes.
In AIMD simulations, inter-atomic interactions are evaluated on-the-fly using elec-tronic structure calculations at every MD step.
KohnSham density functional theory (KS-DFT) with planewave (PW) basis set is widely used for performing AIMDsimulations of periodic condensed matter systems. Thechoice of the exchange-correlation (XC) functionals forthe KS-DFT calculations determines the accuracy ofthe predicted properties. Generalized gradient approx-imation (GGA) type XC functionals are most com-monly used to carry out AIMD using PW KS-DFT. How-ever, these functionals suffer from self interaction error(SIE) where the XC functional erroneously includesthe unphysical self-interaction of electron density withitself. Due to the SIE, the XC functionals tend to over-delocalize the electron density, leading to the underesti-mation of bandgap of solids, reaction barriers, and disso-ciation energies.
SIE can be minimized by using hybrid function-als, where a certain portion of the Hartree-Fock(HF) exchange energy is added to the GGA exchangeenergy.
Hybrid functionals are generally knownto improve the prediction of energetics, structures, elec-tronic properties, chemical reaction barriers and band a) Electronic mail: [email protected] gap of solids.
Also, hybrid functional based AIMDsimulations have been shown to improve the descriptionof the structural and dynamical properties of liquids as well as the accuracy of the computed free energysurfaces.
However, AIMD simulation with hybridfunctionals and PW basis set is extremely time consum-ing due to the high computational cost associated withthe exact exchange energy evaluation. Thus, it is nota common practice to perform hybrid functional basedAIMD simulations for systems containing several hun-dred or more atoms.Various possible strategies have been proposed tospeed-up such calculations, for example, utilization oflocalized orbitals, usage of the multipletime step (MTS) algorithms, use of coordinate-scaling approach, employment of massively parallelalgorithms and other approaches.
Recently, weproposed a robust method for performing efficienthybrid functionals and PW based AIMD, where a MTSintegrator scheme was employed based on the adaptivelycompressed exchange (ACE) operator formalism. Inthe proposed method, the ionic forces were artificiallypartitioned into computationally cheap fast forces andcomputationally costly slow forces with the help of an ap-proximate representation of the exact exchange operator(i.e. the ACE operator). Using the MTS algorithm, thecomputationally cheap fast forces were computed morefrequently as compared to the computationally costlyslow forces, thereby reducing the computational cost ofperforming such AIMD calculations. Additionally, ourmethod has been shown to be efficient and accurate inpredicting the structural and dynamical properties of re-alistic condensed matter systems.It was realized that the construction of the ACE op-erator, which has to be done at least once in everyMD time step, is the computational bottleneck. In thepresent work, we systematically improve the computa-tional speed by employing localized orbitals for the con-struction of the ACE operator. Specifically, the local-ized selected column of the density matrix (SCDM) or-bitals are used to build the ACE operator. The SCDMmethod for obtaining localized orbitals is computation-ally efficient, and the procedure is not iterative (unlikein the case of computing Maximally Localized WannierFunctions). Recently, Carnimeo et al. also employeda methodology to speed-up ACE operator constructionbased on the localized SCDM orbitals. Their implemen-tation in the Quantum ESPRESSO code reported a speed-up of at least 3-4 fold for hybrid functionals and PWbased exact exchange energy calculations.Our approach presented in this paper, termed as s-MTACE, opens up the possibility of performing hy-brid DFT based MD simulations of large condensedmatter systems providing a speed up of 10 fold ormore. After presenting the theory and computationaldetails, benchmark calculations are reported to demon-strate the accuracy and the computational efficiencyof the method. In particular, we have investigatedthe structural and dynamical properties of bulk wa-ter for benchmarking purpose. Subsequently, we usedour approach to perform free energy calculations oftwo important proton transfer reactions, namely methylformate hydrolysis reaction in neutral water employ-ing periodic AIMD simulations (Figure 1), and pro-ton transfer reaction within the active site residues ofclass-C β -lactamase enzyme complexed with cephalothindrug using quantum-mechanical/molecular-mechanical(QM/MM) simulations (Figure 2). Here we compare thefree energies of chemical reactions computed using bothPBE (GGA) and PBE0 (hybrid) density functionals. Ourresults shed light on the performance of the PBE andthe PBE0 functionals in predicting free energy barriersof reactions involving proton transfers in solutions and inenzymes from finite temperature MD simulations. II. THEORY
KS-DFT calculations with hybrid functionals requiresthe application of the exact exchange operator V X = − P N orb j | ψ j ih ψ j | r on each of the KS orbitals | ψ i i at eachstep of the self consistent field (SCF) iterations: V X | ψ i i = − N orb X j | ψ j i D ψ j (cid:12)(cid:12)(cid:12) ( r ) − (cid:12)(cid:12)(cid:12) ψ i E , = − N orb X j v ij ( r ) | ψ j i , i = 1 , ...., N orb (1)with v ij ( r ) = D ψ j (cid:12)(cid:12)(cid:12) ( r ) − (cid:12)(cid:12)(cid:12) ψ i E . (2)Here, N orb is the total number of occupied orbitals.Computation of v ij ( r ) is efficiently done in reciprocal space using Fourier transform. If N G is the totalnumber of PWs, the computational cost for doing Fouriertransform scales as N G log N G by using the Fast Fouriertransform (FFT) algorithm. The total computationaloverhead scales as N N G log N G , as the operation of V X on all the KS orbitals requires N evaluations of v ij ( r ). As a result, such calculations are highly compu-tationally intensive for typical molecular systems of ourinterest with about 100 atoms.Whereas, the recently developed ACE operatorformulation could greatly reduce the computationalcost of such calculations. In the ACE operator formu-lation, the full rank V X operator is approximated bythe low rank ACE operator V ACEX = − P N orb k | P k ih P k | .Here, {| P k i} is the set of ACE projection vectors whichcan be computed by the decomposition of the V X oper-ator; See Appendix A for more details. Construction of {| P k i} requires evaluation of { V X | ψ i i} , which is a com-putationally costly step, because of the N number ofevaluations of v ij ( r ). However, once the V ACEX operatoris constructed, the evaluation of the action of this opera-tor on KS orbitals can be done easily with N numberof inner products as V ACEX | ψ i i = − N orb X k | P k i h P k | ψ i i , i = 1 , ...., N orb . (3)The advantage of the ACE approach is that the cost ofapplying the V ACEX operator on each KS orbitals is muchless as compared to V X operator.In our earlier method, we took advantage of thisparticular feature of the ACE operator to combine withthe MTS scheme for performing efficient hybrid func-tional based AIMD. However, the construction of the V ACEX operator, which has to be done at least once in ev-ery MD time step, is computationally demanding. FromEquation (2), it is clear that v ij ( r ) will be zero for non-overlapping orbital pairs and such pairs of orbitals willnot contribute to the sum in Equation (1). Thus, it ispossible to speed-up the construction of the ACE oper-ator by employing a screened set of localized orbitals. For systems with finite band gap, a unitary transforma-tion can be carried out to localize the KS orbitals in realspace as | φ k i = N orb X i | ψ i i u ik , (4)where u ik ≡ ( U ) ik and U is a unitary matrix. Now, bya screening procedure, where least contributing orbitalsare neglected, we expect to achieve a substantial decreasein the number of orbitals involved in the evaluation of { V X | ψ i i} .In our present study, we employed the SCDMmethod to localize the canonical KS orbitals in realspace. In our computations, we employed the followingcutoff criteria Z d r (cid:12)(cid:12) φ i ( r ) φ ∗ j ( r ) (cid:12)(cid:12) > ρ cut , to screen the orbital pairs entering in Equation (1). Anorbital pair i - j is only considered during the computationof the V ACEX operator if the above criteria is satisfied.The ACE operator computed through the above men-tioned screening scheme is denoted as V s-ACEX (screenedACE operator).Now, similar to our earlier work, we split the indi-vidual ionic force components (in hybrid functional basedAIMD) for a system containing N atoms as F exact K = F s-ACE K + ∆ F K , K = 1 , · · · , N (5)with ∆ F K = (cid:0) F exact K − F s-ACE K (cid:1) . Here, F exact and F s-ACE are the ionic force computed with the full rankexchange operator V X and the low rank V s-ACEX opera-tor constructed (only) at the beginning of SCF cycles,respectively. As the V s-ACEX operator closely resemblesthe V X operator, the differences in the ionic force com-ponents of F exact and F s-ACE are very small. Now wemake the assumption that F s-ACE is the fast force and∆ F is the slow force. We have verified this assumptionduring our benchmark calculations. Finally, we employthe reversible reference system propagator algorithm (r-RESPA) to compute the computationally costly ∆ F forces less frequently as compared to the cheaper F s-ACE forces, thereby speeding-up the calculations. Specifically, F s-ACE is computed at every small time step δt and ∆ F is computed at every n steps (i.e., with a time step∆ t = nδt ). Here, the small time step δt can be cho-sen as per the time scale of fast forces ( F s-ACE ) that arecheaper to compute. III. COMPUTATIONAL DETAILS
All the calculations presented here were carried outemploying a modified version of the
CPMD program where the proposed s-MTACE method has been imple-mented. We used PBE0 and PBE functionals for allthe calculations. The effect of the core electrons were in-corporated using the norm-conserving Troullier-Martintype pseudopotentials, generated for PBE functional.For expanding the wavefunctions in PW basis set, a cut-off energy of 80 Ry was used. We carried out Born-Oppenheimer molecular dynamics (BOMD) simulationsto perform MD simulations at microcanonical ( N V E )and canonical (
N V T ) ensembles. Nos´e–Hoover chainthermostat was employed to perform canonical ensem-ble AIMD simulation at 300 K temperature. The wave-function convergence criteria was set to 10 − a.u. forthe orbital gradients during the SCF iterations. For theinitial guess of the wavefunctions at every MD step, Al-ways Stable Predictor Corrector Extrapolation scheme of order 5 was used. A. Benchmark Calculations: 32 Water System
We carried out benchmark calculations for bulk water,modelled using 32 water molecules taken in a periodicsupercell of dimensions 9.85 ˚A × × ∼ − . Two classes ofsimulations were performed to benchmark the efficiencyand the accuracy of our method: FIG. 1. The formation of the gem-diol intermediate ( P ) fromthe reactant ( R ) during the methyl formate hydrolysis in neu-tral aqueous medium. The attacking water molecule is shownin green color. HSer O N z H z Tyr O h Lys H z NH Lys
HHSNOR R OO H z KYK + Y - N z H z Tyr O h Lys H z NH Lys
HHH z HSer O SNOR R OO (a)(b) Lys315Tyr150Lys67Ser64Cephalothin
FIG. 2. (a) Proton transfer from protonated Lys to Tyr in presence of cephalothin ( K + Y − → KY ). The residues in-volved in the proton transfer are shown in magenta color,and the hydrogen atoms involved are shown in green color.(b) A representative snapshot of the solvated non-covalentlycomplexed drug-substrate (class-C β − lactamase bound tothe cephalothin drug molecule) in K + Y − state. Active siteresidues of class-C β − lactamase and the cephalothin drugmolecule are shown in CPK style. The atoms shown in glossyspheres were treated by QM. Atom color code: C (black), O(red), N (blue), S (yellow), H (white), and capping H (green). (a) VV : Here, we performed BOMD simulations inNVE/NVT ensemble with the standard velocity Verletscheme with a timestep of ∆ t fs. The value of ∆ t usedfor the runs will be specified later in the below sections.(b) MTS-n : Here, we performed BOMD simulationsin NVE/NVT ensemble with s-MTACE scheme with an
FIG. 3. The SCDM-localized orbitals { φ j } (wireframe volu-metric representation) considered in the computation of v ij for a given SCDM-localized orbital φ i (solid-red volumetricrepresentation) while using the cutoff ρ cut = 2 . × − fora periodic system containing 64 water molecules. Atom colorcode: Red (O), White (H). Isovalue used for the volumetricrepresentation is 1 × − a.u. inner time step, δt = 0 . t =0 . × n fs.Additionally, for these MTS-n runs, three differentcutoffs ( ρ cut ) were chosen for constructing V s-ACEX oper-ators: ρ cut = 2 . × − , . × − , and 2 . × − . Forobtaining a realistic picture, we mention in passing thatwhile using ρ cut = 2 . × − , Equation 2 is mostly com-puted over the localized orbitals within the first solvationshell (Figure 3). B. Application: Methyl Formate Hydrolysis
A cubic periodic simulation cell with side length of10.1 ˚A was chosen for modelling methyl formate hy-drolysis in neutral water. The system contained onemethyl formate molecule and 32 water molecules. Weperformed two sets of simulations with the GGA/PBEand hybrid/PBE0 XC functionals. BOMD simulationswere carried out to perform MD simulations in canoni-cal ensemble at T = 300K temperature employing Nos´e-Hoover chain thermostats. For the PBE calculations,the standard velocity Verlet scheme was employed witha timestep of ∆ t = 0 .
48 fs. Whereas, for the PBE0 runs,s-MTACE with δt = 0 .
48 fs and ∆ t = 7 . n = 15)were performed. We used ρ cut = 2 . × − for the screen-ing of the SCDM orbital pairs.Well-sliced metadynamics (WS-MTD) approach wasemployed to compute the free energy surface of themethyl formate hydrolysis reaction. Two collective vari-ables (CVs), s = { s , s } , were chosen to explore the freeenergy surface of the reaction. The distance betweenthe carbonyl carbon atom (C) of the methyl formate andthe oxygen atom (O ) of the attacking water molecule ( d [C–O ]) is chosen as the first CV ( s ), see Figure 1.Sampling along this CV was performed in a controlledmanner using the umbrella sampling like bias potential W h ( s ) = 12 κ h (cid:16) s − d (0) h (cid:17) , h = 1 , ..., M . (6)Here, M is the total number of umbrella windows used. κ h and d (0) h are the restraining force constant and theequilibrium value of the h -th umbrella restraint, respec-tively.The coordination number ( CN ) of the oxygen atom(O ) of the attacking water with all the hydrogen atoms(H w ) of the solvent water molecules is chosen as the sec-ond CV ( s ) CN [O : H w ] = N Hw X i =1
11 + ( d i /d ) , (7)with d = 1 .
30 ˚A. Here, N H w is the total number ofH w atoms and d i is the distance between the O atomand the i -th H w atom. This CV was sampled employ-ing the well tempered metadynamics (WT-MTD) biaspotential V b ( s , t ) = X τ 05] ˚A) for PBE0 and PBE based sim-ulations, respectively. C. Application: Protonation State of Active Site Residuesof Class-C β -Lactamase The equilibrated initial structure for the sol-vated non-covalently complexed drug-substrate (class-C β − lactamase bound to the cephalothin drug molecule)was taken from an earlier work by our group; See Fig-ure 2(b). The protein-drug complex was solvated with13473 TIP3P water molecules in a periodic box of size78 × × 76 ˚A . We used the CPMD/GROMOS in-terface as available in the CPMD program to carry outthe hybrid QM/MM canonical ensemble MD simulations.The protein side chains of Lys , Tyr , Lys , Ser ,and the cephalothin molecule (except the thiophene ring)were treated by QM, and the remaining part of the sys-tem including the solvent molecules were modelled byMM. Capping hydrogen atoms were used whenever theQM/MM boundary cleaves a chemical bond. In to-tal 66 atoms were taken in a QM box with dimensions18 × × 22 ˚A ; see Figure 2(b) for details. The QMpart was treated using KS-DFT and PW. The core elec-trons of all the QM atoms were accounted using norm-conserving Troullier-Martin type pseudopotentials at thelevel of PBE and PW cutoff of 70 Ry was taken. MMpart was treated using parm99 AMBER force-field.The QM/MM interaction was handled using the elec-tronic coupling scheme developed by Laio et al . TheQM charge density interacts explicitly with any MMpoint charge within a range of 15 ˚A, beyond which theinteraction was accounted by considering only the multi-pole expansion of the charge density up to the quadrapoleterm. Nearest neighbor lists were updated every 50 steps,and the long-range interaction cutoff of 20 ˚A was taken.QM/MM BOMD simulations were carried out and thetemperature of the system was maintained at 300K us-ing two separate Nos´e-Hoover chain thermostats for theQM and MM subsystems.We performed two sets of simulations: (a) PBE : thestandard velocity Verlet scheme was employed with atimestep of ∆ t = 0 . 48 fs and the PBE XC functional wasused. (b) PBE0 : the s-MTACE scheme with δt = 0 . 48 fsand ∆ t = 7 . n = 15) were considered andthe hybrid PBE0 XC functionals was used. We took ρ cut = 2 . × − for screening the SCDM-localized or-bital pairs.Umbrella sampling approach was employed to modelthe proton transfer reaction within the active siteresidues of class-C β -lactamase enzyme. In particular,we compute here the free energy barrier for the pro-ton transfer between the Lys N ζ and Tyr O η atoms -0.06-0.030.000.030.06 F x on O a t o m ( au ) F s-ACE D F Time (fs) -0.04-0.020.000.020.04 F x on H a t o m ( au ) (a)(b) FIG. 4. One of the components of F s-ACE and ∆ F on anarbitrarily chosen (a) oxygen and (b) hydrogen atoms for asystem containing 32 water molecules in a periodic box duringAIMD simulation using PBE0 functional. Screening of theSCDMs-localized orbitals was performed with ρ cut = 2 . × − . of class-C β -lactamase in the presence of the substrate;see Figure 2(a). The CN of the Tyr O η atom with allthe three Lys H ζ atoms ( CN [Tyr O η :Lys H ζ ]) wasbiased using the umbrella potential as in Equation (6).The definition of the CN is similar to Equation (7). Atotal of 24 windows were placed in the range 0.10 to 0.90.Structures with CN lying between 0.0 and 0.5 resemble K + Y − , while those between 0.5 and 1.0 are similar to KY state; See Figure 2(a). The details of the umbrellabias potentials are reported in Appendix B2. We gener-ated 24 (35) ps long trajectories per each umbrella win-dow during PBE0 (PBE) based QM/MM simulations.The biased probability distributions obtained from theseindependent simulations were then combined using theWHAM technique to get the Boltzmann reweightedprobability distribution, and hence the free energy. IV. RESULTS AND DISCUSSIONA. Benchmark Calculations To test the accuracy and efficiency of the s-MTACEmethod, we considered a periodic system containing 32water molecules modelling bulk water. In Figures 4(a)and (b), we have shown one component of F s-ACE and∆ F (calculated every n = ∆ t/δt steps) on an arbitrarilychosen oxygen and hydrogen atom, respectively. It is ap-parent that ∆ F is slowly varying as compared to F s-ACE .Additionally, the magnitude of ∆ F is ∼ 100 times smalleras compared to F s-ACE . Thus we conclude that the timescale separation implied in the s-MTACE method is rea-sonable, as seen in the case of MTACE method. It is apparent that the efficiency of the MTS schemescrucially depends on the choice of the parameter n . Inpractise, one has to determine an optimal value of n bymonitoring the drift in total energy. For this purpose, wecompared the total energy fluctuations in velocity Verlet( VV ) and MTS runs ( MTS-n ) with different ∆ t for 32-water system in N V E ensemble. To measure the magni- FIG. 5. Measure of the fluctuations in total energy,log (∆ E ), for different ∆ t values in VV and MTS-n simu-lations using PBE0 functional. ∆ E is calculated using Equa-tion (11) over 5 ps long trajectories. In these runs, thescreening of the SCDM-localized orbitals was performed with ρ cut = 2 . × − . g O - O ( r ) g O - H ( r ) VVMTS-5MTS-7MTS-150 1 2 3 4 5 r (Å) g H - H ( r ) Frequency (cm -1 ) I n t en s i t y ( a r b . un i t ) (a)(b)(c)(d) r cut =2.5×10 -2 FIG. 6. Radial distribution functions (RDFs) from VV , MTS-5 , MTS-7 and MTS-15 simulations of bulk water us-ing PBE0 functional: (a) O-O, (b) O-H, and (c) H-H. (d)Power spectrum of the system computed from these simula-tions is also presented. tude of the total energy ( E ) fluctuations, we computedthe quantity ∆ E = (cid:28)(cid:12)(cid:12)(cid:12)(cid:12) E − h E ih E i (cid:12)(cid:12)(cid:12)(cid:12)(cid:29) , (11)for various values of n . In Figure 5, log (∆ E ) is plot-ted as a function of the timestep ∆ t for VV and MTS-n runs, similar to the analysis done in the earlier works. In the case of VV runs, it is clear that the total energyfluctuations increase with ∆ t . For the MTS-n runs, theinner timestep δt is kept fixed at ∼ t = n δt was varied. In Figure 5, the slope of the curve is smaller for the MTS-n runs as comparedto the VV runs, indicating that the effect of increasingtimestep is less profound in the MTS runs. The sta-bility of the MTS-n runs with n = 15 (i.e. MTS-15 )demonstrates that it is a good choice for production runs.Whereas, MTS-n simulations with n greater than 15were showing substantially high total energy fluctuationsand therefore, they were not considered further. Sameconclusions about the optimal n value were obtainedwhile using the different ρ cut values (see Appendix C).To benchmark the accuracy, stability, and efficiency ofour proposed method, we performed a few sets of N V T simulations with 32-water bulk system. The details ofthe simulation length, average temperature and drift intotal energy for these VV and MTS-n simulations arereported in Appendix D. The accuracy of the method isbenchmarked by comparing the structural and dynam-ical properties of bulk water obtained from VV (whichinvokes no assumption) and MTS-n simulations. In par-ticular, for comparing the structural properties, partialradial distribution functions (RDFs) were computed; SeeFigure 6(a)-(c). It is clear that the location of the peaksand the peak heights of the RDFs from the MTS-n simulations are in excellent agreement with those fromthe VV run. For comparing the dynamical properties,we computed the power spectrum for the same system(shown in Figure 6(d)) by taking the Fourier transformof the velocity auto-correlation function. We have ob-served that the frequencies and the intensities of the spec-tra from VV and MTS-n simulations are in excellentagreement with each other. These results demonstratethat the s-MTACE scheme provides an accurate descrip-tion of structural and dynamical properties. We haveobserved a similar trend for different values of ρ cut ; SeeAppendix E.Now, we study the computational efficiency of ourmethod for which we have compared the average compu-tational time for generating 1 fs trajectory with variousparameter settings on an identical computing platform.The CPU time and the achieved speed-ups for various MTS-n runs as compared to the VV run are reportedin Table I for periodic systems with 32, 64 and 128 watermolecules. Here, we defined the speed-up for a MTS-n run as the ratio between the CPU time per fs in VV ( t VV ) and MTS-n runs ( t MTS-n ):speed-up = t VV t MTS-n . (12)The time taken per fs in a MTS-n run is calculated as: t MTS-n = t forceexact + n t forces-ACE n (cid:18) δt fs (cid:19) . (13)Here, t forceexact and t forces-ACE are the times taken for F exact and F s-ACE force calculations, respectively. In every n MDsteps, the force using the exact exchange operator (i.e. F exact ) is calculated only once, and the force using theapproximate s-ACE operator (i.e. F s-ACE ) is computed n times. The total CPU time for the computation of the F exact and F s-ACE forces are decomposed into variouscontributions in Table II.From Table II, it is clear that by the screening of or-bitals, the computational time required for the construc-tion of the ACE operator has significantly decreased, re-sulting in a net speed up in the force calculation perMD step; see t forces-ACE in Table II for various values of ρ cut and Figure 7. However, we have noticed that the num-ber of SCF cycles required for computing energy/forces( N SCF ) increases with screening, due to poor initial guessof wavefunctions. Due to this reason, the overall speed-up in generating 1 fs long trajectory (on average) is notproportional to n (See Table I). For a 32 water peri-odic system, we could achieve the maximum speed up(of ∼ ρ cut = 2 . × − with MTS-15 . Theperformance was better compared to lower values of ρ cut for the same value of n (Figure 8). As expected, MTS-15 performed better than MTS-7 and MTS-5 . Speed-upsof ∼ 11X and ∼ 13X were observed for 64 and 128 watermolecules systems, respectively. w i t hou t sc r een i ng r c u t = . · − r c u t = . · − r c u t = . · − C P U t i m e ( s ) FIG. 7. Average computing time for the construction of V s-ACEX for various values of ρ cut for a periodic system con-taining 32 water molecules and using 120 identical processors;see also Table II. B. Application: Methyl Formate Hydrolysis The hydrolysis of carboxylic esters is one of the wellstudied reaction in the field of chemistry, biochem-istry and industrial chemistry. Methyl formate hy-drolysis serves as the simplest model for the hydrol-ysis of carboxylic esters and has been well studiedexperimentally and theoretically. In our presentwork, we modelled the methyl formate hydrolysis reac-tion in neutral water. This reaction happens throughproton transfer between solvent and the solute. Neutralhydrolysis of methyl formate leads to the formation offormic acid and methanol. It follows a stepwise mecha-nism, where addition of a solvent water molecule to thecarbonyl group results in a stable gem-diol intermediate(see Figure 1), followed by the decomposition of the in-termediate to the final products. In our work, wefocus only on the elementary step leading to the forma-tion of the gem-diol intermediate. The free energetics 32 water w i t hou t sc r een i ng r c u t = . · − r c u t = . · − r c u t = . · − C P U t i m e pe r f s t r a j e c t o r y ( s ) FIG. 8. Average computational time for generating 1 fs tra-jectory with PBE0 for a system containing 32 water moleculesin a periodic box using the exact V X operator ( VV ), and s-MTACE with n = 15 ( MTS-15 ) for various values of ρ cut .All the computations were performed using identical 120 pro-cessors. The computational time reported here is averagedover 630 MD steps. CN [ O : H w ] D F ( kc a l m o l − ) d [C−O ] (Å) 1 1.5 2 2.5 CN [ O : H w ] D F ( kc a l m o l − ) (a) PBE(b) PBE0PP RR FIG. 9. Free energy surfaces computed using (a) PBE and(b) PBE0 density functionals for the R → P reaction are pre-sented. Contour lines are drawn every 2 kcal mol − . for this reaction step were computed using the WS-MTDtechnique.For reweighting the WT-MTD bias potential, we used t min = 20 ps and t max = 35 ps. For the umbrella win-dows near the transition state region, we took t max as46 (50) ps for the PBE0 (PBE) simulations to obtainbetter statistics. Here, t min was chosen such that thefree energy estimate is independent of the initial config-uration. The time series plots of the CVs for some ofthese umbrella windows are shown in Appendix F. Thereconstructed free energy surfaces are shown in Figure 9.The convergence of free energy barriers as a function of TABLE I. Average computational time for generating 1 fs trajectory ( t CPU ) for periodic systems containing 32, 64 and 128water molecules using PBE0. The averages were calculated over N MD number of MD steps. The achieved speed-ups for various MTS-n runs are compared with the VV run. System size No. of CPU cores Method N MD ρ cut = 2 . × − ρ cut = 1 . × − ρ cut = 2 . × − t CPU (s) speed-up t CPU (s) speed-up t CPU (s) speed-up32 water 120 VV 630 518 1.0 509 1.0 512 1.0 MTS-5 630 127 4.1 129 4.0 136 3.8 MTS-7 630 103 5.0 98 5.2 101 5.1 MTS-15 630 70 7.4 57 8.9 56 9.264 water 160 VV 105 3962 1.0 3969 1.0 3976 1.0 MTS-5 105 856 4.6 936 4.2 1059 3.8 MTS-7 105 656 6.0 692 5.7 792 5.0 MTS-15 105 396 10.0 371 10.7 402 9.9128 water 200 VV MTS-5 15 6400 5.1 7607 4.3 8525 3.8 MTS-7 21 4914 6.6 5562 5.9 6020 5.4 MTS-15 45 2611 12.5 2749 11.9 2998 10.9TABLE II. The decomposition of the total computing time for force calculations for a periodic system containing 32 watermolecules. Timings are in seconds when using identical 120 CPU cores. Here, t Cons-ACE is the computing time for the constructionof V s-ACEX , t s-ACE is the computing time per SCF cycle using V s-ACEX , N SCFs-ACE is the average number of SCF cycles during thecomputation of F s-ACE , t forces-ACE is the total computing time for the calculation of F s-ACE , t exact is the computing time per SCFcycle using V X , N SCFexact is the average number of SCF cycles for F exact force calculation, t forceexact is the total computing time forthe calculation of F exact , t forces-ACE = t Cons-ACE + N SCFs-ACE t s-ACE , and t forceexact = N SCFexact t exact .Method ρ cut t Cons-ACE t s-ACE N SCFs-ACE t forces-ACE t exact N SCFexact t forceexact VV – – – – – 23.6 10.4 245.4 MTS-15 without screening 23.6 0.15 10.1 25.1 23.6 8.0 188.8 MTS-15 . × − MTS-15 . × − MTS-15 . × − F ‡ ) of the reaction R → P (Figure 1) using PBE and PBE0 functionals are comparedwith the experimental estimate.Method ∆ F ‡ (kcal mol − )PBE 18.4PBE0 20.8Experiment t max is reported in Appendix G1. The converged free en-ergy barriers and the experimental estimate are listed inTable III.From the free energy surfaces computed using bothPBE and PBE0, it can be seen that the locations of theminimum corresponding to the reactant state R and theintermediate P are nearly the same. However, the topol-ogy of the free energy surfaces near the transition stateregion has some difference. Additionally, the free en-ergy barriers differ between the two functionals. Theconverged free energy barrier computed using PBE is18.4 kcal/mol, while that using PBE0 is 20.8 kcal mol − . The PBE0 free energy barrier is about 2.4 kcal mol − higher than that of PBE, in line with the trends seenfor the hydrolysis of formamide in water. Notably,the free energy barrier computed from PBE0 is closerto the experimentally “estimated” free energy barrier of28.8 kcal mol − at 298 K for the reaction under neu-tral condition. For simple amides, the free energy barrierfor neutral hydrolysis was reported to be between 21.9 to23.8 kcal mol − . The PBE0 estimate of free energy barrier deviates sub-stantially ( ∼ − ) from the experimental data.This may be due to several reasons. It is likely thatunder neutral aqueous conditions, the rate determiningstep is the subsequent elimination and not the forma-tion of the gem-diol, different to what is observed un-der acidic condition. Finite size effects can also affectthe results, especially since the transition state config-urations have an additional proton dissolved in water.Clearly, a more detailed investigation is warranted in thisdirection to understand the deviation from the experi-mental data. Nonetheless, the observed effect of hybridfunctional on the free energy barrier is the most signif-icant result of this study. We also note in passing thatthe previous computational studies have reported a rangeof free energy barriers for this reaction at 298 K. Gu-naydin et al. reported that autoionization of water isthe rate-determining step and the free energy barrier waspredicted to be 23.8 kcal mol − (which in turn was takenfrom experiment data). Arroyo et al. reported a freeenergy barrier of 28.17 kcal mol − based on their MDsimulations using empirical potentials. C. Application: Proton Transfer Reactions within theActive Site Residues of Class-C β -Lactamase β –Lactam antibiotics are commonly prescribed againstbacterial infections. They inhibit the bacterial cell–wallsynthesizing enzymes known as penicillin binding pro-teins (PBPs). However, the clinical efficacy of theselife saving drugs is progressively deteriorating due to theemergence of drug–resistance in bacteria, primarily asso-ciated with the expression of β –lactamases. Theseenzymes hydrolyze β -lactam antibiotics in an efficientmanner, preventing the drug molecules to react withPBPs. In the past, our group has contributed to theunderstanding of the mechanism of hydrolysis of differ-ent classes of β –lactamases. Here, we focus onclass–C serine β –lactamases that is majorly responsiblefor hospital acquired infections. It was found that theprotonation states of the active site residues Lys andTyr of class-C β -lactamase play a crucial role in deter-mining the hydrolysis mechanism. In the Michaeliscomplex of cephaloathin and the enzyme, two protona-tion states are possible: K + Y − and KY . In K + Y − ,Lys is in the protonated form, while Tyr is in itsdeprotonated form. On the other hand, both Lys andTyr are in their neutral forms in the KY state (Fig-ure 2(a)). Depending on the protonation state, the gen-eral base that activates Ser could differ, rendering dif-ferent acylation mechanisms.In the previous study, the free energy barrier for theproton transfer between Lys and Tyr ( K + Y − → KY ) was found to be small ( ∼ is hydrogen bonded toSer in the KY state and the former residue acts as thegeneral base. Here, we revisited the same problem andcomputed the free energy barrier separating the KY and K + Y − states. It is known that PBE can underestimatethe proton transfer barriers, and thus studying the aboveproblem using a higher level of theory, especially includ-ing the contributions of HF exchange, is crucial. Wecomputed free energies using QM/MM MD with PBEand PBE0 density functionals, and compared their per-formance.Umbrella sampling simulations were carried out andthe reconstructed free energy surface along the coordi-nate CN [Tyr O η : Lys H ζ ] using the two density func-tionals are shown in Figure 10. The convergence of thefree energy barriers as a function of simulation lengthis reported in Appendix G2. Also, the error in thecomputed free energy was calculated using the methodgiven by Hummer et al. For both the functionals,the location of the reactant minima ( K + Y − ) is around CN ∈ [0 . , . 4] (unitless). Whereas, the product minima ( KY ) is located near CN ∈ [0 . , . 9] (unitless). However,we can notice qualitative and quantitative differences infree energy profiles. At the level of PBE, the free en-ergy barrier for the proton transfer in forward direction( K + Y − → KY ) is 1 . − , and the barrier for theproton transfer in the reverse direction ( KY → K + Y − )is 0 . − . On the other hand, at the PBE0level, the barrier for the forward proton transfer is 1 . − , and the barrier for the reverse proton transferis 2 . − . The error in the free energy estimatesis less than 0.3 kcal mol − . Notably, K + Y − is predictedto be the most stable state with PBE whereas, KY is themost stable state with PBE0. In the earlier study us-ing PBE and QM/MM metadynamics simulation, it wasreported that both of the states are equally stable withonly ∼ − barrier for their inter conversion (inthe presence of the substrate). This is consistent with theresults of the present work using PBE. Although, at thePBE0 level, the most stable state is found to be different,the free energy difference between the two protonationstates and the free energy barriers for proton transfer re-actions are small. This implies that quick proton transferbetween the Tyr and Lys residues can take place atambient temperature. Therefore we conclude that class-C β − lactamase catalyzed acylation mechanism of the hy-drolysis of cephalothin, wherein the deprotonated form ofLys activates Ser , is justified at the PBE0 level. K + Y − KY D F ( kc a l / m o l ) CN [Tyr O h :Lys H z ]PBEPBE0 FIG. 10. Free energy as a function of CN [Tyr O η : Lys H ζ ] coordinate computed for theproton transfer reaction ( K + Y − ⇋ KY ) using umbrellasampling simulations at the level of PBE (blue) and PBE0(red) density functionals. Errors in the free energy estimateswere computed following Ref. 101. V. CONCLUSIONS In summary, we presented a new scheme named s-MTACE to perform efficient hybrid functional basedAIMD simulations using PW basis set. In particular,we have shown that screening the localized orbitals andcomputing ACE operator using the screened orbitals is an0efficient way to speed up our earlier method. Bench-mark results show that dynamic and structural prop-erties can be computed accurately using the s-MTACEmethod. We achieved a computational speed-up of ≈ 10 for a periodic 32-water molecules system compared tothe conventional implementation of the hybrid functionalPW-DFT. For a larger system with 128 water molecules,we have achieved a speed-up of about 13.Using our implementation, we show that it is possi-ble to carry out computationally demanding free energycalculations at the level of hybrid density functionals.We demonstrated this by performing free energy calcula-tions in two systems involving proton transfer reactions.First, our implementation was combined with the WS-MTD method to study the methyl formate hydrolysis inneutral aqueous solution. The two dimensional free en-ergy surface for this reaction was computed with PBEand PBE0 functionals. We observed that the free en-ergy barrier computed using PBE is 2.4 kcal/mol lowerthan PBE0, similar to our observations in earlier stud-ies for a different system. By employing our methodwithin the framework of QM/MM, we investigated theproton transfer reaction between the active site residuesin the Michaelis complex of class-C β − lactamase andcephalothin antibiotic. We found that the stability ofprotonation states and the free energy barriers for protontransfer are altered while switching from PBE to PBE0.The general observation of underestimation of protontransfer barriers by PBE compared to PBE0 is consistentwith the detailed benchmarking studies by Adamo andco-workers through structural optimizations of severalgas-phase reactions. As this work enables us to performfree energy calculations at the level of hybrid function-als within AIMD simulations in an efficient manner, wehope that the presented method and the results obtainedfrom our calculations are of great importance to scientistsworking in the area of computational catalysis. ACKNOWLEDGMENTS The authors thank Tobias Kl¨offel and Bernd Meyerfor helpful discussions. Financial support from the Sci-ence and Engineering Research Board (India) under theMATRICS scheme (Ref. No. MTR/2019/000359) andfrom the German Research Foundation (DFG) throughResearch Unit FOR 1878 (funCOS) and CollaborativeResearch Center SFB 953 (project number 182849149)is gratefully acknowledged. SM and VT thank the Uni-versity Grant Commission (UGC), India, and IITK fortheir Ph.D. fellowships. Computational resources wereprovided by the HPC facility (HPC2013) at IITK, theErlangen Regional Computing Center (RRZE) at FAUand SuperMUC-NG (project pn98fa) at Leibniz Super-computing Centre (LRZ). Appendix A: Construction of the Adaptively CompressedExchange Operator According to the ACE operator formalism, the ACEoperator ( V ACEX ) can be constructed through the follow-ing series of simple linear algebra operations. First, theexact exchange operator ( V X ) has to be applied on theset of KS orbitals {| ψ i i} as | W i i = V X | ψ i i , i = 1 , ...., N orb . (A1)Now, ACE formalism defines V ACEX as V ACEX = N orb X i,j | W i i B ij h W j | . (A2)Here, B = M − , and the matrix M has elements M kl = h ψ k | V X | ψ l i . (A3)Now, the Cholesky factorization of − M gives M = − LL T , (A4)where, L is a lower triangular matrix. Then B is com-puted as B = − L − T L − . (A5)Finally, the V ACEX operator can be rewritten as V ACEX = − N orb X k | P k ih P k | , (A6)with {| P k i} as the columns of the matrix P , which isdefined as P = WL − T . (A7)1 Appendix B: Details of the Umbrella Sampling BiasParameters1. Methyl Formate Hydrolysis TABLE B.1. Umbrella sampling parameters: Here, d (0) h is in˚A and κ h is in Hartree Bohr − . h d (0) h κ h h d (0) h κ h 2. Protonation State of Active Site Residues of Class-C β -Lactamase TABLE B.2. Umbrella sampling parameters: Here, d (0) h isunitless and κ h is in Hartree. h d (0) h κ h h d (0) h κ h Appendix C: Comparison of the Drift in Total Energy forVV and MTS Simulations in the NVE Ensemble TABLE C.1. The drift in total energy and log (∆ E ) forvarious VV and MTS-n simulations in NV E ensemble. Theresults of MTS-n simulations with various ρ cut values arereported.Method ρ cut Timestep Drift a log (∆ E )(fs) (au/ps/atom) VV – 0.48 1.8 × − -6.77 VV – 0.96 1.6 × − -6.14 VV – 1.44 2.7 × − -5.76 MTS-3 . × − × − -5.01 MTS-5 . × − × − -5.41 MTS-7 . × − × − -5.54 MTS-15 . × − × − -5.45 MTS-30 . × − × − -3.79 MTS-3 . × − × − -5.00 MTS-5 . × − × − -5.47 MTS-7 . × − × − -5.11 MTS-15 . × − × − -5.32 MTS-30 . × − × − -3.65 MTS-3 . × − × − -4.92 MTS-5 . × − × − -5.16 MTS-7 . × − × − -5.07 MTS-15 . × − × − -4.84 MTS-30 . × − × − -3.58 a Drift was calculated as |h E ( t ) − E (0) i| / number of atoms / time.Here, E ( t ) is the total energy at any time t . Appendix D: Accuracy of VV and MTS Simulations in theNVT Ensemble TABLE D.1. The simulation time length, average tempera-ture ( h T i ) and drift in total energy for various VV and MTS-n (for various values of ρ cut ) runs in NV T ensemble.Method ρ cut Timestep Length h T i Drift a (fs) (ps) (K) (au/ps/atom) VV – 0.48 10 301 1 . × − MTS-5 . × − . × − MTS-7 . × − . × − MTS-15 . × − . × − MTS-5 . × − . × − MTS-7 . × − . × − MTS-15 . × − . × − MTS-5 . × − . × − MTS-7 . × − . × − MTS-15 . × − . × − Drift was calculated as |h E ( t ) − E (0) i| / number of atoms / time.Here, E ( t ) is the total energy at any time t . Appendix E: Comparison of Structural and DynamicalProperties g O - O ( r ) g O - H ( r ) VVMTS-5MTS-7MTS-150 1 2 3 4 5 r (Å) g H - H ( r ) Frequency (cm -1 ) I n t en s i t y ( a r b . un i t ) (a)(b)(c)(d) r cut =2.0×10 -3 FIG. E.1. Radial distribution functions (RDFs) for bulk wa-ter (32 water system) simulation using VV , MTS-5 , MTS-7 and MTS-15 trajectories at the level of PBE0: (a) O-O, (b)O-H, and (c) H-H. (d) Power spectrum of the same systemcomputed from VV , MTS-5 , MTS-7 and MTS-15 trajec-tories are shown. Screening of the SCDMs was performedwith ρ cut = 2 . × − . g O - O ( r ) g O - H ( r ) VVMTS-5MTS-7MTS-150 1 2 3 4 5 r (Å) g H - H ( r ) Frequency (cm -1 ) I n t en s i t y ( a r b . un i t ) (a)(b)(c)(d) r cut =1.0×10 -2 FIG. E.2. Radial distribution functions (RDFs) for bulk wa-ter (32 water system) simulation using VV , MTS-5 , MTS-7 and MTS-15 trajectories at the level of PBE0: (a) O-O, (b)O-H, and (c) H-H. (d) Power spectrum of the same systemcomputed from VV , MTS-5 , MTS-7 and MTS-15 trajec-tories are shown. Screening of the SCDMs was performedwith ρ cut = 1 . × − . Appendix F: Time Series Plots of the CVs for UmbrellaWindows Near the Transition State Region CN CN CN PBE0PBE0.511.52 CN Time (ps) CN FIG. F.1. Fluctuations in the CN [O : H w ] CV with simula-tion time for the umbrella windows near the transition stateregion (from 1.70 to 1.86 ˚A) using PBE and PBE0 functionals. CN CN CN PBE0PBE11.522.5 CN Time (ps) CN FIG. F.2. Fluctuations in the CN [O : H w ] CV with simula-tion time for the umbrella windows near the transition stateregion (from 1.90 to 2.05 ˚A) using PBE and PBE0 functionals. Appendix G: Convergence of the Free Energy Barriers1. Methyl Formate Hydrolysis 18 19 20 21 22 42 44 46 48 50 38 40 42 44 46 D F ‡ ( kc a l / m o l ) t max (ps) t max (ps) PBEPBE0 FIG. G.1. Convergence of the free energy barrier (∆ F ‡ ) forthe reaction R → P with different t max values using PBE andPBE0 density functionals. Here, t min = 20 ps was taken. 2. Protonation State of Active Site Residues of Class-C β -Lactamase D F ‡ ( kc a l / m o l ) t max (ps) t max (ps) KY fi K + Y - (PBE)K + Y - fi KY(PBE)KY fi K + Y - (PBE0)K + Y - fi KY(PBE0) FIG. G.2. Convergence of the free energy barrier (∆ F ‡ ) withdifferent t max values for the reaction KY → K + Y − usingPBE functional (blue circles) and PBE0 functional (red cir-cles); and K + Y − → KY using PBE functional (blue squares)and PBE0 functional (red squares). Here, t min = 15 ps wastaken for PBE functional; and t min = 4 ps was taken for PBE0functional. REFERENCES M. E. Tuckerman, Statistical Mechanics: Theory and MolecularSimulation (Oxford University Press, Oxford, 2010). A. P. Leach, Molecular Modelling. Principles and Applications (Pearson Education Limited, Upper Saddle River, New Jersey,2001). D. Frenkel and B. Smit, Understanding Molecular Simulation:From Algorithms to Applications (Academic Press, San Diego,2002). M. E. Tuckerman and G. J. Martyna, “Understandingmodern molecular dynamics: Techniques and applications,”J. Phys. Chem. B , 159–178 (2000). D. Marx and J. Hutter, Ab Initio Molecular Dynamics: BasicTheory and Advanced Methods (Cambridge University Press,Cambridge, 2009). R. Iftimie, P. Minary, and M. E. Tuckerman, “Ab initiomolecular dynamics: Concepts, recent developments, and futuretrends,” Proc. Natl. Acad. Sci. U.S.A. , 6654–6659 (2005). M. E. Tuckerman, “Ab initio molecular dynamics: ba-sic concepts, current trends and novel applications,”J. Phys.: Condens. Matter , R1297–R1355 (2002). M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D.Joannopoulos, “Iterative minimization techniques for ab initiototal-energy calculations: molecular dynamics and conjugategradients,” Rev. Mod. Phys. , 1045–1097 (1992). A. D. Becke, “Density-functional exchange-energyapproximation with correct asymptotic behavior,”Phys. Rev. A , 3098–3100 (1988). C. Lee, W. Yang, and R. G. Parr, “Development of the Colle-Salvetti correlation-energy formula into a functional of the elec-tron density,” Phys. Rev. B , 785–789 (1988). J. P. Perdew, K. Burke, and M. Ernzerhof,“Generalized gradient approximation made simple,”Phys. Rev. Lett. , 3865–3868 (1996). W. Koch and M. C. Holthausen, A Chemist’s Guide to DensityFunctional Theory (WILEY-VCH, New York, 2001). A. J. Cohen, P. Mori-S´anchez, and W. Yang, “In-sights into current limitations of density functional theory,”Science , 792–794 (2008). J. P. Perdew and A. Zunger, “Self-interaction correction todensity-functional approximations for many-electron systems,”Phys. Rev. B , 5048–5079 (1981). P. Mori-S´anchez, A. J. Cohen, and W. Yang, “Lo-calization and delocalization errors in density func-tional theory and implications for band-gap prediction,”Phys. Rev. Lett. , 146401 (2008). R. M. Martin, Electronic Structure: Basic Theory and PracticalMethods (Cambridge University Press, Cambridge, 2004). A. D. Becke, “Density-functional thermochemistry. III. The roleof exact exchange,” J. Chem. Phys. , 5648–5652 (1993). J. P. Perdew, M. Ernzerhof, and K. Burke, “Rationale for mix-ing exact exchange with density functional approximations,”J. Chem. Phys. , 9982–9985 (1996). J. Heyd, G. E. Scuseria, and M. Ernzerhof, “Hy-brid functionals based on a screened Coulomb potential,”J. Chem. Phys. , 8207–8215 (2003). C. Adamo and V. Barone, “Toward reliable density functionalmethods without adjustable parameters: The PBE0 model,”J. Chem. Phys. , 6158–6170 (1999). C. J. Cramer and D. G. Truhlar, “Density functional the-ory for transition metals and transition metal chemistry,”Phys. Chem. Chem. Phys. , 10757–10816 (2009). B. G. Janesko, T. M. Henderson, and G. E. Scuseria,“Screened hybrid density functionals for solid-state chemistryand physics,” Phys. Chem. Chem. Phys. , 443–454 (2009). Q. Wan, L. Spanu, F. Gygi, and G. Galli, “Electronic structureof aqueous sulfuric acid from first-principles simulations withhybrid functionals,” J. Phys. Chem. Lett. , 2562–2567 (2014). A. J. Cohen, P. Mori-S´anchez, and W. Yang, “Challenges fordensity functional theory,” Chem. Rev. , 289–320 (2012). H. Xiao, J. Tahir-Kheli, and W. A. Goddard III, “Accurateband gaps for semiconductors from density functional theory,”J. Phys. Chem. Lett. , 212–217 (2011). Q. Zhao and H. J. Kulik, “Stable surfaces that bindtoo tightly: Can range-separated hybrids or DFT+Uimprove paradoxical descriptions of surface chemistry?”J. Phys. Chem. Lett. , 5090–5098 (2019). N. Gerrits, E. W. F. Smeets, S. Vuckovic, A. D. Powell,K. Doblhoff-Dier, and G.-J. Kroes, “Density functional theoryfor molecule–metal surface reactions: When does the general-ized gradient approximation get it right, and what to do if itdoes not,” J. Phys. Chem. Lett. , 10552–10560 (2020). B. G. Janesko and G. E. Scuseria, “Hartree-Fock orbitals signif-icantly improve the reaction barrier heights predicted by semilo-cal density functionals,” J. Chem. Phys. , 244112 (2008). G. F. Mangiatordi, E. Br´emond, and C. Adamo, “DFT andproton transfer reactions: A benchmark study on structure andkinetics,” J. Chem. Theory Comput. , 3082–3088 (2012). T. Todorova, A. P. Seitsonen, J. Hutter, I.-F. W. Kuo, and C. J.Mundy, “Molecular dynamics simulation of liquid water: Hybriddensity functionals,” J. Phys. Chem. B , 3685–3691 (2006). C. Zhang, D. Donadio, F. Gygi, and G. Galli,“First principles simulations of the infrared spec-trum of liquid water using hybrid density functionals,”J. Chem. Theory Comput. , 1443–1449 (2011). R. A. DiStasio Jr., B. Santra, Z. Li, X. Wu, and R. Car, “Theindividual and collective effects of exact exchange and disper-sion interactions on the ab initio structure of liquid water,”J. Chem. Phys. , 084502 (2014). B. Santra, R. A. DiStasio Jr., F. Martelli, andR. Car, “Local structure analysis in ab initio liquid water,”Mol. Phys. , 2829–2841 (2015). A. Bankura, B. Santra, R. A. DiStasio Jr., C. W. Swartz, M. L.Klein, and X. Wu, “A systematic study of chloride ion solvationin water using van der Waals inclusive hybrid density functionaltheory,” Mol. Phys. , 2842–2854 (2015). F. Ambrosio, G. Miceli, and A. Pasquarello, “Structural, dy-namical, and electronic properties of liquid water: A hybridfunctional study,” J. Phys. Chem. B , 7456–7470 (2016). E. Sevgen, F. Giberti, H. Sidky, J. K. Whitmer, G. Galli,F. Gygi, and J. J. de Pablo, “Hierarchical coupling of first-principles molecular dynamics with advanced sampling meth-ods,” J. Chem. Theory Comput. , 2881–2888 (2018). S. Mandal, J. Debnath, B. Meyer, and N. N. Nair, “Enhancedsampling and free energy calculations with hybrid functionalsand plane waves for chemical reactions,” J. Chem. Phys. ,144113 (2018). S. Mandal and N. N. Nair, “Efficient computation of freeenergy surfaces of chemical reactions using ab initio molec-ular dynamics with hybrid functionals and plane waves,”J. Comput. Chem. , 1790–1797 (2020). S. Chawla and G. A. Voth, “Exact exchange in ab initiomolecular dynamics: An efficient plane-wave based algorithm,”J. Chem. Phys. , 4697–4700 (1998). X. Wu, A. Selloni, and R. Car, “Order- N implemen-tation of exact exchange in extended insulating systems,”Phys. Rev. B , 085102 (2009). M. Chen, L. Zheng, B. Santra, H.-Y. Ko, R. A. DiStasio Jr.,M. L. Klein, R. Car, and X. Wu, “Hydroxide diffuses slowerthan hydronium in water because its solvated structure inhibitscorrelated proton transfer,” Nat. Chem. , 413–419 (2018). F. Gygi, “Compact representations of Kohn-Sham invariant sub-spaces,” Phys. Rev. Lett. , 166406 (2009). F. Gygi and I. Duchemin, “Efficient computation ofHartree–Fock exchange using recursive subspace bisection,”J. Chem. Theory Comput. , 582–587 (2013). W. Dawson and F. Gygi, “Performance and ac-curacy of recursive subspace bisection for hy-brid DFT calculations in inhomogeneous systems,”J. Chem. Theory Comput. , 4655–4663 (2015). A. P. Gaiduk, C. Zhang, F. Gygi, and G. Galli, “Structural andelectronic properties of aqueous NaCl solutions from ab initiomolecular dynamics simulations with hybrid density function-als,” Chem. Phys. Lett. , 89 – 96 (2014). A. P. Gaiduk, F. Gygi, and G. Galli, “Density and compressibil-ity of liquid water and ice from first-principles simulations withhybrid functionals,” J. Phys. Chem. Lett. , 2902–2908 (2015). H.-Y. Ko, J. Jia, B. Santra, X. Wu, R. Car, andR. A. DiStasio Jr., “Enabling large-scale condensed-phasehybrid density functional theory based ab initio molec-ular dynamics. 1. theory, algorithm, and performance,” J. Chem. Theory Comput. , 3757–3785 (2020). H.-Y. Ko, B. Santra, and R. A. DiStasio Jr., “Enablinglarge-scale condensed-phase hybrid density functional the-ory based ab initio molecular dynamics II: Extensions tothe isobaric-isoenthalpic and isobaric-isothermal ensembles,”(2020), arXiv:2011.07209 [cond-mat.mtrl-sci]. M. Guidon, F. Schiffmann, J. Hutter, and J. VandeVondele,“Ab initio molecular dynamics using hybrid density function-als,” J. Chem. Phys. , 214104 (2008). E. Liberatore, R. Meli, and U. Rothlisberger, “A versatile mul-tiple time step scheme for efficient ab initio molecular dynamicssimulations,” J. Chem. Theory Comput. , 2834–2842 (2018). S. Fatehi and R. P. Steele, “Multiple-time step ab initiomolecular dynamics based on two-electron integral screening,”J. Chem. Theory Comput. , 884–898 (2015). M. P. Bircher and U. Rothlisberger, “Exploiting coordinatescaling relations to accelerate exact exchange calculations,”J. Phys. Chem. Lett. , 3886–3890 (2018). M. P. Bircher and U. Rothlisberger, “From a week toless than a day: Speedup and scaling of coordinate-scaled exact exchange calculations in plane waves,”Comput. Phys. Commun. , 106943 (2020). V. Weber, C. Bekas, T. Laino, A. Curioni, A. Bertsch,and S. Futral, “Shedding light on Lithium/air batteriesusing millions of threads on the BG/Q supercomputer,” in (Phoenix, AZ, USA, 2014) pp. 735–744. I. Duchemin and F. Gygi, “A scalable and accurate al-gorithm for the computation of Hartree–Fock exchange,”Comput. Phys. Commun. , 855 – 860 (2010). N. Varini, D. Ceresoli, L. Martin-Samos, I. Girotto,and C. Cavazzoni, “Enhancement of DFT-calculationsat petascale: Nuclear magnetic resonance, hybrid den-sity functional theory and Car–Parrinello calculations,”Comput. Phys. Commun. , 1827–1833 (2013). T. A. Barnes, T. Kurth, P. Carrier, N. Wichmann,D. Prendergast, P. R. Kent, and J. Deslippe, “Improvedtreatment of exact exchange in Quantum ESPRESSO,”Comput. Phys. Commun. , 52–58 (2017). V. Bolnykh, J. M. H. Olsen, S. Meloni, M. P. Bircher, E. Ip-politi, P. Carloni, and U. Rothlisberger, “Extreme scalabil-ity of DFT-based QM/MM MD simulations using MiMiC,”J. Chem. Theory Comput. , 5601–5613 (2019). J. Vinson, “Faster exact exchange in peri-odic systems using single-precision arithmetic,”J. Chem. Phys. , 204106 (2020). S. Mandal and N. N. Nair, “Speeding-up ab initio molecu-lar dynamics with hybrid functionals using adaptively com-pressed exchange operator based multiple timestepping,”J. Chem. Phys. , 151102 (2019). L. Lin, “Adaptively compressed exchange operator,”J. Chem. Theory Comput. , 2242–2249 (2016). W. Hu, L. Lin, A. S. Banerjee, E. Vecharynski, andC. Yang, “Adaptively compressed exchange operator forlarge-scale hybrid density functional calculations withapplications to the adsorption of water on silicene,”J. Chem. Theory Comput. , 1188–1198 (2017). A. Damle, L. Lin, and L. Ying, “Compressed representationof Kohn–Sham orbitals via selected columns of the density ma-trix,” J. Chem. Theory Comput. , 1463–1469 (2015). I. Carnimeo, S. Baroni, and P. Giannozzi, “Fast hybriddensity-functional computations using plane-wave basis sets,”Electron. Struct. , 015009 (2019). P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B.Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli,M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso,S. de Gironcoli, P. Delugas, R. A. DiStasio Jr., A. Ferretti,A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerst-mann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko,A. Kokalj, E. K¨u¸c¨ukbenli, M. Lazzeri, M. Marsili, N. Marzari,F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la Roza,L. Paulatto, S. Ponc´e, D. Rocca, R. Sabatini, B. Santra,M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thon-hauser, P. Umari, N. Vast, X. Wu, and S. Baroni, “Advanced ca- pabilities for materials modelling with Quantum ESPRESSO,”J. Phys.: Condens. Matter , 465901 (2017). M. Tuckerman, B. J. Berne, and G. J. Martyna,“Reversible multiple time scale molecular dynamics,”J. Chem. Phys. , 1990–2001 (1992). J. Hutter et al., CPMD: An Ab Initio Electronic Structure andMolecular Dynamics Program , see (ac-cessed on December 31, 2020). T. Kl¨offel, G. Mathias, and B. Meyer, “Integrating stateof the art compute, communication, and autotuning strate-gies to multiply the performance of ab initio moleculardynamics on massively parallel multi-core supercomputers,”Comput. Phys. Commun. , 107745 (2021). N. Troullier and J. L. Martins, “Efficient pseudopotentials forplane-wave calculations,” Phys. Rev. B , 1993–2006 (1991). G. J. Martyna, M. L. Klein, and M. Tuckerman, “Nos´e–Hooverchains: The canonical ensemble via continuous dynamics,”J. Chem. Phys. , 2635–2643 (1992). J. Kolafa, “Time-reversible always stable predictor–correctormethod for molecular dynamics of polarizable molecules,”J. Comput. Chem. , 335–342 (2004). S. Awasthi, V. Kapil, and N. N. Nair, “Sampling free energysurfaces as slices by combining umbrella sampling and metady-namics,” J. Comput. Chem. , 1413–1424 (2016). G. M. Torrie and J. P. Valleau, “Monte Carlo free en-ergy estimates using non-Boltzmann sampling: Ap-plication to the sub-critical Lennard-Jones fluid,”Chem. Phys. Lett. , 578 – 581 (1974). A. Barducci, G. Bussi, and M. Parrinello, “Well-temperedmetadynamics: A smoothly converging and tunable free-energymethod,” Phys. Rev. Lett. , 020603 (2008). S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen,and P. A. Kollman, “The weighted histogram analysis methodfor free-energy calculations on biomolecules. I. The method,”J. Comput. Chem. , 1011–1021 (1992). R. Tripathi and N. N. Nair, “Thermodynamic and kinetic sta-bilities of active site protonation states of class–c β -lactamase,”J. Phys. Chem. B , 4741–4753 (2012). T. E. Cheatham III, P. Cieplak, and P. A. Kollman,“A modified version of the Cornell et al. force fieldwith improved sugar pucker phases and helical repeat,”J. Biomol. Struct. Dyn. , 845–862 (1999). A. Laio, J. VandeVondele, and U. Rothlisberger,“A Hamiltonian electrostatic coupling scheme for hy-brid Car-Parrinello molecular dynamics simulations,”J. Chem. Phys. , 6941–6947 (2002). Intel Xeon CPU E5-2670v2 (20 cores/node) with 128 GB ofRAM per node . J. P. Guthrie, “Hydration of carboxylic acids and esters.evaluation of the free energy change for addition of wa-ter to acetic and formic acids and their methyl esters,”J. Am. Chem. Soc. , 6999–7003 (1973). J. F. Marlier, T. G. Frey, J. A. Mallory, and W. W. Cleland,“Multiple isotope effect study of the acid-catalyzed hydrolysisof methyl formate,” J. Org. Chem. , 1737–1744 (2005). O. Jogunola, T. Salmi, K. Er¨anen, J. W¨arn˚a, M. Kan-gas, and J.-P. Mikkola, “Reversible autocatalytic hydrol-ysis of alkyl formate: Kinetic and reactor modeling,”Ind. Eng. Chem. Res. , 4099–4106 (2010). T. D. Vu, A. Seidel-Morgenstern, S. Gr¨uner, and A. Kienle,“Analysis of ester hydrolysis reactions in a chromato-graphic reactor using equilibrium theory and a rate model,”Ind. Eng. Chem. Res. , 9565–9574 (2005). J. P. Guthrie, “Hydration of thioesters. evaluation of thefree-energy changes for the addition of water to somethioesters, rate-equilibrium correlations over very wide rangesin equilibrium constants, and a new mechanistic criterion,”J. Am. Chem. Soc. , 5892–5904 (1978). C.-G. Zhan, D. W. Landry, and R. L. Ornstein,“Reaction pathways and energy barriers for alkaline hy-drolysis of carboxylic acid esters in water studied bya hybrid supermolecule-polarizable continuum approach,”J. Am. Chem. Soc. , 2621–2627 (2000). J. R. Pliego, Jr. and J. M. Riveros, “A theoretical analy-sis of the free-energy profile of the different pathways in thealkaline hydrolysis of methyl formate in aqueous solution,”Chem. Eur. J. , 1945–1953 (2002). J. Pranata, “Ab initio study of the base-catalyzed hydrolysis ofmethyl formate,” J. Phys. Chem. , 1180–1184 (1994). H. Gunaydin and K. N. Houk, “Molecular dynamics pre-diction of the mechanism of ester hydrolysis in water,”J. Am. Chem. Soc. , 15232–15233 (2008). S. Tolosa Arroyo, A. Hidalgo Garc´ıa, M. Moreno Alvero, andJ. A. Sans´on Mart´ın, “Theoretical study of the neutral hydroly-sis of methyl formate via a concerted and stepwise water-assistedmechanism using free-energy curves and molecular dynamicssimulation.” Struct. Chem. , 909–915 (2011). B. P. Callahan, Y. Yuan, and R. Wolfenden, “The burden borneby urease,” J. Am. Chem. Soc. , 10828–10829 (2005). A. L. L. East, “On the hydrolysis mechanisms of amides andpeptides,” Int. J. Chem. Kinet. , 705–709 (2018). E. Y. Klein, T. P. Van Boeckel, E. M. Martinez, S. Pant, S. Gan-dra, S. A. Levin, H. Goossens, and R. Laxminarayan, “Globalincrease and geographic convergence in antibiotic consumptionbetween 2000 and 2015,” Proc. Nat. Acad. Sci. , E3463–E3470 (2018). C. L. Tooke, P. Hinchliffe, E. C. Bragginton, C. K. Colenso,V. H. Hirvonen, Y. Takebayashi, and J. Spencer, “ β -lactamases and β -lactamase inhibitors in the 21st century,”J. Mol. Biol. , 3472 – 3500 (2019). K. Bush, “Bench-to-bedside review: The role of β -lactamasesin antibiotic-resistant Gram-negative infections,” Crit. Care ,224 (2010). C. Bebrone, P. Lassaux, L. Vercheval, J.-S. Sohier, A. Jehaes,E. Sauvage, and M. Galleni, “Current challenges in antimicro-bial chemotherapy,” Drugs , 651–679 (2010). R. Tripathi and N. N. Nair, “Deacylation mechanism and ki-netics of acyl–enzyme complex of class–c β –lactamase andcephalothin,” J. Phys. Chem. B , 2681–2690 (2016). R. Tripathi and N. N. Nair, “Mechanism of acyl–enzyme com-plex formation from the Henry–Michaelis complex of class–c β -lactamases with β –lactam antibiotics,” J. Am. Chem. Soc. ,14679–14690 (2013). S. Awasthi, S. Gupta, R. Tripathi, and N. N. Nair, “Mechanismand kinetics of aztreonam hydrolysis catalyzed by class–c β –lactamase: A temperature-accelerated sliced sampling study,”J. Phys. Chem. B , 4299–4308 (2018). C. K. Das and N. N. Nair, “Elucidating the molecular basis ofavibactam mediated inhibition of class a β –lactamases,” Chem.Eur. J. , 9639–9651 (2020). C. K. Das and N. N. Nair, “Hydrolysis of cephalexin andmeropenem by New Delhi metallo– β –lactamase: the substrateprotonation mechanism is drug dependent,” Phys. Chem. Chem.Phys. , 13111–13121 (2017). F. Zhu and G. Hummer, “Convergence and error estimation infree energy calculations using the weighted histogram analysismethod,” J. Comput. Chem.33