Iterative quantum-classical path integral with dynamically consistent state hopping

As full quantum mechanical calculations require computational effort that scales exponentially with the number of particles, simulations of condensed phase processes often utilize classical trajectories for the majority of degrees of freedom. Classical molecular dynamics simulations with appropriate potential fields are known to yield satisfactory results in many situations. However, classical descriptions are not meaningful in cases of tunneling or nonadiabatic transitions. Fortunately, such strict quantum mechanical effects typically involve just one or only a few degrees of freedom (often represented in terms of discrete quantum mechanical states), suggesting the use of mixed quantum-classical methods. The simplest and oldest quantum-classical method is the Ehrenfest approximation,¹ where each classical trajectory of the environment experiences a force averaged with respect to the wavefunction of the quantum system. Because of this mean field-type assumption, the simple Ehrenfest model is known to lead to incorrect products and branching ratios.² Much effort has been devoted to the development of mixed quantum-classical methods that do not suffer from the shortcoming of the simple Ehrenfest model.^3–8 These methods are powerful and intuitive but include uncontrolled approximations. Recent attempts to develop rigorous quantum-classical methods include the quantum-classical Liouville dynamics formulation^9,10 and the quantum-classical path integral (QCPI).^11–15 A very attractive alternative involves mapping each quantum state on a pair of action-angle variables,^16–19 thus replacing discrete states by continuous degrees of freedom and allowing classical trajectory-based treatments for all degrees of freedom.

In this paper, we present recent developments on the QCPI, a rigorous formulation for calculating the time evolution of a discrete system interacting with a polyatomic environment described in terms of classical trajectories. Unlike common quantum-classical treatments, QCPI is free of approximations, besides the classical trajectory description of the environment. As such, important quantum mechanical interference effects are fully accounted for, leading to correct branching ratios and product distributions. This is achieved by treating the quantum system in terms of Feynman paths,^20,21 which—just like the trajectories that describe the dynamics of the environment—are local, allowing accurate treatment of the system-environment interaction and avoiding the mean field requirement of the original Ehrenfest approximation. The QCPI formulation is rather simple: the solvent experiences a force determined by the state of the system along a particular quantum path, which determines its trajectory. In turn, the solvent supplies a phase to the system, which alters the quantum mechanical amplitude along the given path.

In spite of the exponential proliferation of quantum paths (and, consequently, of classical trajectories) with propagation time, QCPI simulations of realistic systems are possible. The feasibility of such calculations stems from the realization²² that the bulk of the effects of the environment (the classical decoherence¹⁴) can be accounted for by a time-local procedure, which may be chosen to utilize classical trajectories. These effects are built into the single-step system propagator employed in the path integral representation of a dynamical observable,¹³ providing an excellent starting point that allows large time steps, and the remaining quantum decoherence¹⁴ (which is associated with the “back-reaction”¹¹) is captured via phase factors that modify the action of the quantum paths. Further, the exponential growth of classical trajectories is avoided by taking advantage of the memory quenching role of the environment to evaluate the path sum via an iterative procedure¹² similar in spirit to the iterative quasi-adiabatic propagator path integral (i-QuAPI) algorithm.^23,24 The iterative decomposition of the path integral in the case of a general environment of anharmonic, coupled degrees of freedom stems from a cumulant analysis of the influence functional,²⁵ and a general algorithm has been presented.²⁶ Test calculations on model dissipative two-level systems (TLSs) have demonstrated the convergence and favorable scaling of the iterative-QCPI (i-QCPI) methodology.

Still, application of the i-QCPI methodology on many problems of interest requires very demanding calculations. For example, charge transfer processes in solution are associated with a large reorganization energy, which implies very strong interactions between system and solvent particles. Accurate treatment of these strong couplings necessitates small path integral time steps, such that thousands of steps may be required for completion of the reaction. At the same time, sluggish motions lead to memory kernels that can be extremely long-lived, apparently preventing a breakup into short segments. Even with the aid of filtering schemes,^27–30 which eliminate the vast majority of paths whose weights are negligible, i-QCPI calculations can become extremely demanding in these cases. Such regimes can be challenging even to the i-QuAPI methodology, which benefits from the availability of the influence functional in analytic form.³¹ An attractive alternative in the latter case is the blip decomposition,^14,32 which under favorable conditions allows full summation of the Feynman paths (even with hundreds of time steps) without requiring memory truncation. It is clear that additional advances are essential in order to investigate such processes in the absence of a harmonic bath assumption.

