The effect of sampling techniques used in the multiconfigurational Ehrenfest method

In this paper, we compare and contrast basis set sampling techniques recently developed for use in the ab initio multiple cloning method, a direct dynamics extension to the multiconﬁgurational Ehrenfest approach, used recently for the quantum simulation of ultrafast photochemistry. We demonstrate that simultaneous use of basis set cloning and basis function trains can produce results which are converged to the exact quantum result. To demonstrate this, we employ these sampling methods in simulations of quantum dynamics in the spin boson model with a broad range of parameters and compare the results to accurate benchmarks. Published by AIP Publishing. https://doi.org/10.1063/1.5020567


I. INTRODUCTION
Direct Dynamics (DD) methods have become a valuable tool for computational chemistry. They use trajectories to guide the motion of nuclei and quantum ab initio electronic structure methods to calculate potential energies of electronic states, forces, and nonadiabatic coupling matrix elements. The advantage of DD methods is that they can be used without any fit to the potential energy surface (PES) and without any preconditions. As a result, new regions of strong nonadiabatic coupling, for example, which determine new mechanisms of photochemical reactions, can be found from DD numerical experiments. Recently Quantum Direct Dynamics (QDD) techniques have been developed, which quantum mechanically treat not only the electrons but also the nuclei. QDD methods include Ab Initio Multiple Spawning (AIMS), 1,2 Variational Multiconfigurational Gaussians (vMCG), 3,4 and Ab Initio Multiple Cloning (AIMC), an ab initio DD version of the multiconfigurational Ehrenfest method. 5 Quantum DD methods rely not on a single trajectory but on an ensemble of trajectories which serve to guide a basis set of Gaussian coherent states required for the nuclear wave packet propagation. In principle, quantum direct dynamics methods represent a "chemical theory of everything," describing chemical dynamics solely with the help of the time dependent Schrödinger equation; however, in practice, it is difficult to achieve good convergence for systems with many nuclear degrees of freedom. Convergence to the exact quantum result relies on having a good sampling of the Gaussian coherent state basis. Recently we have developed a number of such sampling techniques within the framework of the Ab Initio Multiple Cloning-Multiconfigurational Ehrenfest (AIMC-MCE) approach. 5 The Multiconfigurational Ehrenfest (MCE) method has been shown to be capable of simulating various systems, from calculating the Franck-Condon spectrum of pyrazine 6 to simulating the photodissociation of pyrrole. 7 When using this method, the wavefunction is projected onto a basis of nuclear coherent states |z k and orthogonal electronic states |ϕ r , coupled through a set of amplitudes. The use of coherent states allows the trajectories to be calculated easily as the classical energy of a point in phase space is always known. The propagation equations for the wavefunction parameters are found by use of the variational method and by substitution into the time dependent Schrödinger equation. As such, MCE represents a fully quantum method with no approximations other than the use of a finite basis set which can, in principle, converge to the exact result. When first introduced, the MCE method was tested using the spin boson model, 8 a paradigmatic model which describes two coupled high dimensional potential energy surfaces, which correspond to two different electronic states. The first version of the multiconfigurational Ehrenfest method introduced in the original paper, 8 which will be referred as MCEv1, used an ansatz to describe the wavefunction given by where r are the indices of the electronic states |ϕ r of the system, |z k are the nuclear coherent states which are related to the positions and momenta of the nuclear trajectories, and the a rk are the quantum amplitudes of electron-nuclear basis states |ϕ r z k . In MCEv1, φ k = r (a rk |ϕ r )|z k represents an Ehrenfest configuration, which is the key building block of MCE, and includes the Ehrenfest trajectories |z k (t) and the electronic state amplitudes a rk which are coupled by quantum close coupling equations. As a result of the coupling of these amplitudes in MCEv1, the Ehrenfest trajectories are not independent in this approach. It should also be noted that in MCEv1, normalization is carried out only on the level of the whole wavefunction (1), not on the level of the individual Ehrenfest configurations.
Later another version of the multiconfigurational Ehrenfest method, referred to as MCEv2, was suggested. MCEv2 uses a different ansatz, 6 given by where an additional amplitude D k is introduced. Unlike with MCEv1, in MCEv2, the amplitudes a rk of different coherent states (i.e., those which have different index k) are not coupled. This has the benefit of each Ehrenfest trajectory |z k (t) being independent of other trajectories. The equations for the amplitudes D k then serve to couple the independent basis Ehrenfest configurations φ k . With MCEv2 being based on independent Ehrenfest trajectories, it is much more suitable for ab initio quantum direct dynamics than MCEv1 as the trajectories of basis configurations can be run one-by-one thus representing a significant improvement to the method. MCEv2 has been implemented and used for simulating the dynamics of small organic molecules. 7,9 It is however known that independent Ehrenfest trajectories do not always reproduce the splitting of the wavefunction between two electronic states and therefore they can misguide the basis set. This can become particularly important for dynamics beyond very short time scales.
To address this problem and to adapt the basis to the dynamics better, the multiple cloning procedure and the use of train basis sets were introduced. 10,11 Despite their earlier use in the simulation of small organic molecules, these sampling techniques have not yet been tested against an exact benchmark. In the paper in which MCEv2 was introduced, 6 the MCEv2 approach was tested against benchmark calculations of a 24D model of pyrazine, yielding a well converged spectrum. The dynamics of pyrazine considered however were on very short time scales. In this paper, therefore, we use the spin boson model to gain a better understanding of convergence with the MCEv2 method. The parameters of the spin boson model can be varied to describe several regimes of nonadiabatic coupling and intramolecular energy exchange. Some cases require longer time dynamics than was needed for pyrazine and, as will be shown in this paper, need more sophisticated sampling techniques. Unlike ab initio DD methods where the vast majority of time is spent on electronic structure calculations, the spin boson model relies on a simple analytical form of the potential and the convergence of MCE quantum dynamics to the exact quantum benchmark can be analysed relatively easily. We will also compare MCEv2 with MCEv1, showing that it is a lot harder to converge MCEv2 than MCEv1 and therefore that MCEv2 is not a viable technique for treating analytical models like the spin boson model, where MCEv1 should be used. It will be demonstrated however that the sampling techniques used in the domain of ab initio direct dynamics allow accurate and well converged results to be obtained for the spin boson model with a broad range of parameters using the MCEv2 method, providing additional validation of our ab initio direct dynamics approach. The main cost of our direct dynamics approach comes from ab initio electronic structure calculations, and the cost of MCEv2 simulations, while high in comparison with MCEv1, is still negligible in comparison with that of the electronic structure. As cloning, the main sampling technique used with MCEv2, is very similar to the procedure used in multiple spawning and train basis sets (also called time displaced basis sets) can be used in AIMS 2 the results of this paper can perhaps be transferred to AIMS as well.

