One Million Quantum States of Benzene - The Journal of Physical

Subscriber access provided by NEW YORK MED COLL

Article

One Million Quantum States of Benzene Thomas Halverson, and Bill Poirier J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.5b07868 • Publication Date (Web): 29 Sep 2015 Downloaded from http://pubs.acs.org on October 3, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

One Million Quantum States of Benzene Thomas Halverson, and Bill Poirier∗ Department of Chemistry and Biochemistry, and Department of Physics, Texas Tech University, P.O. Box 41061, Lubbock TX 79409-1061, USA E-mail: [email protected]

Abstract In this study, we compute all of the dynamically relevant vibrational quantum states of benzene, using an “exact” quantum dynamics (EQD) methodology. Benzene (C6 H6 ), in addition to being a very large molecule for EQD (12 atoms, 30 vibrational modes), also has a very large number of vibrational states—around 106 in all, lying within 6500 cm−1 of the ground state. The EQD methodology developed here uses a phase space picture to optimize the truncation of a harmonic oscillator basis—not only with respect to the molecular system of interest, but also the targeted spectral range. By employing several such EQD calculations, targeted to different spectral ranges, a “hybridized” dataset is constructed that provides the most accurate results everywhere. In particular, more than 500,000 states are converged to 15 cm−1 or better.

Introduction In the past few years we have seen a surge in the development of computational methods and hardware that have led to an explosion in the capabilities of computational chemistry. Within the context of exact molecular quantum dynamics (EQD)—e.g., the calculation of quantum vibrational states of molecular nuclei using no approximations except Born-Oppenheimer, with which we concern ourselves here—these developments have greatly extended the domain and size of molecular 1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Thomas Halverson et al.

One Million Quantum States of Benzene

systems that can now be considered. Much attention has been paid to the dimensionality of the systems considered, which until recently was largely restricted to just three-to-four atom molecules. Equally important however, are the other system parameters—e.g., masses and energy—that can also lead to a dramatic increase in the number of “dynamically relevant” quantum states. The chief difficulty stems from the so-called “curse of dimensionality” or exponential scaling that plagues traditional basis set methods: for D degrees of freedom, a naive EQD N × N matrix representation requires N = nD max , where nmax is the number of basis functions per degree of freedom. All competitive state-of-the-art EQD methods attempt to mitigate or circumvent the exponential scaling problem. Different techniques are used, but always these include some sort of customization of the basis set for the problem at hand, designed to minimize the value of N that is required. Thus, multiconfiguration time-dependent Hartree (MCTDH) 1 customizes a truncated direct-product basis (DPB) for the time-evolving wavefunction, at each instant of time. MCTDH can handle tens or even hundreds or thousands of degrees of freedom, and is ideally suited for time-dependent applications. In the time-independent (e.g. vibrational spectrum) context, MCTDH can also be applied— albeit not with the ease or spectroscopic accuracy of other purpose-built EQD techniques, such as those developed by T. Carrington, Jr. 2,3 The latter class of calculation typically combines basis optimization (via pruning, contraction, and/or dimensional combination) to reduce N, together with sparse iterative methods that circumvent direct matrix diagonalization and thereby enable much larger N values to be considered (e.g., N ∼ 108 ). Such methods extend EQD calculations to systems with up to perhaps 8 atoms or so, and also provide extremely high, “spectroscopic” accuracy (e.g., ∼ 0.01 cm−1 ). However, in most implementations, iterative methods only compute a few quantum states at one time, generally near the bottom of the spectrum—when in fact for large D, the number of dynamically relevant states is typically much much larger. For six-to-twelve-atom systems, the number of vibrational states within say, ∼ 6000 cm−1 of

the ground state can be on the order of 104 –106 . 4 Fortunately, the total number, K, of such dynamically relevant states does not increase exponentially with D—though it does increase as a

2


Page 2 of 45

Page 3 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




rather steep power law, 4–6 rendering calculation via iterative methods extremely tedious. What is desired is a methodology that can compute all dynamically relevant states for systems in this size range, from a single (or small number of) calculation(s). In this paper, we explore just such a method, building on earlier basis optimization ideas stemming from a tremendously useful and versatile phase space picture (PSP). 7–15 Indeed, in a series of previous papers, 4–6,15–19 we have used the PSP to design phase-space-localized basis sets (PSLBs, either Weyl-Heisenberg wavelet or momentum-symmetrized Gaussian) that formally defeat exponential scaling in the large N limit. More specifically, these methods achieve perfect asymptotic efficiency—i.e., the ratio of the number of converged states to the total basis size, (K/N), approaches unity as N → ∞. PSLB methods are ideally suited for the combination of large D and very large K, as discussed above. In particular, they have been applied to model systems with up to D = 27 and K ≈ 5000, 19

as well as to the six-atom D = 12, K ≈ 10, 000 acetonitrile application. 4 However, a limitation for large D is that only states in the lowest part of the spectrum are converged to spectroscopic accuracy—with the rest of the (K/N) → 1 states typically converged to ∼ 101 cm−1 . Substantially better accuracy can be achieved throughout the lower part of the spectrum using a harmonic oscillator basis set (HOB), provided an effective HOB truncation strategy is employed. 4 Truncated HOB basis sets have been used quite successfully for many years, e.g. in the EQD M ULTIMODE codes of J. Bowman and coworkers, 20 which rely on “polyad” truncation (Sec. ). Though highly effective at the very bottom of the spectrum where the system is least anharmonic, the polyad approach quickly loses accuracy as the energy increases. In this paper, we develop a much more flexible and accurate “smart” HOB truncation scheme, based on the PSP, which enables us to extend the efficacy of the HOB approach to a point much higher in the spectrum, so that a very large number of states may be obtained. Formally, even PSP HOB does not defeat exponential scaling, and so there must always be a point in the spectrum beyond which PSLB wins out, in terms of efficiency. In other words, below a certain point in the spectrum, PSP HOB provides more accurate energies than PSLB, but for energies above that point, the PSLB energies are more accurate—assuming comparable

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



basis sizes. However, where exactly this changeover point occurs is a moot point in the case of the benzene application of this paper, because of symmetry. Whereas permutation symmetry may be easily incorporated into a PSLB calculation, the same is not generally true for the many point group symmetry operations that characterize a highly-symmetric molecule such as benzene. The HOB, on the other hand, can indeed exploit most of the latter symmetry—thus reducing the whole Hamiltonian matrix down to independent symmetry blocks of a tractable size. So PSP HOB calculations are used exclusively here, because the corresponding PSLB matrices would be too large to diagonalize directly. The resultant symmetry blocks are small enough that they can be directly diagonalized on massively parallel computing clusters, which is the approach adopted here. It should be noted, however, that PSP methods do not require the use of massively parallel computers; it is the PSP approach itself, rather than the number of available cores, that overcomes the curse of dimensionality. Indeed, a D = 15 PSLB calculation was performed over a decade ago on a single-core computer. 18 Even for the current benzene application, for which the potential energy surface (PES) is taken to be a quartic force field, 21 the individual symmetry blocks are still sparse, and could therefore be amenable to an iterative treatment. We do not do so here because a direct diagonalization provides all of the energy eigenvalues at once, as discussed. This is much more convenient, and ultimately far more computationally efficient, if all dynamically relevant states are to be computed; however the trade off for this particular strategy is the need for massively parallel computers. In this paper, we apply the above ideas—massive parallelization, PSP HOB truncation, and symmetry adaptation—to compute all of the K ∼ 106 dynamically relevant vibrational states of benzene (C6 H6 ). With twelve atoms (D = 30), these calculations represent the extreme cutting edge of what is currently possible using “brute force” direct diagonalization methods, applied across very large computing clusters. It is very useful to know the limits of current capabilities— although we stress at the outset that more sophisticated refinements will be introduced in future that will enable still larger calculations to be performed with much greater accuracy on far smaller computers. In any event, for the present calculation, all K ∼ 106 states were computed to an 4


Page 4 of 45

Page 5 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




accuracy of a few tenths to a few tens of cm−1 ’s, using a total basis size of N ∼ 6 × 106 . Symmetry considerations reduce this impossibly large calculation down to eight separate block calculations of around ∼ 700, 000 basis functions each. Even with the reduced basis sizes, these calculations represent some of the largest direct diagonalizations ever attempted for the purposes of actual calculations. Nevertheless—and despite requiring millions of core hours—these calculations were relatively simple to run using our computational software package S WITCH BLADE, for which the user need specify but a single convergence parameter (Sec. ). One of the most important “take home messages” of this work is that a PSP analysis is extremely useful for optimizing a wide variety of basis sets—not just PSLBs. A second, critically important message is that, for large complex molecules like benzene, optimal results are obtained using different basis sets in different parts of the spectrum. We have already mentioned that PSLB is best at high energies, with PSP HOBs performing better at lower energies. This is, of course, due to the increased anharmonicity at higher energies, since a PSLB maintains (and can even improve 6 ) its efficiency in the presence of system anharmonicity/coupling where any HOB necessarily loses efficiency. However, even in the lower part of the spectrum where PSP HOB is more efficient, there are different ways to apply PSP HOB truncation, each with its own optimal spectral range. The beauty of the PSP is that it provides an explanation of exactly which HOB truncation should be applied where, and why. Hence, by appropriately combining multiple forms of PSLB and HOB truncation, we can compute a spectrum across the entire dynamically relevant range that is as accurate as possible. Moreover, the S WITCH BLADE codes provide exactly such flexibility (hence the term “switch”), thereby allowing true “smart” truncation to be easily applied. Since no calculations of this kind, i.e. D ≈ 30, K ≈ 106 , have ever been performed before, we conclude this Introduction by taking a moment to justify the motivation for doing so. In particular, with one million available quantum states, it is tempting to assume that we are operating in the “classical” limit—in which a classical treatment would perform just as well. In fact, this is not at all correct; we are still very much in the quantum regime, owing to the very large number of ways that relatively few quanta can be distributed across D = 30 degrees of freedom. Indeed, with just

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 6 of 45


five vibrational quanta total (i.e., the sum over all normal modes) there are 35!/(30! 5!) = 324, 632 quantum states in all—i.e., on the order of K. The second natural misconception is that, with so few quanta, anharmonicity/coupling may not be very significant, and so a harmonic approximation should work well. This, too, is incorrect. From Ref. 21, a typical cubic force field coefficient for benzene is on the order of 100 cm−1 —giving rise to many transition matrix elements of around, � � � � e.g., 1185 cm−1 (between the �v = 4 and �v = 5 overtones for a single mode). Moreover, there are many hundreds of anharmonic coefficients in all, once again owing to the large value of D.

