When is classical mechanics a reasonable approximation?

1. SMA5233 Particle Methods and Molecular DynamicsLecture 2: Physical Principles and Design Issues of MDA/P Chen Yu ZongTel: 6874-6877Email: phacyz@nus.edu.sghttp://bidd.nus.edu.sgRoom 08-14, level 8, S16 National University of Singapore

2. When is classical mechanics a reasonable approximation? • In Newtonian physics, any particle may possess any one of a continuum of energy values. In quantum physics, the energy is quantized, not continuous. That is, the system can accomodate only certain discrete levels of energy, separated by gaps. • At very low temperatures these gaps are much larger than thermal energy, and the system is confined to one or just a few of the low-energy states. Here, we expect the `discreteness' of the quantum energy landscape to be evident in the system's behavior. • As the temperature is increased, more and more states become thermally accessible, the `discreteness' becomes less and less important, and the system approaches classical behavior.

3. When is classical mechanics a reasonable approximation? • For a harmonic oscillator, the quantized energies are separated by DE=hf, where h is Planck's constant and f is the frequency of harmonic vibration. • Classical behavior is approached at temperatures for which KBT>>hf, where KBis the Boltzmann constant and KBT = 0.596 kcal/mol at 300 K. Setting hf = 0.596 kcal/mol yields f = 6.25/ps, or 209cm-1. So a classical treatment will suffice for motions with characteristic times of a ps or longer at room temperature.

4. Classical Mechanics - in a Nutshell • Building on the work of Galileo and others, Newton unveiled his laws of motion in 1686. According to Newton: • A body remains at rest or in uniform motion (constant velocity - both speed and direction) unless acted on by a net external force. • In response to a net external force, F, a body of mass m accelerates with acceleration a=F/m. • If body i pushes on body j with a force Fij, then body j pushes on body i with a force Fji=-Fij.

5. Classical Mechanics - in a Nutshell • For energy-conserving forces, the net force Fi on particle i is the negative gradient (slope in three dimensions) of the potential energy with respect to particle i's position: Fi= - ▼iV(R) where V(R) represents the potential energy of the system as a function of the positions of all N particles, R. • In three dimensions, ri is the vector of length 3 specifying the position of the atom, and R is the vector of length specifying all coordinates. In the context of simulation, the forces are calculated for energy minimizations and molecular dynamics simulations but are not needed in Monte Carlo simulations. • Classical mechanics is completely deterministic: Given the exact positions and velocities of all particles at a given time, along with the function , one can calculate the future (and past) positions and velocities of all particles at any other time. The evolution of the system's positions and momenta through time is often referred to as a trajectory.

