Non-adiabatic dynamics, where systems non-radiatively transition between electronic states, plays a crucial role in many photo-physical processes, such as fluorescence, phosphorescence, and photoisomerization. Methods for the simulation of non-adiabatic dynamics are typically either numerically impractical, highly complex, or based on approximations which can result in failure for even simple systems. Recently, the Semiclassical Monte Carlo (SCMC) approach was developed in an attempt to combine the accuracy of rigorous semiclassical methods with the efficiency and simplicity of widely used surface hopping methods. However, while SCMC was found to be more efficient than other semiclassical methods, it is not yet as efficient as is needed to be used for large molecular systems. Here, we have developed two new methods: the accelerated-SCMC and the accelerated-SCMC with re-Gaussianization, which reduce the cost of the SCMC algorithm up to two orders of magnitude for certain systems. In most cases shown here, the new procedures are nearly as efficient as the commonly used surface hopping schemes, with little to no loss of accuracy. This implies that these modified SCMC algorithms will be of practical numerical solutions for simulating non-adiabatic dynamics in realistic molecular systems.

## I. INTRODUCTION

Standard molecular dynamics protocols frequently evoke the Born-Oppenheimer approximation, i.e., nuclear motion is significantly slower than electronic. Thus, nuclear dynamics are controlled by the instantaneous configuration of the electrons. Full quantum treatment of the nuclear dynamics, while very useful in small systems,^{1–3} is limited only to a few degrees of freedom (few atoms). For large systems, one usually treats the nuclei classically, while treating the electrons quantum mechanically, the so-called mixed quantum-classical treatment.^{4–7} One can build a potential energy surface (PES), either *on-the-fly* or by sampling nuclear configurations before dynamical simulation. Nuclear motion is determined using Newton’s equations, with forces which depend on this PES. This effectively reduces the non-local (in phase space) Schrödinger equation to a set of localized Newton’s equations. This “adiabatic” approach successfully couples quantum chemistry methods, which are also local in phase space, with classical molecular dynamics and has been effectively used to study ground state properties of large systems.^{8,9} Additionally, dynamics of molecules in the excited state can be similarly simulated using excited state PESs determined by excited state methods.^{10–12}

However, there are a large number of chemical processes, such as non-radiative relaxation processes,^{13–15} photoisomerization,^{16–19} and intersystem crossing,^{20–22} which may involve transitions between electronic states. If the energy gap between coupled PESs approaches the scale set by the inverse time associated with nuclear motion, i.e., the vibrational frequencies, then the Born-Oppenheimer approximation breaks down and “non-adiabatic” transitions can occur.^{23} This causes a change in the electronic state without the absorption or emission of a photon.

It is highly desirable for applications in photochemistry, optically functional materials, solar energy, etc., to efficiently simulate non-adiabatic dynamics in a similar manner as has been well established for adiabatic dynamics. That is, quantum mechanics is used to treat electrons, while classical mechanics is used for nuclei. Whereas for “adiabatic” dynamics, the separation between quantum and classical systems is, more or less, clearly established by the Born-Oppenheimer approximation, for non-adiabatic dynamics, the separation between quantum and classical effects is less clear due to the non-adiabatic coupling (NAC) between nuclei and electrons.^{24}

The simplest mixed quantum-classical method for simulating non-adiabatic dynamics is the mean-field based Ehrenfest method.^{25} The nuclear dynamics are controlled by an average force, based on electronic state populations. These populations are determined by an approximation to the electronic Schrödinger equation. However, the dynamics of many systems are not well described by an average PES. Thus, Tully developed the stochastic fewest-switches surface hopping (FSSH) algorithm,^{26} which uses many trajectories to describe the dynamics of the nuclei on the PES of a single electronic state, with a finite probability to hop to other electronic states over the course of the dynamics. This method is simple to implement and converges relatively quickly. Thus, it has become the most practical tool for simulating non-radiative relaxation dynamics in molecular systems.^{27–34} However, the probability to hop is determined based on *ad hoc* approximations which often lead to erroneous results in simple systems.^{26}