Theory and Background The basic computational approach used herein is very straightforward: we seek to solve the time independent Schrödinger equation by diagonalizing the Hamiltonian matrix, as represented in a (truncated) direct product of HOB functions, i.e., ˆ j �, Hi j = �Φi |H|Φ

(1)

where � � �Φi = |n1 n2 . . . nD � = ΠD |nk �. k=1

(2)

Using only DPB truncation, the above method is extremely inefficient due to exponential scaling, N = nD max , as discussed in Sec. . To mitigate the curse of dimensionality, we must employ some correlated, non-DPB truncation criterion, such as

D

∑ αk nk ≤ Emax.

(3)

k=1

� � Simply stated, we restrict our HOB functions, �n1 , n2 , . . . , nD , in such a way that a weighted sum of the individual mode excitations is less than some maximum value. This greatly reduces N, and can

also greatly improve HOB efficiency—provided that the αk ’s are chosen in a way that is tailored for the specific application. The question then becomes, how best to choose the αk ’s? More 6


Page 7 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




specifically, how should these be chosen to optimize a given calculation within a given desired spectral range? In the past, these coefficients have been chosen in somewhat ad hoc fashion. 2,3,20 To fully understand the problem, a PSP analysis proves invaluable. The basic idea is that the K quantum states of interest—i.e., those that lie below some energy threshold—define some well-defined classical region of phase space. 7,15 Any representational basis meant to compute all K states must, as a necessary condition, incorporate this classical region, within which most of the quantum probability resides. However, if higher accuracy is desired, the representational basis must extend beyond the classical region, in order to also incorporate quantum tunneling probability. Thus, whereas the classical PSP provides an invaluable starting point, a more rigorous analysis requires a quantum PSP treatment—specifically, the Wigner quasi-probability distribution function (PDF) (i.e., the Wigner transform of the density operator). 22–24 A quantum PSP approach has been used in the past to both explain and rectify the origins of exponential scaling, ultimately giving rise to the successful PSLB methods. 16–18 In addition, however, a quantum PSP analysis explains why the PSLB methods provide limited accuracy for large D and K; it is because the tunneling behavior for PSLB’s (after canonical orthogonalization) is not an especially good fit to that of the desired high-lying quantum states. 16 In contrast, for quantum states that are largely harmonic-like, the true tunneling behavior is more closely approximated by an HOB representation. 4 We will address this in some detail in the remainder of this section, and also use the quantum PSP to provide a better understanding of how best to choose the αk values across different spectral domains.

Classical and short-range tunneling behavior To fully understand how truncation affects basis efficiency, we must perform a PSP analysis. However, the high dimensionality of benzene makes it difficult to gain insight directly from the full problem. To make things easier to visualize, we project the 60D phase space of the full Hamiltonian down to a 4D reduced phase space associated just with modes 2 (in-plane anti-symmetric stretch) and 5 (out-of-plane anti-symmetric bend), where the mode labeling is that of Wilson. 25 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 8 of 45


The result is as follows:

h(2) (q2 , p2 , q5 , p5 ) = H(0, 0, q2 , p2 , 0, 0, 0, 0, q5 , p5 , . . . , 0, 0).

(4)

Obviously, there are many mode pairs that could have been chosen, but (2, 5) exhibits a moderate or “typical” amount of anharmonicity and mode coupling for the benzene system. We will first consider two natural choices for the αk coefficients in Eq. (3). For the first choice, all modes are equally weighted, αk = 1; this is called “polyad truncation” (PT). For the second choice, the coefficients are chosen to be the harmonic frequencies of the corresponding normal modes, i.e. αk = ωk ; this is called “frequency truncation” (FT). Of course, these are not the only choices, but we shall find it convenient to regard these two as the opposite ends of a spectrum of possibilities. The quantity of interest is the projection operator for a uniform ensemble of K eigenstates of the D = 2 Hamiltonian, h(2) , i.e. 15

K

ρˆ K = ∑ |Ψi ��Ψi |.

(5)

i=1

We then transform ρˆ K to a phase space function via the Wigner-Weyl prescription, 22–24 yielding the Wigner PDF, ρK (q2 , p2 , q5 , p5 ), for Eq. (5). The position-space projection of this Wigner PDF can be seen in Fig. Figure 1(a). Of course, calculating ρˆ K requires a priori knowledge of the eigenvectors, which we obtained by diagonalizing the Hamiltonian matrix of Eq. (1) for Hˆ = hˆ (2) and {|Φ j �}Nj=1 , using a large HOB. Specifically, a sufficiently large N was chosen such that K = 20 states were all converged to ∼ 0.01 cm−1 or better. The choice K = 20 corresponds roughly to the dynamically relevant range for the h(2) system; from Fig. Figure 1(b), the anharmonicity/coupling is readily apparent. Next, we calculated projection operators for candidate representational HOBs, ρˆ N =

N

∑ |Φ j ��Φ j |,

(6)

j=1

using the two truncation methods, PT and FT, with Emax chosen for each such that N ≈ K. This 8


Page 9 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




is the smallest possible N value with any chance of spanning the desired classical phase space � � region. From Eq. (2), we know that for PT, the set of retained 2D basis functions, �n2 , n5 , corresponds to a “symmetric triangle” of points in (n2 , n5 ) space, which is obviously invariant with respect to coordinate exchange. This coordinate exchange permutation symmetry is evident in the corresponding ρN Wigner PDF, as can be seen in Fig. Figure 2. Indeed, from the figure, it is clear that ρN exhibits not only exchange symmetry but also pure rotational symmetry, with respect to the q2 and q5 coordinates. In any event, nmax is clearly the same for both coordinates. All of this is in stark contrast to the ρK distribution of Fig. Figure 1, which exhibits an oblong shape that is much more extended along q5 than q2 , owing to the smaller ω5 < ω2 frequency. In any event, for N ≈ K, the desired ρK and representational ρN PDFs are not even qualitatively close to each other, and so the PT basis yields very low efficiency. On the other hand, the classical PSP predicts that FT should do a much better job at capturing the desired classical phase space region, when N ≈ K. This is because Eq. (3) effectively becomes the true classical Hamiltonian, apart from anharmonicity/coupling. More specifically, the FT choice provides a higher nmax for q5 than for q2 because of the smaller ω5 value—exactly as it should. Indeed, as indicated in Fig. Figure 3, the resultant FT ρN Wigner PDF is an oval elongated along q5 , that closely resembles the true ρK —at least within the classical region. It is this similarity in shape that results in much higher basis efficiency for FT than for PT, 4,19 in cases when N ≈ K. In any event, this result should not be too surprising, given that PT uses no system-specific information at all to determine its αk ’s, whereas FT at least incorporates harmonic frequencies. However, the FT customization is relatively limited, as it does not take into account any of the mode coupling and anharmonicity, which are also evident in Fig. Figure 1(b). Specifically, these manifest as distortions from a purely oval shape, i.e. the “rightward skew” of ρK . Of course for increasing K, these effects become increasingly pronounced, rendering any HOB approach less and less efficient.

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Long-range tunneling behavior The previous analysis does provide a great deal of insight, but it does not tell the entire story. In particular, PT has been used for quite some time, 20 and can produce excellent results. So how is this to be reconciled with the conclusions of Sec. ? The crucial point is that the latter applies only when N ≈ K—in other words, when one expects all, or nearly all, of the computed states to be “converged.” This requires matching the desired and representational classical phase space regions as closely as possible, as discussed in Sec. . From a quantum PSP analysis, however, we know that the classical phase space region only contains something like 98% of the total ρK probability (at least for the most highly-excited states). 6,16–19 Therefore, although (K/N) may indeed be close to unity in an FT calculation, the accuracy of many of the computed states is typically fairly low—i.e., on the order of tens of cm−1 ’s (or worse, if anharmonicity and mode coupling are substantial). On the other hand, if we do not mind setting N � K, we can achieve far better accuracies for the

K computed states using PT, rather than FT. 4 To understand why this is the case, we must apply the quantum PSP to examine the long-range tunneling behavior of the requisite PDFs—since highly accurate convergence requires the inclusion of deep tunneling probabilities, from well beyond the

classical phase space region. In Fig. Figure 4 we present the long range tunneling of the ρK Wigner PDF for K = 5, for the D = 2 h(2) Hamiltonian of Sec. . Here, in order to properly reflect the higher desired target accuracy, the calculation was converged to ∼ 0.0001 cm−1 or better for each of the K = 5 eigenstates. In Fig. Figure 4, the innermost contour corresponds roughly to the classical phase space region, whereas the larger concentric contours correspond to the rapidly decaying values of ρK in the quantum tunneling region, down to well below numerical precision. Any representational basis that captures the region defined by the largest contour should therefore suffice to compute the lowest K = 5 quantum states numerically “exactly.” The main conclusion to draw from the figure is that for higher accuracies, the contours become increasingly circular, not elliptical. The reason for this is simply understood in terms of the long-range decay of the HOB functions themselves— which scales as the same Gaussian for all normal modes, regardless of excitation, nk , or frequency, 10


Page 10 of 45

Page 11 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




ωk . Of course, the practical upshot is that PT—with its perfectly spherically symmetric ρN contours— becomes the most efficient HOB truncation when N � K. Equivalently, a PT calculation will be more accurate than a corresponding FT calculation of the same size, for states lying at the very bottom of the spectrum. We can see this more explicitly in Figs. Figure 5 and Figure 6—contour plots for the PT and FT ρN Wigner PDFs, respectively, for N = 45, using the same convention as in Fig. Figure 4. From the figures, we see that the PT ρN density is indeed spherically symmetric, whereas the FT ρN remains anisotropic throughout (albeit decreasingly so, with decreasing ρN , as is to be expected). Therefore, if efficiency is not the primary concern, PT will provide the most well converged HOB results at the very bottom of the spectrum. Even for large D, PT can thereby achieve very high accuracies, but only for a relatively small number of quantum states. Conversely, FT is the much more attractive choice for computing extremely large numbers of energetically excited eigenstates with excellent efficiency—but is limited in terms of the accuracy that it can achieve. Of course, there are also a host of “intermediate” options, to be considered next.