The present paper presents a significant advance in this direction. First, classical memory, which is associated with trajectories of the isolated solvent in the absence of forces that depend on the instantaneous state of the system, is incorporated in the system propagators and thus automatically accounted for in the i-QCPI algorithm.¹⁵ Thus, the main challenge is the treatment of the quantum memory, which is a manifestation of the “back-reaction,”¹⁴ i.e., the perturbation of the solvent trajectories by the force exerted by the system along each quantum path. Even though the quantum memory may be extremely long-lived, the effective memory that needs to be treated may be rather short. This is so because quantum memory manifests itself as a phase, and the contributions from its long tail tend to average to zero in the full path sum. We point out that the memory truncation procedure in the original i-QCPI algorithm, which involves randomly selecting a single trajectory branch prior to the onset of full memory treatment,¹² retains some quantum memory along the small subset of the paths that are kept during the period preceding the memory interval, within which such memory is treated explicitly. This treatment implies that the cumulative effects of this memory do not vanish, preventing convergence. To cure this situation, we introduce a modified iterative QCPI algorithm, which is designed to eliminate quantum memory over the propagation interval that precedes the explicit memory length. Further, we introduce a dynamically consistent state hopping (DCSH) procedure that incorporates some of the quantum decoherence effects in the reference propagator. We show that even in the most challenging regimes discussed above, this procedure reduces the need for explicit memory treatment to a small number of time steps, accelerating calculations by orders of magnitude.

Section II gives an overview of the QCPI formulation, with emphasis on the surviving memory effects that are of a quantum mechanical nature. Sections III and IV focus on optimizing trajectory selection in the pre-memory interval, in order to capture as much as possible of this remaining quantum mechanical memory into the system propagators, thus accelerating convergence. Specifically, Section III discusses how spurious memory can be eliminated, in a way analogous to our recent analysis in the special case of a harmonic bath.³³ Section IV presents a probabilistic state hopping procedure for selecting the solvent trajectory branch in accord with instantaneous state populations. Representative calculations on test models are shown in Section V, and some concluding remarks are given in Section VI.

The QCPI formulation combines a fully quantum mechanical path integral description of a few degrees of freedom with a classical trajectory treatment of the environment. As discussed in the Introduction, conventional quantum-classical approaches, which employ a wavefunction description of the quantum degrees of freedom, are inherently tied to Ehrenfest (i.e., mean field) approximations, as a result of the fundamental incompatibility of Newton’s (local) and Schrödinger’s (nonlocal) formulations. No such assumptions need to be made in QCPI because the quantum paths are local in space.

The QCPI algorithm treats the quantum system in terms of a finite number of discrete states in a discrete variable representation (DVR).³⁴ The matrix representation of the system guarantees that the method captures the system’s quantum dynamics exactly and efficiently. The site representation is essential to a position-like path integral description of the time evolution.³⁵ Within this formulation, classical trajectories of the solvent experience a force from the instantaneous position of the system along a particular quantum mechanical path. In turn, each of these trajectories supplies a phase that modifies the quantum mechanical amplitude of the system along the given path. Summing with respect to all discrete paths of the quantum system and integrating over the phase space of the classical variables lead to the exact result.

The Hamiltonian describing the isolated M-state system is written in the form

where $σ n$ represent discrete states or a system-specific DVR,³⁶ i.e., the system position operator is

The total Hamiltonian has the form

where q, p are the coordinates and momenta of the particles comprising the system’s environment (the “solvent”). Observables pertaining to the quantum system are conveniently obtained from the reduced density matrix,

where Δt is the path integral time step. The QCPI expression for the reduced density matrix of the system takes the form

where P is the initial phase space density of the environment and Q is the quantum influence function,^11,12 which contains all the dynamical effects due to the interaction of the solvent with the quantum system, which modify the phase space density of the solvent.

Eq. (2.5) is quite general and allows either a semiclassical^37–39 or a quasiclassical^40,41 treatment of the environment. A forward-backward semiclassical treatment of the solvent dynamics^37–39 leads to separate forward and backward trajectories, while quasiclassical expressions^40,41 employ trajectories that experience the average of the forces exerted by the system at its instantaneous forward and backward configurations. In the quasiclassical form with time-dependent reference propagators,¹³ the quantum influence function is given by the path sum,