Exact quantum dynamics of the nuclei is non-local in electronic energy space and in phase space. The mean field Ehrenfest method fails because it cannot treat this non-locality in the electronic energy space (the nuclei propagates only on one surface, which is a weighted average of all the adiabatic surfaces). This non-locality in electronic state is accounted for in FSSH (by running many trajectories which hop between states), but FSSH still cannot properly treat nuclear dynamics which strongly depend on the electronic state. A FSSH trajectory solves an approximate Schrödinger equation which is a function of its position at a given time. It is only exact in the limit that the dynamics of the trajectory do not depend on the particular electronic state which it occupies. This dependence is important when hops result in a large rescaling of the momentum or there is a large difference in the forces on the different PESs. This results in the often discussed “over-coherence” problem, and many approximations are capable of inserting some decoherence to correct this.^{35–40} However, if dynamics on separate PES are different, but not enough to completely separate the nuclear wavepackets, then two problems arise: (1) the approximate methods for decoherence may not provide the correct decoherence rate and (2) if the wavepackets recombine, there will be interference effects. In FSSH, the assumption that electronic dynamics are determined purely by the nuclear dynamics of the single trajectory will not properly account for this, and it will only be correct in the high kinetic energy limit (when differences in PES are insignificant).

Significant effort has been applied to improving the surface hopping method by including more sophisticated considerations.^{35,38–44} Alternatively, rigorously derived semiclassical methods, such as the quantum-classical Liouville equation, significantly improve on the accuracy of Ehrenfest or FSSH.^{45–50} In the Quantum Classical Liouville Equation (QCLE) method, trajectories corresponding to populations (diagonal elements of density matrix) propagate on single PES. However they can do half-hops, changing into coherence (non-diagonal density matrix element), propagate on an averaged PES, and acquire a phase factor. Then they can do another half-hop to return to being populations. This method does properly account for the phase difference between different “paths” in the dynamics. However, as mentioned before, it is expensive because the number of trajectories required for convergence is of the order of 100 000s. Some approximations can be used to help with this, but they can lower accuracy.^{47} The FSSH, with decoherence and phase corrections, may be considered as an approximate scheme for solving the QCLE, at least in the limit that dynamics on different PESs are similar.^{51} Other semiclassical methods, such as the Meyer–Miller (MM), Stock and Thoss Hamiltonian, and surface hopping Herman-Kluk semiclassical initial value representations (HK-SCIVR), are also accurate but expensive.^{52–56} Symmetric windowing has been shown to somewhat improve convergence rates for the MM model (down to 50 000–100 000 for Tully’s test problems); however, it is essentially a smoothing procedure and can give erroneous results for sharp transitions.^{26,57,58} Additionally, methods based on path-integral formalism, while formally rigorous, can be difficult to converge.^{59–61} Ring-polymer molecular dynamics (RPMD) has been successfully applied to ground state calculations for large systems and has recently been extended to non-adiabatic dynamics.^{62–65} RPMD does a good job of incorporating quantum effects at the single surface level, i.e., tunneling and zero-point energy (ZPE). While this method is quite effective for describing equilibrium statistics, it may become problematic if one seeks information on non-equilibrium properties, such as properties of relaxation dynamics. A lack of phase information (assumption of instant decoherence) may limit the accuracy and applicability in some non-adiabatic systems, including those addressed in this paper. RPMD has been combined with FSSH in an *ad hoc* way, as an attempt to include tunneling and ZPE into FSSH.^{66} However, this does nothing for the fundamental issues of FSSH which we address here.

We recently developed and discussed the semiclassical Monte Carlo (SCMC) method for describing non-adiabatic dynamics.^{67,68} In practice, it is a surface hopping algorithm, like FSSH and QCLE, with correct treatment of phases of trajectories (like QCLE). The advantage over QCLE comes from the fact that it is based on an expansion of the semiclassical wavefunction, instead of the density matrix. This leads to faster convergence rates, on the order of 10 000 trajectories, in addition to avoiding the averaged PES propagation.^{67} However, this is still far from the order of 1000 trajectories that is achieved with FSSH.^{26} In this paper, we discuss new algorithms which we have developed to bring the cost of SCMC down to this level of thousands of trajectories. The first is a numerical reorganization of the algorithm, which provides tremendous cost savings (5–10×), with no loss of information. The second is an approximation that groups of trajectories leaving a region of non-adiabatic coupling can be re-Gaussianized, creating a single trajectory which has the correct parameters. This is similar to the approach of *ab initio* multiple spawning (AIMS).^{69} Both of these new methods show essentially no loss of accuracy from the original SCMC method, at least in the model problems thus considered. This article is organized as follows: in Section II, we rewrite the necessary equations for calculating the SCMC wavefunction. In Section III, we discuss the numerical reorganization of the algorithm, which we call the accelerated semiclassical Monte Carlo or A-SCMC. In Section IV, we describe the A-SCMC procedure coupled with the re-Gaussianization, the A-SCMC-RG. We have included, as the supplementary material,^{72} an extensively detailed description of the algorithms, details of the re-Gaussianization procedure, and descriptions of all model potentials used in the paper.

