9.11.1 Theory

Even in cases where the Born-Oppenheimer separation is valid, solving the electronic Schrödinger equation may only be half the battle. The remainder involves the solution of the nuclear Schrödinger equation for its resulting eigenvalues and eigenfunctions. This half is typically treated by the harmonic approximation at critical points, but anharmonicity, tunneling, and low-frequency (“floppy”) motions can lead to extremely delocalized nuclear distributions, particularly for protons and for non-covalent interactions.

While the Born-Oppenheimer separation allows for a local solution of the electronic problem (in nuclear space), the nuclear half of the Schrödinger equation is entirely non-local and requires the computation of potential energy surfaces over large regions of configuration space. Grid-based methods, therefore, scale exponentially with the number of degrees of freedom, and are quickly rendered useless for all but very small molecules.

For equilibrium thermal distributions, the path integral (PI) formalism provides both an elegant and computationally feasible alternative. The equilibrium partition function can be written as a trace of the thermal, configuration-space density matrix,

The density matrix at inverse temperature $\beta ={(k_{B}T)}^{-1}$ is defined by the last equality. Evaluating the integrals in Eq. (9.63) still requires computing eigenstates of $\widehat{H}$, which is generally intractable. Inserting $N-1$ resolutions of the identity, however, one obtains

Here, the density matrices appear at an inverse temperature $\beta /N$ that corresponds to multiplying the actual temperature $T$ by a factor of $N$.

The high-temperature form of the density matrix can be expressed as

which becomes exact as $T\to \mathrm{\infty }$ (a limit in which quantum mechanics converges to classical mechanics), or in other words as $\beta \to 0$ or $N\to \mathrm{\infty }$. Using $N$ time slices, the partition function is therefore converted into the form

with the implied cyclic condition $x_{N+1}\equiv x_1$. Here, $V(x)$ is the potential function on which the “beads” move, which is the electronic potential generated by Q-Chem.

Equation 9.66 has the form

where the form of the effective potential $V_{\mathrm{eff}}$ is evident from the integrand in Eq. (9.66). Equation (9.67) reveals that the path-integral formulation of the quantum partition function affords a classical configurational integral for the partition function, albeit in an extended-dimensional space The effective potential describes a classical “ring polymer” with $N$ beads, wherein neighboring beads are coupled by harmonic potentials that arise from the quantum nature of the kinetic energy. The exponentially-scaling, non-local nuclear quantum mechanics problem has therefore been mapped onto an entirely classical problem, which is amenable to standard treatments of configuration sampling. These methods typically involve molecular dynamics or Monte Carlo sampling. Importantly, the number of extended degrees of freedom, $N$, is reasonably small when the temperature is not too low: room-temperature systems involving hydrogen atoms typically are converged using roughly $N\approx 30$ beads. Therefore, fully quantum-mechanical nuclear distributions can be obtained at a cost only roughly 30 times that of a classical AIMD simulation. Path integral Monte Carlo (PIMC) is activated by setting JOBTYPE = PIMC.

The single-bead $(N=1$) limit of the equations above is simply classical configuration sampling. When the temperature (controlled by the PIMC_TEMP keyword) is high, or where only heavy atoms are involved, the classical limit is often appropriate. The path integral machinery (with a single “bead”) may be used to perform classical Boltzmann sampling. In this case, the partition function is simply

and this is what is ordinarily done in an AIMD simulation. Use of additional beads incorporates more quantum-mechanical delocalization, at a cost of roughly $N$ times that of the classical AIMD simulation, and this is the primary input variable in a PI simulation. It is controlled by the keyword PIMC_NBEADSPERATOM. The ratio of the inverse temperature to beads ($\beta /N$) dictates convergence with respect to the number of beads, so as the temperature is lowered, a concomitant increase in the number of beads is required.

Integration over configuration space is performed by Metropolis Monte Carlo (MC). The number of MC steps is controlled by the PIMC_MCMAX keyword and should typically be $\gtrsim {10}^5$, depending on the desired level of statistical convergence. A warm-up run, in which the PI ring polymer is allowed to equilibrate without accumulating statistics, can be performed using the PIMC_WARMUP_MCMAX keyword.

As in AIMD simulations, the main results of PIMC jobs in Q-Chem are not in the job output file but are instead output to ($QCSCRATCH/PIMC in the user’s scratch directory, thus PIMC jobs should always be run with the -save option. The output files do contain some useful information, however, including a basic data analysis of the simulation. Average energies (thermodynamic estimator), bond lengths (less than 5 Å), bond length standard deviations and errors are printed at the end of the output file. The $QCSCRATCH/PIMC directory additionally contains the following files:

• •

  BondAves: running average of bond lengths for convergence testing.

• •

  BondBins: normalized distribution of significant bond lengths, binned within 5 standard deviations of the average bond length.

• •

  ChainCarts: human-readable file of configuration coordinates, likely to be used for further, external statistical analysis. This file can get quite large, so be sure to provide enough scratch space!

• •

  ChainView.xyz: Cartesian-formatted file for viewing the ring-polymer sampling in an external visualization program. (The sampling is performed such that the center of mass of the ring polymer system remains centered.)

• •

  Vcorr: potential correlation function for the assessment of statistical correlations in the sampling.

In each of the above files, the first few lines contain a description of how the data are arranged.

One of the unfortunate rites of passage in PIMC usage is the realization of the ramifications of the stiff bead-bead interactions as convergence (with respect to $N$) is approached. Nearing convergence—where quantum mechanical results are correct—the length of statistical correlations grows enormously, and special sampling techniques are required to avoid long (or non-convergent) simulations. Cartesian displacements or normal-mode displacements of the ring polymer lead to this severe stiffening. While both of these naïve sampling schemes are available in Q-Chem, they are not recommended. Rather, the free-particle (harmonic bead-coupling) terms in the path integral action can be sampled directly. Several schemes are available for this purpose. Q-Chem currently adopts the simplest of these options, Levy flights. An $n$-bead segment (with ) of the ring polymer is chosen at random, with the length $n$ controlled by the PIMC_SNIP_LENGTH keyword. Between the endpoints of this segment, a free-particle path is generated by a Levy construction, which exactly samples the free-particle part of the action. Subsequent Metropolis testing of the resulting potential term—for which only the potential on the moved beads is required—then dictates acceptance.

Two measures of the sampling efficiency are provided in the job output file. The lifetime of the potential auto-correlation function $⟨V_0{V}_{\tau }⟩$ is provided in terms of the number of MC steps, $\tau$. This number indicates the number of configurations that are statically correlated. Similarly, the mean-square displacement between MC configurations is also provided. Maximizing this number and/or minimizing the statistical lifetime leads to efficient sampling. Note that the optimally efficient acceptance rate may not be 50% in MC simulations. In Levy flights, the only variable controlling acceptance and sampling efficiency is the length of the snippet. The statistical efficiency can be obtained from relatively short runs, during which the length of the Levy snippet should be optimized by the user.