where U_ref is the time evolution operator for the time-dependent Hamiltonian,

along the particular trajectory and

with ΔV_sol = H − H_ref. The reference Hamiltonian H_ref augments the system by the external time-varying field V_sol generated by a classical trajectory q^ref, p^ref of the solvent; the latter is propagated subject to the forces of the pure solvent that are independent of the state of the system along each particular system path. In the simplest form, the “classical path approximation,” the solvent trajectories are obtained by solving the classical equations of motion with the system coordinate fixed. Each solvent initial condition gives rise to a different reference Hamiltonian. The corresponding propagator is obtained by solving the time-dependent Schrödinger equation in the basis of system states,

Each sequence of system coordinates $s 0 ± , s 1 ± ,…, s N ±$ in Eq. (2.6) specifies a pair of forward-backward paths of the quantum system. Thus, the quantum influence function involves M^2N such path pairs. A classical trajectory of the solvent experiences a force exerted by the system along each of these forward-backward paths. Specifically, in the forward-backward semiclassical formulation,^37,38 the forward and backward system paths give rise to distinct solvent trajectories, while the quasiclassical treatment^40,41 employs trajectories that experience the average of the forces along the forward and backward paths of the system. For simplicity and higher efficiency, we restrict attention to the quasiclassical version of the QCPI expression, where the phase space density in Eq. (2.5) is given by the Wigner transform of the initial density operator, i.e.,

where n is the number of solvent degrees of freedom. (We note that the full evaluation of the n-dimensional Fourier-type integral in Eq. (2.10) represents an extremely challenging task in itself because the oscillatory phase in the integrand leads to a Monte Carlo sign problem. Approximate methods used for evaluation of the Wigner function are mostly based on Gaussian approximations^41–45 or imaginary time semiclassical propagation.⁴⁶ An approximate but quite accurate trajectory-based method for evaluating the Wigner transform of the density operator has recently become available.⁴⁷ This method, which exploits the classical adiabatic theorem, is entirely trajectory-based, thus simple and ideally suited for use within the QCPI framework.) In the quasiclassical version, the M^2N system path pairs give rise to an equal number of solvent trajectories from each initial condition q₀, p₀.

Eq. (2.6) also contains a phase ΔΦ, which is given by the difference in the net forward-backward action along the forced trajectory and the reference trajectory from the same initial condition. In the limit of a single path integral step (N = 1) and for diagonal endpoints (⁠$s 0 + = s 0 −$⁠, $s N + = s N −$⁠), the forward and backward system paths are identical, thus the net forward-backward phase vanishes and Eq. (2.6) reduces to the ensemble-averaged classical path (EACP) approximation.^11,12 It has been shown^12,14 that in the case of a harmonic environment, the EACP approximation captures all the effects associated with stimulated phonon events (the “classical decoherence”), as well as a portion of spontaneous phonon emission (the “quantum decoherence mechanism”). Corrections to the EACP approximation are related to the “back-reaction,” i.e., the deviations of solvent trajectories from the system-independent reference via forces exerted by the system along each quantum path. These corrections enter through the phase factors in Eq. (2.6) and require evaluation of the full path sum, which contains M^2N terms.

The exponential proliferation of trajectories with propagation time makes the full evaluation of the QCPI expression impractical in most cases. This proliferation is the quantum-classical manifestation of nonlocality in the influence functional from the solvent, which is seen most transparently in the analytical Feynman-Vernon expression³¹ for a harmonic bath. However, this nonlocality is related to (multi-) time correlation functions of the environment²⁵ and thus tends to die out eventually, as correlation functions decay irreversibly in condensed media. These ideas led to the development of iterative path integral methods for system-bath Hamiltonians,^48–52 which were initially developed for harmonic environments and recently extended to fermionic baths;^53–55 they have also been formulated for general, nonlinear environments,²⁶ and recently exploited to devise an iterative QCPI methodology.¹² The latter proceeds by regarding two solvent trajectories as identical if they were obtained from the same system path in the recent past, i.e., over a time interval τ_m equal to the “memory length.” By dropping the branches of the system paths that precede this time interval, one maintains a constant number of trajectories as time progresses. If the memory time spans L path integral time steps, i.e., τ_m = LΔt, the i-QCPI algorithm requires the generation and storage of M^2L trajectories from each initial phase space point. Since even in its zeroth-order (EACP) limit the QCPI formulation accounts for the classical part of the decoherence mechanism in its entirety, the relevant value of τ_m is associated with the quantum mechanical part of the decoherence process, i.e., the net phase in Eq. (2.6), which arises from the forces exerted on the environment by the system along each Feynman path (the so-called “back reaction”).

