9.9.4 Thermostats: Sampling the NVT Ensemble

A stochastic, white-noise Langevin thermostat (AIMD_THERMOSTAT = LANGEVIN) combines random “kicks” to the nuclear momenta with a dissipative, friction term. The balance of these two contributions mimics the exchange of energy with a surrounding heat bath. The resulting trajectory, in the long-time sampling limit, generates the correct canonical distribution. The implementation in Q-Chem follows the velocity Verlet formulation of Bussi and Parrinello,¹⁴⁴ which remains a valid propagator for all time steps and thermostat parameters. The thermostat is coupled to each degree of freedom in the simulated system. The MD integration time step (TIME_STEP) should be chosen in the same manner as in an NVE trajectory. The only user-controllable parameter for this thermostat, therefore, is the timescale over which the implied bath influences the trajectory. The AIMD_LANGEVIN_TIMESCALE keyword determines this parameter, in units of femtoseconds. For users who are more accustomed to thinking in terms of friction strength, this parameter is proportional to the inverse friction. A small value of the timescale parameter yields a “tight” thermostat, which strongly maintains the system at the chosen temperature but does not typically allow for rapid configurational flexibility. (Qualitatively, one may think of such simulations as sampling in molasses. This analogy, however, only applies to the thermodynamic sampling properties and does not suggest any electronic role of the solvent!) These small values are generally more appropriate for small systems, where the few degrees of freedom do not rapidly exchange energy and behave may behave in a non-ergodic fashion. Alternatively, large values of the time-scale parameter allow for more flexible configurational sampling, with the tradeoff of more (short-term) deviation from the desired average temperature. These larger values are more appropriate for larger systems since the inherent, microcanonical exchange of energy within the large number of degrees of freedom already tends toward canonical properties. (Think of this regime as sampling in a light, organic solvent.) Importantly, thermodynamic averages in the infinite-sampling limit are completely independent of this time-scale parameter. Instead, the time scale merely controls the efficiency with which the ensemble is explored. If maximum efficiency is desired, the user may externally compute lifetimes from the time correlation function of the desired observable and minimize the lifetime as a function of this timescale parameter. At the end of the trajectory, the average computed temperature is compared to the requested target temperature for validation purposes.

An alternative thermostat approach is also available, namely, the Nosé-Hoover thermostat⁷⁵⁶ (also known as a Nosé-Hoover “chain”), which mimics the role of a surrounding thermal bath by performing a microcanonical ($NVE$) trajectory in an extended phase space. By allowing energy to be exchanged with a chain of fictitious particles that are coupled to the target system, $NVT$ sampling is properly obtained for those degrees of freedom that represent the real system. (Only the target system properties are saved in $QCSCRATCH/AIMD for subsequent analysis and visualization, not the fictitious Nosé-Hoover degrees of freedom.) The implementation in Q-Chem follows that of Martyna,⁷⁵⁶ which augments the original extended-Lagrangian approach of Nosé^{822, 823} and Hoover,⁴⁸¹ using a chain of auxiliary degrees of freedom to restore ergodicity in stiff systems and thus afford the correct $NVT$ ensemble. Unlike the Langevin thermostat, the collection of system and auxiliary chain particles can be propagated in a time-reversible fashion with no need for stochastic perturbations.

Rather than directly setting the masses and force constants of the auxiliary chain particles, the Q-Chem implementation focuses instead, on the time scale of the thermostat, as was the case for the Langevin thermostat described above. The time-scale parameter is controlled by the keyword NOSE_HOOVER_TIMESCALE, given in units of femtoseconds. The only other user-controllable parameter for this function is the length of the Nosé-Hoover chain, which is typically chosen to be 3–6 fictitious particles. Importantly, the version in Q-Chem is currently implemented as a single chain that is coupled to the system, as a whole. Comprehensive thermostatting in which every single degree of freedom is coupled to its own thermostat, which is sometimes used for particularly stiff systems, is not implemented and for such cases the Langevin thermostat is recommended instead. For large and/or fluxional systems, the single-chain Nosé-Hoover approach is appropriate.