## II. THEORY

In Ref. 67, we reported on the perturbative expansion of the time-dependent semiclassical wavefunction, Ψ(*t*) (represented in the basis of adiabatic electronic states), in powers of the non-adiabatic coupling, **d**. This leads to a series of terms which correspond to different numbers of electronic state transitions (hops). Each hop comes with an additional non-adiabatic coupling term, inside another time-ordered integral. Thus, all possible numbers of hops, state combinations, and hopping times are included. Similar expansions are used for solving the QCLE and surface hopping method for non-adiabatic HK-SCIVR.^{48,70}

Here, we simply rewrite the final result. For readability, we assume a two-state system and that at *t* = 0, the wavefunction has a non-zero component for a single electronic basis state, *s*. Generalization to the case of an initially mixed wavefunction and multiple electronic states is absolutely straightforward,

where the Gaussians are defined as

Here, **R**_{C} is the position vector of the Gaussian, **P**_{C} is the momentum vector, $ Z \xaf \xaf C $ is the complex inverse width matrix, and *S _{C}* is the phase (classical action). The dependence on (

**t**,

**s**) means the variable depends on each time integration variable, {

*t*

_{1},

*t*

_{2}, …}, and each state in the sum, {

*s*,

*s*′,

*s*}.

**D**(

**t**,

**s**) is the product of NAC at each integration time. As an example for the integrand in the second term in Eq. (1),

*G*

_{s,s′,s}(

**R**,

*t*

_{1},

*t*

_{2}), we have

and

Here, **F**_{s}(*t*) is the force defined by the position at time *t* and the electronic potential energy surface for *s*, and **v**_{s}(*t*) is the velocity additionally defined by the initial velocity. $\Delta P C s \u2032 , s $ is the momentum shift during a hop from state *s* to *s*′. This shift rescales momentum in the direction of the NAC vector in order to conserve total energy.

Thus, we have an infinite sum of weighted “trajectories” which evolve on their initial potential energy surface,^{7} until at some time they hop to another surface (assuming the trajectory has enough kinetic energy). The momentum is rescaled and the weight of the trajectory is scaled by the NAC for the old and new states at the position of the trajectory at the time of the hop. This continues until the trajectory has gone through all its hops. The trajectory is gaining a phase factor as it propagates, which controls the behavior of the ensemble of trajectories.^{7} Since we cannot use an infinite number of trajectories, we solve these multi-dimensional integrals over hop times by Monte Carlo integration. In our previous publications, we ran an ensemble of trajectories which have a finite probability to hop at a given time, sorted them into different hop orders and used them to sample the hop time integrals.^{67,68} We show similar convergence rates for one and two dimensional model problems; thus, we suspect that scaling with respect to the number of degrees of freedom is dominated by the quantum chemistry method used for dynamics and not the SCMC procedure. In Secs. III and IV, we present a new, more efficient, method for solving Eq. (1). The details of the method are described in the supplementary material.^{72}

## III. A-SCMC

We previously reported on the SCMC method for calculating non-adiabatic dynamics.^{67,68} The method solves Eq. (1) by using independent surface-hopping trajectories as sample points for solving the different integrals by Monte Carlo integration.^{71} The accuracy has been shown for many model systems, including systems with multiple level crossings and up to two dimensions. This accuracy comes with the price of slow convergence, relative to the standard FSSH algorithm. Like FSSH, the SCMC method is based on running independent trajectories and summing their contributions (see Fig. 1(a)). Unlike FSSH, in SCMC, the weight of each trajectory is used for Monte Carlo integration and each trajectory carries the important phase information. From this point, we will refer to this method as the independent trajectory (IT-)SCMC. This algorithm is very simple, trivially parallelizable, and conversion of standard surface hopping code to IT-SCMC is almost trivial. However, this approach results in a significant amount of repetitive calculation. For example, in a two-level single crossing system (Tully problem 1),^{26} one can expect ∼20% of trajectories to not hop. For example, if 25 000 trajectories are run, about 5000 trajectories will be identical (zero hop) trajectories. Additionally, each trajectory only samples the integral which corresponds to the number of undergone hops. However, a trajectory which followed the exact same path, but did not make the last hop, would be, up to that last hop, a perfectly good sample of the lower order integral. The goal of this section is to lay out an alternative approach to generate trajectories, the accelerated (A-)SCMC method, which provides tremendous cost savings by eliminating essentially all redundant calculations present in the IT-SCMC method.