A. The multi-configurational Ehrenfest method
The wavefunction for the MCE method is propagated through the equations of motion for the amplitudes and the centres of the Gaussians which are found by the variational method. 12 A full derivation of the equations of motion is given elsewhere (see, for example, Refs. 6 and 8), and so only a brief outline will be given here. Throughout this section, indices will be represented such that j, k refer to configurations, m will refer to degrees of freedom, r, s will refer to potential energy surfaces, and n will refer to the total number of potential energy surfaces. Furthermore, it should be assumed that = 1.

The coherent state equations
In this paper, we will use the |z notation to denote the Gaussian coherent state which in coordinate representation is given by a Gaussian wave packet, where It is characterised by its position and momentum and in chemistry is often denoted as |q, p (further information about the correspondence between the two systems of notation can be found in Ref. 13). It is a defining property of coherent states that the |z state is an eigenstate of the annihilation operator with eigenvalues z. In the representation used here, the creation and annihilation operators can be given bŷ and so for a single dimension which allows z and z * to be used as the variables of propagation in place of q and p. The above equations are the one-dimensional versions, representing a point in phase space, and can be easily transformed to the multidimensional versions by way of for a system of M dimensions. An ensemble of such coherent states can serve as a basis set for the description of nuclear

Determination of the initial conditions of the system
Typically it is assumed that the wavepacket starts entirely on a given electronic state. For MCEv2, as the electronic states of the quantum system are orthonormal, the initial values for the amplitudes a rj in the MCEv2 anzatz (2) are simply a rj = 1 on the initial PES and a rj = 0 otherwise. The initial values for the amplitudes D j arise from the identity formula where Ω −1 jk is the inverse of the overlap matrix Ω jk = φ j (t) φ k (t) . By applying this identity to the initial wavefunction |Ψ(0) , the values for D j can be found from the resulting system of linear equations. The initial conditions for MCEv1 are identical to those for MCEv2 if the amplitudes are combined, and as such Eq. (10) can be used to find the initial value for a rj in Eq. (1) on the initial PES, while the amplitude remains equal to zero on any other PESs.