Hybrid truncation: the “best of all worlds” The above analysis implies—and preliminary calculations suggest 4 —that for large D, the choice of αk coefficients has a marked effect on the performance of the HOB, in a way that can be tailored remarkably well for specific spectral domains. The closer the αk are to the ωk , the higher we can go in the spectrum and still obtain reasonably well converged results. As all of the αk approach the same value, however, we obtain very high accuracy in the lowest part of the spectrum, but quickly lose even qualitatively correct behavior with increasing excitation energy (with errors typically growing nearly exponentially 19 ). The above facts suggests that a “compromise” between FT and PT might work best in intermediate energy domains—particularly in cases where D is large, and/or when there is a substantial

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 12 of 45


spread in ωk values. One way to explore this idea is to replace Eq. (3) with D

∑ (ωk )µ nk ≤ Emax,

(7)

k=1

where 0 ≤ µ ≤ 1. This allows us to work with the single truncation parameter, µ, rather than with D separate parameters, αk . By comparing Eqs. (3) and (7), it becomes clear that µ = 0 corresponds to PT, whereas µ = 1 corresponds to FT. In addition to these two choices, we also consider µ = 0.5, lying “halfway” between PT and FT. As discussed above, this choice is expected to work best at intermediate energy ranges. As a pragmatic measure, we don’t actually work with µ = 0 when performing “PT” calculations, but rather, we use the small value, µ = 0.1, instead. The reason is that the extremely high symmetry of true PT as applied to a D = 30 calculation would result in very large sudden jumps in N as Emax is increased. This makes it very difficult to “tune” the PT basis size—a property which is desirable, among other reasons, because it enables accurate convergence testing and comparison with the other methods. For example, for µ = 0, D = 30, and 4 < Emax ≤ 5, we have N = 324, 632, but for 5 < Emax ≤ 6, N jumps suddenly to 1,947,792. On the other hand, by breaking the symmetry, µ = 0.1 provides much finer gradations in N, while still providing a truncated basis that is very similar to PT. For benzene for instance, µ = 0.1 and Emax = 180.0 yield N = 315, 061, whereas Emax = 185.0 yields N = 414, 212. Since we can now ensure the similarity of N values from one truncation method to the next, we can now also properly compare the different methods—with respect to both numerical convergence, and accuracy. It is worth noting that all HOB representations as implemented here are rigorously variational—meaning that all computed energy levels converge from above to the exact results, as Emax is increased. Of course, for comparable N values, a lower µ calculation should yield more accurate—and therefore lower—energy levels, Ei , in the lower part of the spectrum, whereas a higher µ calculation will yield lower/more accurate eigenvalues further up in the specµ

trum. In any event, for a given eigenstate, i, the lowest energy level, Ei , computed across different

12


Page 13 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




µ values, is the most accurate. This suggests a natural way to combine the different truncation schemes together, to create a “best of all worlds” dataset: one simply compares a sorted list of computed energy levels for µ

each separate µ calculation; for each eigenstate i, one then selects the lowest of the Ei values as the most accurate computed energy level for that state. We call this the “hybrid truncation” (HT) dataset. In principle, HT does not even require that the respective basis sizes for each calculation, µ

denoted NEmax , be comparable to each other. However, if one wishes to quantify the notion of µ

convergence for the HT dataset, to do this in a consistent way requires that the various NEmax values for the different µ calculations be similar. The HT dataset will be better than that of any specific µ value, throughout the entire spectrum. Generally speaking, one expects the lowest-lying states to come from µ = 0.1 (“PT”), the intermediate states from µ = 0.5, and the high-lying states from µ = 1.0 (FT). However, there are some exceptions to this basic pattern—such as a low-lying overtone for a low-frequency mode with multiple vibrational excitations. PT does not do well for such states, even if Ei is rather low, because it treats all modes equally. In any event, the final hybridized set of eigenvalues will always be more accurate than any of its constituents, since the varational principle ensures that the lowest computed energy is always the most accurate.

Computational Details The PES used for benzene is a full D = 30 quartic force field potential (FFP), developed by N. Handy and coworkers. 21 They used ab initio self-consistent field calculations, with double-zeta plus polarization basis sets, to compute all quadratic, cubic and quartic force constants, omitting coefficients < 1.0 cm−1 . They then used S PECTRO 26,27 to compute anharmonic frequentheory

cies, νk

expt

ues, νk

, via second-order perturbation theory. They then compared with experimental val-

(although benzene has only four infrared-active modes), 28 and discovered substantial expt

disagreement, with |νk

theory

− νk

| > 100 cm−1 in some cases. To partially correct for this, 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



they then defined a new FFP for which all of the new harmonic frequencies were adjusted via expt

ωk → ωk + νk

theory

− νk

(Table Table 1). Finally, they used S PECTRO to recompute the anhar-

monic frequencies, together with a few combination bands. Due to the accuracy limitations of second-order perturbation theory—and the rather large magnitude of the semi-empirical corrections employed—it is not expected that the more accurate spectroscopic calculations conducted here will yield better agreement with experiment. Nevertheless, this FFP serves as a very realistic and challenging benchmark application for our PSP methods. The reduced normal mode coefficients as presented in Ref. 21 can be expanded to the full D = 30 FFP using the Henry and Amat symmetry relations. 29,30 It should be noted that Ref. 21 contains the correct number of reduced coordinate coefficient values needed to properly expand the full FFP, save one. The term, Φ16a 16a 6a 6a is not reported, even though its symmetry partner, Φ16a 16a 6b 6b , is reported. To be on the safe side, and since the numerical value of the reported Φ16a 16a 6b 6b coefficient is only 2.9 cm−1 , we simply ignore both terms. For our own calculations, we computed energy eigenvalues via direct matrix diagonalization on massively parallel computers, using our S WITCH BLADE codes. Although the PSP methods can in principle be applied to arbitrary PESs (Sec. ), S WITCH BLADE is designed to work with quartic FFP’s of arbitrary system dimensionality. All matrix elements are analytically represented, leading to the aforementioned rigorously variational convergence. In the past, S WITCH BLADE has been applied to a range of model applications, and real molecules up to D = 12. 4,5,19 It was originally designed only for momentum-symmetrized Gaussian PSLB’s; however, the recent generalization for truncated HOB’s allows for much more versatility. For the D = 30 benzene application, symmetry adaptation of the basis set is required in order to reduce the basis size N to a tractable size, even for a massively parallel computer. Benzene belongs to the D6h point group, which has 12 singly- and doubly-degenerate irreducible representations (irreps). The HOB representation does not respect the full D6h symmetry of the problem. However, the HOB functions transform very straightforwardly (with ±1 character) under each of the three commuting symmetry operations, iˆ, σˆv , and Cˆ2 . This is easily exploited to allow for an 8–fold 14


Page 14 of 45

Page 15 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




reduction in N (eight independent symmetry blocks), reducing the largest basis size from N = 5, 706, 342 down to N = 719, 588. Direct diagonalization for matrices of even this size are quite challenging, requiring wall times of ∼ 40 hrs despite being run in parallel across 7200 cores. Most calculations were performed on the Texas Advanced Computing Center, Lonestar 4 facility.

Results and Discussion Three sets of calculations were done using the benzene FFP outlined in Ref. 21 and Sec. , corresponding to three different µ values: µ = 0.1; µ = 0.5; µ = 1.0. For each µ value, calculations for three different basis sizes were conducted: small; medium; large. One such calculation was performed for each of the eight symmetry blocks (Sec. ), resulting in a total of 72 separate calculations in all. The normal mode frequencies, ωk , used in Eq. (7), are the adjusted benzene frequencies as presented in Ref. 21 and also Table Table 1. µ

As discussed, Emax values were chosen so that small, medium, and large basis sizes, NEmax , are comparable across µ. The resultant basis sizes as used in each of the 72 calculations can be seen in Table Table 2, along with the irrep decomposition of each symmetry block. Specifically, each symmetry block corresponds to a unique singly-degenerate irrep (eight in all), plus one half of a corresponding doubly degenerate irrep (four in all). The splitting of the doubly-degenerate modes across two blocks (chosen to be adjacent odd-even rows in the table) makes it possible to uniquely identify the D6h irrep for any given state, by comparing computed energy levels across the block pair. All estimates of numerical convergence are accomplished by calculating differences between the large N calculation, and the corresponding medium or small N calculation. This provides two separate convergence estimates for each eigenstate. The presentation of two such convergence values allows for a better understanding of convergence, and also helps to rule out spurious or pseudo-convergence.

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 16 of 45


Low-lying energy levels For the lowest lying 30+ states in each symmetry block (257 total), we used energetic ordering of the ωk , together with symmetry considerations, to assign state labels based on mode excitations. The symmetry characters for the overtone and combination states were determined using the D6h direct product table. This procedure is straightforward—except that care must be taken to exclude antisymmetric combinations when there are multiple excitations of the same doubly degenerate mode. Excluding these forbidden combinations, we deduce the following combination rules for multiple products of identical doubly-degenerate modes: E ⊗ E = A1 ⊕ E1 E1 ⊗ E1 ⊗ E1 = B1 ⊕ B2 ⊕ E1 E2 ⊗ E2 ⊗ E2 = A1 ⊕ A2 ⊕ E1 E ⊗ E ⊗ E ⊗ E = A1 ⊕ E1 ⊕ E1 .

(8)

Note that the 1 and 2 subscripts do not play a roll in even-power products, but do do so when the power is odd. The low-lying states as described above are presented in Tables Table 3–Table 10, for which one table corresponds to each symmetry block. State labels are presented in Column II, in “reduced” (D = 20) normal mode notation [e.g., (6, 7) → 6, (8, 9) → 7, etc.; see Table Table 1]. The corresponding frequencies, νi = Ei − E0 (where E0 is the ground state energy), are presented in Column V. Since PT is expected to be most accurate at the bottom of the spectrum, the fre0.1 . Columns VI and VII present quencies correspond to the large µ = 0.1 calculation—i.e., N14.23 large

∆νi = νi − νi

—i.e., the numerical convergence of the medium and small calculation frequen-

cies, respectively, relative to the large calculation. About half of the frequencies are converged to subwavenumber accuracy, with most of the rest converged to within just a few cm−1 ’s. However, ν’s that result from a large excitation of a single low-frequency mode tend to be less well converged, as anticipated in Sec. . 16


Page 17 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




Along with the states we computed, Tables Table 3–Table 10 also show ν’s from Ref. 21 obtained using second-order perturbation theory (where available), in Column IV. Since identical FFP’s were employed here as in Ref. 21, a direct comparison of Columns IV and V provides an estimate of the error associated with second-order perturbation theory. These errors are generally on the order of 10 cm−1 ’s, although in some cases they are significantly larger. In any event, we report many more energy levels than in Ref. 21, as the latter authors were primarily concerned with just those states for which experimental data is available.

