Multiple-Timestep ab Initio Molecular Dynamics Using an Atomic Basis

Article pubs.acs.org/JPCA

Multiple-Timestep ab Initio Molecular Dynamics Using an Atomic Basis Set Partitioning Ryan P. Steele*

Downloaded by CENTRAL MICHIGAN UNIV on September 10, 2015 | http://pubs.acs.org Publication Date (Web): September 9, 2015 | doi: 10.1021/acs.jpca.5b05850

Department of Chemistry and Henry Eyring Center for Theoretical Chemistry, University of Utah, 315 South 1400 East, Salt Lake City, Utah 84112, United States ABSTRACT: This work describes an approach to accelerate ab initio Born− Oppenheimer molecular dynamics (MD) simulations by exploiting the inherent timescale separation between contributions from different atom-centered Gaussian basis sets. Several MD steps are propagated with a cost-efficient, low-level basis set, after which a dynamical correction accounts for large basis set relaxation effects in a time-reversible fashion. This multiple-timestep scheme is shown to generate valid MD trajectories, on the basis of rigorous testing for water clusters, the methanol dimer, an alanine polypeptide, protonated hydrazine, and the oxidized water dimer. This new approach generates observables that are consistent with those of target basis set trajectories, including MDbased vibrational spectra. This protocol is shown to be valid for Hartree−Fock, density functional theory, and second-order Møller−Plesset perturbation theory approaches. Recommended pairings include 6-31G as a low-level basis set for 6-31G** or 6-311G**, as well as cc-pVDZ as the subset for accurate dynamics with aug-cc-pVTZ. Demonstrated cost savings include factors of 2.6−7.3 on the systems tested and are expected to remain valid across system sizes.

■

INTRODUCTION Multiple-timestep molecular dynamics (MTS-MD) techniques1−21 take advantage of the inherent timescale separation found in molecular potentials. Rapidly fluctuating components of the underlying potential/force are updated at each timestep of an MD simulation, whereas the contribution of slowly varying components is included on a concomitantly longer timescale. In the often-encountered regime in which the slowly varying contributions are also long-ranged (and, therefore, computationally demanding), this MTS-MD partitioning leads to appreciable computational savings. Simulations of condensed-phase and biomolecular systems, for example, update the long-range electrostatic contributions on a timescale that is appreciably longer than the natural MD timestep dictated by the highest-frequency covalent stretches.1,2,7,12,13,16,17,19 Such techniques are now well established and are typically based on the time-reversible r-RESPA algorithm.4 The development of analogous protocols for electronic structure theory-based MD is less clear. Good motivation exists to do so, however, as ab initio MD (AIMD) can properly account for cases of bond-breaking, strong polarization, and charge-transfer effects. The central design challenge for multiple-timestep ab initio molecular dynamics (MTS-AIMD) is the lack of obvious separability in the underlying molecular forces. In electronic structure theory (EST), no partitioning of the potential function into short-/long-range contributions is enforced, nor is the potential forcibly segregated into bonded and nonbonded interactions. This lack of user intervention is, of course, one of the main allures of EST; fundamental chemical interactions (or lack thereof) are assigned on the basis of the principles of quantum mechanics, rather than by the © XXXX American Chemical Society

guiding hand of the user. This lack of separability, however, also poses new challenges for MTS-AIMD and requires rethinking the appropriate partitionings for efficient MTS protocols.22−28 Many of the approaches to molecular electronic structure involve two broad approximations,29 which include (a) the nelectron expansion of the total wavefunction to incorporate electron correlation and (b) the 1-electron expansion of molecular orbitals (MOs) in a computationally tractable basis set, such as atom-centered Gaussian basis functions. Recent work from our group showed that MTS-AIMD is viable when taking advantage of the timescale disparity of the former expansion.28 In particular, Hartree−Fock (HF) can effectively be used as the inner-timestep method with second-order perturbation theory (MP2) as the outer-timestep method in an MTS-AIMD scheme. The main question to be addressed in the present work is whether timescale separability also occurs in the latter approximation (the expansion of MOs in a basis set) and whether this property can be exploited for accelerated AIMD simulations. The accuracy of EST calculations depends critically on the quality of this basis set expansion,29−35 and many now-standard atom-centered basis sets have been developed throughout the history of EST to enable such calculations.36−39 Furthermore, the accuracy of the n-electron expansion is intrinsically tied to Special Issue: Dynamics of Molecular Collisions XXV: Fifty Years of Chemical Reaction Dynamics Received: June 18, 2015 Revised: July 31, 2015

A

DOI: 10.1021/acs.jpca.5b05850 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Table 1. Harmonic Frequencies (cm−1) for the H-Bonded O−H Stretch of the Methanol Dimer, Using the B3LYP Hybrid Density Functional

the quality of the basis set because the inclusion of virtual orbitals requires functions with higher angular momentum and larger spatial extent. This accuracy comes at a steep cost, however. Density functional theory (DFT) and HF calculations, for example, scale with the fourth power of the size of the basis set [O(N4)]. Advanced algorithms that scale linearly with respect to system size40−42 still do not change this scaling with respect to basis set size, and in fact, extended basis sets delay the onset of this linearly scaling regime to even larger chemical systems. Therefore, the development of new approaches that include accurate basis sets without undue computational burden is critical for the continued extension of AIMD to chemically relevant complexes.

ωO−H

ω Δb

3-21G 6-31G cc-pVDZ aug-cc-pVDZ cc-pVTZ aug-cc-pVTZ

3220 3445 3665 3661 3675 3663

1747 1247 103 122 293

The aug-cc-pVTZ basis set is used for the reference frequency, ωref. Differential basis set frequency contribution, ωΔ = (|ωref2 − ωlow2|)1/2, in cm−1. a b

■


basisa

MOTIVATION: TIMESCALE CONSIDERATIONS The anticipated timescale separation of basis set components in AIMD can first be demonstrated with existing computational machinery, and these preliminary tests motivate the subsequent analysis. In particular, the contribution of basis set components to vibrational frequencies and force components may be used to estimate the relative timescales in a dynamics simulation. The contribution of basis sets to the quality of the vibrational frequencies has long been known.43,44 However, basis sets rarely change harmonic frequencies by orders of magnitude or even factors of two. They simply refine the frequencies to an accuracy that is sufficient to provide reliable validation and/or interpretation of experimental spectra assignments, such as the distinguishing of multiple possible isomers of a chemical complex.45−48 These admittedly rigorous experimental demands provide ample motivation for pursuing simulations with large basis sets, but the (relatively) smaller changes in frequencies afforded by extended basis sets makes the prospect for MTS-AIMD appear feasible. In particular, an analysis of the role of basis sets in harmonic frequencies can be devised. The harmonic, H-bonded O−H stretch frequency for the methanol dimer, (CH3OH)2 (Figure 1), is provided for a series of common basis sets, using the

where ωref signifies the reference frequency in the target basis set. The frequency of the differential contribution of the larger basis set is denoted as ωΔ and may, on the basis of the above expression, be defined as ωΔ =

ωref 2 − ω low 2