In the IT-SCMC method, see Fig. 1(a), all trajectories run and then they are sorted by the number of hops and intermediate states to determine which integral in Eq. (1) they sample. In A-SCMC, see Fig. 1(b), the 0th order integral is solved using a single trajectory. The 1st order is solved using numerical quadrature, requiring only a few trajectories branching off of the 0th order trajectory. Higher order integrals are solved by Monte Carlo integration, creating new branches off of the other trajectories one order lower. Due to the time ordering of the integrals in Eq. (1), each hop will occur later than the previous hop, decreasing the average length of the dynamics with increased order of hop. This allows us to easily calculate these higher order integrals, whereas in the IT-SCMC, a large increase in trajectories is required to generate high order trajectories. Each integral is solved one-by-one up to a pre-specified order (8th order except where otherwise specified). Alternatively, convergence with respect to order of hop can be checked before moving on to the next order. Since one-hop trajectories are more expensive (on average) than higher order trajectories, due to their earlier average start time, it is advantageous to reduce their number. Thus, using quadrature instead of MC integration for the one-hop trajectories results in up to a factor of 2 reduction in cost. However, in order to generate random samples for two-hop trajectory initial conditions, we must build splines of the 2-D (first hop and second hop times) surfaces for the time-dependent parameters of the trajectories. We discuss in detail the A-SCMC algorithm in the supplementary material.^{72}

To demonstrate the improved efficiency of the A-SCMC method, we calculate the relative standard deviation (*δ*_{rel}) of the calculated percent transmitted on the lower surface for Tully’s problems 1 and 2. *δ*_{rel} is plotted as a function of the average number of “effective trajectories” (see Fig. 2). The effective trajectory is defined as the total number of calls to the quantum chemistry package (QCP) for the simulation (i.e., force calculations), summed over all trajectories, divided by the number of time steps in a single full length trajectory (a trajectory that starts at the initial time). In an independent trajectory method, the number of effective trajectories is equivalent to the total number of trajectories, since each trajectory is full length. Assuming that the call to the QCP is the most time consuming step in the dynamics (as is the case for *on-the-fly* dynamics of realistic systems), then this is a useful measure of the total computational cost of the simulation. In A-SCMC, most trajectories branch off of a lower order trajectory, requiring less time steps and less calls to the QCP than a full length trajectory. In Figure 2, the initial wave vector is chosen so that transmission on each surface is approximately 50%. For a single crossing (P1), the FSSH and A-SCMC converge at a nearly identical rate. A-SCMC achieves a 5% relative error by ∼500 effective trajectories and 2.5% relative error by ∼2000 effective trajectories. The IT-SCMC, with renormalization, achieved these levels of convergence using 5000 and 25 000 trajectories, respectively.^{67} Like the IT-SCMC result, the double crossing (P2) problem is slower to converge (see inset of Fig. 2). The reasons for this are presented in Ref. 67 and are still applicable here. To give the best assessment of the benefits of A-SCMC, the initial position and length of the simulation are chosen to coincide with the wavepacket just entering the region of non-adiabatic coupling at the beginning and just leaving by the end. If the initial position was decreased, the number of effective trajectories would decrease, as more of the calculation would be considered redundant. Alternatively, if the time was extended so that trajectories leave the region of coupling, the number of effective trajectories would increase. The choice to apply the method only in the region of significant non-adiabatic coupling is also inline with the re-Gaussianizing procedure which will be discussed in Sec. IV.

## IV. A-SCMC-RG

The re-Gaussianizing (RG) procedure is based on the following approximation. After all trajectories (which will be re-Gaussianized together, see below) leave the region of non-adiabatic coupling (i.e., **d** ⋅ **P** ≈ 0 or |**d**| ≈ 0 for all trajectories), one can replace the thousands of Gaussians (SCMC trajectories with complex weights), produced by either SCMC procedure, with a small number of new Gaussians without significant loss of information. If trajectories do not exit the region of non-adiabatic coupling, then the SCMC procedure will continue without interruption. Formally, this is just a basis set reduction. This approximation is reasonable in the same limit that SCMC is applicable, i.e., the momentum of the wavepacket is reasonably high; thus, different hop times do not spread the wavepacket too much and passing through non-adiabatic coupling does not take the wavepacket far from Gaussian. However, the wavepacket may branch into a small number of Gaussians and this is still tractable. In the case where passing through the region of NAC takes the wavepacket very far from Gaussian, it will require a large number of Gaussians in order to describe the on-going dynamics of the true wavefunction. As this case will be very difficult to sample, it is not clear that any method based on classical dynamics and trajectories will provide meaningful results and be able to reasonably converge.