High-lying energy levels Reporting and properly assigning the low-lying ν’s is important, but the true power of the PSP approach is not realized until we consider more highly-excited states. As previously discussed, there are around one million such states that are dynamically relevant for benzene. Table Table 11 presents frequencies for a selection of dynamically relevant excited states, i, ranging from i = 1 to i = 1, 000, 000, across all eight symmetry blocks (Column I). To generate this table, the computed energy levels from all blocks were first concatenated together, and then resorted by energy. Columns II, V, and VIII, list frequencies, νi , as computed using the large calculation for each of the three truncations, µ = 0.1, µ = 0.5, and µ = 1.0, respectively. The remaining columns indicate numerical convergence, of the three medium and small calculations, relative to the corresponding large calculations. Note that here we present convergence of the energy levels—i.e., large

∆Ei = Ei − Ei

—rather than the frequencies. The third row in the table indicates which basis set µ

was used for each calculation, using NEmax notation. Table Table 11 also serves to show how each truncation method performs, in different regions of the spectrum. We find that the numerical convergence patterns are not as “clean” as anticipated in Sec. —nor indeed, as has been observed previously. 4 To be sure, µ = 1.0 (FT) is least well converged near the bottom of the spectrum, and the small µ = 0.1 (PT) calculation is very poorly converged as i → 106 , in accord with expectation. On the other hand, the medium µ = 0.1 calculation appears the most converged near the top of the spectrum. In reality, this is most likely a 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



pseudo-convergence, owing to the (nearly)-perfectly-symmetric PT behavior discussed in Sec. . It might also be due to the limited range of basis sizes considered; even though the large-to-small basis ratio used here is nearly 3×, which is quite large, it may be that for D ≈ 30, even larger ratios are needed. In any event, one solid indication that this is a pseudo-convergence, besides the fact that the small µ = 0.1 errors are so large, is the fact that the computed frequency, e.g. for i = 106 , is several hundred cm−1 ’s larger than for either of the other two µ values. We must be a bit careful in comparing frequencies rather than energy levels across methods, as only the latter are variational. Thus, the fact that µ = 1.0 provides the smallest ν1 value should not be taken as an indication that this is the most accurate frequency value, but rather, that it corresponds to the least well converged ground state. That said, a “lowest” frequency comparison is legitimate for differences that are greater than the worst convergence of the ground state energy— in this case, around 2.5 cm−1 ’s. By this reckoning, one sees that µ = 0.1 is indeed least accurate at the highest energies, and moreover, the intermediate µ = 0.5 calculation becomes most accurate in the intermediate range, between i = 1000 and i = 10, 000. With one million states to choose from, it might in principle be the case that the small sampling considered in Table Table 11 happens not to be representative. We therefore next turn our attention to statistics gathered across all 106 states. Tables Table 12–Table 14 indicate the total numbers of energy levels that are converged to within different error tolerances, for the µ = 0.1, µ = 0.5, and µ = 1.0 calculations, respectively. We see the same basic patterns as described above. For example, µ = 0.1 yields 318 states that are converged to within 4 cm−1 , whereas µ = 1.0 yields only one such state. Conversely, µ = 1.0 yields 652,066 levels to within 15 cm−1 —i.e., nearly 2/3 of the desired 106 states—whereas µ = 0.1 yields only 33,233 such states. Although this is in keeping with the expected pattern, in actuality the latter number should probably be even smaller, owing to the aforementioned pseudo-convergence. This would then also place the performance of the µ = 0.5 calculation more in line with expectations; otherwise, a purely convergence-based assessment does not place the µ = 0.5 as the best choice anywhere in the spectrum! In reality, µ = 0.5 does perform best at intermediate energies—for reasons discussed above, and for more 18


Page 18 of 45

Page 19 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




compelling reasons described in Sec. . There is another pattern evident in Tables Table 12–Table 14 that merits discussion. This is the pronounced “efficiency cliff” in performance, observed also in previous PSP calculations. 6,18,19 It can be described this way: for a small increase in N, and for a specific error tolerance, one observes a massive increase in K. We see this very clearly in Table Table 12, where for ∆E = 40 cm−1 , as N increases modestly from around two million to three million, K increases by nearly one thousand fold, from K = 908 to K = 724, 266. The efficiency cliff can also manifest “vertically” as well as horizontally—meaning that for a single calculation, there is a very large drop in efficiency associated with a modest reduction of the error tolerance. This is most evident in Table Table 14 Column IV, where the modest reduction of ∆E from 40 cm−1 to 15 cm−1 decreases K from 652,066 to just three! The efficiency cliff is therefore both “good” and “bad”; it is what provides huge numbers of converged states, but at the same time, it is what limits us from substantially improving the accuracy of those states beyond a certain point.

Hybrid Truncation (HT) dataset All of the above serves to indicate that where large D and K energy level calculations are concerned, we should: (a) take care when interpreting convergence and accuracy; (b) rely on a range of different truncation methods. Both of these ideas come together in the HT procedure of Sec. , which also introduces a third crucially important element: the variational principle. Because the variational principle enables us to rigorously compare energy levels across truncation methods, it provides the unambiguously most accurate results, while at the same time mitigating the vagaries that individual truncation methods might exhibit at one point or another throughout the spectrum (such as pseudo-convergence). Moreover, with regard to individual calculations, “the more the merrier”; with HT, it can never hurt to add new calculations into the mix, even if they are not very accurate. In the present context, for instance, based on the discussion in Sec. , one might be tempted to exclude the µ = 0.5 calculations from the hybridization. In reality though, not only does it not hurt 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



to do so, it actually helps very substantially. Indeed, hybridization is the most rigorous procedure for determining which method is most useful where—because it does not depend on a convergence comparison between basis sets of two (substantially) different sizes. Table Table 15 shows the number of states that contribute to the HT dataset, as culled from the large calculations for each of the three truncation methods (throwing in the small and medium calculations would not change anything). Here, we see a very clear pattern in accord with expectation: the fewest and presumably mostly lowest-lying states, about 15,000 in all, come from µ = 0.1; an intermediate number, about 400,000, come from µ = 0.5; the vast majority, around 5.3 × 106 , come from µ = 1.0. As also indicated in the table, the relative percentages are reliably consistent even across symmetry blocks. Table Table 16, Column III presents frequencies, νi , for selected eigenstates i, as taken from the HT dataset. Note that these represent the most accurate benzene frequencies presented in this paper. In Columns IV and V, we present energy level convergence errors, ∆EiHT , for the medium and small calculations of the same µ value as that from which the given HT state i came from (indicated in Column II). Defined in this fashion, the HT convergence errors for individual states are the same as those for the constituent truncation methods, and so we do not necessarily observe an enormous reduction. The full advantage of HT may not be completely evident in the convergence data, but it is nevertheless very substantial. This is because HT optimizes for accuracy rather than convergence. In any event, Table Table 17 offers a more complete picture of HT convergence, by providing statistical data across all 106 states, in the same format as Tables Table 12–Table 14. From Table Table 17, it is clear that HT is as well converged as the best µ calculation, throughout the entire energy range, as desired.

Conclusions By combining phase space ideas with massive parallelization, we have computed all one million of the dynamically relevant vibrational quantum states for the 12-atom benzene molecule, to a converged accuracy of ∼ 10−1 –101 cm−1 . Moreover, regarding the most accurate frequencies as 20


Page 20 of 45

Page 21 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




presented in Table Table 16 Column III, we have done so using only 24 calculations in all. Such a feat—i.e., combining large D with very large K in an EQD context—would have been unheard of just a few years ago. Yet, despite the advanced computer hardware used, it is the PSP approach itself that has primarily made this possible. No developments in hardware alone can overcome the EQD exponential scaling problem. In contrast, PSP ideas are what have led to the creation of the only EQD methods (PSLB) formally known to circumvent exponential scaling. 16–18 Here, we have extended the utility of PSP ideas still further, by using them to develop an HOB truncation scheme that is not only customized for individual molecular systems, but also tailored for specific regions of the energetic spectrum. Thus, the PSP has lead to both the creation of new methods (PSLB) and the reinventing of older ones (truncated HOB). Using the PSP, also, we have learned from this study that large D, very large K calculations can benefit from a range of different EQD calculations. When combined with the variational principle, this enables the construction of a “hybrid” dataset that is not only optimally accurate throughout the entire spectral range, but also exhibits more stable numerical convergence. The PSP is so useful because it provides simple, clear, intuitive insight, offering answers to a broad range of potential questions. Over the years, it has been used to optimize discrete variable representations, 31,32 explain exponential scaling, 15 and overcome it. 7,16–18 It not only provides optimal EQD basis set representations, it also offers powerful a priori tools for estimating (and improving) their accuracy and convergence scaling. We anticipate that PSP ideas will continue to provide useful insight in new and unanticipated ways—particularly for large D, very large K applications of the sort considered here. For the present benzene study, we have chosen to use massively parallel computers because such resources are available to us—and because these will be the tools of cutting edge computational science, going forward. Also, massively parallel calculations provide many eigenvalues at once. So for large molecules, for which there are a great many dynamically relevant quantum states, such calculations are extremely convenient, as discussed in Sec. . However, it is more than just this. From Sec. , we learn that the convergence of PSP HOB exhibits a pronounced efficiency

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