(Note that this quantity is not merely the difference in harmonic frequencies.) The values of this quantity are included in Table 1. The results indicate that the differential contribution of the atomic basis sets ranges from 48% of the aug-cc-pVTZ target frequency for STO-3G36 to 3% for cc-pVDZ. Because the necessary MD timestep is tied to the magnitude of the frequency (a common rule-of-thumb is roughly 1/20th of the timescale of the fastest frequency), the ωΔ values in this table suggest that a MTS-AIMD scheme based on basis set partitions should be viable. The role of these basis set contributions can also be analyzed directly in the time-dependent force components of standard AIMD trajectories. Shown in Figure 2 are two Cartesian force components along a 1000-step trajectory (with a 20 au timestep) of the water tetramer, (H2O)4, starting from a nonequilibrium geometry and thermal (300 K) initial momenta. The red trace in these plots corresponds to the HF/6-31++G**36,37 force component throughout the trajectory. The black trace depicts only the differential change in force, relative to HF/6-31G, at the same geometry. The upper panel, which displays a single force component for a hydrogen atom, clearly indicates a timescale separation between the two contributions. The total force oscillates on the timescale of the fast O−H stretch, and although the differential force does contain some low-amplitude, short-timescale coupling to this motion, the majority of its contribution occurs on a longer timescale. The lower panel of Figure 2, corresponding to a force component of an oxygen atom, demonstrates that frequency separation does not occur in all Cartesian components. In this case, an amplitude separation occurs instead, and this (less interesting) case is still amenable to MTS because fewer timesteps are required to integrate such a low-amplitude contribution. To determine whether this behavior is specific to small water clusters, the dynamics of one conformer of the alanine tetrapeptide, (Ala) 4 , has been simulated with HF/6-31G**, and results are depicted in Figure 3 for a Cartesian force component of the H-bonded hydrogen atom. Once again, a timescale (and amplitude) separation between the 6-31G and 6-31G** basis set contributions is observed. Therefore, these examples demonstrate that ample motivation exists to pursue an MTS-AIMD algorithm based on a basis set partitioning. This paper includes a brief discussion of the

Figure 1. Optimized structure of the methanol dimer, (CH3OH)2, using B3LYP/aug-cc-pVTZ.

B3LYP density functional,49,50 in Table 1 as a concrete reference point. The change in frequency throughout the table, compared to the aug-cc-pVTZ39,51 reference, is only 12%, but often this change in frequency would be sufficient to render comparison with experiment invalid. Furthermore, an examination of the frequency changes alone as a function of basis set does not directly suggest the inherent timescales. Under the assumption that the normal modes of a molecule do not appreciably change between two basis sets, the harmonic potential along a given mode may be written as 1 1 1 V (q) = μωref 2q2 = μω low 2q2 + μωΔ2q2 2 2 2 B


Article



Figure 2. Representative force components of hydrogen (upper panel) and oxygen (lower panel) throughout an AIMD trajectory of the water tetramer, (H2O)4, using HF/6-31++G**. The total force component is shown in red; the incremental force difference, compared to HF/6-31G, is shown in black.

standard MD integration schemes, such as Velocity Verlet. Importantly, all tabulated dynamical quantities are still computed at the target level of theory. For MTS-AIMD formulated with basis set partitioning, the total electronic energy can be written as E = Elow + ΔE, where the change in energy simply accounts for remnant relaxation effects in the target basis set, relative to a chosen low-level basis set. In this work, no approximation is made for this relaxation, and the second SCF procedure is fully converged. The corresponding molecular force used in classical MD trajectories is then −∂E/∂x = −∂Elow/∂x − ∂ΔE/∂x; inner timesteps are updated solely according to the low-level force, whereas outer timesteps are subsequently updated according to −∂ΔE/∂x. All of these quantities can be readily obtained by standard computations in quantum chemistry software packages, and the present r-RESPA MTS-AIMD scheme has been implemented in a development version of Q-Chem.52 The implementation makes a global basis set transition within the existing AIMD routines and interfaces with the r-RESPA module as needed. The lone algorithmic improvement beyond a straightforward r-RESPA implementation involves the self-consistent field (SCF) initial guess at the outer timestep. As will be discussed in the Results section, MTS-AIMD contains additional cost considerations compared to a spatially separable force field, for example. In particular, the low-level and high-level energies/ forces must be computed at the geometry of the outer timestep, which adds a potential cost overhead to MTS-AIMD. The approach used in this work is to project the converged MOs from the smaller basis set into the larger basis set’s space, by using algorithms that were originally developed for dual basis set methods.53−60 This projection is rigorous and results in properly idempotent initial density matrices. Furthermore, because the molecular structure does not change between lowand high-level computations at the outer timestep, no timereversibility concerns arise out of this “recycling” of orbitals, as

Figure 3. Representative force components of the strongly hydrogenbonded hydrogen atom throughout an AIMD trajectory of the alanine tetrapeptie, (Ala)4, using B3LYP/6-31G**. The total force component is shown in red; the incremental force difference, compared to B3LYP/ 6-31G, is shown in black. The inset depicts the equilibrium structure of the peptide, as well as the Cartesian force component examined in the figure.

implementation of this approach, as well as several tests of the accuracy and efficiency of this new algorithm.

■

METHODS Full details of multiple-timestep MD methodsin particular, the reversible reference system propagator approach (rRESPA)4can be found elsewhere and will not be repeated here. Briefly, MTS-MD involves a series of nsub MD steps performed at a chosen lower level of theory, followed by a subsequent dynamical correction at the target level of theory. The r-RESPA protocol partitions these contributions to create an approach that is both time-reversible and amenable to C


Article



3. To What Extent Does the Basis Set Dependence of Correlated EST Methods Affect Their Performance? Correlated-electron methodshere taken to include both density functional theory (DFT) and wave-function-based approachesare, in principle, well suited to a basis set partitioning. The appreciable dependence of wavefunctionbased correlation on basis set,31,35,39,51 in particular, requires special focus in the present context. With a water trimer as a test case, the performance of MP2/aug-cc-pVTZ in the context of basis set partitionings is examined, with energy fluctuations as the metric of interest. Furthermore, the recently published28 ability to perform a HF/MP2 MTS splitting suggests that examination of cumulative effects is warranted. Could a small basis set HF calculation, for example, be used as the innertimestep method for large basis set MP2 trajectories, or is this combined change a bridge too far? The water trimer analysis is, therefore, repeated by using no correlation for the inner timesteps. Initial conditions and trajectory lengths are the same as those described previously. The MP2 calculations are performed within the resolution-of-the-identity approximation.57,65−70 4. What Efficiency Gains Should Be Expected? An effective nsub-fold multiple-timestep scheme should, in principle, lead to a computational speedup by a factor of roughly nsub. This simple assertion requires that the inner timesteps are sufficiently efficient to make the outer timestep cost-dominant and that no additional overhead is encountered. The O(N4) scaling with respect to basis set size suggests that the former requirement should be satisfied for reasonably distinct basis set sizes. However, unlike the conclusions drawn for correlationbased partitioning, this relative cost savings should be independent of system size. On the basis of the allowed outer timestep lengths determined for the previous questions, these factors will be assessed in two regimes: small complexes with large basis sets [(H2O)2 with HF/aug-cc-pVTZ] and large complexes with modest basis sets [(H 2 O) 1 6 with HF/6-31G**]. Further analysis on a difficult-to-converge case (the oxidized water dimer, H4O2+) will also be performed in the context of the SCF guess recycling considerations mentioned above, using a density functional which has recently shown promise for this open-shell ion.71 Cost analyses will include average timing comparisons, as well as detailed component timings to assess where savings and overhead occur.

can be the case when using information from previous timesteps.61 By this approach, the overhead of the low-level method at this particular timestep can actually be recycled into an accelerated outer-timestep force calculation. The Introduction addressed the physical motivation for pursuing an MTS-AIMD approach based on basis set partitioning. Here, a series of testable questions is used to guide the methodology for assessing for the new method: 1. Are the Resulting Dynamics Valid AIMD Trajectories? MD trajectories for a pair of water clusters (the tetramer and 16-mer) are examined, with a focus on the magnitude of the total energy fluctuations throughout the trajectory. Although the dynamics are rigorously time-reversible, step-tostep fluctuations still occur, and the magnitude of these fluctuations is a sensitive test of the MD integration algorithm. The metric used here to analyze this quantity is the mean relative absolute energy fluctuation, averaged over the trajectory, and defined as ⟨δ⟩ =

1 NMD

NMD

∑ i=1

|Ei − E0| E0