In practice, the A-SCMC-RG procedure is similar to AIMS.^{69} Both involve a single wavepacket entering a region of non-adiabatic coupling, and multiple (ideally few) wavepackets emerging from the region. Both AIMS and SCMC wavepackets are Heller type Gaussians, though in SCMC, they are “free-thawed” and in AIMS, they are “frozen.”^{7} The real difference is in how those new wavepackets are generated. In AIMS when new wavepackets are “spawned,” their weights are found using the nuclear overlap matrices and their time derivatives. Multiple approximations, such as saddle-point approximations and the finite basis size, are required to make the AIMS feasible for realistic sized systems.^{69} Here, we generate many more packets, but find their weights by the SCMC procedure, and then, when leaving the region of non-adiabatic coupling, reduce the number of wavepackets by re-Gaussianization.

The momentum, position, width, phase, and real weights of the new Gaussians are calculated from the SCMC trajectories. The details of this procedure are described in the supplementary material.^{72} For systems in which coupled regions are far enough apart, this leads to incredible savings in the number of trajectories required. Interestingly, FSSH is likely to fail in systems with multiple well-separated level crossing, due to the incorrect treatment of the phases.

We determine the number of new Gaussians to use based on the number of possible “significantly different” paths for the trajectories. For example, in Tully’s double crossing problem (P2), an incoming wavepacket (on the lower surface) with enough momentum to transmit on either surface has four possible paths: (1) through both crossings without hopping, (2) hopping on the first crossing, then through the second on the upper level, (3) through the first crossing on the lower level, then hopping up the second, and (4) hopping up at the first crossing, then down at the second. A trajectory that hops up in one crossing and then down in the same will not constitute a “significantly different” path, it acts as a correction to the trajectory that does not hop at all. Alternatively, reflection, or in higher dimension, a significant momentum shift, could cause a “significantly different” path. To determine points along a trajectory in which a small change in hop time will correspond to a different path, we check for points in the dynamics where

This can be seen by looking at Eq. (10) of Ref. 67, where the right hand side would either diverge or the direction of the rescaled momentum becomes undefined. Thus, the rescaled momentum cannot be assumed to be a slowly varying function of time before and after passing through this point. Luckily, this is also a condition defining zero probability of hopping, providing a gap between trajectories spawned before and after this point. The specifics of how we implement this grouping and how to calculate the parameters of the new Gaussians are discussed in the supplementary material.^{72} It should be noted that for the case of high momentum and well-separated single-level crossings, each crossing will generate just two new Gaussians, one on each surface. As another example, following this procedure will correctly find three groups of trajectories for Tully problem 3 (extended coupling with reflection) at lower momentum, corresponding to the transmitted wavepacket on the lower surface and the two reflected wavepackets on the upper and lower surfaces.

We demonstrate that the RG procedure does capture the correct parameters for the new Gaussians by looking at three model problems (Figure 3). The first is a version of Tully’s double crossing, but with an extended middle portion (Fig. 4 (upper)). The A-SCMC-RG result reproduces the exact solution very well, whereas the FSSH procedure does not correctly reproduce the phase of the Stueckelberg oscillations at lower momentums. This is more prominent than in Tully’s original double crossing problem due to the increased distance between the crossings. FSSH cannot correctly account for the phase shift that comes with the separation of the wavepackets on the upper and lower surfaces. This phase shift is more drastic when the wavepackets have more of a chance to separate. We can take this even further by going to a system with two avoided crossings which are even further separated. In this case, we see that the A-SCMC-RG continues to accurately reproduce exact results, but the FSSH results are completely inaccurate (Fig. 4 (middle)). To further convince of the accuracy of the A-SCMC-RG result, we look at the wavefunction (real part) after passing through both crossings. The A-SCMC-RG wavefunction clearly reproduces the correct amplitudes and phases (evident by almost the exact same interference profile). The only deviation from the exact result seems to be in the tails of the wavefunction. This is likely due to ignoring the second derivative of the potential surface in the “thawed” Gaussian dynamics used.^{7}

