Iterative blip-summed path integral for quantum dynamics in strongly dissipative environments

In addition to their utility as simplified models of condensed phase phenomena,^1–10 system-bath (also known as generalized spin-boson) Hamiltonians offer useful quantitative descriptions of a variety of processes, such as those in crystalline solids (where the harmonic bath consists of lattice phonons), chemical reactions (where the bath modes correspond to molecular vibrations orthogonal to the reaction path^11,12), and even in some situations where anharmonicity at the microscopic level is quite significant. In the latter case, the harmonic bath¹³ is commonly associated with linear response (Gaussian statistics) and has been shown¹⁴ through a path integral analysis to emerge as the dominant term whenever the process of interest is affected by a large number of degrees of freedom. Recent mean-field calculations in crystals,¹⁵ as well as a rigorous quantum-classical path integral (QCPI) simulation of a charge transfer reaction in solution,¹⁶ corroborated the quantitative nature of effective harmonic bath descriptions for certain processes in strongly anharmonic environments.

Early studies of the dissipative effects from a harmonic bath on barrier crossing and tunneling were based largely on analytical approximations. The first comprehensive study of tunneling in two-level systems (TLSs) employed the non-interacting blip approximation,⁴ which allowed the approximate summation of the path integral with the dissipative effects from the harmonic bath included through the Feynman-Vernon influence functional.¹ Analytical studies of barrier crossing established the turnover^17,18 between weak and strong friction regimes, and energy relaxation was investigated using the Redfield approximations.¹⁹ Early numerical simulations using Monte Carlo sampling^20,21 to evaluate the discretized path integral were able to capture the short-time behavior of dissipative tunneling systems.^22–26 The development of quasi-adiabatic propagator path integral (QuAPI) methods with well behaved propagators,^27,28 optimal grid representations,²⁹ and numerically exact iterative decomposition schemes of the path sum,^30–41 which do not suffer from the Monte Carlo sign problem,^42,43 extended the feasibility of real-time path integral calculations to long times. Several alternative path integral-based and hybrid methods have also been reported in recent years (see, for example, Refs. 44–50). One of these approaches leads to hierarchical equations of motion^46,47 (HEOMs) and requires computational effort that depends strongly on the form of the bath spectral density. Finally, the multiconfiguration time-dependent Hartree (MCTDH) methodology⁵¹ has also been used to propagate weakly dissipative systems at low temperatures.⁵² 

The major task common to the numerical path integral schemes described above is the inclusion of the influence functional to alter the quantum mechanical amplitude of the bare system. When the system-bath coupling is weak, the dissipative effects captured in the influence functional are rather subtle, amounting to a modest modification of the dynamics from purely coherent to underdamped oscillations. On the other hand, strong system-bath coupling leads to a qualitatively different behavior, characterized by monotonic decay. The dissipative effects of the bath are so strong in this regime that the coherent dynamics of the isolated system no longer serves as a reasonable zeroth order description. Accounting for such a dramatic change of the dynamics is a challenging task that requires major computational effort. It is clear that a more efficient strategy would result from the use of an incoherent model as the starting point; rather than adding decoherence to the oscillatory system amplitude, one should start with a completely incoherent approximation and gradually include quantum coherence effects until convergence is achieved. Recent work⁵³ pointed out that such a fully incoherent starting point is readily obtained by forcing the forward and backward branches of the path integral to be identical. In the semiclassical limit, such a procedure forms the basis of forward-backward semiclassical and linearized semiclassical treatments,^54–65 which are quite adequate for systems (e.g., fluids) where quantum coherence effects are naturally quenched.^66–71 

In the fully incoherent limit, the discretized double path sum is replaced by a single path sum where the quantum mechanical amplitude along each forward path is multiplied by its complex conjugate corresponding to the identical backward path, thus no phases survive. Quantum coherence arises from phases contributed by forward-backward path pairs that differ over at least one time step. Such elementary time intervals (of length equal to the path integral time step) have been termed⁴ blips. Thus, the completely incoherent result is the zero-blip limit of the discretized path sum, and corrections are obtained by including path pairs that contain one or more blips.⁵³ 

