In the preceding discussion, it has been assumed that the N particles in the system under consideration could be treated as classical point particles. In many cases, this treatment is justifiable, however, there is a large class of systems for which such an approximation is not valid. In general, systems where hydrogen/proton motion is important, for example, proton transfer processes, often have significant nuclear quantum effects. The problem of treating nuclear quantum effects in a system at finite temperature requires the solution of a quantum statistical mechanical problem. One approach that has been applied with considerable success is based on the Feynman path integral formalism of statistical mechanics [94, 95].
Consider the quantum canonical partition function for a single particle in one spatial dimension. The partition function is given by the trace:
where the trace is carried out in the coordinate basis. Assuming H=T+U, where T is the kinetic energy operator, and U is the potential, the Trotter theorem, Eq. (6.4), allows to be expressed as in the limit . The Trotter theorem expression for is then substituted into Eq. (8.1), and an identity operator in the form of is inserted in between each factor of , yielding
Then, using the fact that
one obtains the final expression for Q as a function of P
where is an effective potential given by
Equation (8.4) is in the form of a configurational partition function for a P-particle system in one dimension subject to a potential . The configurational partition function can also be expressed in a quasi phase space form by recognizing that the prefactor can be written as a product of P uncoupled Gaussian integrals:
In Eq. (8.6), the constant is an overall constant that ensures equality of Eqs. (8.6) and (8.4). In addition, the mass m', being a fictitious mass, is arbitrary, a fact that can be exploited in devising an MD scheme for Eq. (8.6) as will be shown below. As was pointed out by Chandler and Wolynes , Eqs. (8.6) and (8.7) together show that, for finite P, the path integral of a single quantum particle is isomorphic to a classical system of P particles subject with a Hamiltonian given by Eq. (8.7). Inspection of Eq. (8.5) shows that the P particles form a closed polymer chain with nearest neighbor harmonic coupling and are subject to a potential U. The classical isomorphism allows molecular dynamics to be used to simulate a finite-temperature quantum system. The extension of the path integral scheme to N particles in three dimensions is straightforward if it is assumed that the particles obey Boltzmann statistics, i.e., all spin statistics are neglected. In this case, the partition function is
where the classical Hamiltonian is given by
In principle, the equations of motion resulting from Eq. (8.9) could be implemented as a MD procedure, from which quantum equilibrium properties of a system could be computed . A number of well known difficulties arise in a straightforward implementation of MD to the path integral. Primarily, since , the force constant of the harmonic coupling increases as P increases, giving rise to a stiff harmonic interaction and a time scale separation. As was shown by Hall and Berne , this time scale separation gives rise to non-ergodic trajectories that do not sample the available canonical phase space. A solution to this problem was first presented in Ref. . There, it was shown that several elements are needed to devise an efficient MD scheme for path integrals. First, a change of variables that diagonalizes the harmonic coupling is introduced. This has the effect of isolating the various time scales present in the Hamiltonian of Eq. (8.9). The change of variables is linear, having the general form
where the matrix is a constant matrix of unit determinant. Two different choices of the matrix , discussed in Ref.  and , lead to the staging and normal mode transformations. The transformed coordinates are known as staging or normal mode variables. If the change of variables is made in Eq. (8.8), then the corresponding classical Hamiltonian takes the form:
where the s-dependent masses result from the variable transformation. For a staging transformation, the masses are , for , while for the normal mode transformation, the masses are proportional to the normal mode eigenvalues. Thus, it is clear that the fictitious masses should be chosen according to and . In this way, all modes will move on the same time scale, leading to maximally efficient exploration of the configuration space.
In addition to variable transformations, it is necessary to ensure that a canonical phase space is generated. This can be achieved via one of the non-Hamiltonian MD schemes for generating the NVT ensemble. It has been found that maximum efficiency is obtained if each Cartesian direction of each mode variable is coupled to its own thermostat, as was clearly demonstrated in Ref. , and multiple time scale integration techniques are employed .
It is worth mentioning that the path integral MD scheme outlined here has been combined with ab initio MD to yield an ab initio path integral Car-Parrinello method [101, 99]. This allows quantum effects on chemical processes to be studied. More recently, the ab initio path integral scheme has been extended to incorporate approximate quantum dynamical properties  via the so called centroid dynamics method [103, 104]. Finally, the path integral MD scheme has been modified to allow path integral simulations under conditions of constant temperature and pressure to be carried out .