Values of ⟨δ⟩ are computed for standard AIMD trajectories with a series of timesteps, using Velocity Verlet, and are compared to the same metric for MTS-AIMD, using a variable number of inner timesteps of length Δt = 20 au (0.484 fs). Several basis sets, of both the Pople and Dunning varieties, are examined as the low levels of theory to assess the role of the small basis set in dictating the quality of the resulting dynamics. Initial conditions for these constant-energy trajectories are composed of the equilibrium molecular structure and randomly assigned (but consistent) momenta that are selected from a 300 K distribution. Trajectories are propagated for a total of 2.42 ps, which is long enough to observe the making and breaking of hydrogen bonds but not long enough to generate global rearrangement of the cluster. Throughout all of the tests in this work, SCF energies are converged to less than 10−8 au for the maximum value of the DIIS error vector,62,63 and twoelectron integrals are screened according to the Schwarz estimate,64 using a cutoff of 10−12 au. Both of these settings are conservative enough to reliably highlight differences among basis sets and MD integration techniques. 2. Are the Resulting Dynamics Faithful to Those Computed with the Target Level of Theory? A minimum requirement for the replacement of standard AIMD with its MTS analogue is the generation of dynamical observables that are faithful representations of those generated by the target level of theory. In this series of tests, vibrational spectra, generated from the Fourier transform of the velocity autocorrelation function, are used as the observable of interest. Prior to Fourier transformation, the autocorrelation functions are dressed by a standard Tukey (tapered cosine) windowing function. A single water molecule, H2O, and protonated hydrazine, N2H5+, are used as test molecules. To avoid the statistical errors that stem from averaging over trajectories with initial conditions sampled from an ensemble, AIMD and MTSAIMD spectra are compared by using single trajectories that are initiated with identical initial (random, 300 K) conditions. To directly compare spectra, the total length of the simulation is held fixed, regardless of the size of the outer timestep. For H2O, a 16.6 ps trajectory is computed; for N2H5+, the trajectory length is 3.6 ps. Both of these trajectory lengths provide sufficient resolution to distinguish differences in methodology.

■

RESULTS The questions detailed in the preceding section were used to guide the design of a stringent series of tests of MTS-AIMD using basis set partitioning. The results of these tests are provided in this section and are organized according to the preceding questions. 1. Validity of MTS-AIMD Trajectories. Average energy fluctuations for the water tetramer are presented in Figure 4 on a logarithmic scale. The reference Velocity Verlet data, using HF/aug-cc-pVTZ, are shown as a solid black line with points. These reference data show an increase in the size of fluctuations as a function of MD timestep. The 20 au timestep chosen for MTS analysis possesses a 0.14 ppm average fluctuation in the total energy throughout the trajectory. Trajectories with timesteps beyond 60 au were unstable and led to rapid dissociation of covalent bonds. The results for MTS-AIMD, which are plotted as a function of the outer timestep (using 20 au inner timesteps), indicate stable trajectories out to at least an 8-fold outer timestep extension. The decrease in ⟨δ⟩ near 60−80 D


Article



Figure 4. Energy fluctuation data for the water tetramer, (H2O)4, using HF/aug-cc-pVTZ. Shown are energy fluctuation averages for standard Velocity Verlet (solid line with points), as well as corresponding averages for MTS-AIMD using a series of smaller basis sets. In all MTS cases, an inner timestep of Δt = 20 au was used throughout. The data are plotted as a function of the outer timestep.

Figure 5. Energy fluctuation data for the water tetramer, (H2O)4, using B3LYP/aug-cc-pVTZ. Shown are energy fluctuation averages for standard Velocity Verlet (solid line with points), as well as corresponding averages for MTS-AIMD using a series of smaller basis sets. In all MTS cases, an inner timestep of Δt = 20 au was used throughout. The data are plotted as a function of the outer timestep.

au for some basis sets is a reproducible small-system artifact that vanishes for larger test systems. As anticipated, the results are sensitive to the choice of lowlevel basis set. The STO-3G basis set, for example, produces ⟨δ⟩ values that are comparable to reference values, suggesting that little improvement is observed for the STO-3G/ aug-cc-pVTZ pairing. The fact that the longer-timestep trajectories are still stable (over this fairly limited trajectory length) with this pairing is worth noting, however, and is a property that is not fully captured by the ⟨δ⟩ metric alone. The split-valence basis sets all perform better as low-level sets, with a general trend toward lower ⟨δ⟩ with a larger low-level basis set. The cc-pVTZ basis set, for example, can extend the timestep by at least a factor of 8 with no increase in fluctuation size. More reasonable low-level basis sets, such as cc-pVDZ or even 6-31G, also perform well enough to be considered viable MD algorithms. Speedups of factors of 2.5−5 should be anticipated for these pairings when equivalent values of ⟨δ⟩ are considered. The HF and DFT methods depend similarly on the quality of the atomic-orbital basis set (aside from occasional issues arising from the convergence of some exchange-correlation functionals72). Confirmation that similar MTS behavior is observed in DFT is provided in Figure 5, with the same test parameters as in Figure 4. For the basis sets tested, B3LYP MTS-AIMD performance is analogous to results observed for HF. The data in Figure 6 indicate that these properties are also transferable to the regime of larger complexes and more modest basis sets. With the 6-31G** basis set as the target basis for an MD trajectory of (H2O)16, three small low-level basis are examined for MTS-AIMD. All three lead to valid dynamics, although results are, as expected, best for the 6-31G low-level basis. Encouragingly, the results are relatively improved in the 16-water case, compared to the tetramer case examined above. Even 8-fold increases in the outer timestep lead to only modest increases in the energy fluctuation values. To assess whether these results remain valid at higher kinetic energies, such as those encountered in high-temperature and quasiclassical simulations, a similar analysis (not shown in the

Figure 6. Energy fluctuation data for a water cluster, (H2O)16, using HF/6-31G**. Shown are energy fluctuation averages for standard Velocity Verlet, as well as corresponding averages for MTS-AIMD using a series of smaller basis sets. In all MTS cases, an inner timestep of Δt = 20 au was used throughout. The data are plotted as a function of the outer timestep.

plots) was repeated for (H2O)8 at 1000 K using the 6-31G/631G** pairing. For this latter test, the standard MD integration step must naturally decrease (Δt = 10 au) to exhibit similar energy fluctuations to the 300 K tests. However, the MTSAIMD algorithm still performed comparably in this regime and exhibited stable trajectories out to at least 14 outer timesteps. 2. Faithfulness of MTS-AIMD to AIMD Observables. Using vibrational spectra as a sensitive test of the inherent dynamics, the fidelity of MTS-AIMD observables to their AIMD analogues is examined in this section. The O−H stretch portion of the vibrational spectrum of a single water molecule, using HF/6-31G**, is presented in Figure 7. Although spectra are obtainable for this small molecule with nearly any electronic structure method, the lack of low-frequency modes in this system actually makes it a somewhat demanding test system for MTS. The figure shows the results of using STO-3G (left panel) and 6-31G (right panel) as the low-level basis sets. The E


Article



Figure 7. Vibrational spectral densities (O−H stretch region) of a single water molecule using the HF/6-31G** target method. Spectra are computed from the Fourier transform of the velocity autocorrelation function. The left panel depicts results using STO-3G as the low-level basis set; the right panel depicts analogous results using 6-31G. The upper two traces show spectra using standard AIMD with the basis set indicated. All subsequent data, indicated by MTS(n), are computed with MTS-AIMD. Densities have been normalized to the density of the asymmetric stretch.

upper two traces in each panel show standard AIMD data with the indicated basis set. The remaining traces show corresponding MTS-AIMD spectra with up to nine inner timesteps. The MTS(n) nomenclature indicates the number of low-level gradients computed per full timestep, with the convention that MTS(1) is merely a repartitioned equivalent of standard Velocity Verlet in the target basis set. The results in Figure 7 indicate that minimal basis sets, such as STO-3G, are not particularly effective for MTS-AIMD simulations when 6-31G** accuracy is desired. Changes, relative to the reference spectrum, begin to appear by MTS(3) and become catastrophic by MTS(9). In contrast,

the 6-31G basis set serves as a sufficient partner to 6-31G**. Apart from very subtle shifts in relative intensities (spectral densities in these plots), the spectra are identical out to the longest timestep examined here [MTS(9)]. Because the relative size disparity between 6-31G and 6-31G** is still rather large, this low-level basis set is hereafter recommended as an MTS partner to basis sets of the polarized double-ζ variety. The high-frequency (upper panel) and mid-frequency (lower panel) vibrational spectra of the hydrazinium cation are presented in Figure 8. The upper and lower traces in each panel depict reference AIMD data using the 6-311G** and 6-31G basis sets, respectively. These two reference spectra are F