These ideas led to the blip decomposition of the path integral⁵³ for a discrete system coupled to a harmonic bath. The blip series was shown to converge rapidly when the process of interest is characterized by incoherent dynamics. In favorable regimes, converged full-memory results for propagation to a few hundred time steps can be obtained by including only a few blips (and thus an extremely small number of terms relative to the full path sum). In addition, it was shown⁵³ that recasting the double path sum in terms of sum and difference coordinates maps entire segments of the inner between-blip sum on the time-local dynamics of an externally driven Hamiltonian. This mapping allows these portions to be evaluated by ordinary matrix procedures, reducing the computational effort dramatically. Even in the limit where all possible path pairs are included, producing the exact quantum mechanical result, the blip-summed path integral (BSPI) amounts to an exponential reduction of the number of operations compared to an explicit evaluation of the forward-backward path sum.

For a given total time, the contribution of forward-backward path pairs eventually decreases exponentially with the number of blips. Thus, once a sufficient number of blips have been included, the blip series converges rapidly. However, the required number $bmax$ of blips increases with propagation time, such that the inclusion of all contributing path pairs eventually becomes impractical. The evaluation of the path integral for very long times is possible via iterative procedures,^30–41 which exploit the decay of memory, i.e., the decreasing magnitude of the nonlocal influence functional interactions with time separation. If implemented within the BSPI formulation, memory truncation will reduce the value of $bmax$ to the number of blips required for convergence only within the memory length, which usually is much shorter than the total propagation time. Conversely, the implementation of the blip decomposition within the iterative path integral methodology offers an excellent way of filtering out path segments of negligible weight and—most importantly—obviates the explicit enumeration of between-blip state sequences, leading to exponential reduction of the number of operations and required storage. Thus, the iterative blip-summed path integral (i-BSPI) leads to dramatic savings in the simulation of long-time incoherent dynamics of quantum dissipative systems.

Section II reviews the blip decomposition of the discretized path integral for a discrete system coupled to a harmonic bath and describes its iterative implementation. The convergence of the procedure is demonstrated in Section III with application to a strongly dissipative asymmetric TLS. Some concluding remarks, along with extensions of the iterative blip methodology, are given in Section IV.

Throughout this paper it is assumed that the quantum mechanical system can be represented in terms of M spin-type or discrete variable representation^72–75 (DVR) states $|σn⟩$⁠, such that the system position operator is

A grid representation of the system coordinate²⁷ is essential for circumventing the “sign problem” that plagues Monte Carlo methods in real time.⁴³ The use of DVR states stems²⁹ from the need for a compact, minimal-sized basis of position-like states, which are necessary for evaluating the system-bath potential coupling in the factored time evolution operator. Further, truncation of the system’s Hilbert space eliminates the rapid oscillations of the real-time propagator,^28,76 leading to well-behaved integrands. In the chosen basis, the system Hamiltonian has the generic form

The system is linearly coupled to a bath of harmonic oscillators,

and the total Hamiltonian is

The characteristics of the bath pertaining to the dynamics of the system can be specified collectively via the spectral density function³ 

Conversely, the spectral density function can be used to obtain discrete bath parameters, which are useful for treating the bath via classical trajectories. Recent work⁷⁷ has shown that it is possible to map a polyatomic environment on a discrete harmonic bath directly through a time correlation function computed through a molecular dynamics (MD) simulation. System observables are obtained from the reduced density matrix

where the trace is with respect to the bath degrees of freedom. It is frequently assumed⁴ that the initial density operator factorizes and that the isolated bath is initially in equilibrium at a temperature $1/kBβ$⁠, where $kB$ is the Boltzmann constant, i.e.,

where $Hb$ is the Hamiltonian for the uncoupled bath (obtained from Eq. (2.3) by setting $cj=0$⁠) and $Zb$ is the canonical partition function of the bath. Expressing each of the forward and backward time evolution operators in the discretized path integral form⁷⁸ with a time step $Δt=t/N$ and integrating out the harmonic bath degrees of freedom leads to the following form for the reduced density matrix,

Here F is the DVR-discretized²⁹ Feynman-Vernon influence functional,¹ 

The short time propagators in Eq. (2.8) are obtained using the quasi-adiabatic splitting,²⁷ which offers the advantage of an increased time step compared to the Trotter⁷⁹ factorization. The coefficients $ηk′k″$ have been given in previous publications^32,80 in terms of integrals of the spectral density function, and recent work⁷⁷ has shown that they can also be obtained directly from MD data.