In our previous report, we tested the IT-SCMC on a one-dimension 3-level 3-crossing system.^{36,67} While we found that the IT-SCMC produced results in good agreement with the exact solution, the required number of trajectories was extremely high (500 000 trajectories). In Fig. 5, we compare A-SCMC-RG results to exact solution and FSSH. For the majority of initial momentums, the A-SCMC-RG reproduces the exact result very well, using less than 4000 effective trajectories, although with the sharp transitions associated with semiclassical trajectories. However, in the region between *k* = ∼ 11 and 15, a trajectory can hop to the upper surface, but be caught inside the well. The method proves inaccurate in this regime, even predicting some probability to remain on the upper surface at long times. This momentum range is truly in the “non-perturbative” regime, as some trajectories spend infinite time in the NAC region. The SCMC methods, A- and IT-, are most appropriate when trajectories can make clean passes through interaction regions, without the possibility of a high numbers of oscillations inside the region of non-adiabatic coupling.

Interestingly, the IT-SCMC will not predict that any density will remain on the upper surface in this regime, as all trajectories will eventually hop down and leave the interaction region. The A-SCMC will sample the integrals which correspond to density remaining on the upper surface, and very high orders would be required to insure cancelation due to changing sign of the integrals. We still see some deviation from the exact result slightly above *k* = 15.5, this is due to the very large hopping induced spreading of the wavepacket on the upper surface at these low momenta. There is not enough room to RG between the middle and right most crossings. This should not be a problem area, but more trajectories, and possibly higher order in hops, are required to achieve convergence in this area. This shows that there is still some room for “fine-tuning,” the process by which trajectories are allocated and the required order of hops is determined.

## V. CONCLUSION

The SCMC procedure is a high accuracy surface hopping method which correctly accounts for phase interference and decoherence. These phase interference effects arise in systems with multiple level crossings. Correct treatment of the phase is most important when there is significant separation between these crossings, relative to the energy gap and initial momentum. The accelerated SCMC procedure increases efficiency compared to the independent trajectory SCMC by up to an order of magnitude. In a single crossing system, A-SCMC is shown to be as efficient as the standard FSSH method. SCMC trajectories can be used to find parameters which allow one to well approximate the calculated wavefunction by a small number (as low as one per surface) of Gaussians, allowing for swift propagation through the regions of negligible non-adiabatic coupling and better scaling with the number of crossings. These two methods can combine to give a speedup of up to two orders of magnitude. Subsequently, reasonable convergence at ∼1000 effective trajectories, similar to the FSSH method, was achieved across many model problems while maintaining high accuracy. The approximate number of trajectories required to achieve convergence for the models discussed here is summarized in Table I. This places the A-SCMC and A-SCMC-RG algorithms as practical numerical solutions applicable to realistic molecular systems.^{31}

Approximate number of effective trajectories for convergence to ∼1% standard deviation populations . | |||
---|---|---|---|

Problem . | SCMC . | FSSH . | A-SCMC-RG . |

P1 (k = 20) | 20 000 | 2000 | 2 000 |

P2 (k = 27) | 100 000 | 2500 | 20 000 (no RG) |

P2 − gap (k = 32) | 250 000 | 2500 | 3 000 |

P1 × 2 (k = 25) | 300 000 | 2000 | 2 000 |

Model X (k = 20) | >500 000 | 2000 | 3 000 |

Approximate number of effective trajectories for convergence to ∼1% standard deviation populations . | |||
---|---|---|---|

Problem . | SCMC . | FSSH . | A-SCMC-RG . |

P1 (k = 20) | 20 000 | 2000 | 2 000 |

P2 (k = 27) | 100 000 | 2500 | 20 000 (no RG) |

P2 − gap (k = 32) | 250 000 | 2500 | 3 000 |

P1 × 2 (k = 25) | 300 000 | 2000 | 2 000 |

Model X (k = 20) | >500 000 | 2000 | 3 000 |

## Acknowledgments

We acknowledge support of the U.S. Department of Energy through the Los Alamos National Laboratory (LANL) LDRD Program. LANL is operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. DE-AC52- 06NA25396. We acknowledge support of the Center for Nonlinear Studies (CNLS) and the Center for Integrated Nanotechnology (CINT) at LANL. V.G. gratefully acknowledges support from Texas A&M University at Qatar via the NPRP Grant No. 6-021-1-005 from the Qatar National Research Fund.