An efficient and stable hybrid extended Lagrangian/self-consistent field scheme for solving classical mutual induction

Polarizable empirical force fields offer a clear and systematic improvement over current generation fixed charge force fields by including many-body effects that allow for a molecular response to evolving heterogeneous environments.^1–6 However, the primary computational expense of a classical polarization model for large systems resides in the solution of a linear system of equations for the induced dipoles. They may be solved exactly by matrix inversion or Cholesky factorization for small systems¹ or approximately by an iterative self-consistent field (SCF) solution to the induced dipoles for the larger systems encountered in a condensed phase simulation. With conservative methods such as successive over-relaxation (SOR),⁷ the SCF costs for an induced dipole convergence of 10⁻⁵–10⁻⁸ D is up to 10 times the energy and force cost of a fixed charge model, and thus better SCF methods are required to reduce this computational cost in order to maintain effective sampling of thermalized systems with many degrees of freedom.

One choice is to implement an effective iterative scheme such as a pre-conditioned conjugate gradient SCF (CG-SCF)⁸ approach or direct inversion in the iterative subspace (DIIS)^9,10 along with a predictor to accelerate the convergence of the SCF problem. However, predictors use information from previous steps and are time irreversible, leading to an inevitable degradation in energy conservation. Another approach is to replace the SCF step with an extended Lagrangian (EL) formulation to avoid any iterative SCF costs.^2,3,11,12 However, this approach can be plagued with problems of accuracy since EL formulations allow the induced dipoles to fluctuate around an average orientation that does not strictly conform to the true electric field vector. Moreover, this approach can suffer from problems of stability and energy conservation in the context of a molecular dynamics (MD) trajectory that forces the time step to be unacceptably short.

Recently, Niklasson and co-workers^13–16 have introduced a hybrid EL/SCF scheme in the context of Born-Oppenheimer molecular dynamics (BOMD) wherein an extended set of auxiliary electronic degrees of freedom serves as an initial guess of the SCF solver. This allows less strict convergence of the ground state electron density due to the benefits of a time reversible Verlet algorithm that realizes excellent energy conservation. However, when analyzed over timescales that exceed ∼100 fs, it was observed that the accumulated error in the auxiliary electronic degrees led to eventual numerical instability due to numerical noise that increases without loss.¹⁷ Subsequently, it was shown that the introduction of a dissipative term can be effectively parameterized to offer the best compromise between maximizing numerical stability while minimizing the inevitable introduction of time irreversibility in the auxiliary equation of motion.¹⁷ 

In this work, we have extended the hybrid EL/SCF approach to the problem of classical polarization, in which the initial guess for the mutual induction calculation is treated by additional dynamical induced dipole variables evolved via a velocity Verlet scheme.¹⁸ The benefit to such an approach is that we can integrate reversible equations of motion of the auxiliary dipole initial guess at the same large time step as the atomic positions, while maintaining superior energy conservation since the polarization response stays near the Born-Oppenheimer surface at looser convergence levels relative to standard SCF solvers. However, like the original hybrid EL/SCF approach used in BOMD for the electron density matrix, we find that the set of auxiliary induced dipoles also exhibit a similar problem in numerical stability. This is manifested as a continued (and likely boundless) increase in the number of SCF cycles to meet even loose (10⁻¹ D) convergence tolerances over the course of a 1 ns trajectory of the AMOEBA14 water model,¹⁹ and we show that adding dissipation to the auxiliary equation of motion as described in Ref. 17 has a detrimental effect on energy conservation for small numbers of SCF cycles.

We have diagnosed the problem in the hybrid EL/SCF scheme applied to classical polarization as arising from resonances in the equations of motion^20–22 that manifests as a buildup of inertia for the auxiliary dipoles. Although in principle we could address the resonance problem with a smaller integration time step, instead we have formulated a new restrained inertial EL/SCF (iEL/SCF) method that for all intents and purposes controls for the resonance problem, analogous to other isokinetic approaches,²² such that the equations of motion of the auxiliary dipoles remain stable and time reversible. The results on the Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) polarization model show that the iEL/SCF method exhibits excellent energy conservation and thermodynamic and dynamic properties, but at greatly relaxed real dipole convergence tolerances, which reduce the number of SCF cycles relative to standard SCF solvers used for the classical polarization solution. As such, the iEL/SCF scheme clearly offers a better choice for classical mutual induction calculations compared to many EL and SCF and hybrid alternatives and is worthy of investigation for application to BOMD as well.