Equation (2.9) contains nonlocal couplings between non-adjacent time points, which prevent the stepwise evaluation of the path integral by means of simple matrix-vector multiplication,⁸¹ a procedure routinely applied to the propagation of low-dimensional systems in the absence of an influence functional. Fortunately, the magnitude of these couplings decays with time separation (as long as the bath is characterized by a continuous spectral density).³⁰ When the time interval $(k′−k″)Δt$ separating the path integral variables $sk′±$ and $sk″±$ exceeds the time over which the bath correlation function is nonzero, these variables become decoupled and need not be summed simultaneously. Inclusion of the time points that span the bath correlation (or “memory”) time leads to multi-time propagators,^30–32 which can be thought of as higher rank tensors³¹ (with respect to system coordinates) or ordinary matrices³³ (with respect to discrete path segments).

Thus, the iterative quasi-adiabatic propagator path integral (i-QuAPI) methodology offers a simple and efficient way of propagating the reduced density matrix in system-bath Hamiltonians characterized by any type of spectral density, provided the bath-induced memory is not very long (i.e., $Δkmax≲ 12$⁠). This is often the case when the bath consists of a continuum of modes, some of which have frequencies which are high compared to the frequency components of the system. The i-QuAPI methodology efficiently captures the coherent-incoherent transition in the tunneling dynamics of dissipative TLSs, as well as subtle phenomena arising from entanglement or time-dependent fields,^82–87 and has found numerous applications in chemistry, biology, or condensed matter physics. The most challenging regime is that of a low-frequency bath strongly coupled to the system, as many time steps must be included in that case to saturate the time nonlocality. Filtering procedures^33,34,39,40 discard the vast majority of path segments, which enter with exponentially small weights, retaining only those that are statistically significant. Such filtering techniques extend significantly the applicability of the i-QuAPI methodology to multi-level systems coupled to relatively slow solvents. The i-QuAPI methodology has been generalized to a bath in equilibrium with the initially populated (“donor”) state of the discrete system⁸⁸ and to equilibrium (non-product) initial conditions^37,38 and thus can be used to calculate linear or nonlinear time correlation functions.⁸⁹ The methodology has also been extended to fermionic baths.^90,91 Last, it has been shown⁹² that iterative decompositions of the real-time path integral are not restricted to harmonic baths but can be implemented in connection with any condensed phase environment as long as the influence functional of the latter can be computed (either exactly or approximately). Most recently, these ideas have found application within the quantum-classical path integral (QCPI) methodology,^93–97 which offers a rigorous treatment of the interaction between a quantum system and a complex polyatomic environment whose dynamics is captured via classical trajectories.

The role of blips becomes clear by expressing the Feynman-Vernon influence functional¹ in the form

where

are the sum and difference coordinates of each time point. It is clear that the two-time pair $k′≥k″$ contributes to the influence functional only if at least the later of the two time points is a blip, i.e., $sk′+≠sk″−$⁠. Through the real part of the exponent of Eq. (2.10), blips give rise to damping, quenching the coherent nature of the dynamics described by the bare system propagators. This type of “classical” decoherence, which is dominant at high temperature, is associated with stimulated phonon absorption and emission.¹⁰ Interestingly, the memory associated with this classical decoherence is removable through the introduction of auxiliary variables.⁹⁸ (This classical memory removal is exploited in the QCPI formulation,⁹³ where the bath is described in terms of trajectories sampled from a phase space distribution. However, within the Feynman-Vernon framework, which completely eliminates the bath variables, all memory interactions have to be included explicitly.) When the system-bath coupling is weak, forward-backward path pairs with several blips play the most important role in changing the dynamics from a fully coherent to an underdamped oscillatory behavior. However, if the dissipative effects of the bath are strong, the contribution of multi-blip paths is exponentially small. Thus, in strongly dissipative systems, the most important contribution to the path sum comes from path pairs with a relatively (in comparison to the total propagation length) small number of blips.

Rewriting the forward-backward path sum in terms of sum and difference coordinates, Eq. (2.8) becomes

where $sk+,sk−$ are obtained from Eq. (2.11). Here time points with $Δsk≠0$ are the blips, and those with $Δsk=0$ are the “sojourns,” time points where the forward and backward paths coincide. (Note that the inner summation index $s¯k$ at each time point in this expression depends on the outer summation index $Δsk$⁠. For example, in the case of a TLS with coordinate values $σ1=1,σ2=−1$⁠, the sum/difference coordinate pair takes on the values $Δs=±2,s¯=0$ and $Δs=0,s¯=±1$⁠.) Equation (2.12) provides a convenient starting point for a dramatic acceleration of the path sum in two ways:

• By restricting the outer sum to those path pairs that do not exceed a particular number $bmax$ of blips. As discussed above, in strongly dissipative environments convergence is achieved with relatively small values of $bmax≪N$⁠.

• Through the replacement of the inner sum by a series of iterative matrix multiplications. This is possible because all couplings of time points between blips are of nearest neighbor type (i.e., there are no sojourn-sojourn couplings in the influence functional). Even if all possible paths (i.e., $bmax=N+1$⁠) are included, this replacement amounts⁵³ to an exponential reduction of effort compared to Eq. (2.8).

In favorable regimes of incoherent dynamics, the blip decomposition enables the evaluation of the path integral to fairly long times.⁵³ For propagation to N time steps, the number of blip arrangements for a given value of $bmax$ is given by

In addition, for an M-state system there are $M2−M$ possible pairs of forward and backward coordinates at each blip time point $k=0,…,N$⁠. Thus the total number of forward-backward blip configurations is

This number grows rapidly with the number N of path integral time steps. In addition, the number of blips required for convergence generally increases with N. To enable propagation to long times, the blip representation of the path integral needs to be re-cast in an iterative form.

To this end, the influence functional couplings are truncated to the memory length of L time steps. The propagation involves an array of forward-backward path segments over the time points $k+1,…,k+L$⁠, which are connected to path segments spanning the time points $k,…,k+L−1$ through couplings arising from the system propagator and the influence functional. The present algorithm differs from the i-QuAPI methodology in two important aspects as follows:

• The maximum number of blips is restricted to the chosen value of $bmax$⁠, and

• Only the coordinates of blip points (plus one additional time point) are stored.

To begin, one generates all possible placements of $1,2,…,bmax$ blips over L time points. From each of these blip arrangements, one generates g blip configurations $r1,…,rg$ by assigning all possible values of the M system states as the forward and backward coordinates of each blip time point (ensuring the system coordinates of the forward and backward paths are not identical). For example, in the case of a TLS, there are two possible configurations of each blip: the one with forward and backward path coordinates $+1$ and $−1$ and the one with forward and backward path coordinates $−1$ and $+1$⁠, respectively. The number of blip configurations is given by Eq. (2.14) with $N=L−1$⁠, i.e.,

These blip configurations are stored and indexed by the blip locations and their forward-backward coordinate values, as well as the blip configurations to which they connect (as explained later). If all possible blips are included, i.e., if $bmax=L$⁠, the blip configurations correspond to the set of all possible forward-backward path segments that span the memory interval, and $g(L,L)=(M2−M+1)L$⁠. (Note that this number is smaller than the total number $M2L$ of forward-backward paths because sojourn segments are not enumerated.)

As seen from Eq. (2.10), all influence functional interactions connect a blip to an earlier blip or a sojourn. To avoid unnecessary increase of storage, the iterative propagation begins from the last time point $N$ and proceeds toward $t=0$⁠. The propagation involves iteration of an array $Rk(sk±,rj)$ that holds the stored blip configurations $rj,j=1,…,g$ over a memory interval corresponding to the time points $k,…,k+L−1$⁠, as well as the coordinate values $sk+=σ1,…,σM$⁠, $sk−=σ1,…,σM$ of the first time point of this interval (regardless of whether this point is a blip or a sojourn).

The initialization step generates the array $𝐑N−L$⁠. Next, each blip configuration over the time segment $N−L,…,N−1$ is augmented by the addition of the forward-backward path coordinates at the time point N (see Figure 1). The value of the array $𝐑N−L$ is obtained by including all the system propagator, blip-blip, and blip-sojourn interactions that connect these $L+1$ time points and summing over all sojourn points (except for point $sN−L±$⁠) and over the coordinates of point N. The sum with respect to sojourn points is performed from the latest toward the earliest time point via a series of $M×M$ matrix-vector multiplications as described in Ref. 53. This procedure avoids explicit enumeration of sojourn segments, which amounts to dramatic savings. Since the number of sojourns is smaller than L, this task scales sub-linearly with the memory length, allowing propagation even when the memory time is very long.