Time propagation of the wavefunction
In the MCEv2 method, the orthonormal number states make up the quantum (electronic) system states and the bath (nuclear) states are supplied by the set of |z k (t) in the ansatz. The wavefunction therefore is described by a linear combination of several configurations φ k (t) , and the propagation of the wavefunction is carried out by way of the simultaneous propagation of |z k (t) , a rk (t), and D k (t). The time dependence for |z k (t) and the a rk (t) amplitudes are found through applying a variational principle to the single configuration Lagrangian of the system, and the time dependence for D k (t) is given through application to the time dependent Schrödinger equation. Through this, and by formulating the amplitude of the electronic states in terms of the action and a smooth pre-exponential factor such that it can be found that for the MCEv2 methoḋ where the classical action is given by The time evolution of the coherent basis is carried out through applying the variational principle for z * , resulting in Hamilton's equations in the z-notation, where (15) Derivation of the time evolution of the quantum amplitude D k yields the equation where and As such, Eqs. (12)- (14) and (16) make up the propagation scheme for the wavefunction in the MCEv2 method.
The equations for the MCEv1 method differ from these equations as the cross-configurational coupling is in the d rk equations of motion rather than D k . This results in the time propagation being determined by Eqs. (13) and (14) and the MCEv1 equation forḋ rk , given by where the δ 2 H (r,r) term is given by

B. Basis function sampling and cloning
Through applying different sampling techniques in the MCE method, convergence can be greatly improved, and so proper construction of the basis set is extremely important. The most straightforward way to construct the basis of coherent state-based Ehrenfest configurations is to create a "swarm" of coherent states covering the initial wavefunction 14 which is constructed from a Gaussian distribution biased to the initial wavefunction |Ψ 0 = |z 0 , where α c is the compression parameter, an inverse of the width of the Gaussian swarm. This parameter determines the area of phase space covered by the initial wavefunction, with higher compression parameters resulting in a denser basis set. The values of α c are dependent on the size and dimensionality of the basis set, with larger basis sets requiring compression parameters in the region of a few hundred up to a couple of thousand, and smaller basis sets with lower dimensionality requiring compression parameters in the range of tens up to a couple of hundreds. The compression parameter is determined automatically so as to best describe the wavefunction with a finite number of basis functions through examination of the initial norm and recursive adjustment of α c and resampling of the basis set until a value of α c is found which gives a norm close enough to unity. While this swarm is suitable in many situations, there have been various improvements of the sampling in the MCE method which can be used. The first involves using coherent state "trains" for initial sampling 14 which applies a "smoothing" to the propagation of the wavefunction. Another improvement involves basis set cloning 10 which grows the basis set when intersections are encountered. Also an operator or the wavefunction itself can be split into a superposition of sets of Gaussian coherent states so that each of these sets can be propagated separately, termed "bit-by-bit" propagation. In this section, all three of these options are described and discussed briefly, with more comprehensive descriptions in Refs. 10, 11, and 14.

Use of basis function trains to improve MCEv2
Coherent state trains, also known as time displaced basis sets, 2 were used in the context of coupled coherent states in 2008 14 as a way of inserting some "regularity" into a random swarm. The argument is that a random swarm, while improving scalability, necessitates a sacrifice in convergence. At the other end of the scale, a regular grid allows extremely fast convergence but scales exponentially, resulting in high numerical expense for all but the smallest of systems. The ideal compromise lies somewhere between these two extremes, with not total regularity but not a true random swarm either.
In a coherent state train, a small compressed random swarm is generated, and for each of the basis functions in this swarm a line of basis functions is formed in phase space along the path of propagation. This allows the basis set to cover a larger area in phase space than is covered by just a compressed random swarm. Due to the structure of the trains, the process of constructing the initial basis set is somewhat different to that described earlier. As with the construction of a random swarm, the initial wave packet |z 0 is calculated first, and a small random swarm is constructed around this initial wave packet using Eq. (22). Following this, a parameter δ trn is set which is defined to be the time displacement between two adjacent basis function "carriages." The swarm is then propagated backwards and forwards in time, while the configuration is saved after every t = δ trn such that a set of coherent states with single configuration amplitudes are obtained, all of which follow the same path in phase space, hence the term "trains." The set of cross configuration amplitudes D k is then calculated over all the single configuration basis function "carriages" and so the wavefunction is spread out over the length of the train. This procedure differs from that reported in Refs. 10 and 11 as the simulation of this model system does not benefit from the separation of the calculation of the trajectories from the calculation of the cross configuration amplitudes and so the entire wavefunction is propagated as a whole throughout.