To fully understand the issues associated with this quantum memory, consider a quantum system and its reduced density matrix at the time NΔt. We refer to the time interval preceding the time (N − L)Δt as the “pre-memory interval” and to the interval comprising the last time length LΔt as the “memory interval.” According to the ideas described in the previous paragraph, the details of the system-solvent interaction prior to the onset of the full path sum treatment are irrelevant; thus, for each initial phase space point, the original i-QCPI algorithm selects and retains randomly one of the (M²)^N−L system path branches within this pre-memory interval. The solvent configuration resulting from that particular randomly chosen system path is then continued within the memory interval, where all memory effects are treated explicitly by summing the QCPI phases over all possible forward/backward path pairs of the quantum system. Figure 1 illustrates this procedure in the case of a TLS. The coordinates of a few representative forward/backward system path pairs during the memory interval are shown in Fig. 1 as solid red/blue line segments, while the coordinates of the single, randomly chosen path pair within the pre-memory interval are shown as dashed red/blue lines.

It is instructive to analyze this situation in light of our understanding of harmonic bath dynamics. In this case, the memory lies in the double sum form of the influence functional, which for a bath initially at equilibrium has the form

where η_k′k″ are the QuAPI-discretized⁵⁶ influence functional coefficients⁴⁹ derived from the bath response function, whose imaginary part is the phase in the QCPI expression. Truncating the memory to L time steps in the i-QuAPI algorithm implies neglecting the contribution of the sum,

This includes a real exponential,

which is associated with the classical decoherence and which is captured in its entirety in the QCPI reference propagators. As a consequence, i-QCPI often converges with smaller values of L than i-QuAPI. (Note, however, that QCPI requires evaluation of individual bath trajectories and the associated phase space integral, which are avoided in QuAPI.)

The contribution to Eq. (3.2) of the (neglected) imaginary part of the response function,

is the quantum contribution to decoherence. This phase is captured only along a single trajectory branch in the i-QCPI algorithm, which in the original version of the method is chosen randomly for each initial phase space point. The lack of a sum over all branches (which would give the exact, full-memory result) tends to exaggerate this phase in i-QCPI, introducing a spurious phase that does not appear in the i-QuAPI algorithm. This spurious phase leads to an apparent quantum memory, which tends to be longer than that required to converge the QuAPI expression. This situation is similar to that encountered in the case of a shifted bath initial condition.³³ To arrive at an i-QCPI algorithm, which will converge with the same quantum memory length as the i-QuAPI scheme, we need to eliminate this spurious memory. The simplest way of achieving this is to set s = 0 in Eq. (3.4), i.e., to evolve the classical trajectories in the pre-memory interval subject to no force from the system (rather than a force dictated by the state of the system along the randomly chosen forward/backward path branch).

In the case of an initially shifted bath, i.e., if the bath is initially in equilibrium with one of the localized system states (the “reactant”), the above analysis suggests that the bath trajectories should be propagated during the pre-memory interval according to the same Hamiltonian (that of the solvent in equilibrium with the reactant state). This simple procedure is intuitive, at least for the early part of the dynamics, since the solvent is still in a state that closely resembles the state in equilibrium with the initial state of the quantum system.

As time progresses, propagation of the classical trajectories on the potential surface associated with the initial state of the system (the “reactant”) gradually becomes less meaningful, as the solvent shifts and eventually equilibrates with other states (e.g., the “product”). Thus, to further accelerate the convergence, we seek to modify the pre-memory treatment to make it consistent with the ensuing dynamics. We specialize to the physically relevant situation where the solvent is initially equilibrated with a localized state of the quantum system. Motivated by the ideas discussed so far, it is desirable for the solvent trajectories to evolve initially on the reactant potential surface, at later times to be distributed between the potential surfaces of the reactant and the product, and eventually (in the case of an asymmetric reaction) to be governed primarily by the state of the product. Below, we describe a probabilistic scheme for achieving this task. The scheme is based on the value of the quantum influence function $Q ( q 0 , p 0 , s k ± = σ i ; k Δ t ) ≡ Q i i ( k Δ t )$⁠, i.e., the value of the reduced density matrix prior to averaging with respect to the phase space of the solvent. Q_ii is real-valued and tends to be distributed primarily between zero and unity, although the tails of this distribution lie outside of this interval.¹¹ 