cliff. 6,18,19 This means that a small increase in basis size can lead to a dramatic increase in the number of states converged to a given accuracy level—creating a strong impetus to work with as large a basis size, N, as possible. Also with this study, one of our goals was to push the limits of our present, simple S WITCH BLADE codes as far as possible. Thus, we sought the largest D and N values for which we could feasibly perform direct diagonalization on a massively parallel computer. That is not to say that our PSP methods require such specialized hardware, however. On the contrary, for somewhat smaller but still highly challenging systems, 4,5,18 hundreds or thousands of dynamically relevant states can be computed utilizing much more modest computational resources. Alternatively, since the Hamiltonian matrix (and even its symmetry blocks) are sparse, it is also possible to combine PSP HOB truncation with iterative methods—thereby extending the realm of applicability of the latter, while simultaneously rendering benzene-like calculations possible on much smaller computing clusters than those employed here. Moreover, one could then tackle substantially larger basis sizes, N. This is an area of future investigation. As discussed in Sec. , however, the downside to the iterative approach is that relatively few states can be computed at a time. Whether iterative or direct eigensolvers are employed, one current drawback of the PSP methods that must be addressed going forward is the accuracy limitation for high-lying states: for large D, one can evidently attain high convergence accuracy or very high efficiency, (K/N), but not both at the same time. Finding ways to simultaneously improve both is therefore an ongoing area of investigation. In particular, early work done by our group resulted in a methodology for dramatically increasing the accuracy of PSLB’s 6 without compromising efficiency. In that study, phase space region operators (PSRO) were employed in order to “round out” the edges of the classical phase space region, so that ρN would better resemble ρK , even in the tunneling region. Though model applications up to D = 4 proved extremely encouraging, challenges for a large D, very large K implementation remain—particular vis-à-vis massive parallelization. Another, slightly less pressing issue going forward is the generalization for arbitrary (i.e., non-

22


Page 22 of 45

Page 23 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




FFP) PESs. Currently, the S WITCH BLADE codes are restricted to quartic FFP’s—although this does not reflect an inherent limitation of the underlying PSP methods themselves. For PSP HOBs, in particular, Gauss-Hermite quadrature 33 is a natural tool for representing arbitrary PESs; moreover, a similar approach has been developed for PSLBs, in a manner that itself avoids exponential scaling. 34 To date, we have not implemented these techniques into S WITCH BLADE, primarily for two reasons. First, for D ≈ 15–30, there are very few non-FFP PESs in existence. Second, the accuracy limitation of most FFP PESs is on the order of tens of cm−1 ’s—i.e., comparable to the worst convergence error that S WITCH BLADE currently provides. As PES quality and PSP accuracy improve, however, the motivation for doing so will increase. We also note, though, that the introduction of quadrature removes the rigorously variational property that is so beneficial for the HT procedure—although how damaging this would be in practice remains to be seen. To summarize: now that we have seen what a straightforward but “brute force” implementation of S WITCH BLADE can achieve, it is time to address the refinements discussed above, head on.

ACKNOWLEDGEMENTS This work was largely supported by a grant from The Robert A. Welch Foundation (D-1523). In addition, both a research grant (CHE-1012662) and a CRIF MU instrumentation grant (CHE0840493) from the National Science Foundation are acknowledged. We also gratefully acknowledge the following entities for providing access and technical support for their respective computing clusters: the Texas Tech University High Performance Computing Center, for use of the Hrothgar facility; the Texas Tech University Department of Chemistry and Biochemistry, for use of the Robinson cluster; the Texas Advanced Computing Center, for use of the Lonestar 4 facility. Calculations presented in this paper were performed using the S WITCH BLADE parallel code.

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



References (1) Meyer, H.-D.; Gatti, F.; Worth, G. A., Eds. Multidimensional Quantum Dynamics: MCTDH Theory and Applications. Wiley-VCH, Weinheim, 2009. (2) Avila, G.; Carrington, Jr., T., Using nonproduct quadrature grids to solve the vibrational Schrödinger equation in 12D. J. Chem. Phys., 2011, 134, 054126. (3) Avila, G.; Carrington, Jr., T., Solving the vibrational Schrödinger equation using bases pruned to include strongly coupled functions and compatible quadratures. J. Chem. Phys., 2012, 137, 174108. (4) Halverson, T.; Poirier, B., Large scale quantum dynamics calculations: Ten thousand quantum states of acetonitrile. Chem. Phys. Lett., 2015, 624, 37–42. (5) Halverson, T.; Poirier, B., Calculation of exact vibrational spectra for P2 O and CH2 NH using a phase space wavelet basis. J. Chem. Phys., 2014, 140, 204112. (6) Lombardini, R.; Poirier, B., Improving the accuracy of Weyl-Heisenberg wavelet and symmetrized Gaussian representations using customized phase-space-region operators. Phys. Rev. E, 2006, 74, 036705. (7) Davis, M. J.; Heller, E. J., Semiclassical Gaussian basis set method for molecular vibrational wave functions. J. Chem. Phys., 1979, 71, 3383–3395. (8) Wilson, K. G. Generalized Wannier Functions. Cornell University preprint, 1987. (9) Daubechies, I., Orthonormal bases of compactly supported wavelets. Commun. Pure Appl. Math., 1988, 49, 909–996. (10) Daubechies, I., The wavelet transform, time-frequency localization and signal analysis. IEEE Trans. Inf. Theory, 1990, 36, 961–1005.

24


Page 24 of 45

Page 25 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




(11) Daubechies, I.; Jaffard, S.; Journé, J., A simple Wilson orthonormal bases with exponential decay. SIAM J. Math. Anal., 1991, 22, 554–573. (12) Daubechies, I. Orthonormal bases of coherent states: The canonical case and the ax + b group. in Coherent States: Past, Present, and Future, pp 103–117, Singapore, 1994. World Scientific. (13) Klauder, J. R.; Sudarshan, E. C. G. Fundamentals of Quantum Optics. W. A. Benjamin, Inc., New York, 1968. (14) Perelemov, A. Generalized Coherent States and Their Applications. Springer-Verlag, Berlin, 1986. (15) Poirier, B., Algebraically Self-Consistent Quasiclassical Approximation on Phase Space. Found. Phys., 2000, 30, 1191–1226. (16) Poirier, B., Using wavelets to extend quantum dynamics calculations to ten or more degrees of freedom. J. Theo. Comput. Chem., 2003, 2, 65–72. (17) Poirier, B.; Salam, A., Quantum dynamics calculations using symmetrized, orthogonal WeylHeisenberg wavelets with a phase space truncation scheme: II. Construction and optimization. J. Chem. Phys., 2004, 121, 1690–1703. (18) Poirier, B.; Salam, A., Quantum dynamics calculations using symmetrized, orthogonal WeylHeisenberg wavelets with a phase space truncation scheme: III. Representations and calculations. J. Chem. Phys., 2004, 121, 1704–1724. (19) Halverson, T.; Poirier, B., Accurate quantum dynamics calculations using symmetrized Gaussians on a doubly dense Von Neumann lattice. J. Chem. Phys., 2012, 137, 224101. (20) Carter, S.; Bowman, J. M.; Sharma, A. R., The ‘MULTIMODE’ approach to ro-vibrational spectroscopy. AIP Conf. Proc., 2012, 1504, 465–466.

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



(21) Maslen, P. E.; Handy, N. C.; Amos, R. D.; Jayatilaka, D., Higher analytic derivatives. IV. Anharmonic effects in the benzene spectrum. J. Chem. Phys., 1992, 97, 4233–4254. (22) Weyl, H., Quantenmechanik und Gruppentheorie. Z. Phys., 1927, 46, 1–46. (23) Wigner, E., On the Quantum Correction For Thermodynamic Equilibrium. Phys. Rev., 1932, 40, 749–759. (24) Hillery, M.; O’Connell, R. F.; Scully, M. O.; Wigner, E. P., Distribution Functions in Physics: Fundamentals. Phys. Rep., 1984, 106, 121–167. (25) Wilson, Jr., E. B., The Normal Modes and Frequencies of Vibration of the Regular Plane Hexagon Model of the Benzene Molecule. Phys. Rev., 1934, 45, 706–714. (26) Willetts, A.; Handy, N. C.; Green, W. H.; Jayatilaka, D., Anharmonic corrections to vibrational transition intensities. J. Phys. Chem., 1990, 94, 5608–5616. (27) Gaw, J. F.; Willetts, A.; Green, W. H.; Handy, N. C. (JAI Press Inc., New York, NY). volume IB. p 169. 1991. (28) Goodman, L.; Ozkabak, A. G.; Thakur, S. N., A benchmark vibrational potential surface: ground-state benzene. J. Phys. Chem, 1991, 95, 9044–9058. (29) Henry, L.; Amat, G., The cubic anharmonic potential function of polyatomic molecules. J. Mol. Spectrosc., 1960, 5, 319–325. (30) Henry, L.; Amat, G., The quartic anharmonic potential function of polyatomic molecules. J. Mol. Spectrosc., 1965, 15, 168–179. (31) Fattal, E.; Baer, R.; Kosloff, R., Phase space approach for optimizing grid representations: The mapped Fourier method. Phys. Rev. E., 1996, 53, 1217–1232. (32) Poirier, B.; Light, J. C., Phase space optimization of quantum representations: Direct product basis sets. J. Chem. Phys., 1999, 111, 4869–4885. 26


Page 26 of 45

Page 27 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




(33) Press et al, W. H. (Cambridge University Press, Cambridge, England). Second edition, 1992. (34) Lombardini, R.; Poirier, B., Rovibrational spectroscopy calculations of neon dimer using a phase space truncated Weyl-Heisenberg wavelet basis. J. Chem. Phys., 2006, 124, 144107.

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 28 of 45


Table 1: Adjusted normal mode harmonic frequencies, ωk , as used in the quartic FFP developed by N. Handy and coworkers. 21 Column I: full D = 30 normal mode labels, for which each doubly degenerate mode receives two labels. Column II: “reduced” (D = 20) normal mode labels, for which each doubly degenerate mode receives a single label. Mode k Reduced 1 1 2 2 3 3 4 4 5 5 6,7 6 8,9 7 10,11 8 12,13 9 14,15 10 16 11 17 12 18 13 19 14 20 15 21,22 16 23,24 17 25,26 18 27,29 19 29,30 20

ωk (cm−1 ) 1008 3208 1390 718 1011 613 3191 1639 1192 866 686 1024 3172 1318 1167 407 989 1058 1512 3191

Symmetry D6h irrep A1g A1g A2g B2g B2g E2g E2g E2g E2g E1g A2u B1u B1u B2u B2u E2u E2u E1u E1u E1u

µ

Table 2: Basis sizes, NEmax , as truncated in accord with Eq. (7), using three different Emax values for each of three different truncation schemes, µ. The Emax values were adjusted in such a way as µ to ensure that small, medium, and large NEmax values are roughly constant across µ. µ

Symmetry Block D6h irreps 1 A1g ⊕ E2gx 2 A2g ⊕ E2gy 3 B1u ⊕ E1ux 4 B2u ⊕ E1uy 5 B1g ⊕ E1gx 6 B2g ⊕ E1gy 7 A1u ⊕ E2ux 8 A2u ⊕ E2uy total —

0.1 N13.28 250063 248580 248947 248584 239179 239881 239505 239897 1954636

µ = 0.1 0.1 N13.88 372485 370870 371390 371054 361028 362210 361787 362491 2933315

0.1 N14.23 719588 717697 717955 717602 707327 709084 708032 709057 5706342