Use of basis function cloning to improve MCEv2
Multiple cloning is a recent inclusion to the MCE method, having been introduced in a coupled coherent state context in two papers published by Makhov et al. 10,11 It takes inspiration from the multiple spawning method of Martínez and Ben-Nun 15-21 and may be viewed simply as a straightforward and convenient way of performing spawning.
Like other methods based on Ehrenfest dynamics, the multiple cloning MCE method propagates the wavefunction not along the potential energy surface classically but on a quantum average of the potential. It is a problem inherent to all methods based on Ehrenfest dynamics however that in a region of nonadiabatic coupling where the population of an Ehrenfest trajectory is split almost equally across multiple potential energy surfaces with different forces, this average is not a faithful representation of the system, propagating the wavefunction subject to a force which is an average of different forces on each electronic state. To remedy this, when a basis function meets these conditions, that basis function is cloned, with one instance projected onto the first potential energy surface and the other projected onto the second potential energy surface (assuming a two-state system).
If, before a cloning event, a single basis function for a two-state system is given by then after cloning two basis functions will exist, given by The determination of when to clone a basis set is dependent upon the force between the potential energy surfaces, as given by In the spin boson model, the differential of the potential is a constant, and as such the maximum of the breaking force can be determined to be where the single configurational amplitudes are equal for both potential energy surfaces. As such, an appropriate condition for cloning would be when |a 1k a 2k | 2 > 0.249. A further necessary condition would be the limiting of cloning events on the same configuration within an appropriate number of time steps, allowing the basis function to move away from the intersection of the two potential energy surfaces, thus preventing multiple cloning events being applied to a single basis function due to the same intersection. As the basis set increases in size, this means that the wavefunction can be better described in phase space, and in addition to this the wavefunction no longer becomes illdefined in the region immediately after passing through an intersection.

Bit-by-bit propagation of an operator
A final sampling procedure used is that of the so-called "bit-by-bit" propagation, whereby the effect of using a very large basis set can be approximated by running a large number of repeat propagations with smaller basis sets with differing initial conditions, thus splitting the propagation into smaller tasks. This is implemented for the spin boson model through the values for the initial coherent state |z 0 , which are found from the density operator of the bath coherent states. This operator is given as a product of 1D density operators such that where the 1D density operator is given bŷ In practice, this means that the M values for the z (m) 0 coordinates are sampled from a normally distributed random swarm centred at (q, p) = (0, 0) with width σ (m) , where where the thermal parameter β = 1/(k B T ). Through this procedure, the populations of system states can be obtained by a simple averaging over a number of repeat propagations N rpt with different initial coherent states |z 0 . 8,22

C. The spin boson model
The spin boson model 23 is a paradigmatic physical model which at its most basic consists of a two-state (spin 1/2) system linearly coupled to a bosonic bath, and is the most simple model to describe the effect of an environment on constructive and destructive quantum interference, also allowing the investigation of decoherence and dampening on the quantum system. 24 In the spin boson model, the two-state system (with diabatic donor and acceptor states |ϕ 1 and |ϕ 2 ) and harmonic bath use the Hamiltonian where H B and H C are the bath and coupling Hamiltonians, respectively, and where the bias detuning parameter and the tunneling amplitude between states ∆ can both be taken to be constant. It can be reasonably assumed that the latter of these parameters is approximately independent of the vibrational degrees of freedom. 25 The partial Hamiltonians H B and H C can be expressed in terms of the creation and annihilation operators 8,26 as and by consideration of Eq. (7) the bath and coupling Hamiltonians for the spin boson model become and (32) In the above equations, the strength of the coupling between the bath and the two-state system is given by the parameter C (m) . Information about the harmonic bath is encapsulated in the spectral density which can be chosen to model various physical systems such as a solvent, phonons of a solid, or other condensed phase environments. 27 A widely used special case is that of the Ohmic spectral density with an exponential cutoff which has a characteristic low-frequency behavior J(ω) ∝ ω and peaks at a cutoff frequency ω c , which defines the time-scale distribution of the bath dynamics, 28 such that where α K is the Kondo parameter. In this case, the quantum system is damped equally at all frequencies, which is the case for many physical systems. 24 An important consideration when applying trajectorybased methods to the spin boson model is the scheme used to discretise the harmonic bath. In this case, the continuous bath spectral density given in Eq. (34) is discretised to the form of Eq. (33) by way of the relation 29 where ρ(ω) is a density of frequencies satisfying It is determined in Ref. 29 and elsewhere that the precise functional form of ρ(ω) does not affect the final answer, provided a large enough total number of bath modes M is used; however, it can affect the total number of bath modes needed to correctly represent the continuum. In this case, ρ(ω) is taken to be where where ω max is the largest frequency of the bath modes considered, taken to be ω max = 5ω c . Using this discretisation, the equation for the coupling coefficient between the system and the bath C (m) can be found as and the frequencies of the bath can be found as The variation of the parameters ω c , α K , and the thermal parameter β therefore allows a range of different systems to be modelled, each requiring different numbers of degrees of freedom and basis set sizes.

III. RESULTS
A. Using trains, cloning, and bit-by-bit propagation to converge MCEv2