In this work, we develop our approach on the classical polarizable force field AMOEBA.^23,24 In addition to fixed multipole electrostatics, the AMOEBA model provides a consistent treatment of intramolecular polarization and intermolecular polarization and uses a physically motivated Thole damping scheme for local polarization effects to avoid the well-known polarization catastrophe that results when mutually inducible sites polarize each other to infinity at short inter-site separation.^25,26 AMOEBA’s many-body polarization energy, U_pol, is given by

where

and the elements of $μ ⇀ ind$ are defined as

where $μ ⇀ i$ is the inducible dipole at atom site i, α_i is the isotropic polarizability of atom i, T_ij is the rank-two interaction tensor between atoms i and j containing derivatives of 1/r_ij according to the permanent multipole expansion, T_ij′ is the corresponding interaction tensor for inducible dipole-dipole interactions, $M j ( d )$ are the permanent multipole moments; the T and (rank-one) M tensors encompass the 13 permanent multipole moments for the AMOEBA potential (q, μ_x, μ_y, μ_z, Q_xx, Q_xy, Q_xz, Q_yx, Q_yy, Q_yz, Q_zx, Q_zy, Q_zz), and the superscripts (d and p) refer to special scaling factors used for electrostatic interactions in AMOEBA.

The SCF method for AMOEBA implemented in TINKER has up until recently used a conservative SOR method⁷ but has been replaced in the TINKER 7.0 release with a preconditioned CG-SCF method using a predictor,⁸ while other AMOEBA implementations have used direct inversion of the iterative subspace.^9,10,27 These more recent SCF solvers are more efficient iterative methods that converge in fewer steps compared to SOR. In our comparisons to the various mutual induction calculation approaches such as EL^2,11 and hybrid EL/SCF schemes,¹³ we use the default CG-SCF with predictor method in TINKER.

We have implemented an extended Lagrangian formalism^11,12 (directly analogous to the Car-Parrinello approach¹²) for AMOEBA in the TINKER package that treats the polarization degrees of freedom as additional dynamic variables in the system,¹¹ allowing us to integrate them on the same footing as the atomic positions and avoid using self-consistent iteration to obtain polarization near the Born-Oppenheimer surface. The extended Lagrangian for point dipoles is given by

where m_i is the mass of atom i, $r ⇀ i$ is the position of atomic center i, $U ( r ⇀ N , μ ⇀ N )$ is the potential energy from the AMOEBA force field with the only difference being that the dipoles are now dynamically integrated and not iteratively converged. As a result there is a kinetic energy contribution from the dipoles given as the second term on the right hand side of Eq. (2), where m_μ,i is a fictitious dipole mass given in units of ps²/ų. In addition, we also thermostat these extended system degrees of freedom to a very low temperature (∼1 K)^2,3 to maintain the polarization close to the Born-Oppenheimer surface using Nosé-Hoover (NH) thermostats. The complete extended system equations of motion for Nosé-Hoover temperature control on both the atomic centers and induced dipoles using a single Nosé-Hoover thermostat are

where the thermostat “position,” η, couples to the physical system to control the temperature and the thermostat mass, Q, which is related to a characteristic time parameter, τ, by Q = N_fk_BTτ². In these equations “*” denotes quantities associated with polarization degrees of freedom so N_f and T are the number of degrees of freedom and temperature of the atomic coordinates, respectively, and N_f* and T_* are the degrees of freedom and the temperature of polarization degrees of freedom, respectively, associated with a low T_*. It should be noted that $F ⇀ μ ⇀ , i = E ⇀ i − μ ⇀ i α i$⁠, which goes to 0 in the limit of a self-consistent solution of the dipoles, $μ ⇀ i = α i E ⇀ i$⁠.