Article



Figure 9. Energy fluctuation data for the water trimer, (H2O)3, using RI-MP2/aug-cc-pVTZ. Shown are energy fluctuation averages for standard Velocity Verlet (solid line with points), as well as corresponding averages for MTS-AIMD using a series of smaller basis sets (blue lines). Also shown are MTS energy fluctuations when HF is used as the low-level method (red lines). In all MTS cases, an inner timestep of Δt = 20 au was used throughout. The data are plotted as a function of the outer timestep.

The most straightforward application of the present basis set partitioning scheme uses RI-MP2 with smaller basis sets as the low-level method for MTS; these results are shown as blue curves in Figure 9. The trajectories are stable out to (at least) eight inner timesteps, and errors remain comparable to those of the original Velocity Verlet trajectories. Little advantage is found for including the diffuse (“aug”) functions in the lowlevel basis set because cc-pVDZ and aug-cc-pVDZ results are nearly identical. The cc-pVTZ basis performs the best but also would possess the largest low-level overhead. Once again, cc-pVDZ appears to be a sufficiently accurate and efficient compromise. Interestingly, the basis set partitioning scheme outperforms the previously published correlation partitioning scheme for this test system. With HF/aug-cc-pVTZ as the low-level method (i.e., only changing the level of correlation in the same basis set), errors are improved, relative to Velocity Verlet, but these results are notably inferior to all of the truncated-basis results with RI-MP2. Significantly, all of the HF results are actually improved by moving to smaller basis sets. When the correlation-based and basis set-based approaches to MTS are combined, the smallest basis sets at the HF level produce the lowest average errors. The HF/6-31G low-level method performs the best among the mixed-method approaches. Whether this result stems from a fortuitous cancellation of errors or from a more systematic trend, such as overbinding from basis set superposition error, which is consistent with correlation effects, remains to be investigated further. One of the cost considerations examined in ref 28 was the fact that the O(N4) [HF] vs O(N5) [MP2] system-size scaling guarantees that a correlation-based MTS-AIMD method will remain efficient, even for very large systems. However, for large basis sets, this regime may not yet be practical. Therefore, the results of this analysis are particularly encouraging, because it appears that the cost overhead of the inner timesteps may be significantly reduced with use of a severely truncated basis set. 4. Computational Efficiency. The underlying aim of any MTS recipe is to generate faithful molecular dynamics

Figure 8. Vibrational spectral densities of protonated hydrazine, N2H5+, using the HF/6-311G** target method in the high-frequency (upper panel) and mid-frequency (lower panel) active regions of its spectrum. Spectra are computed from the Fourier transform of the velocity autocorrelation function. The upper trace in each panel depicts results from the target 6-311G** basis set; the lower trace depicts results from the low-level 6-31G basis set alone. The middle trace in each panel shows MTS-AIMD spectra using seven additional (eight total) low-level timesteps, denoted as MTS(8).

appreciably different, particularly in the highest-frequency N−H stretch, and no uniform scaling factor43,44 could recover this discrepancy. The middle trace shows the MTS(8) spectrum using the 6-31G/6-311G** pairing. The spectra are nearly indistinguishable on the scale of the plots. In fact, to within the somewhat limited resolution (3 cm−1) of the computed spectra, the highest-frequency N−H stretch occurs at exactly the same frequency in the two methods. Therefore, these examples indicate that MTS-AIMD trajectories faithfully reproduce large basis set dynamical observables, provided that a low-level basis set of at least 6-31G quality is utilized. 3. MTS Effects with Electron Correlation. Results for MTS-AIMD with DFT have already been examined above and were shown to parallel the behavior for Hartree−Fock, as expected. The basis set behavior for wavefunction-based, correlated-electron methods is examined in this section. The average energy fluctuations for a trajectory of the water trimer, using RI-MP2/aug-cc-pVTZ as the target methodology, are depicted in Figure 9. G


Article

The Journal of Physical Chemistry A Table 2. Computational Cost Analysis for (H2O)16 Using HF/6-31G** ⟨timings⟩a basis set

energy

force

total

estimated speedup, Tb

actual speedup

⟨SCF cycles⟩

6-31G MTS(8) MTS(8)c 6-31G**

56.8 55.5/199.8 55.5/181.9 196.6

13.9 13.7/70.4 13.7/70.3 69.9

71.0 104.1 101.8 266.9

5.0 5.0

2.56 2.62

12.2 12.0/12.0 12.0/10.8 12.0

a Serial CPU sec per 1 fs simulation. For MTS entries, low/target values are shown. For MTS energy and force timings, artificial adjustment for a 20 au timestep has been used to enable direct comparison to the target basis set’s value. For total step timings, the full 160 au timestep has been used. b Based on eq 1. cData using projected low-level MO coefficients as target-level SCF guess.

Table 3. Computational Cost Analysis for (H2O)2 Using HF/aug-cc-pVTZ


⟨timings⟩a basis set

energy

force

total step


actual speedup

⟨SCF cycles⟩

cc-pVDZ MTS(8) MTS(8)c aug-cc-pVTZ

2.6 2.7/82.4 2.6/76.5 81.3

0.5 0.5/40.1 0.5/40.2 39.8

3.1 18.5 17.8 121.2

7.7 7.7

6.6 6.8

13.0 13.0/12.0 13.0/10.8 12.0

a Serial CPU sec per 1 fs simulation. For MTS entries, low/target values are shown. For MTS energy and force timings, artificial adjustment for a 20 au timestep has been used to enable direct comparison to the target basis set’s value. For total step timings, the full 160 au timestep has been used. b Based on eq 1. cData using projected low-level MO coefficients as target-level SCF guess.

properties at decreased computational cost. The fidelity of trajectories has been demonstrated above, but the extent of expected computational speedup comprises the final analysis of this work. For molecular potentials based on piecewise-additive force fields, the use of nsub inner timesteps should immediately lead to a factor of ≲nsub reduction in the computational cost of the simulation (neglecting any potential increase in ⟨δ⟩ values). A cost analysis is slightly more nuanced in an AIMD context for three reasons, which are analyzed independently in this section. First, meaningful MTS acceleration requires that the lowlevel method be sufficiently expedient to compute. For a basis set partitioning in MTS-AIMD, this factor is controlled by the O(N4) scaling of SCF energy and force calculations with respect to basis set size N. MP2 theory scales as O(N3)−O(N4) with respect to the same quantity, depending on the system size, but the following analysis will be restricted to SCF-based dynamics. Therefore, the ratio of low- and target-level basis set sizes is the critical metric for this consideration. The previous sections have demonstrated effective pairings for MTS-AIMD, and basis set ratios for these systems provide concrete reference points. In water clusters (or any system with a similar ratio of heavy/light atoms), for example, the 6-31G/6-31G** pairing exhibits a target-to-low basis set size ratio of 1.92. Given the fourth-order scaling of HF and DFT, the low-level method should be 1.924 = 13.7 times faster than the target-level method. More ambitious partitionings in the large basis set regime further enhance this disparity. The cc-pVDZ/ aug-cc-pVTZ pairing exhibits a size ratio of 3.75, which leads to a 198-fold speedup in the low-level contribution. These formal analyses indicate that both pairings are sufficiently disparate to relegate the cost of the low-level steps to a minor overhead, rather than a computational bottleneck. The second cost consideration involves the fact that the lowlevel method’s energy and force must be computed at the same molecular configuration during the last inner timestep and lone outer timestep for each full MD step. In piecewise-additive analytic potentials, this overhead is not present because all pieces are physically distinct. In a basis set-based MTS-AIMD prescription, the fact that two molecular forces must be

computed at the outer timestep is considered as an overhead. The relative cost, T, of a target AIMD timestep, compared to an equivalent MTS-AIMD timestep, can be shown to be T=