Basis Size, NEmax µ = 0.5 0.5 0.5 0.5 N215.0 N223.0 N237.0 245756 364971 716302 244591 363543 714140 244617 363467 714052 244411 363206 713560 242506 361022 709793 243423 362203 711467 243275 361808 711087 243774 362491 712037 1952353 2902771 5702438

28


1.0 N6600.0 259479 258156 257997 257692 257050 258023 258024 258582 2065003

µ = 1.0 1.0 N6800.0 351713 350166 349825 349458 349274 350428 349979 350637 2801480

1.0 N7278.0 716205 713919 712008 711428 710810 712552 713549 714547 5705018

Page 29 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




large

large

large

Table 3: Frequency values, νi = Ei − E0 , for the first symmetry block of benzene, calculated using µ = 0.1 (Column V). State labels (Column II) were arrived at using energetic ordering and symmetry. Column III is the frequency for the assigned state as obtained via straightforward addition of the harmonic frequencies shown in Table Table 1. Column IV lists second-order perturbation theory values, 21 where available. Columns VI and VII indicate the convergence of the large calculated frequencies, ∆νi = νi − νi . Index State i Label 1 ν6 2 2ν16 3 2ν16 4 ν1 5 ν11 + ν16 6 ν9 7 2ν6 8 2ν6 9 2ν11 10 ν16 + ν17 11 ν16 + ν17 12 ν6 + 2ν16 13 ν6 + 2ν16 14 ν6 + 2ν16 15 2ν4 16 ν4 + ν10 17 ν1 + ν6 18 4ν16 19 4ν16 20 4ν16 21 ν8 22 ν11 + ν17 23 ν6 + ν11 + ν16 24 ν6 + ν11 + ν16 25 ν4 + ν5 26 2ν10 27 2ν10 28 ν6 + ν9 29 ν6 + ν9 30 ν1 + 2ν16 31 ν1 + 2ν16 32 3ν6 33 3ν6 34 ν5 + ν10

Block 1: A1g ⊕ E2gx ωi Ref. 21 νi (cm−1 ) −1 −1 0.1 0.1 (cm ) (cm ) N14.23 N13.88 613 608.0 611.4227 1.241 814 — 797.499 0.1989 814 — 798.6001 0.1937 1008 993.0 997.6235 0.6237 1093 — 1060.814 0.5794 1192 1178.0 1180.374 0.3706 1226 — 1221.237 0.1486 1226 — 1221.338 0.15 1372 — 1317.674 1.296 1396 — 1357.361 0.5605 1396 — 1368.903 0.4642 1427 — 1408.702 2.659 1427 — 1409.768 2.604 1427 — 1409.835 2.58 1436 — 1418.584 0.1918 1584 — 1570.846 0.7995 1621 — 1595.525 0.434 1628 — 1599.05 16.05 1628 — 1599.934 16.86 1628 — 1602.644 45.12 1639 1601.0 1614.455 34.89 1675 — 1615.814 38.43 1706 — 1673.92 6.919 1706 — 1673.927 6.954 1729 — 1691.419 0.6959 1732 — 1733.907 1.466 1732 — 1739.31 1.097 1805 — 1788.659 0.4096 1805 — 1788.955 0.3503 1822 — 1800.083 6.449 1822 — 1801.185 5.726 1839 — 1834.905 3.473 1839 — 1835.206 3.559 1877 — 1852.65 1.3

29


0.1 N13.28 0.4796 -0.2221 -0.2735 -0.1515 0.5612 -0.4409 -0.5442 -0.5395 1.528 0.3564 0.2383 6.248 6.151 6.126 -0.5252 0.3279 -0.1629 15.5 16.5 82.98 71.2 75.77 60.21 61.04 47.13 6.021 15.23 -0.2232 -0.3289 25.01 26.24 5.967 6.101 0.7709


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 30 of 45


Table 4: Frequency values for the second symmetry block of benzene calculated using µ = 0.1 truncation (see Table Table 3 for further details). Index i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

State Label ν6 2ν16 ν11 + ν16 ν9 2ν6 ν3 ν16 + ν17 ν16 + ν17 ν6 + 2ν16 ν6 + 2ν16 ν6 + 2ν16 ν4 + ν10 ν1 + ν6 4ν16 4ν16 ν8 ν11 + ν17 ν6 + ν11 + ν16 ν6 + ν11 + ν16 2ν10 ν6 + ν9 ν6 + ν9 ν1 + 2ν16 3ν6 3ν6 ν5 + ν10

Block 2: A2g ⊕ E2gy ωi Ref. 21 νi (cm−1 ) 0.1 0.1 (cm−1 ) (cm−1 ) N14.23 N13.88 613 608.0 611.4227 1.241 814 — 798.6001 0.1937 1093 — 1060.814 0.5794 1192 1178.0 1180.374 0.3706 1226 — 1221.237 0.1486 1390 1350.0 1352.563 0.3641 1396 — 1357.361 0.5605 1396 — 1375.937 0.2826 1427 — 1408.702 2.659 1427 — 1409.736 2.59 1427 — 1409.768 2.604 1584 — 1570.846 0.7995 1621 — 1595.525 0.434 1628 — 1599.934 15.17 1628 — 1602.644 14.15 1639 1601.0 1614.455 34.89 1675 — 1615.814 38.43 1706 — 1673.811 6.997 1706 — 1673.927 6.954 1732 — 1739.31 1.097 1805 — 1787.628 0.3808 1805 — 1788.955 0.3503 1822 — 1801.185 5.726 1839 — 1834.89 3.475 1839 — 1835.206 3.559 1877 — 1852.65 1.3

30


0.1 N13.28 0.4796 -0.2735 0.5612 -0.4409 -0.5442 -0.1882 0.3564 -0.3505 6.248 6.139 6.151 0.3279 -0.1629 14.61 13.79 71.1 69.84 64.74 66. 15.23 -0.2645 -0.3289 23.91 5.968 6.101 0.7709

Page 31 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




Table 5: Frequency values for the third symmetry block of benzene calculated using µ = 0.1 truncation (see Table Table 3 for further details). Index i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

State Label ν12 ν18 ν4 + ν16 ν10 + ν16 ν10 + ν16 ν4 + ν11 ν5 + ν16 ν19 ν10 + ν11 ν6 + ν12 ν5 + ν11 ν6 + ν18 ν6 + ν18 ν4 + ν17 ν4 + ν6 + ν16 ν4 + ν6 + ν16 ν6 + ν15 ν12 + 2ν16 ν12 + 2ν16 ν10 + ν17 ν10 + ν17 2ν16 + ν18 2ν16 + ν18 2ν16 + ν18 ν6 + ν10 + ν16 ν6 + ν10 + ν16 ν6 + ν10 + ν16 ν6 + ν10 + ν16 ν6 + ν14 ν4 + 3ν16 ν4 + 3ν16 ν15 + 2ν16 ν5 + ν17

Block 3: B1u ⊕ E1ux ωi Ref. 21 νi (cm−1 ) −1 −1 0.1 0.1 (cm ) (cm ) N14.23 N13.88 1024 1010.0 1015.64 0.5037 1058 1037.0 1040.98 0.5263 1125 1108.3 1111.178 0.2046 1273 1243.1 1259.015 0.7022 1273 1243.1 1260.279 0.5853 1404 — 1374.838 0.563 1418 1386.8 1377.769 0.5704 1512 1484.0 1496.231 0.2296 1552 1518.7 1525.564 1.393 1637 1618.1 1624.036 0.1759 1697 — 1632.53 1.359 1671 — 1649.967 0.493 1671 — 1650.055 0.4182 1707 1642.2 1672.044 0.5468 1738 — 1722.69 2.898 1738 — 1722.842 2.901 1780 1756.1 1755.745 0.8737 1838 — 1811.882 3.233 1838 — 1812.909 3.138 1855 1810.4 1827.072 1.323 1855 1810.4 1832.062 1.591 1872 — 1832.857 6.086 1872 — 1833.708 6.206 1872 — 1840.386 3.911 1886 — 1871.566 4.139 1886 — 1872.881 7.33 1886 — 1873.693 6.696 1886 — 1873.96 6.662 1931 1915.6 1917.16 8.177 1939 — 1919.232 29.65 1939 — 1923.571 29.18 1981 — 1941.576 23.21 2000 1955.2 1951.647 16.51

31


0.1 N13.28 -0.3097 -0.2779 -0.3797 0.767 0.4829 0.1013 0.3962 -0.3996 1.219 -0.5759 1.077 -0.1047 -0.1394 -0.05071 5.941 5.921 0.4065 6.587 6.507 0.7704 1.065 11.84 12.11 8.052 6.231 10.8 10.3 10.25 8.591 32.87 34.67 56.97 49.08


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 32 of 45


Table 6: Frequency values for the fourth symmetry block of benzene calculated using µ = 0.1 truncation (see Table Table 3 for further details). Index i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

State Label ν18 ν4 + ν16 ν15 ν10 + ν16 ν10 + ν16 ν14 ν5 + ν16 ν19 ν10 + ν11 ν6 + ν12 ν6 + ν18 ν6 + ν18 ν4 + ν17 ν4 + ν6 + ν16 ν4 + ν6 + ν16 ν6 + ν15 ν12 + 2ν16 ν10 + ν17 ν10 + ν17 2ν16 + ν18 2ν16 + ν18 2ν16 + ν18 ν6 + ν10 + ν16 ν6 + ν10 + ν16 ν6 + ν10 + ν16 ν6 + ν10 + ν16 ν6 + ν14 ν4 + 3ν16 ν4 + 3ν16 ν15 + 2ν16 ν15 + 2ν16 ν5 + ν17 ν4 + ν6 + ν11 ν5 + ν6 + ν16

Block 4: B2u ⊕ E1uy ωi Ref. 21 νi (cm−1 ) −1 −1 0.1 0.1 (cm ) (cm ) N14.23 N13.88 1058 1037.0 1040.98 0.5263 1125 1108.3 1111.178 0.2046 1167 1150.0 1147.751 0.4723 1273 1243.1 1260.279 0.5853 1273 1243.1 1260.781 0.407 1318 1309.0 1315.612 0.2943 1418 1386.8 1377.769 0.5704 1512 1484.0 1496.231 0.2296 1552 1518.7 1525.564 1.393 1637 1618.1 1624.036 0.1759 1671 — 1649.085 0.4198 1671 — 1650.055 0.4058 1707 1642.2 1672.044 0.5468 1738 — 1722.842 2.901 1738 — 1722.936 2.913 1780 1756.1 1755.745 0.8737 1838 — 1812.909 3.138 1855 1810.4 1827.072 1.323 1855 1810.4 1832.857 1.429 1872 — 1833.006 5.937 1872 — 1833.858 5.811 1872 — 1840.386 3.911 1886 — 1871.566 4.139 1886 — 1872.881 7.117 1886 — 1873.205 7.006 1886 — 1873.693 6.696 1931 1915.6 1917.16 8.177 1939 — 1919.498 28.69 1939 — 1923.571 25.31 1981 — 1940.87 11.88 1981 — 1941.576 23.21 2000 1955.2 1951.647 16.86 2017 — 1988.645 7.355 2031 — 1991.223 6.859