All results reported here are for pure water systems of 512 molecules, although all methods described should be generalizable to any molecular system; we use the water parameters of the AMOEBA14 water model.¹⁹ All simulations started from a pre-equilibrated box and long-range electrostatics were treated with particle-mesh Ewald²⁸ with a real-space cutoff of 9 Å. The equations of motion of the atomic degrees of freedom were integrated using the velocity Verlet method¹⁸ and the Nosé-Hoover formalism with a fourth-order chain was used for temperature control²⁹ with a τ of 0.1 ps. We used a time step of 1.0 fs for base AMOEBA simulations. For the extended Lagrangian simulations, we used a T_* of 1.0 K and a τ_* of 0.1 ps with 1500 iterations of the thermostat per step; the “mass” associated with inducible dipoles was set to 3.6 × 10⁻⁹ ps²/ų; we explored time steps from 0.25 to 1.0 fs depending on the method used to solve mutual induction.

We have adapted the approach developed by Niklasson and colleagues,^13–16 originally formulated for BOMD, but now extended by us to classical mutual induction calculations. In particular, this original hybrid EL/SCF method introduced an initial guess for SCF calculations by propagating a set of auxiliary electronic degrees of freedom in a time-reversible manner. In its original form, these auxiliary variables corresponded to the electronic ground state density matrix and here, we adapt these to the case of classical polarization formulated as induced dipoles.

In the spirit of Niklasson et al.,^13,15,16 we can define an extended Lagrangian given by Eq. (4) for the induced dipoles

where $μ ⇀ SCF N$ represents the set of all converged real induced dipoles, and we introduce another set of induced dipoles, $μ ⇀ N$⁠, which are the initial guesses to the iterative solution of $μ ⇀ SCF N$⁠. This auxiliary set of induced dipoles is restrained to stay near the true self-consistent values via the final term in Eq. (4) using a harmonic function where m_μ,i and ω are the fictitious mass and a universal frequency that determines the curvature of the harmonic well, respectively. Applying the Lagrangian equation of motion to Eq. (4) in the limit of m_μ,i → 0 yields the equations of motion for atomic centers and induced dipoles,

Eq. (5a) shows that equations of motion for the atomic centers are propagated in the usual way, except that the iterative solution uses an initial guess that is propagated by the auxiliary electronic degrees of freedom in Eq. (5b). We integrate both equations of motion using time-reversible velocity Verlet integration.¹⁸ To determine $μ ⇀ SCF N$⁠, we still use CG-SCF, but now the time reversible auxiliary dipoles serve as an initial guess. We chose ω to be $2 /Δt$⁠, where Δt is the time step,¹³ which we set to 1 fs.

The one drawback of the original hybrid EL/SCF scheme^13–16 is that over longer trajectories, the propagated auxiliary dipoles systematically degrade as a reasonable initial guess for the subsequent SCF steps using reasonable time steps. Figure 1(a) shows that just beyond the 100 fs time scale the original hybrid EL/SCF^13,15,16 requires an increasing number of SCF cycles, eventually reaching up to 5 SCF iterations to meet even a loose criteria of 10⁻¹ D after 1 ns. To put that in perspective, the standard CG-SCF scheme with a predictor requires 5 SCF cycles to reach a convergence of 10⁻⁶ D at all timescales. In fact the hybrid EL/SCF scheme starts at 6 SCF steps at this corresponding 10⁻⁶ D convergence level and increases to 8 SCF cycles over the 1 ns trajectory. In an attempt to alleviate this problem, we considered the use of higher-order symplectic integrators, which use multiple force evaluations per time step. The dotted line in Figure 1(a) gives the results of a 4th-order integrator³⁰ applied to this hybrid EL/SCF scheme at the lowest level of convergence. These results show that while higher-order integration stretches the time scale of the increase of SCF cycles it does not eliminate the problem. The number of SCF cycles would be expected to continue to rise based on the data shown in Figure 1(a).

The reason for this is presented in Figure 1(b) which shows that over the course of the simulation, the ensemble average of the auxiliary dipole velocities $〈 μ ⇀ ̇ i 2 〉$⁠, a “pseudo temperature,” increases continuously throughout the simulation, and this inertia eventually swamps the harmonic restoring force that aims to keep the auxiliary dipoles close to the real, converged dipoles, ultimately leading to instability in the equations of motion Eq. (5b). This seems to be a problem with resonances²⁰ existing in the auxiliary dipole equations of motion that are on a faster time scale than that experienced by the real induced dipoles through the subsequent SCF solver. As Figure 2 shows, the auxiliary dipoles show high frequency behavior compared to their real counterparts owing to the optimal choice for their characteristic frequency, $2 /Δt$⁠, and their direct coupling in the auxiliary potential (Eq. (4)) leads to corruption of the dynamics.