t target tMTS

=

4 nsubNtarget 4 4 + Ntarget nsubNsub

(1)

In the Ntarget ≫ Nsub limit, this ratio reduces to nsub, as anticipated. In a more practical regime, however, this cost ratio is diminished. The extent to which it is diminished is examined in three reference complexes, the data for which are shown in Tables 2−4. The data in these tables are presented in the form of serial CPU time (second) per simulated femtosecond. For the energy and force timings in the MTS entries, a timestep of 20 au is assumed for both inner and outer steps, simply for the sake of direct comparison to the AIMD timings. In contrast, the “Total” timings properly account for the actual timestep used in the timing simulations. Table 2 shows trajectory timing data for (H2O)16 with HF/6-31G** as the target level of theory. With the 6-31G lowlevel basis, the anticipated speedup for MTS(8), when low-level overhead is included, is a factor of 5. The observed speedup in actual timings is about 2.6. This reduced speedup is mainly due to the fact that the low-level SCF/force calculations are only 4− 5 times faster than those with 6-31G**, rather than the relative cost predicted by the O(N4) scaling. This discrepancy is due to a combination of factors, including (a) the fact that the quartic scaling analysis neglects the different cost of integrals involving higher angular momenta and (b) the use of incremental Fock matrix construction64 during the SCF procedure, which affects each basis set differently and results in up to a factor of 5 change in cost of the Fock build by the last SCF steps. Nonetheless, appreciable speedups in this basis set size regime are still observed. In fact, the resulting timestep is not much more costly than the original low-level dynamics method alone; 83% of the cost difference has been recovered by the MTS approach. Timings are more consistent with T estimates in the large basis set size regime. Analogous timing information for the H


Article

The Journal of Physical Chemistry A Table 4. Computational Cost Analysis for (H2O)2+ Using MPW1K/aug-cc-pVTZ ⟨timings⟩a basis set

energy

force

total step


actual speedup

⟨SCF cycles⟩

cc-pVDZ MTS(8)c aug-cc-pVTZ

15.3 15.2/193.6 294.9

1.1 1.1/47.9 47.5

16.4 46.7 342.7

7.7

7.3

24.1 24.5/17.4 26.1

a


Serial CPU sec per 1 fs simulation. For MTS entries, low/target values are shown. For MTS energy and force timings, artificial adjustment for a 20 au timestep has been used to enable direct comparison to the target basis set’s value. For total step timings, the full 160 au timestep has been used. b Based on eq 1. cData using projected low-level MO coefficients as target-level SCF guess.

advantages, beyond the obvious trait of increased accuracy. First, basis set superposition error (BSSE) is reduced. Eliminating this artifact from MD trajectories is challenging, especially for intramolecular BSSE,78−81 and enabling large basis set MD is the most straightforward solution. Second, the MTS method presented in this work retains much of the accurate basis set’s influence on the high-frequency data contained in MD simulations. In contrast to restrained/ constrained dynamics approaches,82,83 which “freeze out” the high-frequency degrees of freedom, the MTS-AIMD method in this work still allows these degrees of freedom to oscillate and, therefore, also includes the frequency shifts associated with molecular interactions involving these degrees of freedom.84−87 As shown above, vibrational spectra are nearly perfect representations of their target-basis analogues. Based on the preceding tests and timings, recommended basis set pairings for HF and DFT include 6-31G for simulations using 6-31G** or 6-311G** as the target basis set, as well as cc-pVDZ as the subset for high-accuracy aug-cc-pVTZ simulations. The recommendation for MP2 simulations with aug-cc-pVTZ is the same. However, the lone test in this work suggests that HF/6-31G may also be an even more efficient MTS partner for MP2/aug-cc-pVTZ without sacrificing much accuracy, although a more detailed examination of this combined approach is warranted. General recommendations for “which functions to remove” in a given basis set are as-yet difficult to assign. In this sense, plane-wave basis sets have a distinct advantage. In the context of extendedLagrangian/Car−Parrinello MD, different Fourier components can be used to efficiently propagate electronic information on differing timescales.23 Such a strategy has solid physical foundation. The analogue of this intuition in the context of atom-centered Gaussian basis sets, particularly when functions of differing angular momenta are involved, is less clear. The recommended pairings above suggest that polarization functions (higher angular momentum) and diffuse functions (larger spatial extent and the atom-centered equivalent of lower-frequency Fourier components) define the most effective timescale splittings. These qualitative splittings also hold true across all systems examined and even across other basis sets. Therefore, these recommendations appear to be robust but remain empirical at present. Core functions have also not been explicitly examined in the present algorithm, although a slow (or even static) evolution of these contributions is expected, on the basis of the widespread success of core potentials88−97 and similarly motivated extended-Lagrangian approaches that treat these contributions separately.98 Such considerations will be the focus of future analyses. Past work in the realm of mixing basis sets53−60 (at the electronic structure theory level) has shown that perturbative corrections to account for basis set relaxation are an effective compromise of accuracy and efficiency. These “dual-basis” or

water dimer and HF/aug-cc-pVTZ is shown in Table 3. The low/target disparity is still below the O(N4) estimate, but the gap is still large enough to lead to more significant savings. A 6.6-fold speedup is observed for eight inner timesteps, suggesting that the low-level overhead, including the low-level energy/force at the outer timestep’s geometry, is (relatively) reduced for this pairing. The final factor which is unique to MTS-AIMD is the ability to circumvent some of this last step’s overhead. The low-level energy computed at this molecular configuration can actually be used to accelerate the outer timestep’s energy calculation. The aforementioned projection of the molecular orbital coefficients from the low-level basis set into the target basis allows for the possibility of reducing the number of SCF cycles in this target basis. Data with and without this modification are shown in Tables 2 and 3. For (H2O)16, the average number of SCF cycles with 6-31G** is reduced from 12.0 to 10.8, which accounts for a modest reduction in total step time. For (H2O)2 and aug-cc-pVTZ, a similar reduction in cycles is observed. To highlight the potential benefit of this approach, however, timings have also been computed for a trajectory of the oxidized water dimer, (H2O)2+. In past studies of these ionized water clusters,71,73,74 the MPW1K density functional75,76 was shown to accurately reproduce many energetic and structural properties. However, a practical observation during these studies was slow and unreliable SCF convergence. (For these studies, a core Hamiltonian guess was found to be most reliable, in combination with a mixed algorithm involving DIIS for early cycles and Geometric Direct Minimization77 (GDM) in later cycles. The choice of this guess is based on empirical performance results, although the fact that this system is both small and charged likely contributes to the success of this particular guess.) As shown in Table 4, the cc-pVDZ and aug-cc-pVTZ basis sets require 24.1 and 26.1 SCF cycles, respectively, for convergence. The low-level “recycled” guess appreciably accelerates the trajectories for this system. The 26.1 SCF cycles at the target level are reduced to 17.4, on average, for the MTS-AIMD trajectory. Because of this additional savings, the overall timestep average cost (a speedup of 7.3) is even closer to the estimate in eq 1. Therefore, although the additional speedups due to this SCF modification are small for easy-to-converge electronic structures, the benefits are significant for more challenging cases and may reasonably be included for all MTS-AIMD calculations.

■

DISCUSSION AND CONCLUSIONS The results of the previous section confirm the most important aspects to be addressed in a MTS-AIMD scheme based on basis set partitioning: The resulting trajectories are valid MD trajectories that faithfully reproduce target-basis observables, and they do so at appreciably reduced cost. The ability to generate trajectories with quality basis sets possesses additional I


Article