The DCSH scheme we developed here proceeds as follows: Trajectories are launched and begin to evolve on the potential surface associated with the initial state of the system, the reactant. The value of Q_ii at the end of each time kΔt step serves to determine which of all solvent trajectories should be retained over the time step [kΔt,  (k + 1)Δt] for use at times (k + L)Δt or later. Specifically, when 0 ≤ Q_ii ≤ 1, the potential surface of state i is selected with probability Q_ii. The potential surface associated with state i is always chosen when Q_ii > 1, and never chosen when Q_ii < 0. The trajectory is then integrated for one time step subject to the force from the quantum system in the chosen state. The solvent trajectory q_ref(t), p_ref(t) determined this way for use as the pre-memory trajectory also serves as the reference trajectory q_ref(t), p_ref(t) that is included in the system propagator.

Thus, the reference trajectory hops between the various states of the quantum system. At early times, when the system is mostly in the reactant state, the reference trajectory tends to evolve subject to the forces on the potential energy surface associated with the reactant, thus the DCSH scheme reverts to the reactant branch selection discussed in Sec. III. As the system population is transferred to the product state, the reference trajectory hops back and forth between reactant and product. Once the system has settled on the product and the solvent equilibrates with respect to that state, the reference trajectory has a high probability of evolving on the product potential surface, in accord with the solvent reorganization that accompanies the dynamics.

The above procedure for determining the solvent reference trajectory is somewhat reminiscent of surface hopping algorithms,^3–5 but the details differ in several important ways. First, our state selection criterion is based strictly on populations. Second, no velocity rescaling to satisfy energy conservation is performed in our method. This is so because the reference trajectory is simply one of all the trajectory branches developing under the force exerted by all system paths, and while these trajectories are continuous in phase space, hopping from state to state leads to sudden energy jumps. Third (and perhaps most important), the probabilistic state hopping scheme is used solely for the purpose of optimizing the reference trajectory, and thus accelerating convergence to the full QCPI result, which is free of approximations.

We illustrate the branching reference QCPI methodology with application to two TLS models interacting with harmonic baths. The TLS Hamiltonian has the usual form

where ħΩ = 1 and $1 , − 1$ are localized states with coordinates σ₁ = 1,   σ₂ = − 1, while the solvent is modeled in terms of a dissipative harmonic bath linearly coupled to the TLS,

The system-bath interaction is captured in the spectral density function,

for which we choose the common Ohmic form

The bath is initially in equilibrium with state 1. Numerically, exact results were obtained using the blip sum method^14,32 (which amounts to a full evaluation of the path sum with the Feynman-Vernon influence functional³¹) and its iterative variant.⁵⁷ In addition, path integral Monte Carlo calculations were performed on the asymmetric model to obtain the exact equilibrium value. The i-QCPI calculations employed the same time step and the convergence with increasing memory length was investigated.

For the first model, we choose a symmetric TLS (ε = 0) coupled to a bath with ω_c = Ω, ξ = 2 at a temperature ħ ω_cβ = 1. Accurate calculations employing the blip sum methodology^14,32 were performed with the time step ΩΔt = 0.125. Full-memory QCPI calculations were also performed with ΩΔt = 0.125 for Ωt < 1 and ΩΔt = 0.25 for Ωt > 1. (The use of a larger time step was possible because of the inclusion of the solvent reference in the QCPI propagators.) In Figure 2, we compare the convergence of QCPI calculations with three different choices of the pre-memory trajectory. These calculations converged with ΩΔt = 0.125. Panel (a) shows results obtained with the original treatment,¹² which chooses the pre-memory trajectory branch at random. This scheme shows the slowest convergence with memory length included. In this case, the L = 0 result is the EACP approximation with the solvent trajectory evolving according to the mean of the reactant and product forces (i.e., s = 0). Panel (b) shows results obtained with the pre-memory treatment suggested in Section III, which eliminates the spurious memory by always choosing the reactant (state 1) branch. Here, the L = 0 result corresponds to the EACP approximation with the solvent trajectory evolving subject to the force from the reactant state. Finally, panel (c) shows the results of the DCSH treatment. At intermediate times, these results are seen to differ more from the exact results than those obtained with the reactant branch because the latter fully removes the spurious memory when the system is still in the reactant state. However, the DCSH scheme performs better at long times, converging most rapidly as a function of memory length. The state sequence associated with five representative DCSH trajectories is shown in Figure 3.