Comparison of the MCEv1 and MCEv2 methods
While the MCEv1 method has shown itself to be capable of simulating high-dimensional model systems and the MCEv2 method has been successful at simulating small organic molecules, 6,7,9,30 there has never been a test on a model system with high dimensionality comparing the two methods on an even footing.
To this end, such a comparison was made using the spin boson model looking initially at a pair of low temperature symmetric wells [ Fig. 1(a)] and at a pair of low temperature asymmetric wells [ Fig. 1(b)] with initial basis sets sampled from compressed random swarms with average compression parameters of α c = 209 for the symmetric wells and α c = 69 for the asymmetric wells. The results are also compared to exact benchmark data calculated using the Multiconfiguration Time Dependent Hartree (MCTDH) method. Using the tunneling energy ∆ as the unit of energy, the symmetric well had parameters ω c /∆ = 2.5, α K = 0.09, β∆ = 5.0, and /∆ = 0 with It can be immediately seen that there is a great disparity between the two results, with the oscillations for the MCEv2 result overemphasised in both cases, while the MCEv1 result follows the MCTDH benchmark almost exactly. Furthermore, it was found that a further increase in the size of the basis set or the number of repetitions did not improve the agreement for the MCEv2 result. When discovered, this result was very surprising as numerically the wavefunctions should be identical and it was thought that the dynamics of the system should not have been affected by the MCEv2 equations.
It is a cause of further surprise that if the simulation is run using a large number of repeat calculations of single Ehrenfest configurations, thus effectively running a simulation using an ensemble of uncoupled basis functions guided by Ehrenfest dynamics, it can be seen that for the symmetric case the MCEv2 population difference is identical [ Fig. 2(a)] and for the asymmetric case the MCEv2 population difference matches for the first few oscillations [ Fig. 2(b)]. This would seem to indicate either that the basis set is behaving as an ensemble of independent non-interacting basis functions due to a loss of coupling or that the basis is behaving as a set of basis functions guided by trajectories which are too similar to each other, thus behaving almost as a single larger basis function. The loss of coupling between the basis functions would, in most cases, indicate the loss of overlap between the coherent states; however, if the guiding trajectories are too similar, this would most likely result in a higher overlap. As such, a comparison of the normalised absolute average overlap, i.e., an average over the absolute values of all elements of the overlap matrix Ω jk , was prepared. As is seen in Fig. 3, the overlap for the MCEv2 simulation is higher than that of the MCEv1 simulation, decaying much slower and not as smoothly. This means that the basis set does not spread to cover as much of an area of phase space when using the MCEv2 equations. It can be seen from this plot that while both methods have the same starting point, the wavefunction spreads out more for the MCEv1 simulation. This is confirmed in Fig. 4 which shows how the basis functions spread in phase space over the course of propagation for both formulations of the MCE method, both starting with the same initial basis   Fig. 1(b).
set. It is clear by comparing the plots that for MCEv1 the wavefuntion spreads out over the course of propagation, as is indicated by the reduction in the relative density of basis functions at the centre of the initial wavefunction, whereas for MCEv2 the relative density of the basis functions at the centre of the wavefunction is still high, having only dropped down to 71% of the initial maximum, compared to a drop to 54% for MCEv1. This effect can be understood by considering that the time propagation equations for the coherent states use the Ehrenfest Hamiltonian H Ehr , which has a dependency on the amplitudes d kr . For the MCEv1 equations, the interconfigurational coupling is contained within these amplitudes and so the coherent states will effectively "push" on each other, spreading the basis functions out to cover a larger area in phase space. As this is not the case for the MCEv2 equations, the coherent states become less spread out and so the basis set cannot adequately describe a sufficient area of phase space to fully account for the quantum mechanics of the system.