The Journal of Physical Chemistry A “jumping”99,100 approaches have been extended to the realm of chemical dynamics60 and have been shown to be an efficient competitor with extended-Lagrangian-based MD. The present MTS scheme exhibits some advantages over these approaches, including the property that any combination of basis sets can, in principle, be mixed instead of requiring rigorous subsets, as well as the fact that gradients in the MTS-AIMD scheme are no more complicated than in standard AIMD. Of course, the two approaches could feasibly be combined, by which the outertimestep forces are computed with dual-basis methods. Finally, the prospects for very long outer timesteps pose the most significant challenge for continued extension of MTSAIMD. Resonance artifacts limit the outer timestep to roughly half the vibrational period of the fastest motion.101−103 Indeed, the upturn in ⟨δ⟩ values for the longest outer timesteps in the plots of this work, even with relatively large low-level basis sets, suggests that the onset of resonance provides a fundamental limit to the computational acceleration afforded by this methodology. However, progress has been made in developing methods to address this limit, including methods involving weakly perturbing thermostats,104 stochastic isokinetic approaches,105 and dynamical partitioning.106 Whether these approaches further extend the maximum outer timestep in the present approach, or whether the inherent timescale limit has already been reached, remains to be investigated. If these resonance limits can be safely surmounted, multitiered approaches to basis set partitioning would also be enabled.

■

Simulations with Long-Range Interactions. Mol. Simul. 1991, 6, 121−142. (8) Watanabe, M.; Karplus, M. Dynamics of Molecules with Internal Degrees of Freedom by Multiple Time-Step Methods. J. Chem. Phys. 1993, 99, 8063−8074. (9) Procacci, P.; Berne, B. J. Computer Simulation of Solid C60 Using Multiple Time-Step Algorithms. J. Chem. Phys. 1994, 101, 2421−2431. (10) Humphreys, D. D.; Friesner, R. A.; Berne, B. J. Simulated Annealing of a Protein in a Continuum Solvent by Multiple-Time-Step Molecular Dynamics. J. Phys. Chem. 1995, 99, 10674−10685. (11) Watanabe, M.; Karplus, M. Simulations of Macromolecules by Multiple Time-Step Methods. J. Phys. Chem. 1995, 99, 5680−5697. (12) Zhou, R.; Berne, B. J. A New Molecular Dynamics Method Combining the Reference System Propagator Algorithm with a Fast Multipole Method for Simulating Proteins and Other Complex Systems. J. Chem. Phys. 1995, 103, 9444−9459. (13) Procacci, P.; Darden, T.; Marchi, M. A Very Fast Molecular Dynamics Method to Simulate Biomolecular Systems with Realistic Electrostatic Interactions. J. Phys. Chem. 1996, 100, 10464−10468. (14) Procacci, P.; Marchi, M. Taming the Ewald Sum in Molecular Dynamics Simulations of Solvated Proteins Via a Multiple Time Step Algorithm. J. Chem. Phys. 1996, 104, 3003−3012. (15) Stuart, S. J.; Zhou, R.; Berne, B. J. Molecular Dynamics with Multiple Time Scales: The Selection of Efficient Reference System Propagators. J. Chem. Phys. 1996, 105, 1426−1436. (16) Schlick, T.; Barth, E.; Mandziuk, M. Biomolecular Dynamics at Long Timesteps:Bridging the Timescale Gap between Simulation and Experimentation. Annu. Rev. Biophys. Biomol. Struct. 1997, 26, 181− 222. (17) Grubmüller, H.; Tavan, P. Multiple Time Step Algorithms for Molecular Dynamics Simulations of Proteins: How Good Are They? J. Comput. Chem. 1998, 19, 1534−1552. (18) Izaguirre, J. A.; Reich, S.; Skeel, R. D. Longer Time Steps for Molecular Dynamics. J. Chem. Phys. 1999, 110, 9853−9864. (19) Sagui, C.; Darden, T. A. Molecular Dynamics Simulations of Biomolecules: Long-Range Electrostatic Effects. Annu. Rev. Biophys. Biomol. Struct. 1999, 28, 155−179. (20) Zhou, R.; Harder, E.; Xu, H.; Berne, B. J. Efficient Multiple Time Step Method for Use with Ewald and Particle Mesh Ewald for Large Biomolecular Systems. J. Chem. Phys. 2001, 115, 2348−2358. (21) Morrone, J. A.; Zhou, R.; Berne, B. J. Molecular Dynamics with Multiple Time Scales: How to Avoid Pitfalls. J. Chem. Theory Comput. 2010, 6, 1798−1804. (22) Luehr, N.; Markland, T. E.; Martínez, T. J. Multiple Time Step Integrators in Ab Initio Molecular Dynamics. J. Chem. Phys. 2014, 140, 084116. (23) Tuckerman, M. E.; Parrinello, M. Integrating the Car–Parrinello Equations. I. Basic Integration Techniques. J. Chem. Phys. 1994, 101, 1302−1315. (24) Gibson, D. A.; Carter, E. A. Time-Reversible Multiple Time Scale Ab Initio Molecular Dynamics. J. Phys. Chem. 1993, 97, 13429− 13434. (25) Anglada, E.; Junquera, J.; Soler, J. M. Efficient Mixed-Force First-Principles Molecular Dynamics. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 68, 055701. (26) Guidon, M.; Schiffmann, F.; Hutter, J.; VandeVondele, J. Ab Initio Molecular Dynamics Using Hybrid Density Functionals. J. Chem. Phys. 2008, 128, 214104. (27) Fatehi, S.; Steele, R. P. Multiple-Time Step Ab Initio Molecular Dynamics Based on Two-Electron Integral Screening. J. Chem. Theory Comput. 2015, 11, 884−898. (28) Steele, R. P. Communication: Multiple-Timestep Ab Initio Molecular Dynamics with Electron Correlation. J. Chem. Phys. 2013, 139, 011102. (29) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry; McGrawHill, Inc.: New York, 1989. (30) Helgaker, T.; Jørgensen, P.; Olsen, J. Molecular ElectronicStructure Theory; John Wiley & Sons Ltd.: England, 2000.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This material is based upon work supported by the National Science Foundation CAREER under CHE-1452596. The support and resources from the Center for High-Performance Computing at the University of Utah are gratefully acknowledged. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575.

■

REFERENCES

(1) Tuckerman, M. E.; Martyna, G. J.; Berne, B. J. Molecular Dynamics Algorithm for Condensed Systems with Multiple Time Scales. J. Chem. Phys. 1990, 93, 1287−1291. (2) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J. Molecular Dynamics Algorithm for Multiple Time Scales: Systems with Long Range Forces. J. Chem. Phys. 1991, 94, 6811−6815. (3) Tuckerman, M. E.; Berne, B. J.; Rossi, A. Molecular Dynamics Algorithm for Multiple Time Scales: Systems with Disparate Masses. J. Chem. Phys. 1991, 94, 1465−1469. (4) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible Multiple Time Scale Molecular Dynamics. J. Chem. Phys. 1992, 97, 1990−2001. (5) Streett, W. B.; Tildesley, D. J.; Saville, G. Multiple Time-Step Methods in Molecular Dynamics. Mol. Phys. 1978, 35, 639−648. (6) Swindoll, R. D.; Haile, J. M. A Multiple Time-Step Method for Molecular Dynamics Simulations of Fluids of Chain Molecules. J. Comput. Phys. 1984, 53, 289−298. (7) Grubmüller, H.; Heller, H.; Windemuth, A.; Schulten, K. Generalized Verlet Algorithm for Efficient Molecular Dynamics J


Article