Niklasson and co-workers attributed this behavior to accumulation of numerical errors throughout a simulation and sought to mitigate it with a Langevin-like scheme¹⁷ that introduces an explicit dissipative force on the motion of the electronic degrees of freedom. The introduction of a dissipative force will inevitably lead to some time irreversibility, and hence, an optimization scheme was introduced in that study to maximize stability and minimize the undesired time irreversibility that will degrade energy conservation.¹⁷ We have implemented the 9th-order version of this scheme and its results are given in Figure 3. Figure 3(a) shows that the dissipative scheme corrects for the increasing number of SCF cycles. However, Figure 3(b) and Table I show that the benefits of small numbers of SCF iterations comes at the cost of unacceptable energy drift at loose levels of convergence due to its time irreversibility. This leads to the conclusion that dissipation of pseudo kinetic energy is important in achieving a stable number of SCF iterations, but that dissipation schemes with acceptable energy conservation do not significantly reduce the number of SCF iterations relative to standard SCF solvers.

Here, we present a different solution to this instability problem by “thermostating” the auxiliary dipoles through modification of their velocities in the time reversible velocity Verlet integration. It is known that corruption of the dynamics due to resonances can be controlled using the isokinetic ensemble,²² although we cannot formally implement an isokinetic scheme due to the m_μ,i → 0 limit that yields Eqs. (5). Since there are no longer any contributions to the total energy from the auxiliary dipoles and thus we cannot formally define their kinetic energy or temperature, we nonetheless show that by rescaling the auxiliary velocities we can execute control on the mean squared velocity (or pseudo temperature), $〈 μ ⇀ ̇ i 2 〉$⁠, that controls the buildup of this inertial pseudo kinetic energy quantity.

For iEL/SCF, we have implemented both a weak coupling Berendsen velocity rescaling scheme³¹ and a fourth-order time-reversible Nosé-Hoover chain (NHC)³² to control for the inertial accumulation. In this context, the distinction between Berendsen and NHC is likely unimportant since the auxiliary dipoles are serving as an initial guess to the SCF equations. Hence, although Berendsen velocity rescaling does not lead to the correct limiting canonical ensemble as does the NHC scheme, we note that the primary property needed is a damping scheme for the auxiliary variables that can be formally proven to be numerically stable, close to the exact solution regardless of convergence level, minimally perturbs time reversibility, and keeps the pseudo kinetic energy of the auxiliary dipoles constant. More sophisticated schemes for isokinetic integrators would perhaps be more desirable,^21,22 but we show that the simple velocity rescaling approach is sufficient in this case, since we are thermostating the auxiliary variables whose only role is to provide an initial guess to the SCF solver of the real inducible dipoles.

We define the Berendsen rescaling factor, α, which scales the velocities propagated by a reversible velocity Verlet integration at each time step in the weak coupling regime

where τ is a rescaling time scale parameter and T_μ is the set pseudo temperature of the auxiliary induced dipoles that corresponds to the desired value of $〈 μ ⇀ ̇ i 2 〉$ and has units of e²Å²/ps². The pseudo temperature chosen for the auxiliary dipoles is chosen to approximately conform to equipartition of energy consistent with a classical harmonic oscillator given the form of the auxiliary dipole potential

We can estimate the maximum auxiliary dipole velocity by using the square of the maximum displacement of the real induced dipole distribution to approximate $〈 ( μ ⇀ SCF , i − μ ⇀ i ) 2 〉$⁠. Using $T μ =〈 μ ⇀ i 2 〉$ and $ω 2 = 2 Δ t 2$ with $〈 ( μ ⇀ SCF , i − μ ⇀ i ) 2 〉∼ ( 0 . 2 e Å )$⁠,² gives a pseudo temperature of ∼10⁵e²Å²/ps² which is what we use here. A discussion of the numerical stability of our method is given in the  Appendix.

