The modular decomposition of the path integral is a linear-scaling, numerically exact algorithm for calculating dynamical properties of extended systems composed of multilevel units with local couplings. In a recent article, we generalized the method to wavefunction propagation in aggregates characterized by non-diagonal couplings between adjacent units. Here, we extend the method to the calculation of reduced density matrices in aggregates where each unit includes an arbitrary number of coupled harmonic bath modes, which may describe intramolecular normal mode vibrations, at finite temperature. The effects of harmonic modes are included through influence functional factors, which involve analytical expressions that we derive. Representative applications to spin arrays described by the Heisenberg Hamiltonian with dissipative interactions and to J-aggregates of perylene bisimide, where all coupled normal modes are treated explicitly, are presented.
For several decades, Feynman’s path integral approach to time-dependent quantum mechanics1,2 has found wide application in all areas of theoretical physics. In combination with Monte Carlo sampling methods,3 its imaginary formulation4 has been used as a numerically exact and efficient simulation tool for computing equilibrium properties of quantum many-body systems obeying Boltzmann or Bose statistics.5,6 Such simulations have led to a wealth of information regarding the structure of low-temperature fluids, superfluids, and the Bose–Einstein condensate.7 Furthermore, the quantum–classical isomorphism8 motivated the introduction of molecular dynamics techniques to equilibrium path integral calculations,9–12 which are more efficient for some types of interactions.
The main obstacle in applying Feynman’s formulation to real-time dynamics is the presence of a highly oscillatory phase, which arises from the complex exponential of the action integral along each quantum path.13 The paradigm of quantum dissipation is based on the interaction of a system with a harmonic bath,14 where a large (often infinite) number of harmonic degrees of freedom can mimic many effects of generic condensed phase environments. For example, the effects of dissipation on tunneling dynamics15 and reaction rate constants16,17 have been thoroughly investigated this way. The effective harmonic bath picture has also been shown to arise as the dominant contribution from a condensed phase environment in the so-called Gaussian (or linear) response limit18 and justifies the parabolic free energy surface model employed in electron transfer theory.19 Recently, rigorous quantum–classical path integral20 (QCPI) simulations of the ferrocene–ferrocenium electron transfer in a solvent with strongly nonlinear intermolecular forces21 established the quantitative accuracy of the effective harmonic bath picture for some processes involving charge transfer. Furthermore, harmonic degrees of freedom arise naturally in the description of intramolecular normal modes, and the reaction path Hamiltonian22 employs exactly this type of system–bath interaction.
A distinct advantage of the path integral approach is that harmonic bath degrees of freedom may be integrated out analytically within this formulation, giving rise to an influence functional23 that captures exactly the effects of the harmonic bath on the dynamics of the quantum system of interest. However, the influence functional is nonlocal in time, and thus, the path integral cannot be evaluated by Markovian matrix-vector procedures. The iterative quasi-adiabatic propagator path integral (i-QuAPI) is based on tensor decompositions of the reduced density matrix,24,25 two-time,26,27 or multi-time correlation functions,28 which capture bath-induced memory effects and allow for propagation to long times. The full i-QuAPI algorithm involves the storage of tensors whose size grows exponentially with the memory length, but path filtering procedures,29,30 coarse graining,31 blip,32,33 and singular value decomposition (SVD)-based matrix product operator (MPO) decompositions34 can dramatically reduce the cost. A recent small matrix disentanglement of the path integral variables35,36 (SMatPI) reduces the required storage to that of ordinary matrices of size equal to that of the bare system, extending the feasibility of numerically exact simulation to multi-state systems.
Unlike path integral-based methods, approaches based on the Schrödinger equation have to treat the bath degrees of freedom explicitly. Among wavefunction-based schemes, the multiconfiguration time-dependent Hartree (MCTDH) method37–39 along with its multilayer40 and effective mode41 extensions has shown impressive success on select system–bath models.42 The main advantage of the MCTDH approach is its ability to also treat anharmonic bath degrees of freedom. However, wavefunction-based methods are naturally suited to zero-temperature properties, and their extension to finite temperature by summing many microcanonical calculations is computationally prohibitive when the environment includes a large number of low-frequency modes. In the special case of a harmonic bath described by a spectral density of the Drude form, the hierarchical equations of motion (HEOM) method43 allows efficient simulation of system–bath dynamics. However, its use in more general situations requires a decomposition of the spectral density into Drude-type components,44 and since the method scales exponentially with respect to the number of such terms, this approach becomes impractical for simulating processes in real, structured environments. Methods based on the density matrix renormalization group45,46 (DMRG) and matrix product states47,48 can be highly efficient for extended systems, but their cost increases steeply with the number of bath degrees of freedom, and thus, these methods are restricted to a small number of vibrational modes. A recent development combining matrix product techniques with a pseudomode representation of the bath49 appears to overcome some of these limitations in cases where the bath spectral density can be efficiently decomposed, but its ability to tackle molecular environments with structured vibrational spectra may be limited.
Recent developments have extended the applicability of real-time path integral methods to more complex situations. In particular, it has been shown that extended systems described in terms of discrete Hamiltonians with sufficiently local interactions can be treated exactly via a modular path integral (MPI) approach,50,51 where the quantum paths of a unit are linked to those of the adjacent unit and summed without encountering a sign problem. The MPI methodology is applicable to systems with a one-dimensional topology, such as chains of spins or molecular aggregates, but can also account for branching or ring structure, and leads to linear scaling with system size. While the structure of the MPI algorithm is simplest in Ising-type models, where the coupling is local in the single-unit basis, we recently extended the methodology to the general case where the coupling is not diagonalizable in any single-particle basis.52 We also introduced an efficient factorization of the MPI linking algorithm that leads to exponential savings.53
In this paper, we develop an MPI methodology suitable for propagation of density matrices, where each of the units includes coupling to a finite-temperature harmonic bath, which might represent intramolecular vibrations that couple to the electronic states of the unit or a generic dissipative environment. We show that the coupled harmonic degrees of freedom can be integrated out analytically, giving rise to influence functional factors that apply separately to each unit, and derive expressions for the coefficients in the discretized influence functional expression.
We emphasize that the development described in this paper is not a method for treating system–bath dynamics. Even in the absence of bath degrees of freedom, calculations on extended quantum systems, such as the quantum Ising and Heisenberg spin models or molecular aggregates of two or more electronic states, are limited by the usual exponential scaling of quantum mechanics. For example, a basis set treatment of a chain containing r spin- units requires storage of a wavefunction with 2r elements, and thus, calculations quickly become impractical when r exceeds 10. On the other hand, the MPI decomposition is characterized by linear scaling with r. Moreover, the inclusion of harmonic bath degrees of freedom is possible with very little additional cost and is just as simple for discrete vibrational modes as it is for a model spectral density. These attributes suggest that the MPI methodology is able to address important problems not amenable to other approaches.
In Sec. II, we describe the aggregate Hamiltonian, where the discrete states of each unit are coupled to harmonic bath modes. For concreteness, we refer to the discrete states of each unit as “electronic states” and to the attached bath degrees of freedom as “vibrational” modes, with the understanding that a unit might be defined in terms of spin states, and the bath degrees of freedom might describe a generic dissipative environment. In Sec. III, we develop the path integral expression and the MPI algorithm for the density matrix and the expectation value of a target operator. In Sec. IV, we describe the procedure for integrating out the harmonic degrees of freedom and derive expressions for the influence functional coefficients. In Sec. V, we apply the methodology to examine the effects of a dissipative bath on the dynamics of Heisenberg spin chains and to simulate excitation energy transfer in J-aggregates of perylene bisimide, where all coupled normal modes of each monomer are treated explicitly. Finally, in Sec. VI, we give some concluding remarks.
II. MOLECULAR AGGREGATE HAMILTONIAN
Consider an array of r units for which the Hamiltonian has the form
where the subscript α = A, B, …, Y, Z labels the monomers. Each unit is described by the Hamiltonian
describes the Mα electronic states |nα⟩. Each of these states is coupled to vibrational degrees of freedom, which are described by the Hamiltonian
is the operator that specifies the coupling to the electronic system. Equation (2.4) is diagonal in the chosen electronic basis. Its diagonal element
describes the vibrational Hamiltonian associated with the electronic state |nα⟩. The states |nα⟩ may be eigenstates of , in which case Eq. (2.4) is the normal mode expansion of the Born–Oppenheimer potential energy surfaces, or may describe diabatic states coupled through nondiagonal terms in Eq. (2.3). The basis states might also describe spin degrees of freedom in magnetic materials or in model Hamiltonians, such as the transverse field Ising chain54–56 or the Heisenberg model,57 which interact with a local phonon or generic dissipative environment.
The molecular units are coupled to their nearest neighbors through the electronic coupling terms
This type of coupling is particularly relevant for modeling the electron-vibrational dynamics in molecular aggregates coupled through electrostatic or other non-covalent interactions. For example, the Frenkel exciton58 coupling form is often used to describe excitation energy transfer in molecular aggregates.
III. MODULAR PATH INTEGRAL REPRESENTATION
In the general case where the electronic coupling is not diagonalizable in a single-unit basis, we construct pairs of adjacent units by partitioning the Hamiltonian as follows:52
etc. We use the letter Z to label the last unit of the chain and Y to label the next-to-last unit, regardless of the number of units. The expectation value of an operator can be obtained from the propagated density matrix,
where ρ(0) is the initial density operator and the trace is with respect to all (electronic and vibrational) degrees of freedom of all units. The variables labeled and indicate the endpoints of the forward and backward propagators.
In a recent paper,52 we described the MPI representation of the electronic propagator, which yields the time evolution of the electronic wavefunction. While the pair structure of the short-time propagators appears to double the number of path integral variables, we showed that the additional auxiliary variable can be eliminated by implementing a swapped-pair propagator sequence. Here, we adopt this structure and extend the MPI algorithm to the density matrix, where the units include vibrational degrees of freedom. We retain the path representation in terms of electronic states only, while the vibrational degrees of freedom appear in operator form and are immediately traced over.
A diagrammatic representation of the MPI linking expressions for the density matrix is shown in Fig. 1. Note that only even-numbered path integral variables remain in the representation of the first unit A, while odd-numbered variables survive in the last unit. To simplify the notation, we have relabeled the variables and using a numbering scheme compatible with the path integral beads of each unit. For the first and last units, we use and . For units B, C, …, Y, the last path integral variable is , so we use the notation .
If ρ(0) is a product of single-unit operators (e.g., all spins of an Ising chain are initially in the “up” state, or a single molecule is electronically excited at t = 0 and the vibrational modes are equilibrated with respect to that state) and if the electronic properties of a particular molecular unit α1 are being probed, the expectation value of the observable operator , where (the identity operator for monomer α) for all α ≠ α1 is given by the following sequence of MPI steps:
where TrA,v denotes the trace with respect to the vibrational degrees of freedom of unit A,
etc. Each of these steps links the paths in the electronic space of a monomer to those of the next unit. The B-to-C and C-to-D steps are performed sequentially until the next-to-last unit is reached. If the chain has an odd number of units, the last step is given by
and the observable is obtained through the sum
The propagators that appear in these expressions are illustrated in Fig. 1. We note that this diagram corresponds to one of the two mirror-image diagrams in the swapped pair MPI expression, which arise from reversing the order of the Trotter splitting. Averaging the results with respect to these two symmetric forms reduces the Trotter error, allowing a larger path integral time step. A factorization of the path linking operations53 leads to a very efficient algorithm that scales as , where Lα is the number of paths of unit α.
This modular, one-monomer-at-a-time decomposition of the path integral is key to ensuring a linear scaling of the computational cost with the number of monomers. Each element of the arrays Rα is a specific sequence of the path variables of monomer α and hence a particular pair of forward and backward paths for this unit. These paths only interact with those of monomer α + 1 with amplitudes that depend on the specifics of the operator and are linked to those through the expressions.
IV. INFLUENCE FUNCTIONAL
Consider the Hamiltonian for a unit of two adjacent monomers. The Hamiltonian for the pair α, α + 1 becomes
We partition this Hamiltonian into pure electronic and electron-vibration components. Noting that and commute, the short-time propagator is given by the following product:
The path integral variables of unit A are separated by time steps 2Δt and are coupled only to those of unit B. Evaluating Eq. (4.2) in the electronic basis according to the swapped-pair representation of the path integral,52 the first two forward propagator elements for the AB pair are
Recall again that are operators in the space of the vibrational degrees of freedom of unit A, with the potential minimum at the coordinate values . Thus, the exponential operators in Eq. (4.3) describe evolution of a multidimensional forced harmonic oscillator, where the time dependence of the force is specified by the particular sequence of electronic states and the corresponding time lengths. Collecting all the short time propagators that correspond to the forward and backward evolution operators and including the vibrational component of the initial density operator , we see that the vibrational modes of unit A augment the electronic amplitude in the path integral by the factor
where Trv denotes the trace with the respect to the vibrational degrees of freedom and is the vibrational partition function for unit A.
The product of operators in Eq. (4.4) describes the time evolution of the vibrational degrees of freedom subject to forces determined by the value of the electron-vibration coupling index sA, which changes at multiples of 2Δt. This familiar form is the time-discretized Feynman–Vernon influence functional,23 augmented by a phase59 that arises from the “counterterms”14,15 in Eq. (2.6). If the bath is initially in thermal equilibrium at the temperature T = 1/kBβ, the influence functional for unit α has the form
where Cα(t′) is the force autocorrelation function of unit α and is the electronic index that specifies the vibrational state along the given path. The correlation function is given by the expression
It is often expressed in terms of the spectral density function,14
To evaluate the influence functional arising from the vibrational modes of a particular unit, we discretize the double time integral, taking advantage of the piecewise constant form of the force. This brings the influence functional to the form of a discrete sum,
where the coefficients are to be determined.
In the case of unit A, the force changes at the end of each double time step, and thus, the discretization of the double integral in Eq. (4.8) is straightforward. For example, it is easy to see that
etc. Performing the time integrals using Eq. (4.6) leads to expressions for the discretized influence functional coefficients.
We now proceed to examine the vibrational influence functional for unit B. Equation (4.3) shows the sequence of forces that arises from the interaction with unit A during the first few time steps. However, because of the splitting of Eq. (3.2), half of the vibrational Hamiltonian of unit B is contained in the propagators for the pair BC. The first two of the forward time propagators are given by the expressions
Combining the vibrational operators that contain B variables in Eqs. (4.3) and (4.11), one sees that the time dependence of the force for unit B is given by the sequence specified in Fig. 2. Using this force sequence to discretize Eq. (4.5), we find
etc. After performing the time integrals, we obtain the discretized influence functional coefficients for the vibrational modes of unit B.
Next, the variables of unit C enter also through the CD propagators. The first of these is split into electronic and vibrational components as follows:
Combining Eqs. (4.11) and (4.13), we observe that the discretization of the force for monomer C is identical to that of monomer B. In fact, it is not hard to see that every non-terminal monomer has the same discretization of the force. Therefore, all non-terminal monomers B, C, D, …, Y have an identical set of discretized influence functional coefficients, which thus need to be evaluated only once.
Finally, we obtain the third and final unique discretization of forces for unit Z. Because unit Z is only coupled to unit Y, the sequence of forces can be completely obtained from the YZ product of propagators. The first two such propagators are
The sequence of forces for unit Z is shown in Fig. 2. Using this force sequence, we obtain the influence functional coefficients for the last unit of the chain. The discretized vibrational influence functional coefficients are given in Tables I–III for all units.
|ηA,00 = ηA,NN|
|ηA,00 = ηA,NN|
|ηα,00 = ηα,N+1,N+1|
|ηα,11 = ηα,NN|
|ηα,10 = ηα,N+1,N|
|ηα,N0 = ηα,N+1,1|
|ηα,00 = ηα,N+1,N+1|
|ηα,11 = ηα,NN|
|ηα,10 = ηα,N+1,N|
|ηα,N0 = ηα,N+1,1|
|ηZ,00 = ηZ,NN|
|ηZ,11 = ηZ,N−1,N−1|
|ηZ,10 = ηZ,N,N−1|
|ηZ,N,1 = ηZ,N−1,0|
|ηZ,00 = ηZ,NN|
|ηZ,11 = ηZ,N−1,N−1|
|ηZ,10 = ηZ,N,N−1|
|ηZ,N,1 = ηZ,N−1,0|
Finally, we note that just as in the case of a single system coupled to a harmonic bath,60 the MPI influence functional coefficients could also be expressed in terms of correlation functions obtained from classical molecular dynamics calculations.
In this section, we illustrate the inclusion of harmonic degrees of freedom on a spin chain, consisting of two-level systems coupled to model bath described by a continuous spectral density, and a molecular aggregate, where we treat explicitly all available normal mode vibrations of each monomer.
A. Heisenberg spin chains
Models of spin arrays exhibit intriguing dynamical effects, as well as quantum phase transitions. They are also of interest to magnetic materials, where external fields can be used to tune the spin tunneling parameter and thus alter the system’s dynamical behavior. Here, we explore the real-time dynamics of one-dimensional Heisenberg spin chains57 with nearest neighbor interactions (see Fig. 3) in the presence of dissipative interactions arising from coupled harmonic bath terms. The total Hamiltonian has the form
where are the Pauli spin matrices for unit α, and are the spin–spin coupling strengths, and is the external field. The dissipative bath is characterized by an Ohmic spectral density,15 given by
where ωc is the cutoff frequency and ξ is the dimensionless Kondo parameter.
Figure 4 shows the dynamics of a Heisenberg spin chain with r = 31 units, where each spin is coupled to its neighboring units, for two sets of spin–spin coupling parameters and . Even in the absence of coupling to a bath, calculating the dynamics of the Heisenberg chain presents a serious challenge. With the chosen value of r = 31 two-level systems, the electronic Hilbert space is equal to 231. Thus, the dissipative Heisenberg chain with 31 units is equivalent to a system of 231 (not just 31) states coupled to a bath. The bath parameters are set to ξ = 0.2, , and the temperature is . We start with all spins aligned in the “up” state and show the probability Pα(t) for one of the edge spins (α = 1) to remain in the “up” state as a function of time, while the rest of the spins are traced over. The results were obtained using 2(N + 1) = 22 path integral variables for each nonterminal spin. The calculation for the last time point took about 77 min for each nonterminal linking step on a single AMD Epyc 7551P 2.0Ghz processor, with a total of approximately 36 h for the 31-unit chain.
In the absence of dissipative interactions, an isolated spin undergoes coherent tunneling oscillations between the “up” and “down” states. As is well known,15 interaction of a single spin with a dissipative bath quenches the coherent dynamics, leading to underdamped oscillations or even monotonic decay to equilibrium. Figure 4 shows that spin–spin couplings are quite effective in quenching the amplitude of the single-spin oscillations, even in the absence of a dissipative harmonic bath. However, the resulting time evolution differs significantly from spin–boson dynamics. Furthermore, the two sets of spin–spin coupling parameters are seen to produce quite different dynamics. The first set of parameters with a smaller Jy coupling leads to an interesting interference effect, which causes an increase in the probability after the first recurrence of diminished amplitude. Interference is almost wiped out in the second case, where a blue shift in the oscillation of the probability is observed. When each spin is also coupled to a dissipative bath, the tunneling oscillations become further damped, and the interference characteristicsobserved with the first set of parameters disappears. Thus, we see that while coupling to neighboring spins generally leads to damping effects similar to those from harmonic dissipative baths, effects of spin–spin interactions are more complex and subtle, and the resulting dynamics depends delicately on the type and strength of coupling between spins.
In Fig. 5, we show the effects of increasing the temperature of the bath for the Heisenberg chain with . As expected, with the increase in available thermal energy, the damping effects become more pronounced, and eventually, the oscillatory pattern almost disappears.
B. Exciton-vibration dynamics in molecular aggregates
As a second example, we apply the procedure described in Secs. II–IV to simulate the excitation energy transfer dynamics in J-aggregates of perylene bisimide (PBI; see Fig. 6), whose photochemical properties have been the subject of the recent electronic structure calculations.61,62 Each PBI molecule has Mα = 2 electronic states, nα = 1, 2, and its electronic Hamiltonian is described by Eq. (2.3) with and eV.62 Each electronic state is coupled to 28 intramolecular normal modes according to Eq. (2.4) with sα,1 = 0 and sα,2 = 2.The neighboring units interact via Frenkel exciton58 coupling terms,
The electronic energies, inter-monomer couplings, normal mode frequencies, and exciton-vibration coupling coefficients were obtained from electronic structure calculations.62 The spectral density is shown in Fig. 6 in terms of Huang–Rhys factors, which are related to the coupling coefficients through the expression
Figure 7 shows the excited state population of the initially excited molecule in a small aggregate of r = 3 PBI molecules and a long aggregate containing r = 25 units at 300 K, with the initial excitation placed on the central unit. In the larger aggregate, a total of 700 discrete normal mode vibrations were treated explicitly via the fully quantum mechanical MPI procedure described in this paper. At the given temperature, the vibrational partition function of each molecular unit is equal to 2.8 × 1012, and thus, explicit treatment of all intramolecular modes via wavefunction-based methods is not practical. The calculation for the longest time point with N = 12 (26 beads) took about 10 min for each nonterminal linking step on a single AMD Epyc 7551P 2.0Ghz processor, with a total of approximately 220 min for the 25-unit aggregate. The coupling to molecular vibrations strongly affects the population dynamics, causing rapid damping of the electronic state population and some smearing of the minor recurrences observed in the long chain.
VI. CONCLUDING REMARKS
The MPI representation of real-time dynamics enables the simulation of extended systems characterized by local interactions with effort that scales linearly with the number of units, allowing the study of long systems that are not amenable to wavefunction-based schemes. As a path integral method, MPI allows the treatment of any number of harmonic bath degrees of freedom that may couple to each unit, at zero or finite temperature, with practically no additional cost. This is possible by obtaining the influence functional for each unit in terms of pre-computed coefficients, which are specific to the particular discretization of the force on the bath arising from the particular splitting of the short-time propagators. In the common case where the coupling between adjacent units is diagonal in the given basis, the influence functional coefficients are identical to those in the QuAPI expression.25,63 When the coupling between units has a general, non-diagonal form, two-unit propagators must be employed with the swapped-pair ordering in order to avoid the unnecessary doubling of path integral variables. This ordering leads to different force sequences and thus different expressions for the discretized influence functional coefficients. In this paper, we have described the procedure for including the effects of harmonic bath modes to the MPI algorithm and derived the expressions required to obtain the influence functional coefficients for the general form of the MPI Hamiltonian. The method can be implemented with a model spectral density or with discrete bath modes representing molecular vibrations.
The methodology we have developed allows the simulation of time-dependent properties in extended systems representing a variety of different physical situations. For example, it can be used to obtain dynamical properties of spin arrays described by quantum Ising or Heisenberg Hamiltonians in which the spins are subject to dissipative effects. Even in the absence of dissipative interactions, these systems continue to present significant challenges to numerically exact methods, as the size of the Hilbert space grows exponentially with chain length. Thus, the MPI methodology can contribute to the understanding of many interesting physics of these systems, particularly the influence of dissipative phonon environments at finite temperature on the dynamics of spin arrays and their behavior relevant to quantum phase transitions. We emphasize that the MPI methodology is not restricted to simple linear arrays of units but can be extended to obtain the dynamics of more complex structures, such as branched chains, cyclic geometries, and dendrimers. We have recently presented MPI calculations of excitation energy transfer in the B850 ring of the light harvesting complex encountered in photosynthetic bacteria.64
Another very interesting situation arises in the dynamics of molecular aggregates, which involve electronically coupled molecular units, each of which includes several tens or hundreds of intramolecular vibrational degrees of freedom. The resulting Hilbert space of these systems can be vast, particularly under ambient conditions, which lead to several populated vibrational states in many of the vibrational modes. The MPI methodology can treat such aggregates fully quantum mechanically. We have presented MPI calculations on excitation energy transfer in J-aggregates of perylene bisimide and (in a separate work64) investigated the energy transfer process in bacteriochlorophyll chains of varying length. Understanding such properties is important for designing structures for efficient solar energy harvest.
Finally, the MPI methodology can be used in conjunction with generalized quantum master equation65–67 (GQME) approaches in order to extend the time length that can be reached solely via the path integral calculation. GQME calculations require input obtained via different methods over the time span of the memory kernel.68–71 If this input is computed via numerically exact methods, then the GQME results at longer times will also be exact. Our recent work using QuAPI and QCPI calculations to build the GQME kernel,72 as well as the SMatPI decomposition35,36 of system–bath dynamics, suggests that the span of the GQME kernel is easily accessible by the time length of converged MPI calculations. Unlike the cases of QuAPI and QCPI, whose iterative decomposition is more efficient than GQME time-extensions,72 the combination of MPI with GQME is very appealing for obtaining numerically exact results on long molecular aggregates over long simulation times. Work along these lines will be reported in future papers by our group.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
This material is based upon work supported by the Air Force Office of Scientific Research under AFOSR Award No. FA9550-18-1-0291. The calculations were performed on a computer cluster funded by the National Science Foundation under Award No. CHE-1665281.