Iteration involves propagating $𝐑k+1$ backwards in time to $𝐑k,k=N−L−1,…,0$⁠. To achieve this, each blip configuration over the time segment $k,…,k+L−1$ must be matched with a blip configuration spanning the time interval $k+1,…,k+L$⁠. This is achieved very efficiently by creating at the beginning of the calculation and storing lists of blip configurations to which each blip configuration “connects.” Such matching configurations are illustrated in Figure 2. The propagation proceeds by including the influence functional interactions from the blips of the segment $k+1,…,k+L$ to the blip or sojourn of point k and the propagator between points k and $k+1$ and summing with respect to all blip configurations in the array $𝐑k+1$ and the coordinates of time point $k+1$⁠, if this is a sojourn.

Finally, once the initial time $(k=0)$ is reached, the desired elements of the reduced density matrix $ρ∼(NΔt)$ are obtained by summing $𝐑0$ with respect to the coordinates of all blip configurations located at $k=1,…,L−1$ and the coordinates $s0±$⁠.

The i-BSPI methodology is illustrated on an asymmetric TLS coupled to a dissipative bath. The TLS Hamiltonian is given by

where $ℏΩ=1$⁠, $ε=5$⁠, and $|1⟩,|−1⟩$ are localized states with coordinates $σ1=1,σ2=−1$⁠, while the spectral density has the common Ohmic form,

with $ξ=4$ and $ωc=2Ω$⁠. The bath is initially in equilibrium with the “donor” state 1 at a temperature $ℏωcβ=0.2$⁠. The chosen model has been shown⁹⁷ to exhibit very long memory. In particular, QCPI calculations with free- or donor-state bath reference⁹⁵ and memory truncation, which automatically account for all classical memory, failed to converge with realistic values of the memory length, indicating that the quantum memory is very long, and convergence was possible only with a dynamically consistent state hopping (DCSH) QCPI scheme⁹⁷ which captures some of the “quantum memory” effects¹⁰ through properly adapted system propagators.

The results of i-BSPI calculations are presented in Figure 3. Because of the strong system-bath coupling, the factorization of the propagator required a small time step $ΩΔt=0.075$⁠. The reduced density matrix was propagated to its equilibrium value, which required $N=600$ time steps. It is seen that the long-time plateau value of the donor population depends sensitively on the memory length included. Convergence to better than 0.01 required a memory length of $7.5Ω−1$⁠, i.e., the inclusion of influence functional couplings spanning $L=100$ path integral time steps. The long-time value of the converged donor state population is in excellent agreement with the thermodynamic value 0.27, which was obtained from independent calculations using the standard path integral Monte Carlo (PIMC) method. The converged $L=100$ i-BSPI results match the DCSH-QCPI results within the Monte Carlo statistical error of the latter.

The i-BSPI calculations converged with just $bmax=3$ blips. For $L=10$⁠, the $bmax=3$ i-BSPI results are identical (to three or more significant figures) to the values obtained with full i-QuAPI calculations with the same memory length. However, the i-BSPI methodology can easily handle much longer memory. With $L=100$⁠, the full i-QuAPI evaluation of the path integral would involve an array of $2200$ elements. While path filtering^33,34,39,40 would reduce this storage substantially, filtered iterative path integral calculations would still need to store path segments of length $L=100$⁠. The explicit enumeration of sojourn coordinates for the filter-retained paths would lead to an array of much larger (and likely unmanageable) size. By contrast, the propagated array in the i-BSPI calculation with $L=100$ stored only about 1.3 × 10⁶ elements.

The i-BSPI methodology presented in this paper offers an extremely efficient way of evaluating the real-time path integral for a discrete system strongly coupled to a harmonic bath that may include fast as well as very slow degrees of freedom. Such regimes are challenging to earlier path integral methods, where the starting point is the bare quantum system, which exhibits coherent behavior. By contrast, the crudest limit of the blip decomposition, the zero-blip term, describes completely incoherent dynamics, which is qualitatively similar to the true behavior of the system in this regime, thus providing a much better starting point. The inclusion of blips gradually adds quantum coherence, until convergence is reached. Because of the strong damping effect of blip-blip interactions, the blip series converges rapidly in strongly dissipative regimes. Just as importantly, the blip decomposition obviates the need for an explicit enumeration of between-blip sojourn path segments, thus eliminating the vast majority of paths, which would necessitate astronomical amounts of storage. For this reason, the BSPI methodology always (even in coherent regimes, where a large number of blips may be required for convergence) leads to exponential acceleration compared to the full evaluation of the discretized path integral with the same time step.