In what follows, we characterize the relative performance of the CG-SCF solver with predictor, the standard EL method, the hybrid EL/SCF scheme with no dissipation, and our new iEL/SCF scheme over a 1.0 ns simulation of the AMOEBA14 polarizable force field for water in the NVE and NVT ensembles. In all cases, we utilize the CG-SCF solver with predictor in which dipoles are converged to 10⁻⁶ D as the gold standard, with the understanding that over longer timescales there will be a noticeable energy drift with the CG-SCF approach even at this relatively tight dipole convergence tolerance.

Figure 4 shows energy conservation in the NVE ensemble for a standard CG-SCF scheme with predictor, a standard EL scheme, and the original EL/SCF scheme with no dissipation. On the nanosecond time scale examined here, we can see a significant energy drift for CG-SCF at a convergence level of 10⁻⁴ D and a severe lack of energy conservation with even looser convergence criteria (Figure 4(a)). Although the energy conservation of the SCF scheme looks stable for convergence levels of 10⁻⁵ D or tighter over the 1 ns trajectory, in fact the energy will eventually drift over longer time scales due to its one-sided convergence and use of a predictor. The standard EL scheme shows poor energy conservation, even for small time steps of 0.25 fs and 0.5 fs for which the EL method is stable (Figure 4(b)); time steps larger than 0.5 fs are numerically unstable.

For the hybrid EL/SCF without dissipation (Figure 4(c)), there is very good energy conservation up to a convergence level of 10⁻³ D, but for looser levels of convergence, the system energy experiences continual adjustment due to the poor quality of the SCF guess generated by the auxiliary dipoles. These deviations do not appear to be divergent and/or large in magnitude compared to the standard SCF and EL methods and eventually settle down at the end of the 1 ns trajectory. However, although there is some improvement in energy conservation for the hybrid EL/SCF without dissipation compared to the standard CG-SCF with predictor and EL schemes at a given level of convergence of the mutual induction, it comes at the cost of increasing numbers of SCF cycles due to the accumulated inertia over time (Figure 1(a)). In fact, the rebounding energy profile at the lowest level of induced dipole convergence for the no dissipation scheme corresponds exactly to a systematic bump up in the number of SCF cycles as the accumulated pseudo kinetic energy in the initial SCF guess gets worse over time (Figure 1(a)). Even so, the hybrid EL/SCF without dissipation outperforms the standard EL and SCF approaches when we fit the energy conservation profile to derive total energy drift rates (Table I).

Figure 5 shows energy conservation for the iEL/SCF using weak coupling Berendsen velocity rescaling and NHC for controlling $〈 μ ⇀ ̇ i 2 〉$⁠, for which mutual induction tolerances <10⁻² D using either thermostating methods show superior energy conservation compared to CG-SGF at tight tolerances >10⁻⁵ D, with the significant additional benefit that fewer SCF steps are required; using the iEL/SCF method, only 4 SCF steps are needed for 10⁻² D and 3 SCF steps for 10⁻¹ D, compared to the 5 SCF steps needed by CG-SCF for a convergence level of 10⁻⁶ D. Both velocity attenuation methods show superior energy conservation compared to all methods over the full range of induced dipole convergence (Table I) and correct for the “rebounding” energy behavior seen in the non-dissipative EL/SCF method at low convergence levels (Figure 4(c)). We note that both thermostating methods begin to show systematic energy drift at 10⁻¹ D levels of convergence. Since the NHC thermostat is formally formulated as a time reversible scheme, we conclude that we have reached an inherent limit to how low we can set the convergence tolerance in the SCF solutions, although the energy drift is at least one if not many orders of magnitude better than any other possible scheme at 3 SCF steps for 10⁻¹ D.

