1 / 56

Exploring Molecular Dynamics in Simulations: Length, Time, Energy Scales & Potential Forces

Molecular Dynamics (MD) history, length scales, time scales, energy scales, hard-sphere potentials, force calculations, data analysis methods, and ensembles in simulations. Understand the principles of MD and key techniques for efficient computations.

carolynmoore
Download Presentation

Exploring Molecular Dynamics in Simulations: Length, Time, Energy Scales & Potential Forces

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Molecular Dynamics Sauro Succi

  2. Molecular Dynamics (MD) began in the mid 50’s with the famous experiments of Alder-Wainwright. They showed violation of Boltzmann’s Molecular Chaos assumption using hundreds of rigid disks. They went on by computing a phase transition in this “virtual-fluid”. Incredible growth ever since: nowadays Multibillion sims can be performed. Still very short of real—life scales, especially in time. History

  3. MD Length Scales Mean interparticle distance: s Dense fluids: s s d s s is the “size” of the molecule, or, more appropriately, the interaction range. In MD we are talking mostly SHORT RANGE:

  4. MD Time Scales Mean interaction time: s Dense fluids: s s d s s The “size” of the molecule, or, more appropriately, the interaction range. In MD we are talking mostly SHORT RANGE:

  5. MD Energy Scales Typical energy in thermal units: 6-12 Lennard-Jones potential s Standard conditions, T=300K

  6. Classical N-body problem Newton equations for Avogadro’s number of molecules… Potential forces 3-body 2-body

  7. Central/Directional Forces Electrostatic forces are aligned (central)

  8. Third Law: Momentum Conservation This implies global momentum conservation:

  9. Energy Conservation

  10. Molecular Dynamics Distinguished Potentials

  11. Hard-Spheres Hard-core, a rigid wall, still very useful

  12. Lennard-Jones Hard-core repulsion (-12) plus Soft-core attraction (-6). Repulsion= electron orbitals Attraction: electrostatic dipoles (van der Waals)

  13. WCA potential Rep only Smooth Only repulsive part is retained, zero at r=r0=(2^1/6)*sigma Smooth at r0. Short-scale structural details are basically Dictated by repulsion. Very useful for colloids.

  14. Long-range Coulomb potential Unscreened Coulomb interactions (one-component charged fluids): A tough cookie: special treatment of boundary conditions (see Lee-Edwards technique)

  15. Running the MD simulation 1. Initialization 2. Time integration 3. Data Analysis

  16. Particles should be placed at a mean distance d above sigma: Initial Configuration: Positions Divide the system in B=(L/rc)^3 boxes and place the particle within the box (random or center, no big deal). The cutoff range rc should be several (3-5) sigmas

  17. Sample from a Gaussian with zero average and variance kT/m Initial Configuration: Velocities Sample from a gaussian (standard from libraries) See also Box-Mueller algorithm Where u1 and u2 are uniform rn in [0,1]

  18. In principle we have “just” to integrate Newton’s equations. It is of utmost importance to secure time-reversibility, i.e. Liouville theorem: Symplectic integrators. In addition, forces should be computed as efficiently as possible because they take most CPU time. Time-evolution

  19. Time marching: position Verlet Euler start-up: 2° order accurate; Symplectic; but velocities are staggered… Velocities:

  20. Velocity-Verlet

  21. VV is symplectic Let Velocity-Verlet: Let q.e.d VV is pretty popular: 2° order, simplectic, time-aligned

  22. Timestep limitation: the gap Energy-conservation demands large safety margins: Molecules entering the hard-core can ruin the simulation: steep growth

  23. This is the most CPU time consuming section of MD simulations, scaling in principle like N*N. However, short-range potentials can be dealt with more economic techniques: Force computation 1. Direct N-body summation 2. Cut-off summation 3. Linked lists

  24. Force Calculation: Direct Summation for (i=1;i<=N-1;i++) { F[i]=0.0; for (j=i+1,j<=N,j++{ xij = x[i]-x[j]; yij = y[i]-y[j]; zij = z[i]-z[j]; rij=pow(xij*xij+yij*yij+zij*zij,0.5) Fij = -a/rij**13 +b/rij**7 F[i]+=Fij; F[j]-=Fij }

  25. Force Calculation: short range But N(N-1)/2 conditional branches… for (i=1;i<=N;i++) { F[i]=0.0; for (j=i+1,j<=N,j++{ xij = x[i]-x[j]; yij = y[i]-y[j]; zij = z[i]-z[j]; rij=pow(xij*xij+yij*yij+zij*zij,0.5) if(rij<=rcut) { Fij = -a/rij**13 +b/rij**7 F[i]+=Fij }}

  26. Linked Lists Each grid cell associates with an ordered list of the particles which interact in that cell For (i=1;i<=N;i++) { g=floor(x[i]/rc); np[g]++; Linkl[g][np[g]]=i; } 6 8 3 1 7 9 5 10 2 4 Link[1]={2,3,5}; Linkl[2]={1,6,9}; Linkl[3]={4,7,8,10};

  27. Force Calculation: Linked List for (i=1;i<=N;i++) { g=floor(x[i]/rc); for (k=1,k<=npc[g],k++{ j = Linkl[g][k]; xij= x[i]-x[j]; Fij = -a/xij**13 +b/xij**7 F[i]+=Fij; F[j]-=Fij }}

  28. Data Analysis

  29. Equlibrium Properties Discard the transient 0<t<t_tra

  30. Given the density: Compute the energy and the temperature: Equation of State Compute the pressure as the diagonal component of the momentum-flux tensor

  31. Correlation Function The correlation function g(r) measures the number of pairs at a mutual distance r. In a homogeneous fluid this depends only on the distance r (Radial Distribution Function). It reveals the short-scale structure of the fluid and determines The virial coefficients of the Equation of State.

  32. Radial Correlation Function of liquids Can be computed by simpy binning the N(N_1)/2 distances r_ij

  33. Ensembles

  34. Micro-canonical: constant energy NVE Ensemble

  35. Canonical NVT Ensemble The boundary exchange heat so that Temperature is kept constant

  36. Grand-Canonical: mass exchange

  37. Thermostats

  38. This defines the kinetic temperature: Isokinetic Ensemble (T=const) The dynamics does not conserve KE, hence one must rescale each velocity to keep the reference Temperature constant: However this does not allow any temperature fluctuation, So it does not sample the canonical NVT ensemble

  39. Mimick a wall collision Choose one particle at random and make it thermal Andersen rescaling Rescale the other N-1 particles (every time-step):

  40. Canonical NVT

  41. Velocity-Verlet-Nose’-Hoover

  42. Constrained-Ensembles Most useful for biochemical applications

  43. Fixed-length bonds Holonomic Constraint: the bond length is constant in time The l-th and m-th particles participate to constraint k

  44. Constrained-Ensembles Under ordinary dynamics the constraint is not fulfilled Add restoring forces which pull the system towards the manifold

  45. SHAKE algorithm 1. Move with ordinary forces to unconstrained positions 2. Apply the constraint force 3. Solve the n nonlin eqs for the n Lagrangian multipliers This is usually done by a nonlinear (Newton-Raphson) iterator

  46. Boundary Conditions

  47. Solid Walls NVE: Particles impinging on the well are reinjected along a random direction and conserved magnitude NVT: Drawn from a Maxwellian at temperature T

  48. Many community codes: CHARMM (Chemistry at Harvard Macromolecular Mechanics) GROMACS (GROeningen Machine for Advanced Chemical Simulations) AMBER (Assisted Model Building with Energy refinement) LAMPS (Large scale Atomic/Molecular Massively Parallel Simulator) MD: state of the Art Special purpose computers ANTON (D. Shaw Research)

  49. Multi-billion MD sim’s

  50. Multi-billion MD sim’s

More Related