Description of QMPFF2. Because QMPFF1 has already been described (15), here we present QMPFF2 by reference to QMPFF1. In both versions of QMPFF, a molecule is represented as a set of interacting atoms, with each atom consisting of a positive point charge (the core) and a diffuse negatively charged electron density (the cloud). This electron density is represented by an isotropic, exponential distribution, which is centered on the core for an isolated atom, but which shifts in an external field, e.g., in the presence of other atoms (see Fig. 5). It is this movement of the electron clouds that provides a natural way to simulate electronic molecular polarizability.

In both versions of QMPFF, nonbonded interactions between atoms consist of four components: electrostatic, exchange, induction, and dispersion. The functional form of the components imitates that of their QM counterparts. Specifically, electrostatics represent the classical Coulomb interaction of point charges and exponential charge densities. Exchange repulsion is a result of cloud–cloud interactions that decay almost exponentially with distance; the precise form was modified from QMPFF1 to QMPFF2, cf. ref. 15 and Eqs. **6** and **7** in *Supporting Text*, which is published as supporting information on the PNAS web site. Induction is simulated by a “spring” attaching each mobile electron cloud to a reference position (Fig. 5). The restraint potential provided by the spring is close to harmonic at small distances, but the stiffness becomes infinite as the distance approaches a limiting value, thus preventing the “polarization catastrophe” (15). For QMPFF1 the induction term is isotropic, whereas for QMPFF2 it is anisotropic (see Eq. **10** in *Supporting Text*) to provide a better description of the molecular polarizability tensor. Dispersion is represented by a term decaying as *r*^{−6} at large distances in QMPFF1, but as a sum of terms decaying as *r*^{−6} and *r*^{−8} in QMPFF2 for greater accuracy, cf. ref. 15 and Eqs. **8** and **9** in *Supporting Text*.

In QMPFF1, only nonbonded interactions were taken into account, i.e., molecules were assumed rigid. In QMPFF2, the bonded interactions are conventionally subdivided into stretching, bending, and torsion terms. The first two terms use a quadratic function, whereas the torsion term uses a threefold cosine function. Also, more atom types are used in QMPFF2 to allow better resolution of related, but electronically distinct, atom types, and to cover other biologically important elements not treated by QMPFF1, cf. ref. 15 and Tables 4 and 5, which are published as supporting information on the PNAS web site. To accommodate the more refined functional forms of the QMPFF2 nonbonded interactions, each atom type requires eight parameters, and each general bond type requires five parameters (two parameters for symmetric bonds). For the bonded interactions, each bond length, bond angle, and torsion angle function uses two, two, and three parameters, respectively, which is comparable with most advanced empirical force fields (2–5). The atom and bond types in a water molecule use a total of 18 nonbonded and 4 bond parameters.

QMPFF parameterization of nonbonded interactions is based on two fundamental principles. (*i*) The parameters of the model are fitted only to data from high-level *ab initio* QM calculations without use of any experimental data [for QMPFF1, the *ab initio* data were calculated at the MP2/aTZ(-hp) level, whereas for QMPFF2 more accurate QM techniques are used, see *Supporting Text*]. (*ii*) Each component of the interaction energy (electrostatic, exchange, induction, and dispersion) is separately fitted to the corresponding component in a decomposition of the QM-derived interaction energies of representative conformations of simple dimers. Such separate fitting prevents errors in various components from compensating for each other and is a crucial requirement for reliable transferability of the force field. In addition, we fit the QM components of dipole vectors and quadrupole and polarizability tensors of single molecules. Notably, no systems with more than two molecules are used in the parameterization. Parameters of bonded interactions, including those in water, are fitted to *ab initio* calculations at the MP2/TZ(-hp) level of the variation of bond lengths, bond angles, and torsion angles about their equilibrium values.

MD Protocol. The QMPFF2 model was incorporated into our in-house MD package, which allows both classical and quantum (path integral) simulations. We performed simulations with an isothermal–isobaric NPT ensemble consisting of 256 flexible water molecules in a cubic box under periodic boundary conditions. Equations of motion were generally integrated for 2 ns with a 1-fs time step by using the velocity Verlet scheme (29) combined with the Nose–Hoover chain thermostat (30) and Berendsen barostat (31). The electron cloud positions were optimized at every time step.

In our PIMD package the staging transform (32) of coordinates is used as well as application of independent thermal bath to every degree of freedom. To evaluate energy and pressure the virial estimators are used. The path integral discretization index *P* for quantum MD was generally chosen to be equal to 4. Certainly simulations with *P* = 4 do not provide the convergence with respect to internal molecular motions (e.g., bond stretching). However, this level is adequate to describe quantum effects in intermolecular motions, as has been verified by a few calculations that were performed with simulation time up to 5 ns and a *P* value of up to 20; for further details see *Supporting Text*.