While evaluation of the full blip sum is often possible for fairly long times, the number of blips required for convergence tends to grow with propagation time, eventually becoming prohibitive. The i-BSPI methodology described in this paper enables propagation to very long times by reducing the required number of blips to those required only within the memory interval. Representative calculations presented in Section III show that convergence is reached rapidly, even in regimes of very long memory. It is thus expected that the i-BSPI methodology will offer an efficient and convenient approach for simulating quantum mechanical processes in sluggish, strongly dissipative environments, which tend to be rather challenging.

While implementation of the blip decomposition always reduces the number of operations, the older i-QuAPI algorithm still offers a better choice in weakly dissipative regimes. This is so for several reasons: First, the number of blips required for convergence tends to be large in such parameter regimes, implying less dramatic gains. Second, the i-QuAPI methodology easily yields results at each time step, while intermediate results cannot be extracted easily in the present version of the i-BSPI algorithm. Last, the i-BSPI code is more complex than i-QuAPI.

The development during the last decade or so of alternative approaches to i-QuAPI has allowed accurate simulations of system-bath dynamics in regimes previously inaccessible. As noted earlier, the QCPI methodology^93–97 (which is designed for treating complex, anharmonic condensed-phase environments) can also be used to treat harmonic bath dynamics by discretizing the spectral density⁹⁹ to produce explicit harmonic oscillators, whose initial conditions are sampled by Monte Carlo. In particular, the cumulative quantum memory form¹⁰⁰ of QCPI provides an alternative and powerful approach that is well suited to strongly dissipative dynamics. Since the influence functional-based i-BSPI avoids the bath discretization and the use of Monte Carlo sampling, it offers higher efficiency. Further, a methodology based on truncated hierarchical equations of motion^46,47 for auxiliary density matrices has emerged as an alternative to path integral/influence functional-based schemes. Unlike the latter, the convergence of the HEOM approach depends critically on the number of exponential terms required to fit the time correlation function of the bath. Thus, while HEOM calculations for the Debye spectral density converge easily, in particular at high temperatures, other spectral density functions can be quite challenging to that method. A stochastic method that is similar in spirit to generalized Langevin dynamics offers yet another alternative for simulating spin-boson dynamics.¹⁰¹ Last, a recently proposed Monte Carlo based “inchworm” method^102,103 appears promising.

Besides its use for simulating dissipative dynamics in harmonic environments, the i-BSPI methodology holds promise for accelerating the convergence of path integral calculations in anharmonic media. Such calculations have recently become possible with the development of the QCPI methodology.^93–97 QCPI offers a rigorous way of describing the interaction of a quantum system with a complex polyatomic environment, which is described through classical trajectories in full atomistic detail. In particular, the local nature of the quantum paths leads to forces on the classical trajectories that are unambiguously determined, obviating the need for Ehrenfest-like averaging over delocalized wavefunctions. The recent simulation¹⁶ of an electron transfer reaction in solution, where 1320 atoms interacting with CHARMM force fields were treated explicitly, showed that such calculations are feasible without any assumptions or uncontrolled approximations. The major cost of QCPI calculations stems from the computation of classical trajectories, whose number exceeds considerably that required in a typical MD simulation. This number can be reduced dramatically by implementing the blip decomposition, as multi-blip path pairs, along with the trajectories they would generate, do not need to be included. (Note that the role of blips in the QCPI formulation is rather subtle, as the solvent effects enter through a pure phase; however, the phases arising from multi-blip configurations will lead to exponentially small contributions upon averaging with respect to trajectory initial conditions.¹⁰) In addition, it is possible to eliminate the path-specific trajectories by treating the quantum decoherence phase within a harmonic approximation.¹⁰⁰ Preliminary work indicates that such an approximation yields excellent results, while allowing QCPI propagation with a single trajectory per initial condition. The i-BSPI decomposition should also be fully applicable, reducing the cost of iterative QCPI calculations to the combined cost of a MD simulation and an i-BSPI propagation. This work is in progress.

This material is based upon work supported by the National Science Foundation under Award No. CHE 13-62826.