Basis set refinements and improvements for the MCEv2 method
In light of the differences between the results given by MCEv1 and MCEv2, modifications to the sampling and propagation must be considered. The first modification was to use coherent state trains. An important consideration is that of the spacing between the basis functions in a particular train, δ trn . The train spacing parameter determines the degree to which the initial basis functions overlap. Too large a spacing would result in a loss of coupling between the basis functions, as this coupling is dependent upon the overlap matrix. As such, care must be taken in choosing the correct spacing between the basis functions and it should be remembered that the optimum spacing parameter will be dependent upon the system being simulated. In some systems, the basis functions may move faster through phase space than in others, resulting in a more rapid reduction of the coupling. It should also be noted that the number of basis functions in the ensemble is of great importance, as the combination of a small δ trn value and a small number of basis functions will result in trains which do not have sufficient size in phase space to properly describe the wavefunction.  Fig. 1(b).
The determination of the best value for δ trn is carried out through a set of preliminary calculations in which the conservation of the norm over the course of propagation is examined. As only the norm is examined, there is no need for comparison against benchmark results at this stage. It was found through this process that for the symmetric wells a train spacing of δ trn = 0.25∆ 1 was sufficient to construct the trains, and an arrangement of 10 trains, each of the 10 basis functions in length totalling N bf = 100 gave the wavefunction a large enough area in phase space. For the asymmetric wells, a smaller train spacing of δ trn = 0.15∆ 1 was needed as larger spacings caused the wavefunction to lose coupling. This is counteracted by an increase in the length of the trains to 20 basis functions in length. The wavefunction uses 10 of these trains for a total of N bf = 200 basis functions.
In Fig. 5(a), it can be seen that the oscillations in the population difference for the symmetric case are reduced in size and the result is much closer to the MCTDH benchmark result with only a slight increase in the size of the oscillations compared to the benchmark. The degree to which the population differences agree at this stage is unsurprising as the symmetric case is considered to be one of the easiest cases of the spin boson model to simulate. Unfortunately, the same cannot be said for the asymmetric case shown in Fig. 5(b), in which the result only follows the MCTDH benchmark for the first oscillation before the oscillations are over-damped and misguided such that the population difference seems offset from the MCTDH benchmark and thus will decay more slowly onto the second electronic state.
Clearly trains alone are not sufficient to fully overcome the differences between the MCEv1 and MCEv2 results. As such, we consider the effect of applying cloning. As the cloning procedure increases the basis set size greatly, this can cause problems with the system requirements of the simulation. As the calculations needed for the propagation of the basis set parameters are dominated by matrix-vector operations, the amount of computational resources required scale with the number of basis functions on the order of (N bf ) 2 . As such, while a simulation using a swarm of N bf = 50 basis functions will take usually less than an hour to complete normally, when cloning is included this can increase the size of the basis set by a factor of up to 2 N cln for N cln cloning events over the course of the simulation, increasing the runtime to a matter of days, and furthermore causing the memory requirements to increase from the range of tens of MB to a few GB. As such, limits must be put on the amount of cloning allowed. Obviously a larger N cln is better as too few cloning events will mean that the simulation cannot benefit properly from the procedure. In symmetric cases of the spin boson model where = 0, we set N cln = 4, as in this case exponential growth of the basis set is encountered. For asymmetric cases where 0, the situation is slightly different, as for this system the wavefunction as a whole is decaying onto the second electronic state and so once the cloned basis functions are placed wholly onto the two states it will only be the function placed on the first electronic state which will experience cloning again. As such, the basis set only grows by a factor of N cln + 1 which allows much more cloning to occur before the size of the basis set becomes unmanageable and so for this case we set N cln = 8.
Figures 6(a) and 6(b) show the degree to which cloning can improve the MCEv2 method for the spin boson model for the symmetric and asymmetric cases, respectively. As with the application of basis set trains, cloning dampens the overlarge oscillations in the population difference from the MCEv2 method, bringing the result closer to the MCTDH benchmark simulations, although not so much as with the use of trains. Also like the application of trains, this dampening is not sufficient to bring the two results into complete agreement for both cases. For the symmetric case, the agreement between the cloned result and the MCTDH benchmark seems to get worse around the second minimum in the population difference. This can be understood by considering the fact that the dynamics of the wavefunction cause each basis function to be cloned around every 0.75∆ 1 -0.8∆ 1 and so the majority of cloning for N cln = 4 would take place in the first 4.5∆ 1 of the simulation. Therefore any deviation caused by insufficient cloning can be expected to be seen after this time. For the asymmetric case, the oscillations are much smaller than for the standard MCEv2 method past the first oscillation but still FIG. 6. Comparisons of the population differences from cloned MCEv2 simulations with those from the uncloned MCEv2 simulations, both using the swarm-type basis set and N rpt = 100 repetitions. These are also compared against those from MCTDH simulations. 31 Symmetric case (a) uses an initial basis set of N bf = 50 basis functions and N cln = 4 cloning events. Asymmetric case (b) uses an initial basis set of N bf = 100 basis functions and N cln = 8 cloning events. too large. It should be noted however that despite the larger oscillations, with cloning the wavefunction seems to decay onto the second electronic state at a similar rate as for the MCTDH benchmark. This is in contrast to the result when using trains, where the wavefunction appears to decay much slower as evidenced by the offset in the population difference. It can also be seen that as with the symmetric case the agreement between the result and the benchmark worsens as the maximum number of cloning events is reached, which for the asymmetric case happens around t = 7.5∆ 1 .
While neither of the improvements considered are sufficient to correct the MCEv2 method alone, a combination of the two methods yields much better results, as shown in Fig. 7. This combination can be referred to as MC-MCE in an analogous way to the ab initio AIMC-MCE method. In Fig. 7(a), it can be seen that the agreement between the population differences from the MC-MCE simulations (i.e., MCEv2 with both cloning and trains) and the MCTDH benchmark calculations 31 is in much better agreement than was seen previously. The agreement is, in fact, almost to the same level as is given by the MCEv1 method seen in Fig. 1(a) with only a slight   FIG. 7. Comparisons of the population differences for cloned MCEv2 simulations using a swarm/train type basis set against those from uncloned swarm-type MCEv2 simulations and those from the MCTDH benchmark calculations. 31  deviation in the population difference right at the end of propagation. The improvement seen in Fig. 7(b) is much better than is seen for either trains or cloning used independently; indeed for the majority of propagation, the agreement is almost complete, with a slight overestimation of the oscillations toward the end of the simulation. This result is extremely encouraging, especially when one considers that the final cloning event occurs around t = 6.75 a.u. and that it is shortly after this that the discrepancies begin. Figure 8 confirms the earlier assertion that with cloning the basis functions are able to spread out more to cover a greater area in phase space, and it can be seen that the decay of the overlap with cloning is more in line with the way in which the basis functions spread when propagated using the MCEv1 equations. Again, however, the decay of the overlap levels out somewhat once cloning has ceased. It can therefore be reasonably expected that if more cloning events were allowed, the agreement with the benchmark calculations would persist for longer. FIG. 8. Comparison of the normalised average overlap between the coherent states for both formulations of the MCE equations with uncloned swarm type basis sets and the cloned swarm-of-trains type basis sets with the MCEv2 equations for the asymmetric case of the spin boson model, using the parameters given for Fig. 7(b).

