Simulating Matter with Molecular Dynamics (Straightforward, Obvious, - PowerPoint PPT Presentation

Theory Potential Break Thermo PBC Verlet Implement Trajectories Simulating Matter with Molecular Dynamics (Straightforward, Obvious, Ridiculously Effective) Rubin H Landau Sally Haerer, Producer-Director Based on A Survey of Computational Physics by Landau, Páez, & Bordeianu with Support from the National Science Foundation Course: Computational Physics II 1 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories What Can You Do with Molecular Dynamic? Problem Will a collection of argon molecules form a liquid and then an ordered solid as the temperature is lowered? Recall introductory chemistry, “ideal” gas Derived PV = nRT via billiard ball collisions off only walls Now nonideal: molecule-molecule interactions Argon: inert gas (closed e shell) ≃ hard spheres 2 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories Molecular Dynamics (MD) Theory MD Overview MD: simulate phys & chem properties: solids, liquids, amorphous, bio MD: F = ma for bulk properties QM = correct description, but ... Bulk: not small- r behaviors QM (DFT): derive V eff ( r ) “HS physics problem from hell” Can’t run 10 23 –10 25 particles Can ∼ 10 6 : proteins; ∼ 10 8 : materials 3 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories Relation MD to Monte Carlo (MC) Thermodynamics waterpermeation.mpg 2 Powerful Simulation Techniques Both: large N particles Both: initial arbitrary, equilibrate MD: microcanonical ensemble: E , V , N fixed MC: canonical ensemble: heat bath, fixed T , N MD: has dynamics ( F = ma ); MC: not, random MD: x i , v i change continuously, calc thermo variables 4 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories Applying Newton’s Laws ( E Determines Visiblitity) Atom-Atom Interactions + + 1st Prin: 18 e - 18 e, Z + V c m d 2 r i dt 2 = F i ( r 0 , . . . , r N − 1 ) (1) ∼ 1000 e-e, e-Z F i = − ∇ r i U ( r 0 , r 1 , . . . ) (2) Ignore internal � Conservative ( E o ), U = u ( r ij ) (3) i < j central V ( r ) N − 1 m d 2 r i du � dt 2 = − (4) r ij = | r i − r j | = r ji dr ij i < j = 0 5 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories Phenomenological Lennard-Jones Potential � Short-Range Repulsion + Long Attraction �� σ � σ � 6 � � 12 u ( r ) = 4 ǫ − r r Lennard-Jones f ( r ) = − du repulsive u(r) dr �� σ � σ � 6 � = 48 ǫ r � 12 − 1 r 2 r 2 r attraction r Attractive r − 6 ǫ ↔ strength Large r van der Waals σ ↔ length scale Weak induced Repulsive r − 12 dipole–dipole e-e small r overlap ⇐ , ⇐ , ⇐⇐ 6 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories Time for A Break 3 3 3 2 2 2 1 1 1 4 4 4 5 5 5 3 3 3 2 2 2 1 1 1 4 4 4 5 5 5 3 3 3 2 2 2 1 1 1 4 4 4 5 5 5 7 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories Calculation of Thermodynamic Variables from MD Large N , Use Statistical Mechanics Equipartition, equilibrium @ T , E = 1 2 k B T / o F Simulation: just translational KE (no rotate, vibrate) � N − 1 � KE = m � v 2 (1) i 2 i = 0 • Time average KE (3 o freedom) ⇒ T = 2 � KE � MD � KE � = N 3 2 k B T (2) 3 k B N • Pressure via Virial theorem, w = <force> ρ = N ρ P = 3 N ( 2 � KE � + w ) , (3) V 8 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories Setting Initial Velocities for Molecules 1D 2D Initial Conditions Should Not Matter Start: velocity distribution characteristic some T ∝ � KE � � = true T since not equilibrium (KE ↔ PE) Randomness: just to speed up calc � 12 1 Random Gaussian = i = 1 r i 12 9 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories Periodic Boundary Conditions (PBC) EG: 1000 particles 10 × 10 × 10 cube 3 3 3 2 2 2 ⇒ 10 3 –8 3 = 488 surface 1 1 1 4 4 4 5 5 5 10 6 particles → 6% 3 3 3 2 2 2 1 1 1 4 4 4 PBC: min surface effects 5 5 5 3 3 3 2 2 2 Replicate box to infinity 1 1 1 4 4 4 5 5 5 Continuous at edges N = 10 3 –10 6 = bulk? Each ∆ t , outside box? Bring image in opposite side: Must have box!  Walls: surface effects x + L x , if x ≤ 0 ,  x ⇒ x − L x , if x > L x ↓ N ⇒ ↑ surface/volume  10 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories Potential Cutoff Computers & Boxes are Finite m i interacts with all m j 3 3 3 2 2 2 1 4 1 4 1 4 5 5 5 ⇒ ∞ interactions 3 3 3 2 2 2 1 1 1 4 4 4 5 5 5 Yet V ( r >> σ ) ≃ 0 3 3 3 2 2 2 1 1 1 4 4 4 5 5 5 V ( r = 3 σ ) ≃ V ( 1 . 13 σ ) / 200 Cut off: u ( r > 2 . 5 σ ) ≡ 0 Lennard-Jones ⇒ only nearest image interaction repulsive u(r) Problem: f = − du ( r cut ) = −∞ dr Small effect since V ( r ) small attraction r 11 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories Verlet & Velocity-Verlet Algorithms: Integrate Eq of Mtn Realistic MD: 3-D Eq of M, 10 10 t’s, 10 3 –10 6 particles rk4 ODE solver good; need quicker, simpler, built-in Verlet: central-difference f ” f i [ r ( t ) , t ] = m (= 1 ) d 2 r i dt 2 ≃ r i ( t + h ) + r i ( t − h ) − 2 r i ( t ) (1) h 2 r i ( t + h ) ≃ 2 r i ( t ) − r i ( t − h ) + h 2 f i ( t ) + O ( h 4 ) ⇒ (2) Efficient: no solve for v ’s; Determine v at end Velocity-Verlet (> stable); FD r & v together: r i ( t ) + h v i ( t ) + h 2 2 f i ( t ) + O ( h 3 ) r i ( t + h ) ≃ (3) v i ( t ) + h a ( t ) + O ( h 2 ) v i ( t + h ) ≃ (4) 12 / 14

Theory Potential Break Thermo PBC Verlet Implement Trajectories 1-D Implementation: PE, KE, E vs Time V , T ⇒ liquid, solid MD.py Energía (J) Start FCC lattice site: V LJ equilibrium 1 FCC particles/cell: 4 N 3 = 32, 108, . . . 2 Start Maxwellian velocity distribution 3 Highlight: ODE, PBC, image, cut off T 4 Energy vs Time ∼ 10 4 –10 5 steps equilibrate (10 − 9 s) for 300 particles 2D box, initialy at 150 k 5 8.00E–19 6.00E–19 Choose largest h that’s stable 6 4.00E–19 2.00E–19 0.00E+00 Compare to Figures Energy (J) 7 0 2E–13 4E–13 6E–13 8E–13 1E–12 1.2E–12 1.4E–12 –2.00E–19 –4.00E–19 Evaluate < E > t , final T –6.00E–19 8 –8.00E–19 –1.00E–18 Relation final and initial T ’s? –1.20E–18 9 –1.40E–18 13 / 14 Time (sec) (568 steps)

Theory Potential Break Thermo PBC Verlet Implement Trajectories Trajectory Analysis (Time To Get To Work) Final Temperature (K) 300 2 200 P 100 1 0 0 0 100 E-19 200 E-19 0.1 0.2 0.3 Initial KE (j) T Output Several Particles’ x & v for Every 100 Steps Start with 1-D, T = 0 @ equilibrium, then 2-D 1 PV = NRT for ideal ( V = 0) gas 2 Increase T , note motion, interactions 3 Create an animation 4 Plot displacements R rms vs T 5 Test for time-reversal invariance 6 14 / 14