Next, we consider in detail how the induced dipole polarization is accounted for in each method within an NVT simulation at 298.0 K. Figure 6 shows the probability density function for the induced dipole magnitude and Figure 7 shows the normalized dipole autocorrelation function at short times (<1 ps) for the AMOEBA14 water model; Figures S1 and S2 in the supplementary material provide the probability density of the in-plane and out-of-plane induced dipole angle as well.³³ We see that the CG-SCF with predictor method shows significant deviations from the gold standard (10⁻⁶ D) starting at convergence levels of the induced dipoles at 10⁻³ D. Whereas the standard EL method shows good agreement with the SCF converged result for time steps of 0.25 fs and 0.5 fs, it fails completely at larger time steps due to numerical instability. By contrast, the hybrid EL/SCF scheme without dissipation reproduces correct polarization probability distributions and the dipole autocorrelation well, even at 10⁻¹ D, although again this comes at the cost of increasing numbers of SCF cycles. The iEL/SCF approach developed here also reproduces all polarization properties well at loose convergence and uses the fewest number of SCF cycles to do so. All of the methods show some degradation in the range of polarization properties at a given level of convergence when utilized in the NVE ensemble (Figures S3–S6 and Table S1 in the supplementary material³³), which is not unexpected since thermostating of the real system variables can mask underlying numerical problems in integrators and/or poorly converged energy and forces. However, the hybrid EL/SCF schemes are superior to the standard EL and CG-SCF at any dipole convergence tolerance.

Table II reports the average potential energy and molecular dipole which were calculated from a single NVT simulation at 298 K, as well as the diffusion coefficient, calculated by taking independent snapshots from the 298.0 K NVT simulation and calculating the molecular mean squared displacement in the NVE ensemble. The SCF scheme gives correct average potential energies and molecular dipoles down to a relatively loose convergence of 10⁻² D, a benefit of thermostating in regards thermodynamic quantities, but gives incorrect diffusion coefficients starting at a much tighter tolerance of 10⁻⁴ D. For the EL method, even using a smaller time step than the other two approaches, we see more significant differences between the average potential energies and molecular dipoles compared to the best converged SCF result and poor results for the diffusion coefficient. Across all convergence levels, all EL/SCF schemes reproduce the thermodynamic and kinetic data quite well, even producing a diffusion coefficient within error bars of the gold standard at much looser convergence levels. Although ultimately the computational cost of the increasing number of SCF cycles is early evidence for eventual numerical instability in the non-dissipative EL/SCF, the iEL/SCF method retains excellent property performance with a stable algorithm that conserves energy at loose convergence and minimizes the number of SCF cycles relative to the standard CG-SCF with predictor method.

We have presented a new adaptation of a hybrid EL/SCF scheme applied to induced dipole polarization in classical simulations that overcomes numerical instability problems that were also observed in the equations of motion for the auxiliary electronic degrees in BOMD simulations.¹⁷ However instead of using a dissipative force that can compromise time reversibility of the auxiliary induced dipole initial guess, we have introduced the use of thermostats applied to the auxiliary dipole velocities to control for the accumulation of a pseudo kinetic energy due to resonances. This approach simultaneously preserves energy conservation and numerical stability at loose real dipole convergence values of 10⁻¹ D, thereby requiring fewer SCF cycles to describe polarization and system properties accurately when compared to a typical CG-SCF convergence of 10⁻⁵–10⁻⁶ D. For the AMOEBA14 water model, this also translates to an overall decrease in the required number of SCF iterations per time step from 5 to 3, resulting in some computational savings.

In the original explicit dissipative schemes explored in Ref. 17, a flexible but complicated optimization scheme was introduced to parameterize the functional form of the needed dissipation but suffers from severe time irreversibility problems at convergence tolerances that yield the most benefit in number of SCF iterations. We have shown that the actual problem arises from the fact that the dynamics of the auxiliary dipoles evolve on a faster time scale than the time evolution of the real induced dipoles and that resonances arise due to their coupling through the potential energy term in Eq. (4). This leads to a buildup of inertia in the auxiliary dipoles that can be solved by either reducing the time step for their equation of motion or by controlling it through pseudo temperature control as we have done here. Our iEL/SCF approach is a simple velocity attenuation scheme that effectively removes the effect of these resonances while still maintaining time-reversibility (within that allowed by the numerical integration scheme) depending on the thermostat method chosen. We are thus able to have better energy conservation and can reproduce molecular properties, even at very loose levels of convergence with only a small number of SCF iterations required. While the iEL/SCF method has been demonstrated on a classical inducible dipole model AMOEBA (and is available in the TINKER program), it may also be useful in the context of ab initio molecular dynamics.

We thank the National Science Foundation Grant No. CHE-1363320 for support of this work. We would like to thank Mark Tuckerman for helpful discussions and enjoyable collaboration.