For the second model, we choose an asymmetric TLS (ε = 5 ħΩ) coupled to a bath with ω_c = 2Ω, ξ = 4 at a temperature ħ ω_cβ = 0.2. The bath reorganization energy is 16 ħΩ, and this strongly dissipative environment gives rise to a slow decay to the equilibrium value $ρ 11 eq =0.27$⁠. Path integral calculations with the Feynman-Vernon influence functional converged with a time step ΩΔt = 0.075. A combination of full and iterative blip sum calculations with memory L = 100 generated results up to Ωt ≃ 50, i.e., these calculations employed N ≃ 650 path integral time steps.

In Figure 4, we compare the convergence of QCPI calculations on this asymmetric, strongly dissipative TLS with the same three choices of the pre-memory trajectory. These calculations converged with ΩΔt = 0.125. As in the case of the symmetric TLS, the random branch criterion employed in the calculations shown in the left panel leads to the slowest convergence with memory length, and the population is seen to decrease very slowly toward the correct value. The L = 0, mean-surface EACP result leads to equal reactant and product populations (⁠$ρ ̃ 11 =0.5$⁠) at long times. This is consistent with the understanding that the free-bath trajectories capture only the classical decoherence, which coincides with the infinite temperature limit. The reactant (state 1) branch criterion used in the calculations displayed in the middle panel leads to much faster convergence. However, the reactant branch EACP result also predicts equal reactant and product populations at long times. Last, the DCSH treatment again leads to much faster convergence than the other pre-memory branch schemes. The corresponding L = 0 results, obtained from the DCSH trajectories, lead to long-time results that favor the product state. The state sequence associated with five representative DCSH trajectories is shown in Figure 5.

The main computational expense in QCPI calculations stems from the branching of classical trajectories that follow the multiple paths of the quantum system. In the iterative version of the algorithm, each phase space point branches into M^2L trajectories within a memory interval that spans L time steps. To maintain a constant number of trajectories as the propagation progresses, one must retain only one of the trajectory branches prior to the onset of the memory time. A judicious choice of the retained branch is essential for reaching converged results with small values of L.

The first and simplest selection criterion we discussed in Section III aims at eliminating spurious memory. We showed that this is achieved by choosing the branch associated with the reactant state. The spurious memory removal achieved this way parallels the treatment of an initially shifted bath discussed in a recent paper by our group.

As time progresses, it is clear that the branch associated with the reactant is no longer optimal, thus the pre-memory trajectory branch needs to hop between states. We have introduced a probabilistic criterion for these hops, which is based on the instantaneous value of the state population along the particular trajectory. The criterion is designed to select primarily the branch associated with the most populated state and to hop back and forth when multiple states are substantially populated. The chosen pre-memory branch also serves as a physically motivated reference trajectory, whose dynamics is incorporated in the reference propagator. Numerical tests on symmetric and asymmetric two-level systems in contact with strongly dissipative, sluggish harmonic environments show that the DCSH selection of the pre-memory branch and solvent reference leads to a substantial acceleration of convergence as a function of memory length, and thus, to exponential reduction of the number of trajectories. We emphasize that the DCSH criterion is not unique; however, this criterion is used exclusively for accelerating convergence and does not constitute an approximation. It is clear that with sufficiently long memory and small time step, the method converges to the full, rigorous QCPI result, regardless of one’s choice of pre-memory branch and solvent reference.

The DCSH-accelerated QCPI methodology will find application to many interesting processes in chemical and biological systems. We have already obtained converged results in a fully atomistic simulation of the ferrocene-ferrocenium charge transfer in solution.⁵⁸ 

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