Download
1 / 56
580 likes | 630 Views
Molecular Dynamics. Sauro Succi. 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
E N D
Molecular Dynamics Sauro Succi
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
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:
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:
MD Energy Scales Typical energy in thermal units: 6-12 Lennard-Jones potential s Standard conditions, T=300K
Classical N-body problem Newton equations for Avogadro’s number of molecules… Potential forces 3-body 2-body
Central/Directional Forces Electrostatic forces are aligned (central)
Third Law: Momentum Conservation This implies global momentum conservation:
Molecular Dynamics Distinguished Potentials
Hard-Spheres Hard-core, a rigid wall, still very useful
Lennard-Jones Hard-core repulsion (-12) plus Soft-core attraction (-6). Repulsion= electron orbitals Attraction: electrostatic dipoles (van der Waals)
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.
Long-range Coulomb potential Unscreened Coulomb interactions (one-component charged fluids): A tough cookie: special treatment of boundary conditions (see Lee-Edwards technique)
Running the MD simulation 1. Initialization 2. Time integration 3. Data Analysis
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
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]
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
Time marching: position Verlet Euler start-up: 2° order accurate; Symplectic; but velocities are staggered… Velocities:
VV is symplectic Let Velocity-Verlet: Let q.e.d VV is pretty popular: 2° order, simplectic, time-aligned
Timestep limitation: the gap Energy-conservation demands large safety margins: Molecules entering the hard-core can ruin the simulation: steep growth
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
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 }
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 }}
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};
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 }}
Equlibrium Properties Discard the transient 0<t<t_tra
Given the density: Compute the energy and the temperature: Equation of State Compute the pressure as the diagonal component of the momentum-flux tensor
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.
Radial Correlation Function of liquids Can be computed by simpy binning the N(N_1)/2 distances r_ij
Micro-canonical: constant energy NVE Ensemble
Canonical NVT Ensemble The boundary exchange heat so that Temperature is kept constant
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
Mimick a wall collision Choose one particle at random and make it thermal Andersen rescaling Rescale the other N-1 particles (every time-step):
Constrained-Ensembles Most useful for biochemical applications
Fixed-length bonds Holonomic Constraint: the bond length is constant in time The l-th and m-th particles participate to constraint k
Constrained-Ensembles Under ordinary dynamics the constraint is not fulfilled Add restoring forces which pull the system towards the manifold
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
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
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)
