Canonical ensemble, Thermostats Lukáš Sukeník March 13, 2017 Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 1 / 31 Canonical ensemble Ensemble in statistical mechanics Represents the possible states of a mechanical system In thermal equilibrium with a heat bath(fixed temperature) The system can exchange energy with the heat bath The states of the system will differ in total energy Describe Boltzmann distribution Principal variables Thermodynamic - absolute temperature (symbol: T) Determine the probability distribution of states Mechanical - number of particles (N), system’s volume (V) Influence the nature of the system’s internal states Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 2 / 31 Thermostats in Molecular dynamics Newton equations → NVE Reality → NpT In vast majority of system → little difference between NVT and NpT Thermostats Modulate the temperature of a system Ensure that the average temperature of a system is the desired one For NVT - Couple the system to a heat bath that imposes desired T Deterministic: Velocity rescale Nose-Hoover Berendsen Stochastic: Langevin Anderson Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 3 / 31 Anderson thermostat Coupling by stochastic collisions Act periodically on a random particle Instantaneous event Between stochastic collisions, system evolves in NVE Simulation Select two parameters T and f T - desired temperature of the system f - frequency of stochastic collisions, strength of coupling to heat bath The simulations proceeds in NVE until a stochastic collision Particle suffering a ”collision” Given a random Momentum from a Boltzmann distribution at T repeat with frequency f Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 4 / 31 Anderson thermostat Newtonian dynamics + stochastic collisions Turns MD simulation into a Markov process Canonical distribution in phase space is invariant under repeated collisions Anderson’s algorithm generates canonical distribution If Markov chain irreducible and aperiodic Disadvantages Algorithm randomly decorrelates velocities Dynamics is not physical Can’t measure dynamical properties Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 5 / 31 Langevin thermostat Motion of large particles through a continuum of smaller particles Langevin equation ¨r = φ − γ˙r + σξ φ - force from positions of particles similar term to ∂V ∂qi damping force γ˙r σξ - Random force The smaller particles move with kinetic energy Give random nudges to large particles Fluctuation-dissipation relation σ2 = 2γmi kB T To recover the canonical ensemble distribution Langevin thermostat Use Langevin equation Assume that atoms being simulated are embedded in a sea of much smaller fictional particles Instances of solute-solvent systems Solute is desired Solvent is non-interesting Solvent influences Solute via random collisions and a frictional force Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 6 / 31 Langevin thermostat Langevin thermostat At each time step ∆t the Langevin thermostat changes the equation of motion so that the change in momenta is ∆pi = (∂φ(q) ∂qi − γpi + δp)∆t γpi - damp the momenta δp - Gaussian distributed random number represents the thermal kicks from the small particles Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 7 / 31 Langevin thermostat Advantages Fewer computations per time step since we eliminate many atoms Include them implicitly by stochastic terms Relatively large time step ∆t - different fastest frequency motions, slower degree of freedom Disadvantages Excluded volume effects of solvent not included Not trivial to implement drag force for non-spherical particles Solute-solvent system, solvent molecules must be small compared to the smallest molecules explicitly considered Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 8 / 31 Berendsen thermostat Main problem of velocity-rescaling method Does NOT allow T fluctuations as in NVT Berendsen thermostat Weak coupling method to external heat bath Corrects deviations of actual T to T0 by multiplying the velocities by a factor λ Allows the temperature fluctuations Tries to minimize local disturbances (like stochastic thermostat does) while keeping the global effects unchanged Proportional time-rescaling Velocities scaled at each time step Rate of change of T is proportional to the difference in temperature Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 9 / 31 Berendsen thermostat Proportional time-rescaling dT dt = 1 τ (T0 − T) τ coupling parameter Exponential decay of the system towards the desired temperature This lead to a modification of the momenta pi → λpi Rescaling parameter λ λ2 = 1 + ∆t τT (T0 T − 1) Note: Velocity rescale λ2 = T0 T Drawbacks Cannot be mapped onto a specific thermodynamic ensemble Interpolation between the canonical and microcanonical ensemble ∆t = τT , fluctuations of Ekin vanish and phase space distribution reduces to NVT τT → ∞, corresponds to an isolated system (NVE) Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 10 / 31 Nose-hoover thermostat Idea Simulate a system which in the NVT ensemble Introduce a fictitious dynamical variable = friction Friction slows down or accelerates particles Measure kinetic energy and energy given by bipartition theorem 1 2 KB T per degree of freedom Scale velocities of particles so that we have desired T Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 11 / 31 Nose-hoover thermostat Varibles ri , - positions vi = ˙ri = dri dt , - velocities pi = mi · ri , - momentum ˙pi = mi ˙vi = mi ai , - force Dynamical equations ˙ri = pi mi ˙pi = fi − pi · ζ(t) fi = −∂V (q) ∂qi ζ physical meaning friction, changes with time ˙ζ = 1 Q N i=1 mi · v2 i 2 − 3N+1 2 kBT Q determines the relaxation of the dynamics of the friction, heat-bath mass, large Q denotes weak coupling T denotes the target temperature Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 12 / 31 Harmonic oscilator, r(0) = 0, p(0) = 1, ζ Q (0) = 1 or 10 Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 13 / 31 Nose-hoover thermostat Wrong behavior - Solution Need invariant probability distribution Nose-Hoover chain method ˙pi = fi − pi · ζ1(t) ˙ζ1 = 1 Q1 N i=1 mi · v2 i 2 − 3N+1 2 kBT + p2ζ2 ˙ζj = 1 Qj Qj−1ζ2 j−1 − 1 2kBT + pj+1ζj+1 Thermostat masses affect dynamics in achieving canonical distribution Large masses - microcanonical distribution (NVE) Small masses - fluctuations of the momenta greatly inhibited Q1 = 3N+1 2 kB T/ω2 Qj = 1 2 kB T/ω2 Allows the thermostats to be in approximate resonance with both the system variables, which are assumed to have fundamental frequency ω, and each other Mass choice in chain method less critical than in single metllod Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 14 / 31 Harmonic oscilator, r(0) = 0, p(0) = 1, ζ Q (0) = 1 Chain dynamics (M=2) Distribution good approximation to NVT Dynamics fills phase space Changes in the initial conditions dont have an appreciable effect on the results Choice of thermostat mass is not critical Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 15 / 31 The Lyapunov exponent The Lyapunov exponent Measure of the degree of chaos present in a dynamical system More chaotic the dynamics of a system → the more quickly it fills phase space Calculation Systems containing M = 1-15 thermostats Wide variety of initial conditions (Q= 1) Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 16 / 31 Lyapunov exponent of harmonic oscillator as a function of the number of thermostats, Q=1 Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 17 / 31 Nose-hoover thermostat - Conclusion Thermostating the extended variable Stiff complex systems (proteins) Difficult to start near equilibrium Large unphysical oscillations in T may develop Additional thermostats will effectively damp such oscillations More stable simulations Summary Very good approximation to the canonical ensemble even in pathological cases Wide application Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 18 / 31 Thermostat Artifacts in Replica Exchange Molecular Dynamics Simulations Replica exchange molecular dynamics (REMD) simulations Enhance the conformational sampling of molecular dynamics Several “replicas” simulated in parallel at different T At regular intervals - attempts to exchange replicas to increase conformational sampling efficiency at lower T After accepted exchange - Particle velocities: Reassigned from a Maxwell-Boltzmann distribution at T Old velocities are scaled (T1/T2)( 1 2 ) and vice versa Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 19 / 31 REMD Berendsen Thermostat Dont produce correct NVT Detailed balance condition of replica exchange is not satisfied Protein folding in water Helix-forming peptide with a weak-coupling Berendsen thermostat Conformational equilibrium is altered Folded state is overpopulated by about 10 % at low T Underpopulated at high T Enthalpy of folding deviates by almost 3 kcal/mol Non-canonical ensembles with narrowed potential energy fluctuations Artificially bias toward replica exchanges between low-energy folded structures at high T and high-energy unfolded structures at low T Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 20 / 31 Folding probabilities Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 21 / 31 REMD Folding probabilities For NVT REMD does not affect folding populations Berendsen - REMD alters the relative populations Folded states become overpopulated at low T Underpopulated at high T Replica exchange molecular dynamics (REMD) simulations Thermostats producing incorrect canonical ensemble REMD distort the configuration-space distributions For REMD use only thermostats correctly representing NVT Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 22 / 31 Simple Quantitative Tests to Validate Sampling from Thermodynamic Ensembles Some aspects of molecular distributions can be checked directly NVE Total energy is conserved with statistically zero drift NVT Potential energy independent on particle momenta (except in systems with magnetic forces) Kinetic energy will follow the Maxwell–Boltzmann distribution Consistency of distribution, estimated using standard statistical methods Average kinetic energy corresponding to the desired T NPT Proper average instantaneous pressure computed from the virial and kinetic energy Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 23 / 31 Difficult tests Proper distribution for the potential energy Proper distribution for total energy of an arbitrary simulation system Many possible distributions which have the correct average temperature or pressure but do not satisfy the proper Boltzmann probability distributions for our specific ensemble of interest Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 24 / 31 Basis of the ensemble validation techniques Thermodynamic ensembles all have similar probability distributions with respect to macroscopic intensive parameters and microstates P(x|β) ∝ exp(βH(p, r)) canonical where P(x|y) indicates the probability of a microstate determined by variable(s) x given a macroscopic parameter(s) y The probability density of observing a specific energy in the canonical ensemble P(E|β) = Q(β)−1 Ω(E)exp(−βE) Ω(E) density of states, Q canonical partition function No knowledge of Ω(E) distribution, cant identify proper distribution Ratio of distributions from two simulations performed at different T Two different β, other parameters the same Unknown Ω(E) cancels lnP(E|β1) P(E|β2) = [β2A2 − β1A1] − [β2 − β1] E which is of the linear form α0 + α1E. Note that linear coefficient α1 = − [β2 − β1] is independent of the (unknown in general) Helmholtz free energies A2 and A1. Helmholtz free energy A = −β−1 lnQ Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 25 / 31 General and easy Can be derived for any of the standard thermodynamic ensembles Require only energy(NVT), volume(NpT) and particle numbers(µVT) NVT Bin Total energies Distributions must be sufficiently close together Statistically well-defined probabilities at overlapping values of E Fit the ratio of the histogram probabilities to a line in overlap region Slope deviates from −(β2 − β1) → not canonical distribution Shortcomings Necessary test for canonical distribution Not sufficient - not a direct test of ergodicity No info whether states of same energy sampled with equal probability No info whether there are states that are not sampled Can be trapped in a portion of allowed phase space Additional tests of convergence or ergodicity required Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 26 / 31 Example Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 27 / 31 Ensemble validation of water simulations 900 TIP3P water molecules in NVT Nosé – Hoover algorithm Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 28 / 31 Lennard-Jones system a - Berendsen thermostat Slope of energy ratios 7 times higher than expected Low β (high T) simulation over-samples that particular kinetic energy Incorrect, overly narrow kinetic energy distribution b - Nosé – Hoover thermostat Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 29 / 31 Conclusion Validity checks Molecular distributions characterized by Boltzmann distributions Easily check for consistency with the intended ensemble Robust and general method Regardless of the details of a simulation Require only 2 simulations with differing external parameters such as T, p or µ Cancel out system-dependent properties(densities of states) Result in linear relationship between the distributions of extensive quantities (energy, volume, enthalpy and number of particles) Slope of the relationship completely determined by the intensive variables set by user Necessary, but not sufficient condition Ergodicity Full sampling of phase space Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 30 / 31 The End Lukáš Sukeník (MU) NVT, Thermostats March 13, 2017 31 / 31