B. Performance of MCEv2 for spin boson model with further parameter sets
In light of the success of MCEv2 with the combination of trains and cloning (referred to as MC-MCE) when applied to the low temperature symmetric and asymmetric cases of the spin boson model, a selection of more challenging cases was considered. As with the previous cases, results are compared against MCTDH benchmark results. A full accounting of the parameters used is given in Table I. The first of these more challenging cases was that of a pair of high temperature symmetric wells. This case considers a system with only 15 degrees of freedom and a high temperature with low coupling between the electronic states.
Considering the low coupling, and the ratio of degrees of freedom to number of basis functions, previous investigations into the spin boson model 31,33 would indicate that this would be a system which should be relatively easy to simulate, although due to the high temperature a large number of repeats are needed to achieve convergence. This is indeed the case and convergence appeared to have been reached by 2000 repeats (although it should be noted that to be consistent with Ref. 8 10 000 repeat calculations were performed for unmodified MCEv2), and it can be seen in Fig. 9 that the MCEv2 result does not deviate by an unreasonable amount from the MCTDH result. The size of the oscillations is however larger for MCEv2, and it can be seen that the introduction of the multiple cloning procedure does serve to dampen them somewhat, with the result matching the MCTDH result past around t = 1.75∆ 1 . This could conceivably be improved with the use of longer trains; however, the fact that this system seemed to benefit from much more closely packed trains than later cases would raise doubt on this supposition.
We next consider a pair of cases that model decoherence and illustrate the importance of having adequate dimensionality for a system. If the dimensionality of this system is too low, then an unphysical ringing is observed after the wavefunction becomes decoherent as can be seen in Fig. 10(a). As seen before, the oscillations in the unmodified MCEv2 result are much too large; however, these are dampened by the use of the multiple cloning procedure, bringing the result much closer to the MCTDH benchmark data, especially in the period before the decoherent section in both Figs. 10(a) and 10(b). Unfortunately, cloning stops around 1450 steps in of 10 000 (around t = 7.25∆ 1 ), which is just before the decoherent section, and so the results deviate after then. This is most noticeable in the unphysical ringing section of Fig. 10(a), while in Fig. 10(b) this is only noticeable between t ≈ 10∆ 1 and t ≈ 15∆ 1 as the wavefunction is transitioning from oscillating to being fully decoherent, after which time the MC-MCE result once again agrees with the MCTDH result. We next consider two limiting cases of the MC-MCE methods. First we consider localised wavefunctions which use very high dimensionality, high cutoff frequencies, strong system/bath coupling, and very low (near zero) temperatures to ensure that the wavefunction remains in the initial electronic state, which can be seen in Fig. 11. While MCEv1 performed well in this situation, 8 MCEv2 displays high frequency oscillations where no oscillations are expected, and also oscillates around a value for the population difference which is higher than that expected from the benchmark result. This behavior is echoed in the MC-MCE result, which only differs from the MCEv2 result in the size of the oscillations. This however is to be expected as in a localised wavefunction the condition of |a 1k a 2k | > 0.249 which triggers cloning will never be reached as by definition there will never be significant population in both electronic states. As such, the only improvement seen is the smoothing from the use of coherent state trains. It should however be pointed out that despite the oscillations, the wavefunction remains localised in the initial electronic state and does not show any signs of decaying.
The second limiting case models tunneling between symmetric wells at low temperature with fairly strong system/bath coupling and can be seen in Fig. 12. This is another case where MCEv1 was able to reproduce the MCTDH result successfully; 8 however, MCEv2 is unsuccessful in this regard, instead localising in the initial electronic state. When the same FIG. 12. Comparison of the MC-MCE and MCEv2 population differences for tunneling between a pair of low temperature symmetric wells with fairly strong system/bath coupling against a MCTDH benchmark calculation, 31 using the parameters ω c /∆ = 7.5, α K = 0.6, β∆ = 5.0, and /∆ = 0 with M = 60 degrees of freedom. For the MCEv2 result, N bf = 50 basis functions and N rpt = 100 repetitions were used, and for the MC-MCE results an initial basis set of 10 trains with 10 basis functions each was used for a total N bf,init = 100 with N cln = 4. (a) uses simple cloning and cloning threshold |a 1k a 2k | > 0.249 with δ trn = 0.05∆ 1 averaged over N rpt = 50 repetitions, and (b) uses an initial blind clone and cloning threshold |a 1k a 2k | > 0.1 with δ trn = 0.15∆ 1 averaged over N rpt = 200 repetitions. procedure is used to generate the initial basis set as for the other cases, the MC-MCE result seen in Fig. 12(a) is virtually unchanged from the MCEv2 result also localising in the initial electronic state, and so as with Fig. 11 the conditions for stimulating cloning are never met. To overcome this, two modifications are made. The first is the lowering of the cloning threshold such that cloning is stimulated when |a 1k a 2k | > 0.1. The second modification involves performing a "blind" cloning (i.e., not dependent upon the value of the amplitudes) at t = 0, whereby for each initial basis function, starting in the first electronic state, ψ k (t) = 0(0|ϕ 1 + 1|ϕ 2 )|z k (t) .
This allows the wavefunction to transfer population into the second electronic state on the other side of a potential barrier during propagation. Figure 12(b) shows the results of including these modifications. To carry out this calculation, it was found that the wavefunction benefited from an increase of the train spacing parameter to δ trn = 0.15∆ 1 and more repetitions were required to achieve convergence compared to that used for the results shown in Fig. 12(a). Through these modifications, it can be seen that the MC-MCE result now agrees much more with the MCTDH benchmark, although still not an exact agreement. A final case can be considered that of a pair of symmetric wells at high temperature with strong system/bath coupling, shown in Fig. 13. As with Fig. 9, the higher temperature necessitates more repeat calculations to achieve convergence than for lower temperatures. Interestingly this case out of all those considered does not require the addition of trains and cloning to the method; indeed, the results indicate that both the MC-MCE and the unmodified MCEv2 results have a similar level of agreement with the MCTDH benchmark. It is likely that due to the combination of high temperature and strong system/bath coupling, this case can be solved semi-classically and so has no need for the measures taken to prevent the MCEv2 calculations from guiding the basis functions in such a way that there is an inadequate spreading of the wavefunction in phase space.

