Q-Chem 5.1 User’s Manual

10.8 Ab initio Path Integrals

10.8.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,

  \begin{equation} \label{eq:primitive_ Z} Z = \mbox{tr}\bigl (e^{-\beta \hat{H}}\bigr ) =\int dx \; \bigl \langle x \bigl | e^{-\beta \hat{H}} \bigr | x \bigr \rangle =\int dx \; \rho (x,x;\beta ) \; . \end{equation}   (10.19)

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

  \begin{equation} Z = \int dx^{}_1 \int dx^{}_2 \cdots \int dx^{}_ N\ \rho \left(x_1,x_2;\frac{\beta }{N}\right) \rho \left(x_2,x_3;\frac{\beta }{N}\right) \cdots \rho \left(x_ N,x_1;\frac{\beta }{N}\right) \; . \end{equation}   (10.20)

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

  \begin{equation} \rho \left(x,x’;\frac{\beta }{N}\right) = \left(\frac{mN}{2\pi \beta \hbar ^2}\right)^{\! 1/2} \exp \left\{ -\left(\frac{mN}{2\beta \hbar ^2}\right) (x-x’)^2 -\left(\frac{\beta }{2N}\right)\Bigl [V(x)+V(x’) \Bigr ] \right\} \end{equation}   (10.21)

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

  \begin{equation} \label{eq:PI_ partition_ fn} Z = \left(\frac{mN}{2\pi \beta \hbar ^2}\right)^{\! N/2} \int dx^{}_1 \int dx^{}_2 \cdots \int dx^{}_ N \; \exp \left\{ -\frac{\beta }{N} \left[ \frac{mN^2}{2\beta ^2\hbar ^2}\sum _{i=1}^ N \left(x_ i - x_{i+1}\right)^2 + \sum _{i=1}^ N V(x_ i)\right] \right\} \; , \end{equation}   (10.22)

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  has the form

  \begin{equation} \label{eq:PI_ equiv_ to_ classical} Z \propto \int e^{-\beta V_{\rm eff}} \; , \end{equation}   (10.23)

where the form of the effective potential $V_{\rm eff}$ is evident from the integrand in Eq. eq:PI_partition_fn. Equation eq:PI_equiv_to_classical 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

  \begin{equation} Z = \int dx\ e^{-\beta V(x)} \end{equation}   (10.24)

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:

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 naive 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 $n<N$) 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 $\langle V_0V_{\tau }\rangle $ 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.

10.8.2 Job Control and Examples

PIMC_NBEADSPERATOM

Number of path integral time slices (“beads”) used on each atom of a PIMC simulation.


TYPE:

INTEGER


DEFAULT:

None.


OPTIONS:

1

Perform classical Boltzmann sampling.

$>$1

Perform quantum-mechanical path integral sampling.


RECOMMENDATION:

This variable controls the inherent convergence of the path integral simulation. The one-bead limit represents classical sampling and the infinite-bead limit represents exact quantum-mechanical sampling. Using 32 beads is reasonably converged for room-temperature simulations of molecular systems.


PIMC_TEMP

Temperature, in Kelvin (K), of path integral simulations.


TYPE:

INTEGER


DEFAULT:

None.


OPTIONS:
 

User-specified number of Kelvin for PIMC or classical MC simulations.


RECOMMENDATION:

None.


PIMC_MCMAX

Number of Monte Carlo steps to sample.


TYPE:

INTEGER


DEFAULT:

None.


OPTIONS:
 

User-specified number of steps to sample.


RECOMMENDATION:

This variable dictates the statistical convergence of MC/PIMC simulations. For converged simulations at least $10^5$ steps is recommended.


PIMC_WARMUP_MCMAX

Number of Monte Carlo steps to sample during an equilibration period of MC/PIMC simulations.


TYPE:

INTEGER


DEFAULT:

None.


OPTIONS:
 

User-specified number of steps to sample.


RECOMMENDATION:

Use this variable to equilibrate the molecule/ring polymer before collecting production statistics. Usually a short run of roughly 10% of PIMC_MCMAX is sufficient.


PIMC_MOVETYPE

Selects the type of displacements used in MC/PIMC simulations.


TYPE:

INTEGER


DEFAULT:

0


OPTIONS:

0

Cartesian displacements of all beads, with occasional (1%) center-of-mass moves.

1

Normal-mode displacements of all modes, with occasional (1%) center-of-mass moves.

2

Levy flights without center-of-mass moves.


RECOMMENDATION:

Except for classical sampling (MC) or small bead-number quantum sampling (PIMC), Levy flights should be used. For Cartesian and normal-mode moves, the maximum displacement is adjusted during the warm-up run to the desired acceptance rate (controlled by PIMC_ACCEPT_RATE). For Levy flights, the acceptance is solely controlled by PIMC_SNIP_LENGTH.


PIMC_ACCEPT_RATE

Acceptance rate for MC/PIMC simulations when Cartesian or normal-mode displacements are used.


TYPE:

INTEGER


DEFAULT:

None


OPTIONS:

$0 < n < 100$

User-specified rate, given as a whole-number percentage.


RECOMMENDATION:

Choose acceptance rate to maximize sampling efficiency, which is typically signified by the mean-square displacement (printed in the job output). Note that the maximum displacement is adjusted during the warm-up run to achieve roughly this acceptance rate.


PIMC_SNIP_LENGTH

Number of “beads” to use in the Levy flight movement of the ring polymer.


TYPE:

INTEGER


DEFAULT:

None


OPTIONS:

$3 \leq n \leq \mbox{{\small PIMC\_ NBEADSPERATOM}}$

User-specified length of snippet.


RECOMMENDATION:

Choose the snip length to maximize sampling efficiency. The efficiency can be estimated by the mean-square displacement between configurations, printed at the end of the output file. This efficiency will typically, however, be a trade-off between the mean-square displacement (length of statistical correlations) and the number of beads moved. Only the moved beads require recomputing the potential, i.e., a call to Q-Chem for the electronic energy. (Note that the endpoints of the snippet remain fixed during a single move, so $n-2$ beads are actually moved for a snip length of $n$. For 1 or 2 beads in the simulation, Cartesian moves should be used instead.)


Example 10.239  Path integral Monte Carlo simulation of H$_2$ at room temperature

$molecule
   0  1
   H
   H 1 0.75
$end

$rem
   JOBTYPE              pimc
   METHOD               hf   
   BASIS                sto-3g
   PIMC_TEMP            298
   PIMC_NBEADSPERATOM   32
   PIMC_WARMUP_MCMAX    10000     !Equilibration run
   PIMC_MCMAX           100000    !Production run
   PIMC_MOVETYPE        2         !Levy flights
   PIMC_SNIP_LENGTH     10        !Moves 8 beads per MC step (10-endpts) 
$end

Example 10.240  Classical Monte Carlo simulation of a water molecule at 500K

$molecule
   0  1
   H
   O 1 1.0
   H 2 1.0 1 104.5 
$end

$rem
   JOBTYPE              pimc
   METHOD               rimp2
   BASIS                cc-pvdz
   AUX_BASIS            rimp2-cc-pvdz
   PIMC_TEMP            500
   PIMC_NBEADSPERATOM   1       !1 bead is classical sampling
   PIMC_WARMUP_MCMAX    10000   !Equilibration run
   PIMC_MCMAX           100000  !Production run
   PIMC_MOVETYPE        0       !Cartesian displacements (ok for 1 bead)
   PIMC_ACCEPT_RATE     40      !During warm-up, adjusts step size to 40% acceptance
$end