32


0.1 N13.28 -0.2779 -0.3797 -0.3226 0.4735 -0.01892 -0.4229 0.3962 -0.3996 1.219 -0.5759 -0.1776 -0.1919 -0.05071 5.921 5.921 0.4065 6.507 0.7704 0.8406 11.69 11.7 8.052 6.231 10.73 10.48 10.3 8.591 32.61 33.77 17.37 56.97 49.28 12.46 43.71

Page 33 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




Table 7: Frequency values for the fifth symmetry block of benzene calculated using µ = 0.1 truncation (see Table Table 3 for further details). Index i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

State Label ν10 ν4 + ν6 ν12 + ν16 ν16 + ν18 ν16 + ν18 ν6 + ν10 ν6 + ν10 ν4 + 2ν16 ν15 + ν16 ν5 + ν6 ν10 + 2ν16 ν10 + 2ν16 ν10 + 2ν16 ν14 + ν16 ν11 + ν18 ν4 + ν11 + ν16 ν5 + 2.ν16 ν11 + ν15 ν1 + ν10 ν4 + ν9 ν16 + ν19 ν16 + ν19 ν4 + 2ν6 ν10 + ν11 + ν16 ν10 + ν11 + ν16 ν12 + ν17 ν11 + ν14 ν6 + ν12 + ν16 ν6 + ν12 + ν16

Block 5: B1g ⊕ E1gx ωi Ref. 21 νi (cm−1 ) 0.1 0.1 (cm−1 ) (cm−1 ) N14.23 N13.88 866 847.0 868.9106 1.193 1331 — 1319.55 0.1728 1431 — 1412.568 0.2031 1465 — 1435.124 0.5187 1465 — 1440.259 0.3492 1479 — 1477.829 0.5258 1479 — 1477.857 0.5405 1532 — 1513.589 3.199 1574 — 1541.753 0.9968 1624. — 1594.471 0.4994 1680 — 1652.792 5.525 1680 — 1654.487 7.068 1680 — 1657.841 6.937 1725 — 1706.712 1.063 1744 — 1715.345 1.983 1811 — 1773.376 7.173 1825 — 1775.103 6.783 1853 — 1814.526 1.763 1874 — 1864.531 1.123 1910. — 1886.809 0.3736 1919 — 1891.143 0.4686 1919 — 1893.704 0.7949 1944 — 1918.898 9.941 1959 — 1919.119 14.17 1959 — 1932.498 4.292 2013 — 1976.261 0.5489 2004 — 1978.52 3.836 2044 — 2001.149 0.998 2044 — 2002.968 0.9453

33


0.1 N13.28 0.4176 -0.5377 -0.4628 0.1307 -0.2709 -0.00129 0.01738 6.113 0.8674 -0.04323 8.152 12.02 11.85 0.706 2.846 11.52 11.02 1.451 0.6467 -0.3275 -0.1496 0.2913 12.6 18.45 5.95 -0.1406 4.562 0.432 0.3693


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 34 of 45


Table 8: Frequency values for the sixth symmetry block of benzene calculated using µ = 0.1 truncation (see Table Table 3 for further details). Index i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

State Label ν4 ν10 ν5 ν4 + ν6 ν12 + ν16 ν16 + ν18 ν16 + ν18 ν6 + ν10 ν6 + ν10 ν4 + 2ν16 ν4 + 2ν16 ν15 + ν16 ν5 + ν6 ν10 + 2ν16 ν10 + 2ν16 ν10 + 2ν16 ν11 + ν12 ν1 + ν4 ν14 + ν16 ν11 + ν18 ν4 + ν11 + ν16 ν5 + 2.ν16 ν5 + 2.ν16 ν1 + ν10 ν4 + ν9 ν16 + ν19 ν16 + ν19 ν4 + 2ν6 ν4 + 2ν6 ν10 + ν11 + ν16 ν10 + ν11 + ν16 ν12 + ν17 ν1 + ν5 ν6 + ν12 + ν16 ν6 + ν12 + ν16

Block 6: B2g ⊕ E1gy ωi Ref. 21 νi (cm−1 ) −1 −1 0.1 0.1 (cm ) (cm ) N14.23 N13.88 718 707.0 710.7318 1.013 866 847.0 868.9106 1.193 1011 990.0 985.8294 0.7636 1331 — 1319.55 0.1728 1431 — 1412.568 0.2031 1465 — 1435.124 0.5187 1465 — 1435.197 0.5004 1479 — 1477.829 0.5258 1479 — 1477.98 0.5059 1532 — 1512.279 3.314 1532 — 1513.589 3.199 1574 — 1541.753 0.9968 1624. — 1594.471 0.4994 1680 — 1652.792 5.525 1680 — 1653.244 7.942 1680 — 1654.487 7.068 1710 — 1679.784 0.5971 1726 — 1704.671 0.5737 1725 — 1706.712 1.063 1744 — 1715.345 1.983 1811 — 1773.376 7.173 1825 — 1774.597 6.64 1825 — 1775.102 6.783 1874 — 1864.531 1.123 1910 — 1886.809 0.3736 1919 — 1893.704 0.7949 1919 — 1894.277 0.7543 1944 — 1914.807 15.8 1944 — 1918.898 14.39 1959 — 1932.498 4.292 1959 — 1932.652 4.337 2013 — 1976.261 0.5489 2019 — 1980.391 1.055 2044 — 2000.981 1.096 2044 — 2001.149 0.998

34


0.1 N13.28 0.2382 0.4176 -0.03369 -0.5377 -0.4628 0.1307 0.1253 -0.00129 -0.02907 6.251 6.113 0.8674 -0.04323 8.152 13.01 12.02 -0.01024 -0.0298 0.706 2.846 11.52 11.05 11.02 0.6467 -0.3275 0.2913 0.228 19.98 18.67 5.95 5.997 -0.1406 0.5496 0.5538 0.432

Page 35 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




Table 9: Frequency values for the seventh symmetry block of benzene calculated using µ = 0.1 truncation (see Table Table 3 for further details). Index i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

State Label ν16 ν17 ν6 + ν16 ν6 + ν16 3ν16 3ν16 ν6 + ν11 ν1 + ν16 ν11 + 2.ν16 ν9 + ν16 ν9 + ν16 ν6 + ν17 ν6 + ν17 2ν6 + ν16 2ν6 + ν16 2ν6 + ν16 ν4 + ν18 2ν11 + ν16 ν3 + ν16 2ν16 + ν17 2ν16 + ν17 2ν16 + ν17 ν6 + 3ν16 ν6 + 3ν16 ν6 + 3ν16 ν6 + 3ν16 2ν4 + ν16 ν9 + ν11 ν4 + ν15 ν10 + ν12 2ν6 + ν11

Block 7: A1u ⊕ E2ux ωi Ref. 21 νi (cm−1 ) −1 −1 0.1 0.1 0.1 (cm ) (cm ) N14.23 N13.88 N13.28 407 398.0 399.4554 2.215 1.637 989 967.0 964.0127 0.8331 0.03661 1020 1006.4 1009.5 0.1543 -0.4344 1020 1006.4 1009.609 0.1545 -0.4301 1221 — 1195.588 3.747 7.246 1221 — 1197.886 3.658 6.937 1299 — 1277.755 0.4754 0.03375 1415 — 1396.802 0.6489 0.3054 1500 — 1456.891 8.104 12.97 1599 — 1570.89 0.4956 -0.02153 1599 — 1572.486 0.4904 -0.05011 1602 — 1577.013 0.333 -0.3081 1602 — 1577.892 0.4404 -0.09615 1633 — 1621.718 2.788 6.25 1633 — 1621.838 2.783 6.205 1633 — 1622.092 2.826 6.254 1776 — 1707.016 16.43 22.03 1779 — 1744.988 0.536 -0.005982 1797 — 1748.287 2.824 4.196 1803 — 1755.565 6.336 10.82 1803 — 1766.872 5.869 9.716 1803 — 1776.905 3.393 6.442 1834 — 1811.556 16.14 18.39 1834 — 1811.721 36.53 36.02 1834 — 1813.19 41.41 40.91 1834 — 1813.441 48.38 67.42 1843 — 1824.234 37.78 77.74 1878 — 1847.195 17.68 59.25 1885 — 1853.7 11.44 54.05 1890 — 1880.929 0.6125 71.19 1912 — 1891.913 6.562 60.49

35



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Table 10: Frequency values for the eighth symmetry block of benzene calculated using µ = 0.1 truncation (see Table Table 3 for further details). Index i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

State Label ν16 ν11 ν17 ν6 + ν16 ν6 + ν16 3ν16 3ν16 ν6 + ν11 ν1 + ν16 ν11 + 2ν16 ν11 + 2ν16 ν9 + ν16 ν9 + ν16 ν6 + ν17 ν6 + ν17 2ν6 + ν16 2ν6 + ν16 2ν6 + ν16 ν1 + ν11 ν4 + ν18 ν4 + ν12 2ν11 + ν16 ν3 + ν16 2ν16 + ν17 2ν16 + ν17 2ν16 + ν17 ν6 + 3ν16 ν6 + 3ν16 ν6 + 3ν16 ν6 + 3ν16 2ν4 + ν16 ν9 + ν11 ν10 + ν12 2ν6 + ν11 2ν6 + ν11

Block 8: A2u ⊕ E2uy ωi Ref. 21 −1 0.1 (cm ) (cm−1 ) N14.23 407 398.0 399.4554 686 673.8 666.9294 989 967.0 964.0127 1020 1006.4 1009.351 1020 1006.4 1009.500 1221 — 1195.588 1221 — 1197.557 1299 — 1277.755 1415 — 1396.802 1500 — 1455.491 1500 — 1456.891 1599 — 1572.486 1599 — 1572.916 1602 — 1577.892 1602 — 1578.656 1633 — 1621.718 1633 — 1621.844 1633 — 1622.092 1694 1665.7 1663.975 1776 — 1707.016 1742 1717.2 1722.742 1779 — 1744.988 1797 — 1748.287 1803 — 1754.645 1803 — 1766.872 1803 — 1776.905 1834 — 1811.395 1834 — 1811.556 1834 — 1813.190 1834 — 1813.441 1843 — 1824.234 1878 — 1847.195 1890 — 1880.929 1912 — 1891.913 1912 — 1892.034

