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,

(9.18) |

The density matrix at inverse temperature is defined by the last equality. Evaluating the integrals in Eq. eq:primitive_Z still requires computing eigenstates of , which is generally intractable. Inserting resolutions of the identity, however, one obtains

(9.19) |

Here, the density matrices appear at an inverse temperature that corresponds to multiplying the actual temperature by a factor of .

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

(9.20) |

which becomes exact as (a limit in which quantum mechanics converges to classical mechanics), or in other words as or . Using time slices, the partition function is therefore converted into the form

(9.21) |

with the implied cyclic condition . Here, is the potential function on which the “beads” move, which is the electronic potential generated by Q-Chem.

(9.22) |

where the form of the effective potential 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 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, , is reasonably small when the temperature is not too low: room-temperature systems involving hydrogen atoms typically are converged using roughly 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 JOBTYP to PIMC.

The single-bead ) 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

(9.23) |

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 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 () 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 , depending on the desired level of statistical convergence. A warmup 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 ) 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 utilizes the simplest of these options, Levy flights. An -bead segment (with ) of the ring polymer is chosen at random, with the length 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 autocorrelation function is provided in terms of the number of MC steps, . 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.

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 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 utilized. For Cartesian and normal-mode moves, the maximum displacement is adjusted during the warmup 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 utilized.

TYPE:

INTEGER

DEFAULT:

None

OPTIONS:

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 warmup 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:

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 beads are actually moved for a snip length of . For 1 or 2 beads in the simulation, Cartesian moves should be used instead.)

**Example 9.216** Path integral Monte Carlo simulation of H 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 9.217** 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 warmup, adjusts step size to 40% acceptance $end