Following the numerical stability analysis of Niklasson and colleagues,¹⁷ we have developed a Verlet-like recursion for the Berendsen rescaling for analysis purposes. Adding two Taylor series expansions of the dipole “position,”

and

to obtain

Using a scaled velocity from the velocity Verlet recursion

we can substitute for $μ ⇀ ̇ n$ in Eq. (A3) to obtain

Finally an approximation for $μ ⇀ ̇ n + 1$ is given by a finite difference

and Eq. (5b) is substituted for $μ ⇀ ̈ n + 1$ to obtain the recursion

where $μ ⇀ n = μ ⇀ i ( t 0 + n Δ t )$ and α_n = α(t₀ + nΔt) are the instantaneous velocity rescaling factor. Note that when there is no pseudo temperature control (α_n = 1), we recover the Verlet recursion, as expected. To analyze the numerical stability of the velocity rescaling scheme, Eq. (A7) has a characteristic equation given by

where λ corresponds to the roots of the characteristic equation and γ is the largest eigenvalue of the iterative response matrix acting on the difference between $μ ⇀ n$ and the true solution of the induced dipoles. Thus, Eq. (A8) allows us to explore various levels of convergence of the solution of the auxiliary dipoles over the range −1 < γ < 1, where in the limit γ → 0 corresponds to the exact solution, while the magnitude of the largest root, $λ m a x$⁠, will determine the stability of the recursion, with $λ max >1$ giving exponentially growing solutions, $λ max <1$ giving exponentially decaying solutions, with an exactly stable solution occurring at $λ max =1$⁠. We note that α≠1 corresponds to increasing time irreversibility in the equations of motion, which should be avoided when possible. Time-reversibility is preserved in our iEL/SCF method when α = 1, whereas a related coupling parameter that multiplies the dissipation force, α_diss, discussed in Ref. 17 remains time reversible when α_diss = 0.

Similarly, we can show that for a Nosé-Hoover chain thermostat (taking into account rescaling that happens twice per time step now) that the Verlet-like recursion is given by

where α_n and α_n−1 now represent the velocity scaling factors that come from a time-reversible Nosé-Hoover method at the end and beginning of the n-th recursion. These two α-values mask the complexity of the Nosé-Hoover chains which are also a function of additional extended system variables; see the work of Martyna et al. for details.³² Again, we can construct the characteristic equation of Eq. (A9) for analysis

Figure 8(a) gives $λ max$ as a function of the SCF convergence, γ, for various values of α (assuming α_n = α_n−1 for the NHC case). Intuitively, we find that for α > 1, then $λ max >1$ and the equations of motion are unstable which would correspond to an accumulation in the pseudo kinetic energy as observed in Figure 1(b). For α < 1, the increasing dissipation will realize stable solutions $( λ max < 1 )$ but at the expense of time-reversibility as α decreases. Thus, both the Berendsen weak coupling and NHC iEL/SCF schemes have the desirable property that the equations of motion can be made stable under incomplete SCF convergence in the full γ interval for α values that are close to that needed for time reversibility. We would like to note that Eqs. (A8) and (A10) are presented for analysis purposes only and the method is truly implemented using a velocity Verlet scheme with thermostat action being applied at the appropriate points within such a scheme.

Figure 8(b) shows the simulated trajectory of the rescaling parameter α during the course of our weak coupling Berendsen velocity rescaling as well as NH scheme which is shown to range from 0.9997 to 1.0003 with an average of ∼0.999 99 such that we are essentially close to the exact time–reversible solution. While velocity rescaling using weak coupling Berendsen formally breaks the time-reversibility of the integration scheme, an α-value so close to 1 corresponds to only a slight disturbance of this reversibility while dissipating the integration error that causes divergence in $〈 μ ⇀ ̇ i 2 〉$ due to resonances. In any event, errors in time reversibility are formally circumvented through use of NHC thermostats, also shown in Figure 8(b), although, at least for our test system of bulk water, the practical differences are largely unimportant. Thus, our diagnosis of the problem in the original hybrid EL/SCF scheme arises from resonances in the auxiliary equations of motion that can be controlled by a simple velocity rescaling scheme that prevents the accumulation of a pseudo kinetic energy for these degrees of freedom.