6. Statistical Mechanics - Calculating Equilibrium Averages • According to statistical mechanics, the probability that a given state with energy E is occupied in equilibrium at constant particle number N, volume V, and temperature T (constant , the `canonical' ensemble) is proportional to exp{-E/ KBT}, the `Boltzmann factor.' Probability ~ exp{-E/ KBT} • The equilibrium value of any observable O is therefore obtained by averaging over all states accessible to the system, weighting each state by this factor.

7. Statistical Mechanics - Calculating Equilibrium Averages • Classically, a microstate is specified by the positions and velocities (momenta) of all particles, each of which can take on any value. The averaging over states in the classical limit is done by integrating over these continuous variables: • where the integrals are over all phase space (positions and momenta _) for the N particles in 3 dimensions. • When all forces (the potential energy V) and the observable O are velocity-independent, the momentum integrals can be factored and canceled: • where is the total kinetic energy, and E=K+V. As a result, Monte Carlo simulations compare V's, not E's.

8. What is molecular mechanics? • Molecular Mechanics (MM) finds the geometry that corresponds to a minimum energy for the system - a process known as energy minimization. • A molecular system will generally exhibit numerous minima, each corresponding to a feasible conformation. Each minimum will have a characteristic energy, which can be computed. The lowest energy, or global minimum, will correspond to the most likely conformation.

9. Energy Minimization Torsion angles Are 4-body Non-bonded pair Bonds Are 2-body Angles Are 3-body U is a function of the conformation C of the protein. The problem of “minimizing U” can be stated as finding C such that U(C) is minimum.

10. Units in force fields Most force fields use the AKMA (Angstrom – Kcal – Mol – Atomic Mass Unit) unit system:

11. Energy Minimization Local minima: A conformation X is a local minimum if there exists a domain D in the neighborhood of X such that for all Y≠X in D: U(X) <U(Y) Global minimum: A conformation X is a global minimum if U(X) <U(Y) for all conformations Y ≠X

12. Molecular Modeling: Energy minimization: Force guided approach: Initialize: Change atom position: xi -> xi+dxi Compute potential energy change: V -> V +dV Determine next movement: Force: Fxi=-dV/dxi; Fyi=-dV/dyi; Fzi=-dV/dzi Atom displacement: dxi=C*Fxi New position: xi=xi+dxi

13. Notes on computing the energy Bonded interactions are local, and therefore their computation has a linear computational complexity (O(N), where N is the number of atoms in the molecule considered. The direct computation of the non bonded interactions involve all pairs of atoms and has a quadratic complexity (O(N2)). This is usually prohibitive for large molecules. Reducing the computing time: use of cutoff in UNB

14. Cutoff schemes for faster energy computation wij: weights(0< wij <1). Can be used to exclude bonded terms, or to scale some interactions (usually 1-4) S(r) : cutoff function. Three types: 1) Truncation: b

15. Cutoff schemes for faster energy computation 2. Switching a b with 3. Shifting or b

16. The minimizers Some definitions: Gradient: The gradient of a smooth function f with continuous first and second derivatives is defined as: Hessian The n x n symmetric matrix of second derivatives, H(x), is called the Hessian:

17. The minimizers Minimization of a multi-variable function is usually an iterative process, in which updates of the state variable x are computed using the gradient, and in some (favorable) cases the Hessian. Iterations are stopped either when the maximum number of steps (user’s input) is reached, or when the gradient norm is below a given threshold. Steepest descent (SD): The simplest iteration scheme consists of following the “steepest descent” direction: Usually, SD methods leads to improvement quickly, but then exhibit slow progress toward a solution. They are commonly recommended for initial minimization iterations, when the starting function and gradient-norm values are very large. (a sets the minimum along the line defined by the gradient)

18. The minimizers Conjugate gradients (CG): In each step of conjugate gradient methods, a search vector pk is defined by a recursive formula: The corresponding new position is found by line minimization along pk: the CG methods differ in their definition of b: - Fletcher-Reeves: - Polak-Ribiere - Hestenes-Stiefel

19. The minimizers Newton’s methods: Newton’s method is a popular iterative method for finding the 0 of a one-dimensional function: x3 x2 x1 x0 It can be adapted to the minimization of a one –dimensional function, in which case the iteration formula is: The equivalent iterative scheme for multivariate functions is based on: Several implementations of Newton’s method exist, that avoid computing the full Hessian matrix: quasi-Newton, truncated Newton, “adopted-basis Newton-Raphson” (ABNR),…

20. What is molecular dynamics simulation? • Molecular Dynamics (MD) applies the laws of classical mechanics to compute the motion of the particles in a molecular system. • Simulation that shows how the atoms in the system move with time • Typically on the nanosecond timescale • Atoms are treated like hard balls, and their motions are described by Newton’s laws.

21. What is molecular dynamics simulation? Molecular dynamics simulations generate information at the microscopic level, including atomic positions and velocities. The conversion of this microscopic information to macroscopic observables such as pressure, energy, heat capacities, etc., requires statistical mechanics. MD Simulation Microscopic information Statistical Mechanics Macroscopic properties

22. Characteristic protein motions Periodic (harmonic) Random (stochastic)

23. Molecular Dynamics Simulation Molecule: (classical) N-particle system Newtonian equations of motion: with Integrate numerically via the „leapfrog“ scheme: with Δt  1fs! (equivalent to the Verlet algorithm)

24. How do you run a MD simulation? • Get the initial configuration From x-ray crystallography or NMR spectroscopy (PDB) • Assign initial velocities At thermal equilibrium, the expected value of the kinetic energy of the system at temperature T is: This can be obtained by assigning the velocity components vi from a random Gaussian distribution with mean 0 and standard deviation (kBT/mi):

25. How do you run a MD simulation? For each time step: Compute the force on each atom: Solve Newton’s 2nd law of motion for each atom, to get new coordinates and velocities Store coordinates Stop X: cartesian vector of the system M diagonal mass matrix .. means second order differentiation with respect to time Newton’s equation cannot be solved analytically: Use stepwise numerical integration

26. What the integration algorithm does • Advance the system by a small time step Dt during which forces are considered constant • Recalculate forces and velocities • Repeat the process If Dt is small enough, solution is a reasonable approximation

27. Choosing a time step Dt • Not too short so that conformations are efficiently sampled • Not too long to prevent wild fluctuations or system ‘blow-up’ • An order of magnitude less than the fastest motion is ideal • Usually bond stretching is the fastest motion: C-H is ~10fs so use time step of 1fs • Not interested in these motions? Constrain these bonds and double the time step

28. Molecular Dynamics Ensembles • The method discussed above is appropriate for the micro-canonical ensemble: constant N (# of particles) V (volume) and ET (total energy = E + Ekin) • In some cases, it might be more appropriate to simulate under constant Temperature (T) or constant Pressure (P): • Canonical ensemble: NVT • Isothermal-isobaric: NPT • Constant pressure and enthalpy: NPH How do you simulate at constant temperature and pressure?

29. Common Ensembles

30. MD in different thermodynamic ensemble • MD naturally samples the NVE (microcanonical ensemble): this is not very useful for most applications. • More useful ensembles include: • Canonical (NVT): requires the particles to interact with a thermostat • Isobaric-Isothermal (NPT): requires the particles to interact with a thermostat and a baristat

31. Temperature • Thermodynamics definition • Measuring T

32. Canonical (NVT) • Several possible approaches • Berendsen Method (velocities are rescaled deterministically after each step so that the system is forced towards the desired temperature) • Stochastic collisions (Anderson: a random particle’s velocity is reassigned by random selection from the Maxwell-Boltzmann distribution at set intervals) • Extended Lagrangians (Nose and Nose-Hoover: the thermal reservoir is considered an integral part of the system and it is represented by an additional degree of freedom: set of deterministic equations that give sampling in NVT)

33. system Simulating at constant T: the Berendsen scheme • Bath supplies or removes heat from the system as appropriate • Exponentially scale the velocities at each time step by the factor : where  determines how strong the bath influences the system Heat bath T: “kinetic” temperature Berendsen et al. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684 (1984)

34. system Simulating at constant P: Berendsen scheme • Couple the system to a pressure bath • Exponentially scale the volume of the simulation box at each time step by a factor : where  : isothermal compressibility P : coupling constant pressure bath where u : volume xi: position of particle i Fi : force on particle i Berendsen et al. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684 (1984)

35. Nose-Hoover Thermostat Introduce a new coordinate, s, representing reservoir in the Lagrangian Q=effective mass associated with s g is a parameter we will fix later With momenta:

36. With p’=p/s Nose-Hoover Thermostat • Extended-system Hamiltonian: • it can be shown that the molecular positions and momenta follow a canonical (NVT) distribution if g = 3N+1

37. Nose-Hoover Thermostat Real and virtual variables are related by: • s can be interpreted as a scaling factor of the time step • treal = tsim/s • since s varies during the simulation, each “true” time step is of varying length

38. Scaled-variables equations of motion constant simulation Dt fluctuating real Dt Real-variables equation of motion (here, I have Dropped the prime) Nose-Hoover Thermostat From the Hamiltonian, can derive the equations of motion for the virtual variables p, r and t.

39. So equations of motions (dropping the prime): (g=3N) IMPLEMENTATION Simplify expression (Nose-Hoover) by introducing a thermodynamic Friction coefficient:

40. An experiment is usually made on a macroscopic sample that contains an extremely large number of atoms or molecules sampling an enormous number of conformations. In statistical mechanics, averages corresponding to experimental observables are defined in terms of ensemble averages;

41. AVERAGES Thermodynamic (Ensemble) average: average over all points in phase space at a Single time Dynamic average: average over a single point in phase space at All times Ergodic hypothesis: for infinitely long trajectory, thermo average= Dynamic average A large number of observations made on a single system at N arbitrary instants in time have the same statistical properties as observing N arbitrarily chosen systems at the same time.

42. Ensemble (thermodynamic) average of observable A(rN,pN): To get a thermo average, need to know the probability of finding the system at each point (state) in phase space. Probability density Partition function

43. In an MD simulation, we determine a time average of A, which is expressed as wheretis the simulation time, M is the number of time steps in the simulation and A(pN,rN) is the instantaneous value of A.

44. THERMODYNAMIC PROPERTIES FROM MD SIMULATIONS Thermodynamic (equilibrium) averages can be calculated via time averaging: Examples: Average potential energy Average kinetic energy

45. Other quantities obtained from MD: Heat Capacity Expressed in terms of fluctuations in the energy

46. Other Quantities • Root Mean Square Deviation (for instance, over Ca atoms) • Radius of Gyration

47. From PDB orXtal Relieve local stresses, Relax bond length/angle distortions velocities (vx, vy,vz) are randomly assigned according to a Maxwellian distribution, P(vi)dv = (m/2pikBT)1/2 exp[ -mv2/2kBT ]dv to overcome “hot spots” due to Initial high velocities T(t) = (1/kB(3N-n)) Si (mivi2) Scale velocities by (T'/T(t))1/2 to get new temperatureT

48. Types of MD Systems Empirical Potentials • Empirical potentials either ignore quantum mechanical effects, or attempt to capture quantum effects in a limited way through entirely empirical equations. Parameters in the potential are fitted against known physical properties of the atoms being simulated, such as elastic constants and lattice parameters. • Empirical potentials can further be subcategorized into pair potentials, and many-body potentials • Pair potentials: the total potential energy of a system can be calculated from the sum of energy contributions from pairs of atoms. • Many-body potentials: the potential energy cannot be found by a sum over pairs of atoms. • Tersoff Potential involves a sum over groups of three atoms, with the angles between the atoms being an important factor in the potential. • Tight-Binding Second Moment Approximation (TBSMA), the electron density of states in the region of an atom is calculated from a sum of contributions from surrounding atoms, and the potential energy contribution is then a function of this sum.

49. Types of MD Systems Semi-Empirical Potentials • Semi-empirical potentials make use of the matrix representation from quantum mechanics. However, the values of the matrix elements are found through empirical formulae that estimate the degree of overlap of specific atomic orbitals. The matrix is then diagonalized to determine the occupancy of the different atomic orbitals, and empirical formulae are used once again to determine the energy contributions of the orbitals. • There are a wide variety of semi-empirical potentials, known as tight-binding potentials, which vary according to the atoms being modeled.

50. Types of MD Systems Ab-initio (First Principles) Methods • Ab-initio (first principles) methods use full quantum-mechanical formula to calculate the potential energy of a system of atoms or molecules. This calculation is usually made "locally", i.e., for nuclei in the close neighborhood of the reaction coordinate. Although various approximations may be used, these are based on theoretical considerations, not on empirical fitting. Ab-Initio produce a large amount of information that is not available from the empirical methods, such as density of states information. • Perhaps the most famous ab-initio package is the Car-Parrinello Molecular Dynamics (CPMD) package based on the density functional theory..