36


νi (cm−1 ) 0.1 0.1 N13.88 N13.28 2.215 1.637 2.354 1.633 0.8331 0.03661 0.1543 -0.4368 0.1543 -0.4344 3.747 7.246 3.704 7.02 0.4754 0.03375 0.6489 0.3054 8.338 13.31 8.104 12.97 0.4904 -0.05011 0.4695 -0.0832 0.4404 -0.09615 0.4416 -0.1146 2.788 6.25 2.785 6.208 2.826 6.254 1.133 0.7216 15.92 15.16 0.7007 6.306 0.536 -0.005982 2.824 4.196 6.64 11.14 5.869 9.716 3.393 6.442 16.3 18.55 36.69 36.19 48.46 67.67 48.38 88.53 40.64 77.97 17.94 60.55 0.6125 27.01 6.562 59.95 6.639 60.08

Page 36 of 45

Page 37 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




large

large

large

Table 11: Comparison of selected frequencies, νi = Ei − E0 , across all eight symmetry blocks of benzene, calculated using µ = 0.1 (Column II), µ = 0.5 (Column V), and µ = 1.0 (Column VIII). The selected states are indicated in energetic order, in Column I. The remaining large columns indicate the convergence of the calculated energy levels, ∆Ei = Ei − Ei . Index i 1 4 7 10 40 70 100 400 700 1000 4000 7000 10000 40000 70000 100000 300000 700000 1000000

0.1 N14.23

399.455 611.423 797.499 868.911 1260.781 1455.491 1599.050 2112.064 2361.765 2513.700 3178.729 3462.713 3651.131 4426.778 4760.956 4982.088 5750.939 6486.772 6856.688

µ = 0.1 0.1 0.1 N13.88 N13.28 2.336 2.619 1.362 1.462 0.320 0.761 1.314 1.400 0.528 0.964 8.460 14.289 16.176 16.478 8.189 46.248 2.303 18.979 11.088 43.624 26.040 56.635 46.806 69.823 41.122 64.949 15.343 50.652 14.108 79.757 23.657 128.802 35.487 439.590 45.343 1105.212 35.439 1650.181

Frequency, νi (cm−1 ) µ = 0.5 0.5 0.5 0.5 N237.0 N223.0 N215.0 400.355 0.546 2.450 611.400 2.515 2.688 797.195 0.914 2.290 869.589 2.608 2.983 1261.452 3.389 55.772 1457.577 54.975 56.639 1617.696 27.934 29.497 2142.069 4.208 5.812 2383.363 10.167 25.192 2529.586 40.731 59.241 3170.714 48.122 62.044 3456.399 43.820 57.962 3648.295 41.617 53.822 4397.544 17.645 29.395 4710.033 24.803 43.205 4916.992 30.431 49.208 5600.084 40.238 65.211 6168.371 70.913 139.941 6427.678 110.588 271.758

1.0 N7278.0

397.065 607.672 793.369 870.514 1272.950 1472.884 1603.272 2131.168 2391.455 2546.907 3190.850 3462.089 3634.691 4377.892 4694.238 4897.594 5548.106 6094.962 6338.072

µ = 1.0 1.0 N6800.0 1.655 5.107 8.771 9.712 9.437 9.086 12.724 41.138 59.276 27.290 7.132 13.832 26.715 16.146 9.778 9.151 25.048 62.954 57.618

Table 12: Number of converged energy levels, K, across all eight symmetry blocks of benzene, calculated using µ = 0.1, for two different basis sizes (Columns II and IV). Efficiencies, (K/N), are also indicated (Columns III and V). An energy level Ei is considered “converged” if ∆Ei = large Ei − Ei is less than the accuracy threshold listed in Column I. Convergence Accuracy (cm−1 ) 150 40 15 4 0.15

0.1 N13.28 K K/N 112563 0.058 908 0.00046 316 0.00016 137 0.000070 0 0.0

37

0.1 N13.88

K 1410487 724266 33233 318 1


K/N 0.48 0.25 0.011 0.00011 0.0000034

1.0 N6600.0 5.556 8.535 11.688 14.818 11.483 19.336 30.400 48.429 67.785 30.133 16.625 25.814 37.856 18.903 13.593 15.842 36.939 70.482 62.023


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Page 38 of 45


Table 13: Number of converged energy levels, K, across all eight symmetry blocks of benzene, calculated using µ = 0.5, for two different basis sizes (Columns II and IV). Efficiencies, (K/N), are also indicated (Columns III and V). An energy level Ei is considered “converged” if ∆Ei = large Ei − Ei is less than the accuracy threshold listed in Column I. 0.5 N215.0 K K/N 728011 0.37 47635 0.024 321 0.00016 26 0.000013 0 0.0

Convergence Accuracy (cm−1 ) 150 40 15 4 0.15

0.5 N223.0

K K/N 1295900 0.45 268581 0.093 5953 0.0020 116 0.000040 0 0.0

Table 14: Number of converged energy levels, K, across all eight symmetry blocks of benzene, calculated using µ = 1.0, for two different basis sizes (Columns II and IV). Efficiencies, (K/N), are also indicated (Columns III and V). An energy level Ei is considered “converged” if ∆Ei = large Ei − Ei is less than the accuracy threshold listed in Column I. Convergence Accuracy (cm−1 ) 150 40 15 4 0.15

1.0 N6600.0 K K/N 1983268 0.96 974760 0.47 48410 0.023 1 0.0000048 0 0.0

1.0 N6800.0 K K/N 2688080 0.96 1843956 0.66 652066 0.23 3 0.000011 0 0.0

Table 15: Total number of energy levels, N µ , contributed to the HT dataset from µ = 0.1 (Column IV), µ = 0.5 (Column V), and µ = 1.0 (Column VI), for each of the eight symmetry blocks of benzene. The total HT dataset size, N HT is indicated in Column III for each block. Totals across all eight symmetry blocks are listed in the last row. Symmetry Block D6h irreps 1 A1g ⊕ E2gx 2 A2g ⊕ E2gy 3 B1u ⊕ E1ux 4 B2u ⊕ E1uy 5 B1g ⊕ E1gx 6 B2g ⊕ E1gy 7 A1u ⊕ E2ux 8 A2u ⊕ E2uy total —

N HT 716205 713919 712008 711428 707327 709084 708032 709057 5687060 38

µ = 0.1 1566 1498 2244 2204 2171 2215 1613 1548 15059

Nµ µ = 0.5 49511 49053 48405 47852 43976 43993 48174 48357 379321


µ = 1.0 665128 663368 661359 661372 661180 662876 658245 659152 5292680

Page 39 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




Table 16: Selected frequencies, νiHT = EiHT − E0HT , across all eight symmetry blocks of benzene, from the HT dataset (Column III). The selected states are indicated in energetic order, in Column I. Column II indicates which µ calculation the energy level comes from. The remaining columns µ indicate the convergence of the calculated energy levels, ∆EiHT = Ei − EiHT , with Column IV for the medium calculation, and Column V for the small calculation. Index Truncation i µ 1 0.1 4 0.1 7 0.1 10 0.1 40 0.1 70 0.1 100 0.1 400 0.1 700 0.1 1000 0.1 4000 0.5 7000 0.1 10000 0.1 40000 0.5 70000 0.5 100000 0.5 300000 0.5 700000 1.0 1000000 1.0

Frequency, νiHT (cm−1 ) HT HT HT Nlarge Nmedium Nsmall 399.455 2.336 2.619 611.423 1.362 1.462 797.499 0.320 0.761 868.911 1.314 1.400 1260.781 0.528 0.964 1455.491 8.460 14.289 1599.050 16.176 16.478 2112.064 8.189 46.248 2361.765 2.303 18.979 2513.700 11.088 43.624 3174.798 46.883 60.979 3460.260 41.366 63.069 3650.437 46.649 64.415 4401.744 17.645 29.395 4714.234 24.803 43.205 4921.193 30.431 49.208 5604.285 40.238 65.211 6159.057 62.954 70.482 6402.168 57.618 62.023

Table 17: Number of converged energy levels, K, across all eight symmetry blocks of benzene, calculated using the HT dataset, for two different basis sizes (Columns II and IV). Efficiencies, (K/N), are also indicated (Columns III and V). An energy level Ei is considered “converged” if µ ∆EiHT = Ei − EiHT is less than the accuracy threshold listed in Column I. Convergence Accuracy (cm−1 ) 150 40 15 4 0.15

HT Nsmall

K 1937837 682272 318 137 0

K/N 0.99 0.35 0.00016 0.000070 0.0

39

HT Nmedium K K/N 2688080 0.96 1835631 0.66 534638 0.19 316 0.0000036 0 0.0



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Figure 1: Position space projection of the ρK Wigner PDF for the density operator of the K = 20 lowest lying eigenstates of the D = 2 h(2) Hamiltonian of Sec. : (a) 3D plot; (b) contour plot.

40


Page 40 of 45

Page 41 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




Figure 2: Position space projection of the ρN Wigner PDF for the density operator of the N = 21 D = 2 representational HOB functions, truncated using PT: (a) 3D plot; (b) contour plot.

Figure 3: Position space projection of the ρN Wigner PDF for the density operator of the N = 20 D = 2 representational HOB functions, truncated using FT: (a) 3D plot; (b) contour plot.

41



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Figure 4: Contour plot of the position space projection of the ρK Wigner PDF for the density operator of the K = 5 lowest lying eigenstates of the D = 2 h(2) Hamiltonian of Sec. . The tunneling region is emphasized.

Figure 5: Contour plot of the position space projection of the ρN Wigner PDF for the density operator of the N = 45 D = 2 representational HOB functions, truncated using PT.

42


Page 42 of 45

Page 43 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




Figure 6: Contour plot of the position space projection of the ρN Wigner PDF for the density operator of the N = 45 D = 2 representational HOB functions, truncated using FT.

43



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 44 of 45

Page 45 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC graphic 324x305mm (72 x 72 DPI)


One Million Quantum States of Benzene - The Journal of Physical

Recommend Documents