Coherent State Mapping Ring-Polymer Molecular Dynamics for Non-Adiabatic quantum propagations
CCoherent State Mapping Ring Polymer Molecular Dynamics forNon-Adiabatic Quantum Propagations
Sutirtha Chowdhury and Pengfei Huo a) Department of Chemistry, University of Rochester, 120 Trustee Road, Rochester, New York 14627,United States
We introduce the coherent state mapping ring polymer molecular dynamics (CS-RPMD), a new method thataccurately describes electronic non-adiabatic dynamics with explicit nuclear quantization. This new approachis derived by using coherent state mapping representation for the electronic degrees of freedom (DOF) and thering-polymer path-integral representation for the nuclear DOF. CS-RPMD Hamiltonian does not contain anyinter-bead coupling term in the state-dependent potential and correctly describes electronic Rabi oscillations.Classical equation of motion is used to sample initial configurations and propagate the trajectories from theCS-RPMD Hamiltonian. At the time equals to zero, the quantum Boltzmann distribution (QBD) is recoveredby reweighting the sampled distribution with an additional phase factor. In a special limit that there is onebead for mapping variables and multiple beads for nuclei, CS-RPMD satisfies detailed balance and preservesan approximate QBD. Numerical tests of this method with a two-state model system show a very goodagreement with exact quantum results over a broad range of electronic couplings.
I. INTRODUCTION
Accurately simulating quantum dynamics effects, in-cluding non-adiabatic electronic transitions and nuclearquantum effects in large-scale condensed phase systems,is one of the central challenges in modern theoreticalchemistry. Direct simulations of the exact quantum dy-namics in these systems remains to be computationallydemanding. It is thus ideal to develop trajectory-basedapproximate methods that scale linearly with respectto the nuclear degrees of freedom (DOF), while at thesame time, accurately describe electronic non-adiabaticdynamics and nuclear quantum effects.Mixed quantum-classical (MQC) and semi-classical(SC) dynamics approaches have already been proved aspromising methods that can accurately describe elec-tronic non-adiabatic transitions. The widely used MQCmethods include Ehrenfest dynamics, surface-hopping(SH) dynamics, and mixed quantum-classical Liouville(MQCL) equation.
The commonly used SC methodsinclude semi-classical initial-value representation (SC-IVR) path-integral methods and linearized path-integral dynamics.
All of these methods use classicaltrajectories to propagate the nuclear DOF, thus signifi-cantly reduce computational cost. However, the classicaldescription of the nuclear dynamics causes inconsisten-cies between quantum and classical mechanics in MQC-based methods, and cannot preserve the quantum ini-tial distribution such as Wigner distribution used in SC-based methods. These deficiencies can lead to problemssuch as the breakdown of detailed balance or zero-point energy leakage.
Imaginary-time path-integral approaches, includingcentroid molecular dynamics (CMD) and ring poly-mer molecular dynamics (RPMD) have been suc- a) Electronic mail: [email protected] cessfully developed and applied to investigate nuclearquantum effects and electronic non-adiabatic dynamicsin large-scale simulations. In particular, RPMD whichresembles classical MD in an extended phase space, pro-vides a convenient approach to compute quantum cor-relation functions and rate constants. In these meth-ods, nuclear quantum statistics are captured with theimaginary-time path-integral formalism, leading to aring-polymer classical isomorphism that describes quan-tum Boltzmann distribution (QBD) in the extended clas-sical phase space. The classical evolution preserves theQBD captured by ring-polymer Hamiltonian due to thesymplectic nature of classical dynamics, and will befree of the zero-point energy leakage problem. Despiteits success in describing quantum effects in condensedphase, RPMD approach is limited to one-electron non-adiabatic dynamics or nuclear quantization, as wellas the lack of the real-time electronic and nuclear coher-ence effects. Recent efforts have been focused on developing RPMDapproaches with electronic-state representation, with avision to accurately describes electronic dynamics andat the same time, preserve QBD. Unfortunately, suchmethods are still missing in the current literature. Forexample, Meanfield RPMD (MF-RPMD) approach preserves QBD; kinetically-constrained RPMD (KC-RPMD) preserves an approximated distribution thatis close to QBD. However, they cannot properly de-scribe electronic coherence because they do not containexplicit electronic state information. Mapping-variableRPMD (MV-RPMD) approach does employ explicitelectronic state variables and preserves the exact QBD,but it cannot accurately capture Rabi oscillations in abare two-state system. On the other hand, mappingCMD, ring polymer surface-hopping (RPSH), ringpolymer Ehrenfest Dynamics, and non-adiabatic map-ping RPMD (NRPMD) are promising methods to pro-vide explicit and accurate electronic dynamics. However,these approaches usually lack rigorous derivations and a r X i v : . [ phy s i c s . c h e m - ph ] N ov they do not preserve detailed balance in general.In this paper, we rigorously derive a state-dependentring-polymer Hamiltonian. Based on that, we developa new RPMD approach, coherent state mapping RPMD(CS-RPMD). Using Meyer-Miller-Stock-Thoss rep-resentation in the coherent state basis, we introduce con-tinuous fictitious phase space variables (mapping vari-ables) to represent the discrete electronic states. Apply-ing the usual path-integral technique, we derive theCS-RPMD partition function expression that providesexact QBD through the extended phase space descrip-tion. Initial distributions in the CS-RPMD are sampledfrom the classical dynamics of the CS-RPMD Hamil-tonian. By using coherent state basis for the mappingvariables, CS-RPMD Hamiltonian does not contain anyinter-bead coupling terms in the state-dependent map-ping potential, leading to an accurate description of theelectronic dynamics and correct Rabi oscillations. In theadiabatic limit and state-independent limit, CS-RPMDreduces back to the regular RPMD, just like any state-dependent RPMD approach, and thus rigorouslypreserves QBD under this limit. In a special case thatthere is only one bead for the mapping variables butstill multiple beads for the nuclear DOF, we can rigor-ously prove that CS-RPMD satisfies detailed balance andpreserves an approximate QBD. While the NRPMD ap-proach assumes a Hamiltonian that closely resembles CS-RPMD Hamiltonian, the current work demonstrates arigorous way to derive this Hamiltonian with a partitionfunction that provides exact QBD, providing a solid the-oretical foundation. II. THEORY
We start with expressing the total Hamiltonian oper-ator of the system as followsˆ H = ˆ T + ˆ V + ˆ H e = ˆ P M + V ( ˆR ) + L (cid:88) n,m =1 V nm ( ˆR ) | n (cid:105)(cid:104) m | , (1)where ˆ T is the nuclear kinetic energy operator, ˆ P is thenuclear momentum operator, M is the nuclear mass. V ( ˆR ) is the state-independent potential operator, andˆ H e = (cid:80) nm V nm ( ˆR ) | n (cid:105)(cid:104) m | is the state-dependent poten-tial operator (electronic part of the Hamiltonian) with L total diabatic electronic states.To derive the CS-RPMD Hamiltonian, we startfrom the canonical partition function defined as Z =Tr en [ e − β ˆ H ], where Tr en = Tr e Tr n represents the traceover both electronic and nuclear DOFs, β = 1 /k B T isthe reciprocal temperature, and ˆ H is the total Hamil-tonian operator defined in Eqn. 1. The partition func-tion can be exactly evaluated as Z = Tr en (cid:81) Nα =1 [ e − β N ˆ H ],with α as the imaginary-time (bead) index, and a highereffective temperature defined as β N = β/N . Fur-ther splitting the Boltzmann operator by trotter expan- sion under the infinite bead limit N → ∞ gives Z =lim N →∞ Tr en (cid:81) Nα =1 [ e − β N ( ˆ T + ˆ V ) e − β N ˆ H e ]. Inserting N copies of the resolution of identity I R = (cid:82) dR α | R α (cid:105)(cid:104) R α | and I P = (cid:82) dP α | P α (cid:105)(cid:104) P α | , and explicitly performing thetrace over the nuclear DOF based on the standard path-integral technique, we have Z = lim N →∞ (cid:90) d { P α } d { R α } e − β N H rp Tr e N (cid:89) α =1 [ e − β N ˆ H e ( R α ) ] , (2)with (cid:82) d { P α } d { R α } = (cid:81) Nα =1 (cid:82) d P α d R α . Here, the ring-polymer Hamiltonian H rp is expressed as follows H rp = N (cid:88) α =1 P α M + V ( R α ) + M β N (cid:126) ( R α − R α − ) , (3)and the state-independent potential operator ˆ H e ( R α ) = (cid:80) n,m V nm ( R α ) | n (cid:105)(cid:104) m | parametrically depends upon α th bead’s nuclear position R α .The above partition function is a common expres-sion and the starting point for all state-dependentRPMD approaches and path-integral Monte-Carlo (PIMC) methods. The only difference amongthese approaches arises from the treatment of the elec-tronic potential term Tr e (cid:81) Nα =1 [ e − β N ˆ H e ( R α ) ]. For exam-ple, in the mean-field RPMD approach, the elec-tronic potential is obtained from a weighted average ofring-polymer in different electronic configurations; in theKC-RPMD approach, the potential is obtained fromthe averaged ring-polymer kink configurations; in theMV-RPMD approach, the electronic states are explic-itly described with mapping variables in the Wignerrepresentation; in the NRPMD approach, the elec-tronic states are described with mapping variables inboth position and momentum bases. II.1. Mapping representation for electronic states
We use Meyer-Miller-Stock-Thoss (MMST) map-ping representation to transform the discrete electronicstates into continuous variables. Based on this rep-resentation, L diabatic electronic states are mappedonto L harmonic oscillators’ ground and first excitedstates with the following relation | n (cid:105) → | ... n ... L (cid:105) =ˆ a † n | ... n ... L (cid:105) . Here, | n (cid:105) is the diabatic state, and | ... n ... L (cid:105) is the singly excited oscillator (SEO) statewith L − n th oscillator in its first excited state. Thus, MMSTformulation provides the following mapping relation | n (cid:105)(cid:104) m | → ˆ a † n ˆ a m , with ˆ a † n = 1 / √ (cid:126) (ˆ q n − iˆ p n ) and ˆ a m =1 / √ (cid:126) (ˆ q m + iˆ p m ) as the creation and annihilation op-erators for harmonic oscillator. With MMST mappingrepresentation, the state-dependent potential operator inEqn. 1 is transformed to (cid:88) n,m V nm ( R α ) | n (cid:105)(cid:104) m | → (cid:88) n,m V nm ( R α )ˆ a † n ˆ a m . (4)Using the above mapping relation, we can rewrite thepartition function as Z = lim N →∞ (cid:90) d { P α } d { R α } e − β N H rp (5) × Tr e N (cid:89) α =1 (cid:104) e − β N (cid:80) nm V nm ( R α )ˆ a † n ˆ a m (cid:105) . To proceed, we need to choose a convenient basis to eval-uate the operators ˆ a † n and ˆ a m inside Z . Recall that coher-ent states | p , q (cid:105) = | p q , ..., p n q n , ...p L q L (cid:105) are the eigen-states of the creation and annihilation operators, withthe following eigen equationsˆ a m | p , q (cid:105) = ( q m + ip m ) √ (cid:126) | p , q (cid:105) ; (cid:104) p , q | ˆ a † n = (cid:104) p , q | ( q n − ip n ) √ (cid:126) , (6)where q ≡ { q , ...q n , ...q L } and p ≡ { p , ...p n , ...p L } .The overlap between the coherent state basis and thediabatic basis can be expressed as follows (cid:104) p , q | n (cid:105) = (cid:104) p , q | ... n ... L (cid:105) = ( q n − ip n ) √ (cid:126) e − ( q T q + p T p ) / (cid:126) (7) (cid:104) m | p , q (cid:105) = (cid:104) ... m ... L | p , q (cid:105) = ( q m + ip m ) √ (cid:126) e − ( q T q + p T p ) / (cid:126) II.2. Derivation of the CS-RPMD Hamiltonian
Expanding the exponential of electronic Hamiltonianoperator in Eqn. 5 up to the linear order of β N , under thelimit that β N →
0, we obtain an equivalent expression asfollows Z = lim N →∞ (cid:90) d { P α } d { R α } e − β N H rp (8) × Tr e N (cid:89) α =1 (cid:34) − β N (cid:88) nm V nm ( R α )ˆ a † n ˆ a m + O ( β N ) (cid:35) To proceed, recall the commutation relationship be-tween the creation and annihilation operators ˆ a † n ˆ a m =ˆ a m ˆ a † n − δ nm . Using this relation, Eqn. 8 becomes Z = lim N →∞ (cid:90) d { P α } d { R α } e − β N H rp (9) × Tr e N (cid:89) α =1 (cid:34) − β N (cid:88) nm V nm ( R α )(ˆ a m ˆ a † n − δ nm ) + O ( β N ) (cid:35) . Now by inserting N copies of the reso-lution of identity for coherent state I p , q =(1 / π (cid:126) ) L (cid:82) d p α d q α | p α , q α (cid:105)(cid:104) p α , q α | in Eqn. 9, andleaving out the higher order terms O ( β N ) under β N → Z ∝ lim N →∞ (cid:90) d { P α } d { R α } e − β N H rp (cid:90) d { p α } d { q α } (10) × Tr e N (cid:89) α =1 (cid:20) | p α q α (cid:105)(cid:104) p α q α | − β N (cid:88) nm V nm ( R α ) (cid:18) ˆ a m | p α q α (cid:105)(cid:104) p α q α | ˆ a † n − δ nm | p α q α (cid:105)(cid:104) p α q α | (cid:19)(cid:21) , with d { p α } d { q α } = (cid:81) Nα =1 d p α d q α Further applying Eqn. 6 to evaluate ˆ a m and ˆ a † n , setting (cid:126) = 1 from now on, and inserting the diabatic projectionoperator P = (cid:80) n | n (cid:105)(cid:104) n | in-between the coherent statebasis to ensure correct projection onto the finite subspaceof SEOs, we obtain the following expression Z ∝ lim N →∞ (cid:90) d { P α } d { R α } e − β N H rp (cid:90) d { p α } d { q α } (11) × N (cid:89) α =1 (cid:104) p α , q α | (cid:88) n | n (cid:105)(cid:104) n | p α + , q α + (cid:105)× (cid:40) − β N (cid:88) nm V nm ( R α ) (cid:20)
12 [ q α + i p α ] m [ q α − i p α ] n − δ nm (cid:21)(cid:41) . Based on a similar derivation procedure developed forthe real-time propagator, now we express the third lineof the above expression back to the full exponential fac-tor, and explicitly evaluate the overlap between the co-herent state basis and the diabatic basis. We arrive at thefinal expression of the coherent state partition function Z cs as the central result of this paper Z cs ∝ lim N →∞ (cid:90) d { P α } d { R α } e − β N H rp (cid:90) d { p α } d { q α } (12) × N (cid:89) α =1
12 ( q α − i p α ) T ( q α + + i p α + ) e − ( q T α q α + p T α p α ) × e − β N (cid:80) nm V nm ( R α ) [ ([ q α ] m [ q α ] n +[ p α ] m [ p α ] n ) − δ nm ] . The above coherent state partition function can bewritten into more compact form as follows Z cs ∝ lim N →∞ (cid:90) d { P α } d { R α } (cid:90) d { p α } d { q α } Γ e − β N H cs , (13)with the weighting factor Γ = N (cid:89) α =1
12 ( q α − i p α ) T ( q α +1 + i p α +1 ) e − ( q T α q α + p T α p α ) , (14)and the CS-RPMD Hamiltonian H cs = N (cid:88) α =1 P α M + V ( R α ) + M β N (cid:126) ( R α − R α − ) (15)+ (cid:88) nm V nm ( R α ) (cid:20)
12 ([ q α ] m [ q α ] n + [ p α ] m [ p α ] n ) − δ nm (cid:21) . The first line in H cs corresponds to the nuclear ring-polymer part of the Hamiltonian, and the second linecorresponds to the non-adiabatic part of the Hamilto-nian which describes electrons-nuclei interactions. Notethat H cs does not contain a potential that couples twoadjacent mapping beads. Similar feature has also beenproposed in the NRPMD Hamiltonian, which ensuresto capture the correct electronic Rabi oscillations. Here,we derived H cs with this feature. To be specific, when theelectronic DOF is decoupled from nuclear DOF, such that V nm ( R α ) = V nm , the non-adiabatic part of H cs becomesMMST Hamiltonian with bead-averaged initial condi-tions, which gives exact frequency for electronic Rabioscillations. Meanwhile, MV-RPMD does containinter-bead coupling for mapping DOF and cannot cap-ture the correct Rabi oscillations. We would like to emphasize several other key fea-tures of CS-RPMD approach. First, CS-RPMD will re-duce back to regular adiabatic RPMD with decoupledelectrons-nuclei limit, including state-independent limit V nm = 0 ( n (cid:54) = m ) and the adiabatic limit βV nm (cid:29) β → it/ (cid:126) , coherent statepartition function Z cs in Eqn. 12 becomes a coherentstate real-time propagator used in the Forward-BackwardMQCL (FB-MQCL) and the Partial-Linearized Den-sity Matrix (PLDM) approach. In FB-MQCL, the Γ term appears as the overlap between two coherent statebases from two consecutive real-time propagators. We further emphasize that the possible mappingRPMD Hamiltonian expression is not unique, for in-stance, we have also obtained the MV-RPMD Hamilto-nian in the coherent state representation as presented inAppendix A. The reason for this non uniqueness is thatthe quantum partition function in Eqn. 2 is an integralof the function of Hamiltonian, thus it is mathematicallypossible to find different Hamiltonian (integrand) thatgives the same quantum partition function (integral).
II.3. CS-RPMD Time Correlation Function
With the derived partition function (Eqn. 13) and H cs (Eqn. 15), we propose to use them to compute the Kubo-transformed time correlation function for operators ˆ A and ˆ B ˜ C AB ( t ) = 1 Z β (cid:90) β Tr[ e − ( β − λ ) ˆ H ˆ Ae − λ ˆ H e i ˆ Ht/ (cid:126) ˆ Be − i ˆ Ht/ (cid:126) ] dλ. (16)Similar to the original RPMD approach, we propose thatthe CS-RPMD correlation function C AB ( t ) = 1 Z cs (cid:90) d { R α } d { P α } d { q α } d { p α } (17) × Γ (0) e − β N ˆ H cs ¯ A (0) ¯ B ( t ) , is an approximate Kubo-transformed time correlationfunction, where ¯ A (0) and ¯ B ( t ) are the bead-averagedestimators for the corresponding operators. Note that Γ (0) is Eqn. 14 evaluated with the initial mapping vari-ables at t=0 , and ¯ B ( t ) is evaluated with the classicaltrajectories generated from the Hamiltonian H cs . In thispaper, we are interested in two types of auto-correlationfunctions: (1) nuclear position auto-correlation func-tion C RR ( t ) where ˆ A = ˆ B = ˆ R , thus ¯ A = ¯ B = N (cid:80) Nα =1 R α , and (2) population auto-correlation func-tion C nn ( t ) where ˆ A = ˆ B = | n (cid:105)(cid:104) n | with the correspond-ing estimator¯ A = ¯ B = 1 N N (cid:88) α =1 [ q α − i p α ] n [ q α +1 + i p α +1 ] n ( q α − i p α ) T ( q α +1 + i p α +1 ) . (18)The CS-RPMD time correlation function can be com-puted by sampling initial configurations with NVT tra-jectories generated from H cs , then propagating the dy-namics with same Hamiltonian H cs to evaluate ¯ B ( t ).Each trajectory is weighted by an initial weighting factor Γ (Eqn. 14). Note that both Γ as well as the projectionoperator estimator are complex. We thus use their com-plex values to accumulate the time correlation function.We obtain results of C AB ( t ) with zero complex valueswithin numerical errors.The novelty of the CS-RPMD formalism is thatthrough the classical evolution of H cs , the initial configu-ration governed by e − β N H cs is preserved for the ensembleof trajectories. Thus, CS-RPMD provides a stable propa-gation scheme for the dynamics and avoids initial config-uration leakage problems. NRPMD approach, on theother hand, uses one Hamiltonian (that closely resemblesMV-RPMD Hamiltonian) to sample the initial configu-rations and then another Hamiltonian (that closely re-sembles CS-RPMD Hamiltonian) to propagate dynamics,and thus might encounter initial configuration leakageproblem. Similarly, in commonly used linearized path-integral approaches (classical Wigner methods), the classical nuclear propagation cannot preserve Wignerinitial distribution, and might cause zero-point energyleakage problem. Just like any imaginary-time path-integral method, CS-RPMD uses classical dynamics inthe extended phase space to quantize nuclei rather thaninitially enforces ZPE through Wigner distribution, thussignificantly alleviating the ZPE leakage problem encoun-tered in the classical Wigner methods. However, CS-RPMD does not preserve QBD in gen-eral, despite the fact that the partition function inEqn. 13 exactly describes the QBD for any quantumstatistics calculations (under the N → ∞ limit for bothnuclear and mapping variables). As can be clearly seenin Eqn. 13, CS-RPMD partition function requires an ad-ditional weighting factor Γ to be multiplied with e − β N H cs in order to recover QBD. This can be easily accomplishedfor quantum statistics calculations. On the other hand,for any time correlation function calculation, the pre-factor Γ (Eqn. 14) is not a constant during the classicalevolution governed by H cs . Thus, CS-RPMD time cor-relation function calculations which only use Γ (0) as an initial weighting factor, does not preserve QBD in gen-eral, especially at the longer time. As a consequence,CS-RPMD in general does not preserve ZPE associatedwith QBD, except at time t = 0, or under the adiabaticor state-independent limit. Despite this deficiency, ournumerical results demonstrate that accurate time cor-relation function can still be obtained. In the state-independent or the adiabatic limit, on the other hand,CS-RPMD reduces back to regular RPMD and preservesQBD, just like any state-dependent RPMD approach. For a special case where there is only one beadfor the mapping DOF (such that all mapping beadscollapse into one), but still multiple beads for thenuclear DOF, Γ = (cid:80) n ( q n + p n ) e − ( q T q + p T p ) = (cid:104) pq | (cid:80) n | n (cid:105)(cid:104) n | pq (cid:105) is indeed the integral of motionof H cs . Thus in this special case, CS-RPMD pre-serves detailed balance, such that C AB ( t ) = C BA ( − t ), with an approximate QBD generated fromone mapping bead and N nuclear beads. Underthis limit, CS-RPMD Hamiltonian in Eqn. 15 becomes H cs = (cid:80) Nα =1 P α M + V ( R α ) + M β N (cid:126) ( R α − R α − ) + (cid:80) nm V nm ( R α ) (cid:2) ( q m q n + p m p n ) − δ nm (cid:3) , which is thering-polymer nuclei under the coherent state MMSTpotential. Note that the projection operator P = (cid:80) n | n (cid:105)(cid:104) n | in the Γ expression constrains the electronicmapping variables within the physical SEO subspace. This gives the partition function Tr[ Γ e − βH cs ] (with Tr ≡ (cid:82) dRdqdp ) rather than simply assuming a classical par-tition function Tr[ e − βH cs ].In addition, CS-RPMD preserves Rabi oscillations aswell as detailed balance under this special one-bead map-ping limit. However, we need to emphasize that by us-ing only one bead for the mapping variables, one cannotfully recover the exact QBD even at t=0. In order toobtain the exact QBD (see Eqn. 10), we explicitly re-quire β N → N → ∞ ) for both electronic andnuclear DOFs in our derivations. Nevertheless, our nu-merical results (provided in Appendix C) suggest thatthe quantum statistics obtained from this approximatedpartition function is close to the exact QBD, and thetime correlation function is close to both exact resultsand the CS-RPMD calculation with multiple beads for allDOFs. We want to further emphasize that this limitingcase has already proven to be useful in recently developedapproaches, such as Ehrenfest RPMD and RPSH which assume one bead for the electronic DOF and mul-tiple beads for the nuclear DOF.Finally, due to the presence of MMST mapping Hamil-tonian in H cs , extra caution is still needed. It is wellknown that this Hamiltonian could exhibit the “invertedpotential” problem when ([ q α ] n + [ p α ] n ) − <
0, aswell as the ZPE leakage problem associated with themapping DOF due to the explicitly incorporated ZPEterm, i.e. , − (cid:80) n V nn ( R α ). These intrinsic deficiencies as-sociated with the MMST Hamiltonian can be addressed and refined with the recent theoretical developments,such as applying semi-classical approximation to treatthe mapping DOF, using ZPE corrections for themapping variables, or using new forms of mappingHamiltonians. II.4. Connections and Differences with Other Methods
CS-RPMD Hamiltonian is derived with the coherentstate basis for mapping variables, and it is closely re-lated to two recently developed state-dependent RPMDformalisms, MV-RPMD and NRPMD. Here we wantto clarify the connections and differences between CS-RPMD and these two methods.First, MV-RPMD Hamiltonian does contain inter-beadcouplings for the mapping DOF, and cannot correctlycapture the electronic Rabi oscillations in a bare two statesystem. In addition, the inter-bead couplings for map-ping DOF in MV-RPMD might cause unphysical oscilla-tion frequency even in the nuclear auto-correlation func-tion calculations, as demonstrated in the result section.We further emphasize that the possible form of the state-dependent Hamiltonian that gives exact QBD for quan-tum statistics calculations is not unique. For instance,we can also obtain the MV-RPMD Hamiltonian in thecoherent state representation, as presented in AppendixA. On the other hand, MV-RPMD approach does pre-serve QBD throughout its dynamical propagation, whichis an appealing feature.Second, the NRPMD approach uses a proposedHamiltonian that closely resembles H cs (Eqn. 15) topropagate dynamics, thus provides accurate electronicRabi oscillations. However, the initial phase space pointsin this method are not sampled from the same Hamil-tonian, but from another one that contains inter-beadcoupling terms in the mapping potential. By doing so,each trajectory in the NRPMD approach might move outof the original phase space distribution due to using twodifferent Hamiltonians. Thus, NRPMD might encounterinitial configuration leakage problem. In the coherentstate basis context, this corresponds to a method thatuses the coherent state MV-RPMD Hamiltonian H cmv (Eqn. 21 in Appendix A) to sample initial configurations,and then uses CS-RPMD Hamiltonian, H cs (Eqn. 15)to propagate trajectories. Besides that, NRPMD alsorequires two estimators for projection operator in or-der to efficiently compute the population correlationfunction. CS-RPMD approach, on the other hand,uses only one estimator (Eqn. 18).Compared to these two recently developed mappingRPMD approaches, the merits of CS-RPMD are (1) pre-serving the electronic Rabi oscillations that MV-RPMDfails to describe, through a rigorously derived H cs thatdoes not contain any inter-bead coupling terms for map-ping DOF, (2) sampling and propagating trajectoriesfrom one derived Hamiltonian, as opposed to NRPMDapproach that uses two Hamiltonians without rigorousderivations, and (3) preserving detailed balance with anapproximate QBD in a special case which contains onlyone mapping bead.However, we must admit that the numerical results ob-tained with CS-RPMD (for the current model systems)are not significantly different than those obtained fromNRPMD. This suggests that both H cs and the Hamilto-nian proposed in NRPMD (which closely resembles H cs )provide accurate dynamics. While the the NRPMD ap-proach simply assumes a Hamiltonian of this form, thecurrent work demonstrates how to rigorously derive thisHamiltonian from a partition function that provides ex-act QBD. III. RESULTS AND DISCUSSION.
To test the accuracy of CS-RPMD, we adapt a com-monly used model system that contains one nuclear co-ordinate and two electronic states ˆ H = ˆ P M + 12 M ω ˆ R + (cid:20) (cid:15) + k ˆ R ∆∆ − (cid:15) − k ˆ R (cid:21) , (19)where ∆ is the electronic coupling, k is the vibronic cou-pling, and 2 (cid:15) is the energy bias between the two diabaticstates. In this paper, we choose a reduced unit systemsuch that M = (cid:126) = 1 and ω = β = k = 1.Table I presents the parameters for all of the modelsystems used in this paper. In particular, Model I and Vare in the adiabatic regime, where ∆ (cid:29) β − ; Model IIand III are in the non-adiabatic regime, where ∆ (cid:28) β − ;Model IV and VI are in the intermediate regime, where∆ ∼ β − . Model III and VI are asymmetric cases withfinite diabatic energy bias 2 (cid:15) , and the rest of the modelsystems are symmetric cases with (cid:15) =0. I II III IV V VI (cid:15)
CS-RPMD correlation functions are computed fromEqn. 17. To evaluate the ensemble average, a total num-ber of 10 initial configurations are sampled from a 2 × a.u. long CS-RPMD NVT trajectory, thermostattedby resampling the nuclear velocities from the Maxwell-Boltzmann distribution at every 0 . / ∆ a.u. Each con-figuration is then further equilibrated with NVE CS-RPMD propagation for another 200 a.u., before beingused to propagate and accumulate the CS-RPMD cor-relation function. All of the correlation functions con-verge at N = 8 or fewer beads. Numerical exact resultsare obtained from discrete variable representation (DVR)calculations. Figure 1 presents nuclear position auto-correlationfunction computed from CS-RPMD (black), the mean- -2.0-1.00.01.02.0 0 5 10 15 20 C RR ( t ) t (a.u.) C RR ( t ) t (a.u.) -2.0-1.00.01.02.0 0 5 10 15 20 C RR ( t ) t (a.u.) (a) (b) (c) (d) C RR ( t ) C RR ( t ) t (a.u.) t (a.u.) -2.0-1.00.01.02.0 0 5 10 15 20 C RR ( t ) t (a.u.) FIG. 1. The Kubo-transformed nuclear position auto-correlation function for model I-IV obtained from CS-RPMD(black solid), the mean-field RPMD (blue dash), and nu-merical exact results (red dot). Results for Model I (sym-metric, adiabatic) are in panel (a), Model II (symmetric,non-adiabatic) are in panel (b), Model III (asymmetric, non-adiabatic) are in panel (c), and Model IV (symmetric, inter-mediate) are in panel (d). field RPMD approach (blue) with expression pro-vided in Appendix B, and numerical exact method (red)for Models I-IV. Model I in Fig. 1a is in the adiabaticregime. In this case, CS-RPMD goes back to the stan-dard RPMD, and agrees with the exact result due to thenear Harmonic adiabatic potential. Mean-field RPMDin this case also gives the same exact result thus notshown here. Model II in Fig. 1b is in the non-adiabaticregime. This is the most challenging case and themost relevant regime for non-adiabatic electron transfer and proton-coupled electron transfer reactions. In thisregime, mean field RPMD starts to break down even at avery short time. CS-RPMD, on the other hand, performsreasonably well compared to exact DVR calculations atthe longer time. Models III and IV are in the interme-diate regime, with asymmetric (Fig. 1c) and symmetric(Fig. 1d) diabatic bias 2 (cid:15) . In this regime, both CS-RPMDand mean field RPMD behave reasonably well.We have also tested the special case when there isonly one bead for mapping DOF, with the results pro-vided in Appendix C. With N = 8 beads for the nu-clear DOF and only one bead for the mapping DOF,we found good agreements for these correlation functionscompared to the results obtained with multiple beads forall DOFs, suggesting a promising practical strategy forpreserving detailed balance as well as providing accuratenon-adiabatic dynamics.Figure 2 presents the comparison between CS-RPMD,NRPMD, and MV-RPMD. In Figure 2(a), Model I isused to illustrate the initial distribution leakage prob-lem that NRPMD encounters. Here we present the re- -2.0-1.00.01.02.0 0 5 10 15 20 C RR ( t ) t (a.u.) -1.00.01.02.0 0 5 10 15 20 C RR ( t ) t (a.u.) C RR ( t ) C RR ( t ) t (a.u.) (a) (b) FIG. 2. Comparison between CS-RPMD with NRPMD and MV-RPMD. (a) The Kubo-transformed nuclear posi-tion auto-correlation function for Model I obtained from CS-RPMD (black), NRPMD (blue), and numerical exact results(red). In this adiabatic case, MV-RPMD agrees with numeri-cal exact result, and thus not shown. (b) Results for Model IIobtained from CS-RPMD (black), MV-RPMD (green), andnumerical exact results (red). In this non-adiabatic case,NRPMD result is close to CS-RPMD, and thus not shown.Magnified plots that correspond to the square regions in bothpanels are provided on the right hand side. sults obtained from NRPMD (blue), CS-RPMD (black)and numerical exact method (red). In this adiabatic testcase, all of the state-dependent RPMD approaches areexpected to reduce to the regular RPMD approach, andreproduce the exact result due to the nearly harmonicadiabatic potential in this model. This is indeed the casefor MF-RPMD , MV-RPMD , and CS-RPMD, all ofwhich can fully recover the exact result under this limit.However, as can be seen from Figure 2(a), the resultfrom NRPMD approach starts to deviate from the exactone at longer times. This is due to the fact that in theNRPMD approach, the Hamiltonian used to propagatetrajectories will not preserve the initial distribution forthe ensemble of trajectories that is sampled from anotherHamiltonian. A similar situation is also encountered inthe coherent state version of NRPMD (result not shown)that samples initial configurations with H cmv (Eqn. 21in Appendix A) and propagates trajectories with H cs (Eqn. 15). CS-RPMD provides accurate long-time dy-namics for this adiabatic test case, by sampling and prop-agating trajectories with the same Hamiltonian H cs andthus preserving the phase space distribution for the en-semble of trajectories governed by H cs throughout thedynamical propagation.In Figure 2(b), Model II is used to provide the com-parison between MV-RPMD and CS-RPMD. Here wepresent the results obtained from MV-RPMD (green),CS-RPMD (black) and numerical exact method (red). As can be seen, the correlation function obtained from MV-RPMD starts to oscillate with a different frequencycompared to the quantum result at a longer time. Thismight happen because the inter-bead couplings for map-ping DOF start to contaminate the physical frequencyof the system. A similar situation is also encountered inthe coherent state MV-RPMD (result not shown) thatsamples and propagates trajectories with H cmv (Eqn. 21in Appendix A). CS-RPMD on the other hand, preservesthe correct oscillation frequency in the correlation func-tion, due to the character of H cs that does not contain theinter-bead coupling for mapping DOF. NRPMD providesa similar result (not shown) compared to CS-RPMDin this non-adiabatic test case. Interestingly, NRPMDand MV-RPMD start to show problematic behaviors atthe opposite limit of the electronic coupling, where CS-RPMD provides reliable results across a broad range ofparameters as demonstrated in Figure 1.Figure 3 presents the nuclear position and the elec-tronic population auto-correlation functions computedfrom CS-RPMD (black) and numerical exact method(red) for models IV-VI. Accurately describing electronicinterference effects (Rabi oscillations) are essential fornon-adiabatic dynamics simulations. Again, CS-RPMDagrees with exact results in the adiabatic regime formodel V presented in Fig. 3a, and provides reasonablygood results for the model systems in the intermediateregimes presented in Fig. 3c-f. NRPMD approach (re-sults not shown) gives nearly identical results for thesecorrelation functions. MV-RPMD on the other hand,cannot correctly capture the electronic oscillations inthese population auto-correlation functions, due to thecontamination of the true electronic Rabi oscillationswith the inter-beads couplings in the mapping ring-polymer Hamiltonian. IV. CONCLUSION.
In this paper, we present CS-RPMD approach, a newstate-dependent RPMD method that can accurately de-scribe electronic non-adiabatic dynamics and nuclearquantum effects. With MMST mapping representation inthe coherent state basis for the electronic DOF and regu-lar path-integral representation for the nuclear DOF, wederive the CS-RPMD Hamiltonian, which closely resem-bles the proposed NRPMD Hamiltonian.
In a spe-cial case where there is only one mapping bead (and stillmultiple nuclear beads), we can rigorously prove thatCS-RPMD preserves detailed balance with an approxi-mate quantum Boltzmann distribution. Numerical re-sults from model systems demonstrate the accuracy ofthis approach across a broad range of electronic couplingregimes.Compared to recently developed state-dependentRPMD approaches , CS-RPMD provides further ap-pealing features, including preserving the electronic Rabioscillations that MV-RPMD fails to describe. In addi- C ( t ) t (a.u.) C ( t ) t (a.u.) -1.000.001.002.00 0 1 2 3 4 5 C RR ( t ) t (a.u.) -1.000.001.002.00 0 1 2 3 4 5 C RR ( t ) t (a.u.) -1.000.001.002.00 0 1 2 3 4 5 C RR ( t ) t (a.u.) t (a.u.) C ( t ) C RR ( t ) (a) (c) (e) (b) (d) (f) C RR ( t ) C RR ( t ) C ( t ) C ( t ) t (a.u.) C ( t ) t (a.u.) FIG. 3. The Kubo-transformed nuclear position and the elec-tronic population auto-correlation functions for model IV-VIobtained from CS-RPMD (black) and numerical exact DVRmethod (red). Results for Model V (symmetric, adiabatic) arein panel (a) and (b), model VI (asymmetric, intermediate) arein panel (c) and (d), and model IV (symmetric, intermediate)are in panel (e) and (f). tion, compared to NRPMD approach which simply as-sumes a Hamiltonian of this form , the current workdemonstrates how to derive this Hamiltonian from a par-tition function that contains exact QBD, providing amore solid theoretical foundation for such methods.The equivalence between RPMD and Kubo-transformed time correlation function has been originallyproposed and recently proved through Matsubaradynamics. We envision a similar rigorous proof of the relation between the CS-RPMD and Kubo-transformed time correlation function, with a coherentstate mapping Matsubara dynamics framework infuture. In addition, we envision to develop a methodthat exactly preserves both the quantum Boltzmanndistribution and the electronic Rabi oscillations, andat the same time, scales linearly with the nuclear DOFand friendly with electronic DOF. Finally, practicaldirections will focus on applying CS-RPMD approachto simulate the coupled electron and proton transferreactions in large-scale condensed phase systems.
Acknowledgement
This work was supported by the University of Rochesterstartup funds. Computing resources were provided bythe Center for Integrated Research Computing (CIRC)at the University of Rochester. We appreciates valuablediscussions with Dr. Tim Hele and Profs. Tom Miller,Jeremy Richardson, and Nandini Ananth. We appreciatethe critical and thorough comments from both reviewers.
Appendix A: Derivation of the coherent-state MV-RPMD Hamiltonian.
Here, we obtain MV-RPMD Hamiltonian in the coherentstate mapping representation. Note that MV-RPMD isnot a new method and has been derived in the Wignerbasis. We start from the partition function expressionin Eqn. 2, perform the trace over the electronic DOFwith coherent state mapping basis, further insert theprojection operator P = (cid:80) n | n (cid:105)(cid:104) n | , and arrive at thefollowing expression Z ∝ lim N →∞ (cid:90) d { P α } d { R α } e − β N H rp (cid:90) d { p α } d { q α }× N (cid:89) α =1 (cid:104) p α , q α | (cid:88) n | n (cid:105)(cid:104) n | e − β N ˆ H e ( R α ) | (cid:88) m | m (cid:105)(cid:104) m | p α + , q α + (cid:105) . The overlap between the diabatic basis and coherentstate basis are explicitly evaluated, and the final expres-sion for coherent state MV-RPMD partition function is Z cmv ∝ lim N →∞ (cid:90) d { R α } d { P α } d { q α } d { p α } sgn(Θ) e − β N H cmv (20)with the coherent state MV-RPMD Hamiltonian H cmv defined as H cmv = N (cid:88) α =1 (cid:18) P α M + V ( R α ) + M β N (cid:126) ( R α − R α + ) (cid:19) + Nβ N (cid:88) α =1
12 ( q T α q α + p T α p α ) − Nβ ln | Θ | . (21)Here, Θ = Re (cid:81) Nα =1 (cid:80) nm [ q α − i p α ] n M nm [ q α +1 + i p α +1 ] m , with M nm ( R α ) = (cid:104) n | e − β N V nm ( R α ) | m (cid:105) . Theseresults have been derived in the Wigner representationof the mapping variables. Appendix B: Meanfield RPMD Hamiltonian.
Here, we derive Mean-field RPMD Hamiltonian withthe coherent state mapping representation. Note thatMF-RPMD is not a new method and has been derivedwithout using mapping representation. Our startingpoint is the CMV-RPMD partition function expressionin Eqn. 20. The mean-field (MF) approximation isapplied by integrating over the electronic mappingvariables to obtain an effective potential for the nuclearDOF. Analytically evaluating the Gaussian integral over d { q α } d { p α } in Eqn. 20, we arrived at the MF-RPMDpartition function expression Z MF ∝ lim N →∞ (cid:90) d { P α } d { R α } e − β N H MF sgn(Θ (cid:48) ) , (22)where H MF = H rp − Nβ ln | Θ (cid:48) | with Θ (cid:48) =Tr e (cid:104)(cid:81) Nα =1 M ( R α ) (cid:105) and H rp is the regular ring-polymer Hamiltonian defined in Eqn. 3. Appendix C: Results for one mapping beadand multiple nuclear beads. -2.0-1.00.01.02.0 0 5 10 15 20 C RR ( t ) t (a.u.) -2.0-1.00.01.02.0 0 5 10 15 20 C RR ( t ) t (a.u.) C RR ( t ) t (a.u.) -2.0-1.00.01.02.0 0 5 10 15 20 C RR ( t ) t (a.u.) (a)(c) C RR ( t ) C RR ( t ) t((a.u.) C RR ( t ) (b)t((a.u.) (d) C RR ( t ) -1.00.01.02.0 0 1 2 3 4 5 C RR ( t ) t (a.u.) C RR ( t ) C ( t ) t (a.u.) C ( t ) t((a.u.) t((a.u.)(e) (f) FIG. 4. CS-RPMD correlation functions obtained withone mapping bead and eight nuclear beads. The Kubo-transformed nuclear position auto-correlation functions formodel I-IV are presented in (a)-(d), and the Kubo-transformed nuclear position and population auto-correlationfunctions for model V are presented in (e)-(f). Results areobtained from CS-RPMD (black) and numerical exact DVRmethod (red). For population auto-correlation function pre-sented in (f), CS-RPMD results are obtained with one map-ping bead and 1 (blue line), 8 (black line), and 16 (green line)nuclear beads.
Here, we present the CS-RPMD results with one map-ping bead and eight nuclear beads. Under this speciallimit, CS-RPMD rigorously preserves detailed balanceas Γ becomes an integral of motion of H cs . However,CS-RPMD is not able to recover the exact QBD evenat t=0, due to the fact that there is only one bead formapping DOF and Eqn. 10 becomes a rough approxima-tion. As can be clearly seen, the CS-RPMD correlationfunctions start to deviate from the exact result evenat t=0. Nevertheless, the numerical results of thesecorrelation functions show good agreement to thosepresented in the main text, suggesting that the nuclearquantization is the most important factor for achievingan accurate result in these test cases. Thus for thecalculations presented here, we expect that the MMSTHamiltonian (which appears under this one mappingbead limit) is accurate for describing non-adiabaticeffects, and nuclear quantization with ring-polymershould be able to capture the nuclear quantum effects,leading to a reasonable numerical accuracy. S. C. Althorpe, N. Ananth, G. Angulo, R. D. Astumian, V. Beni- wal, J. Blumberger, P. G. Bolhuis, B. Ensing, D. R. Glowacki,S. Habershon, S. Hammes-Schiffer, T. J. H. Hele, N. Makri, D. E.Manolopoulos, L. K. McKemmish, T. F. M. III, W. H. Miller,A. J. M. andTatiana Nekipelova, E. Pollak, J. O. Richardson,M. Richter, P. R. Chowdhury, D. Shalashilin, and R. Szabla,Faraday Discuss , 311 (2016). J. Tully, J. Chem. Phys. , 1061 (1990). J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, andN. Bellonzi, Annu. Rev. Phys. Chem. , 387 (2016). L. Wang, A. Akimov, and O. V. Prezhdo, J. Phys. Chem. Lett , 2100 (2016). D. M. Kernan, G. Ciccotti, and R. Kapral, J. Phys. Chem. B. , 424 (2008). H. Kim, A. Nassimi, and R. Kapral, J. Chem. Phys , 084102(2008). C. Hsieh and R. Kapral, J. Chem. Phys. , 22A507 (2012). C. Hsieh and R. Kapral, J. Chem. Phys. , 134110 (2013). R. Kapral, J. Chem. Phys. , 134110 (2013). W.H.Miller, J. Phys. Chem. A. , 2942 (2001). W.H.Miller, J. Phys. Chem. B. , 1405 (2009). X. Sun, H. Wang, and W. H. Miller, J. Chem. Phys (1998). Q. Shi and E. Geva, J. Phys. Chem. A. (2004). P. Huo and D. F. Coker, J. Chem. Phys. , 201101 (2011). J. C. Tully, J. Phys. Chem. Lett. , 22A301 (2012). J. Liu and W. H. Miller, J. Phys.: Condens. Matter , 073201(2015). P. V. Parandekar and J. C. Tully, J. Chem. Theo. Comp. , 229(2006). J. R. Schmidt, P. V. Parandekar, and J. C. Tully, J. Chem. Phys. , 044104 (2008). U. Muller and G. Stock, J. Chem. Phys. , 77 (1999). S. Habershon and D. E. Manolopoulos, J. Chem. Phys. ,244518 (2009). J. Cao and G. Voth, J. Chem. Phys. , 5106 (1994). Jang and G. Voth, J. Chem. Phys. , 2371 (1999). S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F.Miller, Annu. Rev. Phys. Chem. , 124105 (2013). I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. , 3368(2004). A. R. Menzeleev, N. Ananth, and T. F. Miller, J. Chem. Phys. , 074106 (2011). T. J. H. Hele, An electronically non-adiabatic generalization ofring polymer molecular dynamics, MChem thesis, Exeter College,University of Oxford (2011). J. R. Duke and N. Ananth, Faraday Discuss , 253 (2016). A. R. Menzeleev, F. Bell, and T. F. Miller, J. Chem. Phys. ,064103 (2014). N. Ananth, J. Chem. Phys. , 124102 (2013). J.-L. Liao and G. A. Voth, J. Phys. Chem. B. , 8449 (2002). P. Shushkov, R. Li, and J. C. Tully, J. Chem. Phys. , 22A549(2012). F. A. Shakib and P. Huo, J. Phys. Chem. Lett. , 3073 (2017). T. Yoshikawa and T. Takayanagi, Chem. Phys. Letts. , 1(2013). J. O. Richardson and M. Thoss, J. Chem. Phys. , 031102(2013). H. D. Meyer and W. H. Miller, J. Chem. Phys , 3214 (1979). G. Stock and M. Thoss, Phys. Rev. Lett. , 578 (1997). U. M¨uller and G. Stock, J. Chem. Phys. , 77 (1999). R. P. Feynman and A. R. Hibbs, Quantum Mechanics and PathIntegrals, Dover Publications, INC. (1965). B. J. Berne and D. Thirumalai, Annu. Rev. Phys. Chem. , 401(1986). D. M. Ceperley, Rev. Mod. Phys. , 279 (1995). D. Chandler and P. G. Wolynes, J. Chem. Phys. , 4078 (1981). M. H. Alexander, Chem. Phys. Lett. , 436 (2001). J. R. Schmidt and J. C. Tully, J. Chem. Phys. , 094103(2007). N. Ananth and T. F. Miller, J. Chem. Phys , 234103 (2010). J. Lu and Z. Zhou, J. Chem. Phys , 154110 (2017). T. J. H. Hele and N. Ananth, Faraday Discuss. , 269 (2016). P. Huo, T. F. Miller, and D. F. Coker, J. Chem. Phys. ,151103 (2013). T. J. H. Hele, M. J. Willatt, A. Muolo, and S. C. Althorpe, J.Chem. Phys. , 134103 (2015). S. Bonella and D. Coker, J. Chem. Phys. , 7778 (2001). S. Bonella and D. Coker, J. Chem. Phys. , 4370 (2003). W. H. Miller and S. J. Cotton, Faraday Discuss. , 9 (2016). S. J. Cotton and W. H. Miller, J. Phys. Chem. A. , 12138(2015). J. Liu, J. Chem. Phys. , 204105 (2016). J. O.Richardson, P. Meyer, M.-O. Pleinert, and M. Thoss, Chem.Phys. , 124 (2017). D. T. Colbert and W. H. Miller, J. Chem. Phys. , 1982 (1992). J. S. Kretchmer and T. F. Miller, J. Chem. Phys. , 134109(2013). B. J. Braams and D. E. Manolopoulos, J. Chem. Phys. ,124105 (2006). T. J. H. Hele, M. J. Willatt, A. Muolo, and S. C. Althorpe, J.Chem. Phys.142