The Journal of Physical Chemistry A (31) Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-Set Convergence of Correlated Calculations on Water. J. Chem. Phys. 1997, 106, 9639−9646. (32) Helgaker, T.; Gauss, J.; Jørgensen, P.; Olsen, J. The Prediction of Molecular Equilibrium Structures by the Standard Electronic Wave Functions. J. Chem. Phys. 1997, 106, 6430−6440. (33) Bak, K. L.; Jørgensen, P.; Olsen, J.; Helgaker, T.; Klopper, W. Accuracy of Atomization Energies and Reaction Enthalpies in Standard and Extrapolated Electronic Wave Function/Basis Set Calculations. J. Chem. Phys. 2000, 112, 9229−9242. (34) Feller, D. Application of Systematic Sequences of Wave Functions to the Water Dimer. J. Chem. Phys. 1992, 96, 6104−6114. (35) Feller, D. The Use of Systematic Sequences of Wave Functions for Estimating the Complete Basis Set, Full Configuration Interaction Limit in Water. J. Chem. Phys. 1993, 98, 7059−7071. (36) Hehre, W. J.; Stewart, R. F.; Pople, J. A. Self-Consistent Molecular-Orbital Methods. I. Use of Gaussian Expansions of SlaterType Atomic Orbitals. J. Chem. Phys. 1969, 51, 2657−2664. (37) Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theoret. Chim. Acta 1973, 28, 213−222. (38) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. SelfConsistent Molecular Orbital Methods. XX. A Basis Set for Correlated Wave Functions. J. Chem. Phys. 1980, 72, 650−654. (39) Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (40) White, C. A.; Johnson, B. G.; Gill, P. M. W.; Head-Gordon, M. The Continuous Fast Multipole Method. Chem. Phys. Lett. 1994, 230, 8−16. (41) Schwegler, E.; Challacombe, M. Linear Scaling Computation of the Hartree–Fock Exchange Matrix. J. Chem. Phys. 1996, 105, 2726− 2734. (42) Strain, M. C.; Scuseria, G. E.; Frisch, M. J. Achieving Linear Scaling for the Electronic Quantum Coulomb Problem. Science 1996, 271, 51−53. (43) Scott, A. P.; Radom, L. Harmonic Vibrational Frequencies: An Evaluation of Hartree−Fock, Møller−Plesset, Quadratic Configuration Interaction, Density Functional Theory, and Semiempirical Scale Factors. J. Phys. Chem. 1996, 100, 16502−16513. (44) Sinha, P.; Boesch, S. E.; Gu, C.; Wheeler, R. A.; Wilson, A. K. Harmonic Vibrational Frequencies: Scaling Factors for HF, B3LYP, and MP2 Methods in Combination with Correlation Consistent Basis Sets. J. Phys. Chem. A 2004, 108, 9213−9217. (45) Relph, R. A.; Guasco, T. L.; Elliott, B. M.; Kamrath, M. Z.; McCoy, A. B.; Steele, R. P.; Schofield, D. P.; Jordan, K. D.; Viggiano, A. A.; Ferguson, E. E.; et al. How the Shape of an H-Bonded Network Controls Proton-Coupled Water Activation in HONO Formation. Science 2010, 327, 308−312. (46) Marsh, B. M.; Duffy, E. M.; Soukup, M. T.; Zhou, J.; Garand, E. Intramolecular Hydrogen Bonding Motifs in Deprotonated Glycine Peptides by Cryogenic Ion Infrared Spectroscopy. J. Phys. Chem. A 2014, 118, 3906−3912. (47) Marsh, B. M.; Zhou, J.; Garand, E. Vibrational Spectroscopy of Small Hydrated CuOH+ Clusters. J. Phys. Chem. A 2014, 118, 2063− 2071. (48) Yang, B.; Wu, R. R.; Berden, G.; Oomens, J.; Rodgers, M. T. Infrared Multiple Photon Dissociation Action Spectroscopy of ProtonBound Dimers of Cytosine and Modified Cytosines: Effects of Modifications on Gas-Phase Conformations. J. Phys. Chem. B 2013, 117, 14191−14201. (49) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (50) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627.

(51) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron Affinities of the First-Row Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 1992, 96, 6796−6806. (52) Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T. B.; Slipchenko, L. V.; Levchenko, S. V.; O’Neill, D. P.; et al. Advances in Methods and Algorithms in a Modern Quantum Chemistry Program Package. Phys. Chem. Chem. Phys. 2006, 8, 3172−3191. (53) Liang; Head-Gordon, M. Approaching the Basis Set Limit in Density Functional Theory Calculations Using Dual Basis Sets without Diagonalization†. J. Phys. Chem. A 2004, 108, 3206−3210. (54) Wolinski, K.; Pulay, P. Second-Order Møller−Plesset Calculations with Dual Basis Sets. J. Chem. Phys. 2003, 118, 9497− 9503. (55) Steele, R. P.; DiStasio, J. R. A.; Shao, Y.; Kong, J.; Head-Gordon, M. Dual-Basis Second-Order Møller-Plesset Perturbation Theory: A Reduced-Cost Reference for Correlation Calculations. J. Chem. Phys. 2006, 125, 074108−11. (56) Steele, R. P.; Shao, Y.; DiStasio, R. A.; Head-Gordon, M. DualBasis Analytic Gradients. 1. Self-Consistent Field Theory. J. Phys. Chem. A 2006, 110, 13915−13922. (57) Distasio, R. A.; Steele, R. P.; Rhee, Y. M.; Shao, Y.; HeadGordon, M. An Improved Algorithm for Analytical Gradient Evaluation in Resolution-of-the-Identity Second-Order Møller-Plesset Perturbation Theory: Application to Alanine Tetrapeptide Conformational Analysis. J. Comput. Chem. 2007, 28, 839−856. (58) Steele, R. P.; Head-Gordon, M. Dual-Basis Self-Consistent Field Methods: 6-31g* Calculations with a Minimal 6-4G Primary Basis. Mol. Phys. 2007, 105, 2455−2473. (59) Steele, R. P.; DiStasio, R. A.; Head-Gordon, M. Non-Covalent Interactions with Dual-Basis Methods: Pairings for Augmented Basis Sets. J. Chem. Theory Comput. 2009, 5, 1560−1572. (60) Steele, R. P.; Head-Gordon, M.; Tully, J. C. Ab Initio Molecular Dynamics with Dual Basis Set Methods. J. Phys. Chem. A 2010, 114, 11853−11860. (61) Herbert, J. M.; Head-Gordon, M. Accelerated, EnergyConserving Born-Oppenheimer Molecular Dynamics Via Fock Matrix Extrapolation. Phys. Chem. Chem. Phys. 2005, 7, 3269−3275. (62) Pulay, P. Convergence Acceleration of Iterative Sequences. The Case of SCF Iteration. Chem. Phys. Lett. 1980, 73, 393−398. (63) Pulay, P. Improved SCF Convergence Acceleration. J. Comput. Chem. 1982, 3, 556−560. (64) Häser, M.; Ahlrichs, R. Improvements on the Direct SCF Method. J. Comput. Chem. 1989, 10, 104−111. (65) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Use of Approximate Integrals in Ab Initio Theory. An Application in MP2 Energy Calculations. Chem. Phys. Lett. 1993, 208, 359−363. (66) Vahtras, O.; Almlö f , J.; Feyereisen, M. W. Integral Approximations for LCAO-SCF Calculations. Chem. Phys. Lett. 1993, 213, 514−518. (67) Eichkorn, K.; Treutler, O.; Ö hm, H.; Häser, M.; Ahlrichs, R. Auxiliary Basis Sets to Approximate Coulomb Potentials. Chem. Phys. Lett. 1995, 240, 283−290. (68) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Auxiliary Basis Sets for Main Row Atoms and Transition Metals and Their Use to Approximate Coulomb Potentials. Theor. Chem. Acc. 1997, 97, 119−124. (69) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: Optimized Auxiliary Basis Sets and Demonstration of Efficiency. Chem. Phys. Lett. 1998, 294, 143−152. (70) Jung, Y.; Sodt, A.; Gill, P. M. W.; Head-Gordon, M. Auxiliary Basis Expansions for Large-Scale Electronic Structure Calculations. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6692−6697. (71) Herr, J. D.; Talbot, J.; Steele, R. P. Structural Progression in Clusters of Ionized Water, (H2O)n=1−5+. J. Phys. Chem. A 2015, 119, 752−766. (72) Mardirossian, N.; Head-Gordon, M. Characterizing and Understanding the Remarkably Slow Basis Set Convergence of Several K


Article