IV. CONCLUSIONS
This article has investigated the influence of sampling on the convergence of the MCEv2 method, which recently has been used extensively for direct dynamics. It has been seen that for systems which cannot be accurately simulated with the unmodified MCEv2 method, the combination of using basis function trains when sampling the initial basis set and employing a basis set cloning procedure during propagation with the MCEv2 method, referred to as MC-MCE can generate results with an excellent level of agreement to benchmark calculations. It has also been shown that this combination is necessary, as using just one of these modifications is not sufficient to correct the disagreement entirely. While the computational cost of the MC-MCE method with these modifications is high compared to that of MCEv1 for the spin boson model due to the increased size of the basis set and a corresponding increase in computational time on the order of (N bf ) 2 , it has been shown that comparative levels of accuracy are achievable. The higher computational cost of MC-MCE is still negligible in comparison to the cost of electronic structure calculations when using the direct dynamics AIMC-MCE method. It has been shown here that MC-MCE yields good results for all cases of the spin boson model except those where quantum tunneling or low temperature effects are the most pronounced (shown in Figs. 11 and 12). Conversely, the most classical case shown in Fig. 13 is the most natural for MCEv2 even without the extra sampling procedures used by MC-MCE.
This result is of great significance as it confirms the validity of the methods applied by Makhov et al. 10,11 in their simulations of ultrafast processes in small organic molecules. In such simulations, due to the structure of the MCEv2 equations, the single configurational equations can be run separately to the calculations for the cross-configurational amplitudes [Eq. (16)]. It is the single configurational equations which contain computationally expensive electronic structure calculations and so separating these from the rest of the simulation can be beneficial. If information is saved at each time step, this can speed up calculations significantly, as each individual configuration propagated can be made into a basis function train at no extra cost, meaning that the train basis set is effortlessly constructed. Cloning also only occurs on this level, meaning that the cross-configurational calculations [Eq. (16)] can be run independently of the modifications discussed here. In confirming the validity of the modified MCEv2 method, the original motivation for creating the second formulation of these equations is finally and fully realised.