(CTG journal of Chemical Theory and Computation. © Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX-XXX pubs.acs.org/JCTC Anomalous Effects of Velocity Rescaling Algorithms: The Flying Ice Cube Effect Revisited Efrem Braun/® Seyed Mohamad Moosavi,*® and Berend Smit*'^'*® ^Department of Chemical and Biomolecular Engineering, University of California, Berkeley, Berkeley, California 94720, United States ^Laboratory of Molecular Simulation (LSMO), Institut des Sciences et Ingenierie Chimiques (ISIC), Valais, Ecole Polytechnique Federale de Lausanne (EPFL), Rue de l'lndustrie 17, CH-1951 Sion, Switzerland O Supporting Information 1500 - O OJ 1000 - CL "> _ ~D § B 500 - S3 o> rc CO 4-1 u 0 - 1990 2020 ABSTRACT: The flying ice cube effect is a molecular dynamics simulation artifact in which the use of velocity rescaling thermostats sometimes causes violation of the equipartition theorem, affecting both structural and dynamic properties. The reason for this artifact and the conditions under which it occurs have not been fully understood. Since the flying ice cube effect was first demonstrated, a new velocity rescaling algorithm (the CSVR thermostat) has been developed and become popular without its effects on the equipartition theorem being truly known. Meanwhile, the use of simple velocity rescaling and Berendsen (weak coupling) thermostat algorithms has not abated but has actually continued to grow. Here, we have calculated the partitioning of the kinetic energy between translational, rotational, and vibrational modes in simulations of diatomic molecules to explicitly determine whether the equipartition theorem is violated under different thermostats and while rescaling velocities to different kinetic energy distributions. We have found that the underlying cause of the flying ice cube effect is a violation of balance leading to systematic redistributions of kinetic energy under simple velocity rescaling and the Berendsen thermostat. When velocities are instead rescaled to the canonical ensemble's kinetic energy distribution, as is done with the CSVR thermostat, the equipartition theorem is not violated, and we show that the CSVR thermostat satisfies detailed balance. The critical necessity for molecular dynamics practitioners to abandon the use of popular yet incorrect velocity rescaling algorithms is underscored with an example demonstrating that the main result of a highly cited study is entirely due to artifacts resulting from the study's use of the Berendsen thermostat. 2000 2010 Year 1. INTRODUCTION By integrating the classical Newtonian equations of motion, molecular dynamics (MD) simulations naturally sample the microcanonical (NVE) ensemble due to conservation laws.1'2 For comparison with experiment, it is often desirable to sample constant-temperature ensembles such as the canonical (NVT) or isothermal—isobaric (NPT) ensembles. In analogy with experiment, these ensembles could be generated by sampling a subspace of a much larger microcanonical system that serves as a heat bath, but such an approach is usually too computationally expensive to implement in practice. Instead, various fhermo-statting algorithms are typically applied to change the Hamiltonian dynamics in a manner such that the intended ensemble is sampled. Many such algorithms have been proposed, and some of the more well-known choices include the following: • Simple velocity rescaling, pioneered by Woodcock3 for thermal equilibration, rescales the velocities of all particles at the end of each time step (it can also be conducted with a less frequent time rescaling period) by a factor X to achieve a target instantaneous temperature: X Ik Ý'1 I ^target \ with ^target = 2ndoffcbTtargey where ndof is the number of degrees of freedom in the system. • The Gaussian thermostat supplements Newton's second law with a force intended to keep the kinetic energy constant:4-p, = —VU, — op„ where a is a Lagrange multiplier determined using Gauss' principle of least constraint to be a = (J^i^rPi/ m,)/(l£iP?/m,)- • Langevin dynamics supplements Newton's second law with terms describing Brownian motion: p, = — VU, — /p, + tj, where y represents a factional dissipative force and t]{t,T,Y,m^) is a stochastic term representing random collisions. • The Berendsen thermostat (going forward, we will refer to this as the weak coupling (WC) thermostat) takes the Langevin equation, removes the stochastic term, and modifies the frictional dissipative force to yield similar temperature time Received: May 8, 2018 Published: August 3, 2018 ACS Publications ©XXXX American Chemical Society A DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX Journal of Chemical Theory and Computation Article with the stochastic term present: where K>, j^DOF^B^targef dependence p. = -VLf - yp. In practice, this is implemented as a smoother version of the simple velocity rescaling technique in which the velocities of all particles are rescaled at the end of each time step by a factor X with X 1 + At / -Kmga A TT\ K ) 1/2 where Tr represents a time damping constant; if it is set equal to the time step, the algorithm recovers simple velocity rescaling, and as the time damping constant approaches infinity, the algorithm recovers conventional microcanonical dynamics. • The canonical sampling through velocity rescaling (CSVR) thermostat is a velocity rescaling algorithm in which the velocities of all particles are rescaled at the end of each time step by a factor X designed such that the kinetic energy exhibits the distribution of the canonical ensemble.9'10 To this end, X —— I ™g" I , where is stochastically drawn from the probability density function P(Xarget) cx i^^/^'e"^-. This algorithm can be adjusted to yield a smoother evolution in a similar manner as the weak coupling thermostat smoothes simple velocity rescaling.9 • The Nose—Hoover thermostat extends the classical Lagrangian to include the additional coordinate s and its time- derivative: u + ln.2 where Q is the effective "mass" associated with s and L is set by the number of degrees of freedom. A single Nose—Hoover thermostat may be used, or chains of thermostats may be implemented to improve ergodicity and to take into account additional conservation laws.13 There exist numerous additional thermostats (e.g., the Andersen thermostat14), and small changes can be made to the listed thermostats, such as implementing the originally global Nose—Hoover thermostat in a local "massive" manner by pairing a separate Nose—Hoover thermostat to each degree of freedom.15 The reader is referred to a noncomprehensive list of reviews and textbooks for additional information.1'16-18 Simple velocity rescaling and the Gaussian thermostat aim to sample the isokinetic ensemble (NVK). However, they are often presented as equivalent to the canonical ensemble with respect to position-dependent equilibrium properties with justification for this based on the argument that the configurational part of the isokinetic ensemble's partition function is exactly equal to that of the canonical ensemble.5'19-22 Meanwhile, the weak coupling thermostat does not correspond to a known ensemble but is rather supposed to sample a configurational phase space intermediate to the canonical and microcanonical ensem-bles.8'23'24 In the 1990s, it was found that the simple velocity rescaling and weak coupling thermostat algorithms introduce an artifact:25'26 the "flying ice cube effect", as coined by Harvey et al.,26 describes a violation of the equipartition theorem observed when using these algorithms in which kinetic energy drains from high-frequency modes such as bond stretching into low-frequency modes such as center of mass (COM) translation. This was shown to affect systems' structural, thermodynamic, and dynamic properties. As it can be proven that the equipartition theorem holds in the canonical, microcanonical, and isokinetic ensembles (see Supporting Information),27 31 a simulation exhibiting the flying ice cube effect is not ergodically sampling any of these ensembles in configurational or momentum phase space. Nonetheless, simple velocity rescaling and the weak coupling thermostat continue to be commonly used17'32 with Cooke and Schmidler32 stating, "By far the most commonly used algorithm for constant temperature MD of biomolecules is the Berendsen heat bath, due to its ease of implementation and availability in standard software packages." Use of the weak coupling thermostat can be approximated by tracking citations of its canonical reference,8 which have continued to grow over time (Figure l). 1200 1000 BOO 600 400 200 - Berendsen et al. Bussi etal. 1935 1990 1995 (2000 Harvey et al. _Year 2005 2010 2015 Figure 1. Citations of Berendsen et al.8 and Bussi et al.9 over time. Data provided by Web of Science, extracted on May 4, 2018. Some technical aspects of the flying ice cube effect are as of yet still unclear. Since Harvey et al.,26 there has been continued discussion about whether the flying ice cube effect may occur with other thermostats.33'34 The CSVR thermostat rescales velocities to yield the canonical ensemble's distribution of kinetic energies similar to how simple velocity scaling yields the isokinetic ensemble's distribution of kinetic energies and the weak coupling thermostat yields a kinetic energy distribution intermediate to the two ensembles. If all velocity rescaling algorithms always lead to the flying ice cube effect, then it maybe suspected that the same flying ice cube artifact occurs when using the CSVR thermostat, which would be worrisome because the CSVR thermostat has been quickly adopted into widespread use (Figure l). In addition, because the Gaussian thermostat has been shown to be similar to simple velocity rescaling,36 it may be suspected that the Gaussian thermostat exhibits the artifact as well. Given the widespread use of these algorithms in MD simulations, more understanding is warranted, and we will show that neither the CSVR thermostat nor the Gaussian thermostat bring about the flying ice cube effect. In the present work, we refer to the "flying ice cube effect" as the term was originally used to describe the violation of the equipartition theorem as caused by velocity rescaling procedures.26 Multiple manifestations of this effect are possible, such as the accumulation of kinetic energy into translational or vibrational degrees of freedom (as we will demonstrate) or the development of temperature gradients,37'38 all of which we refer to as the flying ice cube effect. Meanwhile, other MD simulation DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX 99 Journal of Chemical Theory and Computation Article methods that fail to conserve energy in the microcanonical ensemble can also bring about equipartition theorem violations, including approximate treatment of long-range electrostatic interactions, certain multiple time step algorithms, constraining molecular geometries with too loose of a tolerance, not updating neighbor lists frequently enough, and using too large of a time step.33'39'40 In some cases, these issues are also referred to as flying ice cube effects, but these are not related to the artifact with which we are concerned.41 43 In this work, we have revisited the simple model system of united-atom diatomic ethane molecules that Harvey et al.26 first used to illustrate the flying ice cube effect. By explicitly calculating the partitioning of kinetic energies between trans-lational, rotational, and vibrational degrees of freedom, we are able to determine which thermostats and conditions lead to the violation of equipartition as well as the manner and degree to which they do so. We go on to rationalize these findings by illustrating how simple velocity rescaling violates balance while the CSVR thermostat satisfies detailed balance. We end by illustrating some severe errors that are directly caused by these subtleties related to thermostatting. 2. SIMULATION DETAILS Diatomic ethane molecule simulations were conducted with the open-source LAMMPS code.44 LAMMPS input scripts are available in the Supporting Information 45 Except where stated otherwise, the simulations consisted of cubic simulation boxes with periodic boundary conditions (PBC) set up by placing the ethane molecules on a simple cubic lattice, equilibrated with a Langevin thermostat for at least 50 ns, switched to the target thermostat for at least a further 50 ns of equilibration, and finally ran with the target thermostat for at least 50 ns of production. We verified that all simulations were conducted for sufficient time periods for the energies to equilibrate and be well sampled. For the simulations in which the COM linear momentum was fixed to zero (stated in the figure captions), the system's linear momentum was zeroed every time step followed by a rescaling of velocities to maintain the same total kinetic energy as before the zeroing had occurred to prevent energy leakage. The equations of motion were integrated with a standard Velocity Verlet algorithm using half-step velocity calculations. The time step used was 0.5 fs, which was found to give adequate energy conservation in the microcanonical ensemble. Thermostat parameters were as follows except where stated otherwise. Simple velocity rescaling was done every time step. The Nose—Hoover chain consisted of three thermostats. The weak coupling, Nose—Hoover, and CSVR thermostats were used with time damping constants (zT) of 100 fs, and the Nose-Hoover thermostat used effective thermostat masses of = NDOfkBTz^ and QJ>1 = fcBTTy.13 When doing simulations in the microcanonical ensemble, the total energy was set such that a simulation temperature equal to the canonical ensemble simulations' target temperature was achieved. The target simulation temperature was set to 350 K, well above the critical temperature of ethane 46 Kinetic energies of each diatomic molecule were partitioned into translational, rotational, and vibrational kinetic energies, as shown in the Supporting Information. In all figures that plot kinetic energies, the error bars shown represent ±1 standard error of the mean. This was calculated by dividing the data from the production timesteps into 20 consecutive blocks, averaging the data for each block and computing the standard error over the 20 data values.1 Error bars are not shown when they would be smaller than the symbols or the line widths. Bonded parameters for the united-atom ethane molecule were taken from Harvey et al.26 (harmonic bond potential U(r) = k(r - r0)2 with r0 = 1.54 A and k = 240 kcal mol"1 A"2), and nonbonded parameters were taken from Martin and Siepmann46 (Lennard-Jones potential with e = 0.195 kcal mol-1, a = 3.75 A, truncated and shifted at 14 A, and no charges). Details on the simulations of benzene in MOF-5 can be found in the Supporting Information. 3. RESULTS AND DISCUSSION 3.1. Examining Equipartition under Different Thermostats. It is instructive to reconsider the simple case previously examined by Harvey et al.,26 that of a single ethane molecule moving in one-dimensional space along its bond axis. In the microcanonical ensemble under perfect energy conservation, the translational kinetic energy will remain constant at its set initial energy, and the vibrational kinetic energy will oscillate. In the canonical ensemble, equipartition states that the translational and vibrational degree of freedom should each have an average kinetic energy of ^kBT. As expected, the Langevin thermostat satisfies the equipartition theorem (see Figure 2). In agreement 0j u 0.010 - ■c u 0.005 - I en IT3 0J < 0.000 j-JJJJJJJMJ-1.1 I i-lu-LLS—i Langevin CSVR Simple WC VR Figure 2. Partitioning of the kinetic energies obtained from one-dimensional MD simulations of a single ethane molecule using various thermostats. Both atoms were given a starting velocity of 100 m s_1 along the same direction as the bond vector. For the thermostats shown, the same energy partitionings were observed regardless of initial bond length and initial COM momentum. The microcanonical, Nose-Hoover thermostat, and Gaussian thermostat results are not shown here because we found that the energy partitionings are dependent on the initial conditions, indicative of these thermostats' well-known lacks of ergodicity that are more manifest for small systems.1,5'13'17'36'47~49 with the work of Harvey et al., we find that simple velocity rescaling and the weak coupling thermostat bring about a violation of equipartition in the kinetic degrees of freedom with all kinetic energy flowing to translational motion, in the plainest illustration of the flying ice cube effect. We find that the CSVR thermostat correctly partitions the energies. We next consider the more complex case of a large number of ethane molecules interacting in three dimensions with anharmonic Lennard-Jones potentials. Each diatomic ethane molecule now has three translational modes, two rotational C DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX Journal of Chemical Theory and Computation Article o 10 > ,. 1—I—I—l_l—I—I—I I I—1—1—I—I—l_l—LJ—I—!_l—J_l—l_l—LJ—I—I I i I—I-- _NVE_Langevin NM-chain_CSVR_Gaussian Simple VR_WC Figure 3. Partitioning of the kinetic energies obtained from MD simulations of 50 ethane molecules in a 30 A cubic simulation box using various thermostats. In all simulations shown, the COM momentum was fixed to zero. Simple Velocity Repealing Weak Coupling Thermostat Figure 4. Partitioning of the kinetic energies obtained from MD simulations performed under the same conditions as in Figure 3 but changing the time step, using (left) simple velocity rescaling and (right) the weak coupling thermostat with the time damping constant maintained at 100 fs. Lines are a guide to the eye. modes, and one vibrational mode, so the equipartition theorem states that these modes' kinetic energies should be equal to ■jfcgT, ^kBT, and ~^kBT, respectively, with a correction of ■|fcBT/Nmolecs to the translational kinetic energy in cases where the COM momentum is constrained. In Figure 3, we show that the Langevin, Nose—Hoover, CSVR, and Gaussian thermostats all exhibit correctly equipartitioned energies as does the microcanonical ensemble. As in the case of the single ethane molecule in one dimension, the simple velocity rescaling and weak coupling thermostat algorithms lead to a violation of equipartition with translational and rotational modes having too much kinetic energy and vibrational modes having too little. 3.2. Equivalence of Simple Velocity Rescaling and the Gaussian Thermostat. Because the thermostatting under simple velocity rescaling does not take place within the equations of motion, this ad hoc temperature control algorithm was initially difficult to investigate theoretically, and its validity was considered questionable.5'6 The algorithm's use was justified on the basis of empirical arguments, such as that simple velocity rescaling and the Gaussian thermostat give similar static and dynamic properties for the Lennard-Jones fluid.19 It was eventually proven that simple velocity rescaling is analytically equivalent to the Gaussian thermostat within an error of (9(timestep) when the velocity rescaling time period is set equal to the time step, which gave support for the legitimacy of using simple velocity rescaling to sample the isokinetic ensemble. However, we have shown that the Gaussian thermostat exhibits correct energy equipartitioning, whereas simple velocity rescaling does not. We prove in the Supporting Information that the isokinetic ensemble should satisfy the equipartition theorem. Thus, it is clear that simple velocity rescaling does not actually sample the isokinetic ensemble. The equivalence of simple velocity rescaling and the Gaussian thermostat under small timesteps leads to the expectation that the flying ice cube effect will gradually disappear under simple velocity rescaling as the time step is decreased. We demonstrate confirmation of this expectation in Figure 4. However, Figure 4 shows that the time step needs to be reduced by over 3 orders of magnitude from typical simulation timesteps before the flying ice cube effect is no longer discerned. Of course, such a decreased time step requires an equivalent 3 orders of magnitude increase in CPU time; if the time step between integrations is so small, the forces on the particles should not need to be recalculated every time step, and thus one could envision implementing a multiple-time-step algorithm to mitigate the increase in CPU time. We also note that, under the weak coupling thermostat, lowering the time step does not correct the energy partitioning. 3.3. Violation of Balance Causes the Flying Ice Cube Effect. The mechanism underlying the flying ice cube effect can DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX Journal of Chemical Theory and Computation Article be elucidated graphically for the first test case we examined, that of a single ethane molecule. In Figure 5, we show this system's phase space, putting translational kinetic energy on the x-zxis and vibrational kinetic energy on the y-axis. Figure 5. Kinetic phase space of a single ethane molecule moving in one-dimensional space along its bond axis under simple velocity rescaling. = fcBTtarget, Ktrans = Ujnx + m2) ("^^T^T^) ' and Kvjb = — I "v"1 )(vlx — vlx)2. Solid lines show a particular path in 2 \ mi + m2 / ' ' phase space between labeled points, referred to in the text. Dotted lines are guides useful to understanding the velocity rescaling moves. Dashed lines show the boundaries of phase space accessible by any sequence of MD and velocity rescalings from lines AB, CD, and AG with the accessible phase spaces shaded. During microcanonical MD, the system can only explore phase space on a vertical line betweeny = 0 andy = Umax because a constant total energy and translational kinetic energy is maintained with energy exchanges only allowed between vibrational kinetic energy and potential energy. Consider a MD simulation initially on such a vertical line in phase space, AB. Under simple velocity rescaling, if a rescaling move is conducted at point B, the system will move to point C; this occurs because the translational and vibrational energies are both scaled by the same factor 1 such that their sum is equal to the target kinetic energy, moving the system to the intersection of the lines y = —x and the target isokinetic line (y = — x + t). Because points B and C have the same configuration with zero potential energy, MD will now explore line CD. Let us examine whether we can reach point B by rescaling from line CD back to a line with the same translational energy of line AB. With a single rescaling, we would need to rescale from point E to point F. From point F, MD will explore phase space on line AG, where the lengths of lines FG and CE are equal with both representing the stored potential energy of the system prior to the rescaling. Obviously, line EF must have a smaller slope than line BC; accordingly, yG will necessarily be smaller than yB. Hence, with a single velocity rescaling, point B cannot be reached. Multiple velocity rescalings from line CD allows us to reach a point with greater vibrational kinetic energy than point G. However, all phase space reachable by any number of velocity rescalings from line CD is bounded by the red dashed line in Figure 5 (see Supporting Information for derivation). Continuing to rescale will continue to shrink the volume of accessible phase space, as rescaling from lines AB to CD to AG lowers the boundary from the blue to the red to the green dashed lines; eventually, accessible phase space will be confined only to the point with all kinetic energy in the translational mode. Notably, the decrease in accessible phase space becomes smaller as velocity rescaling occurs closer to the isokinetic line. In a simulation, this occurs when the time step between velocity rescalings is reduced. This explains why the flying ice cube effect is reduced under simple velocity rescaling by decreasing the time step (Figure 4). 3.3.1. Monte Carlo Perspective. We can view the combination of MD and velocity scaling moves as a Monte Carlo simulation. Hence, our previous example shows that simple velocity rescaling violates the condition of balance.1'50 In contrast, the CSVR thermostat can explicitly be proven to sample the desired distribution by considering the condition of detailed balance. Let us assume that we do a large and random number of MD steps between velocity rescaling moves. We define A as the set of all configurations of the system with a total energy EA. The flow of configurations from set A to set B is given by K(A —> B) = P(EA) 'Z'Z'Z P"IEA) '1 ?\ '2 ?l S(E(i[, p") - -EA)a(r", p" -> r\, p2)S(E(rl p2) - EB) (l) where rj,pj is the configuration with position vector rj, momentum vector pj, p(rJ,pJlEA) is the probability of finding the configuration rj,pj from all configurations with energy EA during MD, and o(rJ,pJ —> r^p-J) is the a priori probability to velocity rescale from configuration rj,pj to configuration r^pj Recognizing that velocity rescaling does not alter positions K(A -» B) = P(EA) X E EP(r"> Pi|£a) r" pi p\ S(E(rn, p") - EA)a(r", p" — r", p"2)S(E(rn, p"2) - EB) (2) Next, recognizing that velocity rescaling can only give one configuration in momentum space with E(r",p2) = EB from starting configuration r",pj and that the acceptance probabilities only involve the kinetic energy K(A -» B) = P(EA) 2 X p"|£a) r" p" S(E(rn, p") - EA)a(K = EA - U(r") -» EB - U(r")) (3) where a(K = EA — U(r") —> EB — U(r")) is the a priori probability to velocity rescale to the configuration having kinetic energy K = EB — U(r") given that we start with a configuration having kinetic energy K = EA — U(r"). Then, recognizing that momentum and position are decoupled, i.e., the number of possible states in momentum space only depends on the total kinetic energy but does not depend on the details of the potential energy surface and that each of these possible states in momentum space are equally likely E DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX Journal of Chemical Theory and Computation Article Kinetic energy (eV per ethane molecule) COM Fixed COM Free £ 0.06 - 10" 101 10r 103 10* 10s 10' _Mdof_ 10° 101 10* 10! 101 10s 10' 10T Figure 6. (top) Probability density function of kinetic energies following P(K) = ,-/*kNdof«- ., -, where 6 is chosen such that the average kinetic /rNDOF/2r(-NDOF/2)' ^ e : 300, and fi0 = (ks X 350 K)_1. (bottom) Partitioning of the kinetic energy (temperature) is the same for all choices of NDOF via fi = /? —^ ND( ndof,0 energies obtained from MD simulations of 50 ethane molecules in a 30 A cubic simulation box using the CSVR thermostat modified such that the target distribution of kinetic energies was set to those shown in the top part of the figure for the proper NDOF0 value, (bottom left) Here, the COM momentum was fixed at 0, and NDOF0 was set to 297. (bottom right) Here, the COM momentum was not fixed after the Langevin thermostat equilibration, allowing the COM momentum to drift, and NDOF0 was set to 300. Lines are a guide to the eye. K(A -» B) = P(EA) 2 a>(EA - U(r"))p(r", p"LEA) r a(K = EA- U(r") -» EB - U(r")) (4) where cl)(K) is the number of configurations in momentum space for a given kinetic energy K (equivalent to the ideal gas microcanonical partition function). Finally, by making the a(K : I7(r") I7(r")) a(K : I7(r") I7(r")) e-^a>(EB ~ U(r")) e ^MEA ~ U(r")) U(r"))A .F/2-l -i!(EA-U(rn)) (EA - U(r")f .F/2-l (6) substitutions p(rn,pn\EA)=Q,N$/EA and P(EA) in which we used the known expression for the ideal gas microcanonical partition function. Eq 6 is satisfied by the CSVR thermostat, which rescales velocities to the target kinetic energy distribution given by the gamma distribution K(A — B) X co{EA - l/(r")) P{K) e-/XKNDOF/2-l e-/XKNDOF/2-l f°° dKK1 Jo ho*n-i-i)K ■f/2 r(NDOF/2) (7) a{K = EA- U(r") -» EB - U(r")) (5) The two flows, K(A —> B) and K(B —> A), are equal if we impose as condition for the a priori probabilities Hence, the CSVR thermostat satisfies detailed balance. 3.3.2. Velocity Rescaling to Other Kinetic Energy Distributions. We have seen that simple velocity rescaling violates balance and brings about the flying ice cube effect, whereas the CSVR thermostat satisfies detailed balance and does not exhibit the artifact. One key difference between these algorithms is that DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX Journal of Chemical Theory and Computation Article Supercritical Saturated Liquid 150 100 50 1 2 ■ MVE NH-c haln Gaussian Simple W WC RDF=1 — NVE --NH chain ---Gaussian ...... Simple VP. - vic .... RDF=1 r (A) Saturated Liquid a o.o2 - en NDOF0) and broader (^dof < ^dof o) than the canonical distribution. In the limit of NDOF —> oo, this method closely approximates simple velocity rescaling or the weak coupling thermostat depending on the time damping constant used. The energy partitionings that resulted from setting these target kinetic energy distributions are shown for simulations in the bottom of Figure 6. It can be seen that, with sharper distributions, the flying ice cube effect is observed with more kinetic energy partitioned in low-frequency modes and less in high-frequency modes. Interestingly, the opposite effect is observed with broader distributions with more kinetic energy partitioned in high-frequency modes and less in low-frequency modes. When the COM momentum is not constrained to zero, a more drastic effect is observed such that rotational kinetic energy decreases both with decreasing NDOF as energy flows to the higher-frequency vibrational modes and with increasing NDOF as almost all energy flows to the lower-frequency translational modes. Only at the canonical kinetic energy distribution (NDOF = 297 and NDOF = 300 for the constrained and not-constrained COM momentum simulations, respectively) is proper equipartitioning observed. 3.4. Sampling Configurational Degrees of Freedom. So far, we have exclusively used kinetic degrees of freedom to show that the simple velocity rescaling and weak coupling thermostat algorithms cause the violation of equipartition. These methods are sometimes used only to sample configurational degrees of freedom, justified on the grounds that the isokinetic ensemble samples the same configurational phase space as the canonical ensemble.5'19-22 Because we have proven that the violation of equipartition is incommensurate with sampling the isokinetic ensemble, it follows that this justification is invalid. We now wish to show this explicitly. To do so, we will DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX Journal of Chemical Theory and Computation Article examine the radial distribution function (RDF), which is solely dependent on configurational degrees of freedom. In Figure 7 (top left), we show the RDFs of the supercritical ethane simulations whose kinetic energy partitionings are shown in Figure 3. The Nose—Hoover, CSVR, Langevin, and Gaussian thermostat simulations exhibit identical RDFs, but the simple velocity rescaling and weak coupling thermostat simulations show a subtly different RDF. Although the difference is slight, it is sufficient to demonstrably disprove the claims that simple velocity rescaling samples the same configurational phase space as the canonical ensemble and that the weak coupling thermostat samples a configurational phase space intermediate between the canonical and microcanonical ensembles.23'24 We next turn to saturated liquid phase ethane simulations for which we show RDFs under various thermostats in Figure 7 (top right). The Nose—Hoover, Langevin, CSVR, and Gaussian thermostats all give identical results typical of a simple diatomic liquid.51 The simple velocity rescaling algorithm once again shows a subtle difference, but the weak coupling thermostat shows a very different RDF more reminiscent of the solid phase than the liquid phase,51 and visualization of the weak coupling thermostat system shows that the ethane molecules have indeed packed into a volume smaller than available in the simulation box Examination of the kinetic energy partitionings in Figure 7 (bottom) shows that most of the kinetic energy is in vibrational modes, which is unexpected because that is the opposite of the usual flying ice cube result. The weak coupling thermostat's results are heavily dependent on the choice of time damping constant with the RDF indicating a solid-like phase for time damping constants from approximately 10 to 150 fs (Figure S4). This effect of intermediate time damping constants giving larger deviations than small or large ones has been seen before in simulations of bulk water, where the effect was attributed to the intermediate time constant matching a characteristic time scale on which dynamical correlations are most pronounced.52 It thus appears that the weak coupling thermostat is susceptible to resonance artifacts, which we have also observed under simple velocity rescaling (Figure Si). 3.5. Contemporary Use of the Simple Velocity Rescaling and Weak Coupling Thermostat Algorithms. Ours is not the first publication to warn against the use of simple velocity rescaling and the weak coupling thermostat.2 '3 '53 Nonetheless, as we have stated, these algorithms continue to be widely used (Figure l). As we have just shown, for some systems the improper velocity rescaling algorithms may not greatly affect the system properties, and there are a slew of studies in which these thermostats are tested for specific systems with some showing artifacts and others showing indistinguishabil-ity. ' ' ' However, slight changes to a system could introduce artifacts in an unpredictable fashion. Rather than testing for the correctness of simple velocity rescaling or the weak coupling thermostat in every specific system, we advocate for the cessation of their use. We find no reason to use simple velocity rescaling or the weak coupling thermostat instead of the CSVR thermostat given their similar ease of implementation, likely similar speeds of equilibration,58 and our study's finding that the CSVR thermostat does not lead to the flying ice cube effect. As a case study on the dangers of continuing to use these thermostat algorithms, we examine a highly cited study in depth, the replication of which initially led us to examine the flying ice cube phenomenon. In 2007, a flexible force field intended for use with MOF-5 was parametrized,59 and it was shortly thereafter used to study the confined transport of guest molecules within the framework. The authors were able to replicate the experimental diffusion coefficient of confined benzene, but they found that this replicability only held when the metal—organic framework (MOF) was allowed to be flexible; when the MOF atoms were held rigid, the benzene diffusion coefficient increased by an order of magnitude. The conclusions of this manuscript are often evoked to question the validity of the rigid framework assumption that is commonly used in many MOF molecular simulation studies. The findings continue to be accepted because it is known that the effect of framework flexibility on guest diffusion is complex,61 though surprise has been expressed62 because a rigid lattice more typically leads to a decrease in the diffusion coefficient for tight fitting molecules.61 In addition, using a different flexible force field for MOF-5,63 it was found that flexibility had little effect on the diffusion coefficient, increasing it by less than a factor of 1.5.64 As the reader now anticipates, Amirjalayer et al.60 used the weak coupling thermostat, which was the default option in the Tinker simulation package at the time (the default has since been changed to the CSVR thermostat).65 As we show in Figure 8, the result of Amirjalayer et al.60 was completely an artifact of 10" 10" ♦ □ Amirjalayer et al., rigid Amirjalayer et al., flexible Weak Coupling, rigid Weak Coupling, flexible Nose-Hoover, rigid Nose-Hoover, flexible 7 4 2.6 2.8 3.0 3.2 3.4 iooo/r (K-1) 3.6 3,8 4.0 Figure 8. Self-diffusion coefficient of benzene in MOF-5 at a loading of 10 molecules per unit cell as a function of inverse temperature. Data are shown for flexible and rigid frameworks and using the weak coupling and Nose—Hoover chain thermostats (use of the CSVR thermostat gives diffusion coefficients that are statistically indistinguishable from use of the Nose—Hoover thermostat). With the weak coupling thermostat, it appears that the framework flexibility has a large effect on the calculated diffusion coefficient, replicating the main finding of Amirjalayer et al.60 However, it is seen that this result is a flying ice cube artifact as no flexibility effect is seen with the Nose—Hoover thermostat. Error bars represent ±1 standard error of the mean using block averaging1 and are not shown for the data from Amirjalayer et al.60 or if they would be smaller than the symbol size. the weak coupling thermostat. Using the same force field, no dependence of the benzene diffusion coefficient on the framework flexibility is observed when a Nose—Hoover or CSVR thermostat is used. Apparently, when the weak coupling thermostat is thermostated to fewer degrees of freedom during rigid framework simulations, the flying ice cube effect becomes more noticeable and kinetic energy is drawn into the translational modes of the guest benzene molecules, accounting DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX Journal of Chemical Theory and Computation Article for the result observed by Amirjalayer et al. We also found that changing the time damping constant of the weak coupling thermostat had a large effect on the diffusion coefficient (Figure S5). As an aside, it is now known that bulk-like vapor and liquid phases of benzene exist in MOF-5 below a critical temperature.66 It is actually improper to calculate the diffusion coefficient at a loading that is within the vapor—liquid phase envelope, e.g., 3— 67 molecules per unit cell at 300 K in this system,66 because there is not a single homogeneous phase present under these conditions. Here, we are not attempting to calculate correct diffusion coefficients of benzene in MOF-5 but rather to compare results with the prior work of Amirjalayer et al.,60 who conducted the simulations at a loading of 10 molecules per unit cell. The importance of framework flexibility on the simulated diffusion coefficient is expected to be independent of the choice of loading. Other errors, varying in severity, are likely present in many of the thousands of studies that have used simple velocity rescaling or the weak coupling thermostat. Occasionally, one of these errors is explicitly pointed out,67'68 but negative replications are not commonly published,69 so the extent to which these articles contain data contaminated by the flying ice cube artifact cannot be estimated. 4. CONCLUDING REMARKS In this work, we have shown that rescaling velocities to a noncanonical distribution of kinetic energies, as is done with the simple velocity rescaling and weak coupling (Berendsen) thermostat algorithms, causes the flying ice cube effect whereby the equipartition theorem is violated. Thus, simple velocity rescaling does not sample the isokinetic ensemble, and the weak coupling thermostat does not sample a configurational phase space intermediate between the canonical and microcanonical ensembles; justifications for their use do not hold. The flying ice cube effect is brought about by a violation of balance causing systematic redistributions of kinetic energy; this violation is lessened as the time step between simple velocity rescalings is decreased, eventually making simple velocity rescaling equivalent to the Gaussian thermostat. Equipartition violation is completely avoided when velocities are rescaled to the canonical distribution of kinetic energies, as is done under the CSVR thermostat, because detailed balance is obeyed. We have identified several simulation parameters that affect the conspicuousness of the flying ice cube effect under simple velocity rescaling and the weak coupling thermostat, including the time step, the thermostat's coupling strength, the frequency of collisions within the simulation (e.g., with a wall), and the system size (see Supporting Information for further investigations of these parameters). However, most of these parameters cannot be adjusted in a manner that eliminates the flying ice cube effect without making simulations prohibitively expensive for relevant systems of contemporary interest. Another reason not to attempt to tune these simulation parameters to allow the use of incorrect thermostatting algorithms is the existence of additional resonance artifacts that occur when the thermostat coupling strengths are set to particular values that are difficult to predict a priori. Finally, we have demonstrated several severe simulation artifacts that the flying ice cube effect can bring about to the system's structural and dynamic properties. These include incorrect RDFs, phase properties, and diffusion coefficients. We have highlighted one case in which the flying ice cube effect has been wholly responsible for the main finding of a highly cited study. Many more such cases are likely present in the literature. We strongly advocate for discontinuing use of the simple velocity rescaling and weak coupling thermostat algorithms in all MD simulations for both equilibration and production cycles. The results of past studies that have used these two algorithms should be treated with caution unless they are shown to be replicable with a more reliable thermostat. In situations where velocity rescaling methods are desirable, such as for fast equilibration of a system,70 the CSVR thermostat should be used instead. ■ ASSOCIATED CONTENT O Supporting Information The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00446. Proof that equipartition should be obeyed in the isokinetic ensemble, details on how we calculated the partitioned kinetic energies, simulations details for the benzene in the MOF-5 system, derivation of the formula for the dashed lines in Figure 5, investigation of simulation parameters that affect the conspicuousness of the flying ice cube effect, and supplementary figures (PDF) LAMMPS and Tinker input files (ZIP) ■ AUTHOR INFORMATION Corresponding Author *E-mail: berend.smit(2>epfl.ch. ORCID® Efrem Braun: 0000-0001-5379-7031 Seyed Mohamad Moosavi: 0000-0002-0357-5729 Berend Smit: 0000-0003-4653-8562 Notes The authors declare no competing financial interest. ■ ACKNOWLEDGMENTS We thank the responders on the LAMMPS mailing list for useful discussion and for giving advice regarding the LAMMPS source code (Axel Kohlmeyer, Steven J. Plimpton, and Aidan P. Thompson were particularly helpful) and Sai Sanigepalli for helping to implement the Tinker simulations. Special thanks go to Rochus Schmid for insightful discussion on the roles of thermostatting and for providing assistance in implementing the Tinker simulations. This research was supported as part of the Center for Gas Separations Relevant to Clean Energy Technologies, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award DE-SC0001015. S.M.M. was supported by the Deutsche Forschungsgemeinschaft (DFG, priority program SPP 1570). This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. ■ REFERENCES (1) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Elsevier Science, 2002. (2) Leimkuhler, B.; Matthews, C. In Molecular Dynamics with Deterministic and Stochastic Numerical Methods; Antman, S. S., Holmes, P., Greengard, L., Eds.; Interdisciplinary Applied Mathematics 39; Springer, 2015. I DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX Journal of Chemical Theory and Computation Article (3) Woodcock, L. V. Isothermal Molecular Dynamics Calculations for Liquid Salts. Chem. Phys. Lett. 1971, 10, 257-261. (4) Evans, D.J.; Hoover, W. G.; Failor, B. H; Moran, B.; Ladd,A.J. C. Nonequilibrium Molecular Dynamics via Gauss's Principle of Least Constraint. Phys. Rev. A: At, Mol, Opt. Phys. 1983, 28, 1016-1021. (5) Nose, S. A. Unified Formulation of the Constant Temperature Molecular Dynamics Methods. /. Chem. Phys. 1984, 81, 511-519. (6) Evans, D. J.; Morriss, G. P. Statistical Mechanics of Nonequilibrium Liquids; Theoretical Chemistry Monograph Series; Academic Press: London, 1990. (7) Schneider, T.; Stoll, E. Molecular-Dynamics Study of a Three-Dimensional One-Component Model for Distortive Phase Transitions. Phys. Rev. B: Condens. Matter Mater. Phys. 1978, 17, 1302-1322. (8) Berendsen, H. J. C; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. /. Chem. Phys. 1984, 81, 3684-3690. (9) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. /. Chem. Phys. 2007, 126, 014101. (10) Heyes, D. M. Molecular Dynamics at Constant Pressure and Temperature. Chem. Phys. 1983, 82, 285-301. (11) Nose, S. A. Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255-268. (12) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At, Mol, Opt. Phys. 1985, 31, 1695-1697. (13) Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nose-Hoover Chains: The Canonical Ensemble via Continuous Dynamics. /. Chem. Phys. 1992, 97, 2635-2643. (14) Andersen, H. C. Molecular Dynamics Simulations at Constant Pressure and/or Temperature. /. Chem. Phys. 1980, 72, 2384—2393. (15) Tobias, D. J.; Martyna, G. J.; Klein, M. L. Molecular Dynamics Simulations of a Protein in the Canonical Ensemble. /. Phys. Chem. 1993, 97, 12959-12966. (16) Morriss, G. P.; Dettmann, C. P. Thermostats: Analysis and Application. Chaos 1998, 8, 321-336. (17) Hiinenberger, P. H. In Advanced Computer Simulation; Holm, C, Kremer, K., Eds.; Advances in Polymer Science; Springer, 2005; Vol. 173, pp 104-149. (18) Tuckerman, M. Statistical Mechanics: Theory and Molecular Simulation; Oxford Graduate Texts; Oxford University Press: Oxford, 2010. (19) Haile, J. M.; Gupta, S. Extensions of the Molecular Dynamics Simulation Method. II. Isothermal Systems. /. Chem. Phys. 1983, 79, 3067-3076. (20) Evans, D. J.; Morriss, G. The Isothermal/Isobaric Molecular Dynamics Ensemble. Phys. Lett. A 1983, 98, 433-436. (21) Minary, P.; Martyna, G. J.; Tuckerman, M. E. Algorithms and Novel Applications Based on the Isokinetic Ensemble. I. Biophysical and Path Integral Molecular Dynamics. /. Chem. Phys. 2003, 118, 2510-2526. (22) Collins, P.; Ezra, G. S.; Wiggins, S. Phase Space Structure and Dynamics for the Hamiltonian Isokinetic Thermostat. /. Chem. Phys. 2010,133, 014105. (23) Morishita, T. Fluctuation Formulas in Molecular-Dynamics Simulations with the Weak Coupling Heat Bath. /. Chem. Phys. 2000, 113, 2976-2982. (24) Morishita, T. Generalized Coupling to a Heat Bath: Extension of the Gaussian Isokinetic Dynamics and Effect of Time Scaling. /. Chem. Phys. 2003, 119, 7075-7082. (25) Lemak, A. S.; Balabaev, N. K. On the Berendsen Thermostat. Mol. Simul. 1994, 13, 177-187. (26) Harvey, S. C; Tan, R. K.-Z.; Cheatham, T. E. The Flying Ice Cube: Velocity Rescaling in Molecular Dynamics Leads to Violation of Energy Equipartition./. Comput. Chem. 1998, 19, 726—740. (27) Callen, H. B. Thermodynamics and an Introduction to Thermo-statistics, 2nd ed.; John Wiley & Sons, 1985. (28) Cagin, T.; Ray, J. R. Isothermal Molecular-Dynamics Ensembles. Phys. Rev. A: At, Mol, Opt. Phys. 1988, 37, 4510-4513. (29) Shirts, R. B.; Burt, S. R; Johnson, A. M. Periodic Boundary Condition Induced Breakdown of the Equipartition Principle and Other Kinetic Effects of Finite Sample Size in Classical Hard-Sphere Molecular Dynamics Simulation./. Chem. Phys. 2006, 125, 164102. (30) Uline, M. J.; Siderius, D. W.; Corti, D. S. On the Generalized Equipartition Theorem in Molecular Dynamics Ensembles and the Microcanonical Thermodynamics of Small Systems. /. Chem. Phys. 2008, 128, 124301. (31) Siboni, N. H.; Raabe, D.; Varnik, F. Maintaining the Equipartition Theorem in Small Heterogeneous Molecular Dynamics Ensembles. Phys. Rev. E 2013, 87, 030101. (32) Cooke, B.; Schmidler, S. C. Preserving the Boltzmann Ensemble in Replica-Exchange Molecular Dynamics. /. Chem. Phys. 2008, 129, 164112. (33) Lingenheil, M.; Denschlag, R; Reichold, R; Tavan, P. The "Hot-Solvent/Cold-Solute" Problem Revisited. /. Chem. Theory Comput. 2008, 4, 1293-1306. (34) Goga, N; Rzepiela, A. J.; de Vries, A. H; Marrink, S. J.; Berendsen, H. J. C. Efficient Algorithms for Langevin and DPD Dynamics. /. Chem. Theory Comput. 2012, 8, 3637—3649. (35) Basconi, J. E.; Shirts, M. R. Effects of Temperature Control Algorithms on Transport Properties and Kinetics in Molecular Dynamics Simulations./. Chem. Theory Comput. 2013, 9, 2887—2899. (36) Nose, S. Constant Temperature Molecular Dynamics Methods. Prog. Theor. Phys. Suppl. 1991,103, 1-46. (37) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. /. Chem. Phys. 1995, 103, 4613-4621. (38) Mor, A.; Ziv, G.; Levy, Y. Simulations of Proteins with Inhomogeneous Degrees of Freedom: The Effect of Thermostats. /. Comput. Chem. 2008, 29, 1992-1998. (39) Chiu, S.; Clark, M.; Subramaniam, S.; Jakobsson, E. Collective Motion Artifacts Arising in Long-Duration Molecular Dynamics Simulations./. Comput. Chem. 2000, 21, 121—131. (40) Eastwood, M. P.; Stafford, K A; Lippert, R. A; Jensen, M. 0.; Maragakis, P.; Predescu, C; Dror, R. O; Shaw, D. E. Equipartition and the Calculation of Temperature in Biomolecular Simulations./. Chem. Theory Comput. 2010, 6, 2045-2058. (41) Sagui, C; Darden, T. A. Molecular Dynamics Simulations of Biomolecules: Long-Range Electrostatic Effects. Annu. Rev. Biophys. Biomol. Struct. 1999, 28, 155-179. (42) Wagner, J. R; Balaraman, G. S.; Niesen, M. J. M.; Larsen, A. B.; Jain, A.; Vaidehi, N. Advanced Techniques for Constrained Internal Coordinate Molecular Dynamics./. Comput. Chem. 2013,34,904—914. (43) Yan, L.; Sun, C; Liu, H. Opposite Phenomenon to the Flying Ice Cube in Molecular Dynamics Simulations of Flexible TIP3P Water. Adv. Manuf. 2013, 1, 160-165. (44) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics./. Comput. Phys. 1995, 117, 1—19. (45) We used the Nov 11, 2016 release of LAMMPS to conduct our simulations. The Gaussian thermostat was not implemented in LAMMPS, so we wrote an extension that integrates the equations of motion given by Minary et al.21 This extension was later incorporated into the LAMMPS code and made publicly available starting with the Jan 6, 2017 update as part of the "fix nvk" command. (46) Martin, G. M.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-AIkanes. /. Phys. Chem. B 1998, 102, 2569-2577. (47) Toxvaerd, S.; Olsen, O. H. Molecular Dynamics at Constant Temperature. Phys. Scr. 1990,1990, 98-101. (48) Tuckerman, M. E.; Liu, Y.; Ciccotti, G.; Martyna, G. J. Non-Hamiltonian Molecular Dynamics: Generalizing Hamiltonian Phase Space Principles to Non-Hamiltonian Systems. /. Chem. Phys. 2001, 115, 1678-1702. (49) Hess, S. Construction and Test of Thermostats and Twirlers for Molecular Rotations. Z. Naturforsch., A: Phys. Sci. 2003, 58, 377-391. (50) Manousiouthakis, V. I.; Deem, M. W. Strict Detailed Balance Is Unnecessary in Monte Carlo Simulation. /. Chem. Phys. 1999, 110, 2753-2756. (51) Chandler,D. Introduction to Modern Statistical Mechanics; Oxford University Press, 1987. J DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX Journal of Chemical Theory and Computation Article (52) Mudi, A.; Chakravarty, C. Effect of the Berendsen Thermostat on the Dynamical Properties of Water. Mol. Phys. 2004, 102, 681-685. (53) Shirts, M. R. Simple Quantitative Tests to Validate Sampling from Thermodynamic Ensembles. /. Chem. Theory Comput. 2013, 9, 909-926. (54) Cheng A.; Merz, K. M. Application of the Nose-Hoover Chain Algorithm to the Study of Protein Dynamics./. Phys. Chem. 1996,100, 1927-1937. (55) Mark, P.; Nilsson, L. Structure and Dynamics of Liquid Water with Different Long-Range Interaction Truncation and Temperature Control Methods in Molecular Dynamics Simulations. /. Comput. Chem. 2002, 23, 1211-1219. (56) Rosta, E.; Buchete, N. V.; Hummer, G. Thermostat Artifacts in Replica Exchange Molecular Dynamics Simulations. /. Chem. Theory Comput. 2009, 5, 1393-1399. (57) Spill, Y. G.; Pasquali, S.; Derreumaux, P. Impact of Thermostats on Folding and Aggregation Properties of Peptides Using the Optimized Potential for Efficient Structure Prediction Coarse-Grained Model./. Chem. Theory Comput. 2011, 7, 1502-1510. (58) Bussi, G.; Parrinello, M. Stochastic Thermostats: Comparison of Local and Global Schemes. Comput. Phys. Commun. 2008,179, 26—29. (59) Tafipolsky, M.; Amirjalayer, S.; Schmid, R. Ab Initio Parametrized MM3 Force Field for the Metal-Organic Framework MOF-5. /. Comput. Chem. 2007, 28, 1169-1176. (60) Amirjalayer, S.; Tafipolsky, M.; Schmid, R. Molecular Dynamics Simulation of Benzene Diffusion in MOF-5: Importance of Lattice Dynamics. Angew. Chem., Int. Ed. 2007, 46, 463—466. (61) Smit, B.; Maesen, T. L. M. Molecular Simulations of Zeolites: Adsorption, Diffusion, and Shape Selectivity. Chem. Rev. 2008, 108, 4125-4184. (62) Seehamart, K.; Nanok, T.; Krishna, R.; van Baten, J. M.; Remsungnen, T.; Fritzsche, S. A Molecular Dynamics Investigation of the Influence of Framework Flexibility on Self-Diffusivity of Ethane in Zn(tbip) Frameworks. Microporous Mesoporous Mater. 2009,125, 97— 100. (63) Dubbeldam, D.; Walton, K. S.; Ellis, D. E.; Snurr, R. Q. Exceptional Negative Thermal Expansion in Isoreticular Metal-Organic Frameworks. Angew. Chem., Int. Ed. 2007, 46, 4496-4499. (64) Ford, D. C; Dubbeldam, D.; Snurr, R. Q. The Effect of Framework Flexibility on Diffusion of Small Molecules in the Metal-Organic Framework IRMOF-1. diffusion-fundamentals.org 2009,11, 1— 8. (65) Ponder, J. W.; Richards, F. M. An Efficient Newton-Like Method for Molecular Mechanics Energy Minimization of Large Molecules. /. Comput. Chem. 1987, 8, 1016-1024. (66) Braun, E.; Chen, J. J.; Schnell, S. K.; Lin, L.-C; Reimer, J. A.; Smit, B. Nanoporous Materials Can Tune the Critical Point of a Pure Substance. Angew. Chem., Int. Ed. 2015, 54, 14349-14352. (67) Leyssale, J.-M.; Vignoles, G. L. Molecular Dynamics Evidences of the Full Graphitization of a Nanodiamond Annealed at 1500 K. Chem. Phys. Lett. 2008, 454, 299-304. (68) Wong-ekkabut, J.; Miettinen, M. S.; Dias, C; Karttunen, M. Static Charges Cannot Drive a Continuous Flow of Water Molecules Through a Carbon Nanotube. Nat. Nanotechnol. 2010, 5, 555—557. (69) Baker, M. 1,500 Scientists Lift the Lid on Reproducibility. Nature 2016, 533, 452-454. (70) Hu, Y.; Sinnott, S. B. Constant Temperature Molecular Dynamics Simulations of Energetic Particle-Solid Collisions: Comparison of Temperature Control Methods. /. Comput. Phys. 2004, 200, 251-266. K DOI: 10.1021/acs.jctc.8b00446 J. Chem. Theory Comput. XXXX, XXX, XXX-XXX