The Journal of Physical Chemistry A Minnesota Density Functionals for Intermolecular Interaction Energies. J. Chem. Theory Comput. 2013, 9, 4453−4461. (73) Mizuse, K.; Kuo, J.-L.; Fujii, A. Structural Trends of Ionized Water Networks: Infrared Spectroscopy of Watercluster Radical Cations (H2O)n+ (n = 3−11). Chem. Sci. 2011, 2, 868−876. (74) Mizuse, K.; Fujii, A. Characterization of a Solvent-Separated Ion-Radical Pair in Cationized Water Networks: Infrared Photodissociation and Ar-Attachment Experiments for Water Cluster Radical Cations (H2O)n+ (n = 3−8). J. Phys. Chem. A 2013, 117, 929−938. (75) Lynch, B. J.; Fast, P. L.; Harris, M.; Truhlar, D. G. Adiabatic Connection for Kinetics. J. Phys. Chem. A 2000, 104, 4811−4815. (76) Adamo, C.; Barone, V. Exchange Functionals with Improved Long-Range Behavior and Adiabatic Connection Methods without Adjustable Parameters: The MPW and MPW1PW Models. J. Chem. Phys. 1998, 108, 664−675. (77) Van Voorhis, T.; Head-Gordon, M. A Geometric Approach to Direct Minimization. Mol. Phys. 2002, 100, 1713−1721. (78) Jensen, F. The Magnitude of Intramolecular Basis Set Superposition Error. Chem. Phys. Lett. 1996, 261, 633−636. (79) Brandenburg, J. G.; Alessio, M.; Civalleri, B.; Peintinger, M. F.; Bredow, T.; Grimme, S. Geometrical Correction for the Inter- and Intramolecular Basis Set Superposition Error in Periodic Density Functional Theory Calculations. J. Phys. Chem. A 2013, 117, 9282− 9292. (80) Asturiol, D.; Duran, M.; Salvador, P. Intramolecular Basis Set Superposition Error Effects on the Planarity of Benzene and Other Aromatic Molecules: A Solution to the Problem. J. Chem. Phys. 2008, 128, 144108. (81) Valdés, H.; Klusák, V.; Pitoňaḱ , M.; Exner, O.; Starý, I.; Hobza, P.; Rulíšek, L. Evaluation of the Intramolecular Basis Set Superposition Error in the Calculations of Larger Molecules: [N]Helicenes and PheGly-Phe Tripeptide. J. Comput. Chem. 2008, 29, 861−870. (82) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (83) Andersen, H. C. Rattle: A “Velocity” Version of the Shake Algorithm for Molecular Dynamics Calculations. J. Comput. Phys. 1983, 52, 24−34. (84) Tainter, C. J.; Ni, Y.; Shi, L.; Skinner, J. L. Hydrogen Bonding and Oh-Stretch Spectroscopy in Water: Hexamer (Cage), Liquid Surface, Liquid, and Ice. J. Phys. Chem. Lett. 2013, 4, 12−17. (85) Tainter, C. J.; Shi, L.; Skinner, J. L. Structure and Oh-Stretch Spectroscopy of Low- and High-Density Amorphous Ices. J. Chem. Phys. 2014, 140, 134503. (86) Baron, R.; Setny, P.; Paesani, F. Water Structure, Dynamics, and Spectral Signatures: Changes Upon Model Cavity−Ligand Recognition. J. Phys. Chem. B 2012, 116, 13774−13780. (87) Paesani, F. Hydrogen Bond Dynamics in Heavy Water Studied with Quantum Dynamical Simulations. Phys. Chem. Chem. Phys. 2011, 13, 19865−19875. (88) Kahn, L. R.; Goddard, W. A. Ab Initio Effective Potentials for Use in Molecular Calculations. J. Chem. Phys. 1972, 56, 2685−2701. (89) Stevens, W. J.; Basch, H.; Krauss, M. Compact Effective Potentials and Efficient Shared-Exponent Basis Sets for the First- and Second-Row Atoms. J. Chem. Phys. 1984, 81, 6026−6033. (90) Hay, P. J.; Wadt, W. R. Ab Initio Effective Core Potentials for Molecular Calculations. Potentials for K to Au Including the Outermost Core Orbitals. J. Chem. Phys. 1985, 82, 299−310. (91) Hay, P. J.; Wadt, W. R. Ab Initio Effective Core Potentials for Molecular Calculations. Potentials for the Transition Metal Atoms Sc to Hg. J. Chem. Phys. 1985, 82, 270−283. (92) Hurley, M. M.; Pacios, L. F.; Christiansen, P. A.; Ross, R. B.; Ermler, W. C. Abinitio Relativistic Effective Potentials with Spin-Orbit Operators. Ii. K through Kr. J. Chem. Phys. 1986, 84, 6840−6853. (93) Dolg, M.; Wedig, U.; Stoll, H.; Preuss, H. Energy-Adjusted Abinitio Pseudopotentials for the First Row Transition Elements. J. Chem. Phys. 1987, 86, 866−872.

(94) LaJohn, L. A.; Christiansen, P. A.; Ross, R. B.; Atashroo, T.; Ermler, W. C. Abinitio Relativistic Effective Potentials with Spin− Orbit Operators. III. Rb through Xe. J. Chem. Phys. 1987, 87, 2812− 2824. (95) Andrae, D.; Häußermann, U.; Dolg, M.; Stoll, H.; Preuß, H. Energy-Adjusted ab Initio Pseudopotentials for the Second and Third Row Transition Elements. Theor. Chim. Acta 1990, 77, 123−141. (96) Ross, R. B.; Powers, J. M.; Atashroo, T.; Ermler, W. C.; LaJohn, L. A.; Christiansen, P. A. Abinitio Relativistic Effective Potentials with Spin−Orbit Operators. IV. Cs through Rn. J. Chem. Phys. 1990, 93, 6654−6670. (97) Stevens, W. J.; Krauss, M.; Basch, H.; Jasien, P. G. Relativistic Compact Effective Potentials and Efficient, Shared-Exponent Basis Sets for the Third-, Fourth-, and Fifth-Row Atoms. Can. J. Chem. 1992, 70, 612−630. (98) Iyengar, S. S.; Schlegel, H. B.; Millam, J. M.; A. Voth, G.; Scuseria, G. E.; Frisch, M. J. Ab Initio Molecular Dynamics: Propagating the Density Matrix with Gaussian Orbitals. II. Generalizations Based on Mass-Weighting, Idempotency, Energy Conservation and Choice of Initial Conditions. J. Chem. Phys. 2001, 115, 10291−10302. (99) Deng, J.; Gilbert, A. T. B.; Gill, P. M. W. Density Functional Triple Jumping. Phys. Chem. Chem. Phys. 2010, 12, 10759−10765. (100) Deng, J.; Gilbert, A. T. B.; Gill, P. M. W. Approaching the Hartree−Fock Limit by Perturbative Methods. J. Chem. Phys. 2009, 130, 231101. (101) Littell, T. R.; Skeel, R. D.; Zhang, M. Error Analysis of Symplectic Multiple Time Stepping. SIAM Journal on Numerical Analysis 1997, 34, 1792−1807. (102) Bishop, T. C.; Skeel, R. D.; Schulten, K. Difficulties with Multiple Time Stepping and Fast Multipole Algorithm in Molecular Dynamics. J. Comput. Chem. 1997, 18, 1785−1791. (103) Ma, Q.; Izaguirre, J. A.; Skeel, R. D. Verlet-I/R-Respa/Impulse Is Limited by Nonlinear Instabilities. SIAM Journal on Scientific Computing 2003, 24, 1951−1973. (104) Morrone, J. A.; Markland, T. E.; Ceriotti, M.; Berne, B. J. Efficient Multiple Time Scale Molecular Dynamics: Using Colored Noise Thermostats to Stabilize Resonances. J. Chem. Phys. 2011, 134, 014103. (105) Leimkuhler, B.; Margul, D. T.; Tuckerman, M. E. Stochastic, Resonance-Free Multiple Time-Step Algorithm for Molecular Dynamics with Very Large Time Steps. Mol. Phys. 2013, 111, 3579− 3594. (106) Chin, S. A. Dynamical Multiple-Time Stepping Methods for Overcoming Resonance Instabilities. J. Chem. Phys. 2004, 120, 8−13.

L


Multiple-Timestep ab Initio Molecular Dynamics Using an Atomic Basis

Recommend Documents