Algorithmic Techniques for Atomistic Modeling Hasan Metin Aktulga - PowerPoint PPT Presentation

Algorithmic Techniques for Atomistic Modeling Hasan Metin Aktulga Dept of Computer Science Purdue University haktulga@cs.purdue.edu May 14, 2009 ayg: chage green to a darker color

Outline 1. Introduction: MD Methods & Reactive Force Fields 2. Our Implementation of Reactive Force Fields (ReaxFF) • General structure of a ReaxFF implementation • Implementation details of major components • Additional features 3. Applications using the Serial Version • Hexane simulations: Validation and performance analysis • Water-Silica Surface Interaction • Si/Ge Nanobar Strain Rates 4. Ongoing & Future Work 5. Parallelization of ReaxFF • Challenges • Our solution approaches 6. Current Userbase

Introduction: Atomistic Modeling Methods • Ab-initio methods: Atomistics with electronic degrees of freedom ↑ – Hartree-Fock(HF) methods → misses E corr accuracy – Post-HF methods complexity ∗ E corr incorporated at great computational expense detail comp. ∗ semi-empirical methods to the rescue demand – DFT-based methods ∗ completely different approach but similar derivations to HF theory ∗ bigger systems, longer simulation times made possible ∗ CPMD is the most popular example approx’s sys size • Classical MD methods: sim. frame – many approximations: no electronic d.o.f. ↓ – electronic effects are mimiced through parametrizations – static bonds, no reactions! – systems at nanoscale, simulation times upto hundreds of ns • Continuum Mechanics: Macroscale systems – additional approximations, no atomistic detail! – modeled using PDE’s.

Reactive Force Fields (ReaxFF): Bridging the Gap ReaxFF Classical MD reactive non-reactive advantages dynamic bonds static bonds fixed charges in general polarization with QEq (except for polarizable FF) dynamic interaction lists static lists challenges complex & costly much simpler energy & force enegies,forces formulas frequent update of charges static charges (expensive!) shorter timesteps ( ≈ 0 . 25 fs) longer timesteps ( 1 − 10 fs) Interesting applications: large system with reactions and charge-transfer Simulation of fuel cells, silica crack propogation, corrosion of silica in water, etc.

General Flowchart of Conventional MD Programs read system conf, control files ❄ initialize the simulation ❄ generate neighbor lists ❪ ❄ compute energy, forces ❄ evolve the system ❄ output energy, trajectroy

ReaxFF Flowchart ✲ read geo, control, ffield initializations compute forces ✙ ✲ generate neighbors compute bonds ✻ ❘ bonded forces ❘ update charges (qeq) ❫ ✲ nonbonded forces ✛ evolve the system ayg: What about valence corrections?

Implementation: Neighbor Generation • 3 different neighbor lists: – near nbrs for bonded forces → bond cut ≈ 4-5 ˚ A, full matrix stored – hbond list for hydrogen bonds → hbond cut ≈ 6-7.5 ˚ A, only for H – far nbrs for non-bonded forces → nonb cut ≈ 10 ˚ A, upper-half only • Bin atoms into 3D grid cells – grid cell dims ≈ 1 2 nonb cut • Verlet lists with delayed re-neighboring not implemented → little benefit • Compressed adjacency list representation . . . a n +1 . . . a n − 1 a n ✙ ✢ ❫ . . . . . . nbrs n − 1 nbrs n +1 nbrs n

Implementation: Computing Forces and Potentials • Bonded Interactions:Similar to classical MD but accounts for dynamic bonds – precursor: bond orders – lone-pair energy, over/undercoordination energies – bond energy – valence energy (with penalty & 3-body conjugation corrections) – dihedral energy (with 4-body conjugation correction) • Hydrogen Bonds – precursor: bond orders – H covalently bonded to X and interacting with Z – can be considered a bonded interaction • Non-bonded Interactions – precursor: charge equilibration (QEq) – electrostatic (Coulomb) energy, van der Waals energy = E bond + E lp + E over + E under + E val + E pen + E 3 conj E system + E tors + E 4 conj + E H − bond + E vdW + E Coulomb

Bonded Interaction: Bond Orders • Prior to bonded forces, compute bond orders based on the new near nbrs – bond list: subset of near nbrs, stored in the same way – uncorrected bond orders and derivatives – store both bo ( i, j ) and bo ( j, i ) → efficient construction of angles, dihedrals – compute bo ( i, j ) only if i < j , otherwise bo ( i, j ) =bo ( j, i ) BO σ ′ ij + BO π ′ ij + BO ππ ′ BO ′ = ij ij ! bα 3 2 rij BO α ′ ij ( rij ) = exp 4 aα 5 r 0 α BO ′ ij · f 1(∆ ′ i, ∆ ′ j ) · f 4(∆ ′ i, BO ′ ij ) · f 5(∆ ′ j, BO ′ ij ) where ∆ ′ BO ij i is the valency of atom i . =

Bonded Interaction: Bond Energy ” pbe 2 ”o n “ “ − D σ e · BO σ BO σ E bond = ij · exp 1 − p be 1 ij − D π e · BO π ij − D ππ · BO ππ e ij • the stronger the bond, the lower the associated energy • sweep over the bond list • compute bond energy between i , j only if i < j

Bonded Interaction: Lone-Pair & Over/UnderCoordination • Lone-pair energy – ∆ lp = n lp opt − n lp i i – energy associated with unpaired electrons of an atom → zero for a fully coordinated atom – single-body interaction → just sweep over atom list • Over/undercoordination energy – ideal # of bonds = # of valence electrons X – ∆ i = bo ( i, j ) − V al i j ∈ nbrs ( i ) – actual # of bonds > ideal # of bonds ( ∆ i > 0 ) → over-coordination – actual # of bonds < ideal # of bonds ( ∆ i < 0 ) → under-coordination – actual # of bonds = ideal # of bonds ( ∆ i = 0 ) → no over/undercoordination energy – functionals of ∆ i and ∆ j ’s → just sweep over atom list

Bonded Interaction: Valence Angle Energy E val = f 7( BO ij, pval 3 , pval 4) · f 7( BO jk, pval 3 , pval 4) · f 8(∆ j, pval 5 , pval 6 , pval 7) · ” 2 ff« „  “ pval 1 − pval 1 · exp − pval 2 · Θ0 − Θ ijk • Θ 0 is the ideal angle, Θ ijk is the actual angle – the closer the Θ ijk to Θ 0 , the lower the energy • f 7 ( BO ij , . . . ) and f 7 ( BO jk , . . . ) ensure E val → 0 as BO ij → 0 or BO jk → 0 • {∀ x, y ∈ bo list i | x < y } that meet certain criteria, compute the energy of < x, i, y • {∀ x, y ∈ bo list i | x > y } copy < y, i, x into < x, i, y (for dihedrals!) • E pen and E 3conj account for corrections in special cases . . . a i a i +1 . . . atom list bond list ✾ q . . . . . . . . . . . . bo i,j 0 bo i,jn bo i,jlast bonds i +1 3body list ✰ ❥ q ◆ . . . . . . . . . . . . < j 0 , i, x < jlast, i, x < jn, i, x < x, i + 1 , y ✛ ✲ ✛ ✲ x ∈ bo list i x, y ∈ bo list i +1

Bonded Interaction: Dihedral Energy 1 = 2 · f 10 ( BO ij , BO jk , BO kl , p tor 2 , 1) · sin Θ ijk · sin Θ jkl · V 123 ( ω ijkl ) E tors • Dihedral angle ω ijkl is the angle between planes defined by positions of i, j, k and j, k, l • f 10 ( BO ij , BO jk , BO kl , . . . ) ensure E tors vanishes smoothly as any of these bonds dissociate • sin Θ ijk and sin Θ jkl ensure that E tors → 0 as Θ ijk → 0 or Θ jkl → 0 • {∀ i, j, k, l ∈ atoms | j < k, < i, j, k ∈ 3body list jk , < j, k, l ∈ 3body list kj } compute the energy associated with ω ijkl • weak but very important in determining the 3D structures • no higher order interactions → no storage necessary

Hydrogen Bonds 8 9 0 @ r 0 1 „ Θ XHZ « + rHZ E XHZ = phb 1 · f 7( BO XH, phb 2 , 1) · sin 4 < = hb · exp : − phb 3 · − 2 A r 0 2 rHZ ; hb • Constraints of a hydrogen bond: – Middle atom must be H – X, Z must be one of N, O, P, F – X − H covalently bonded, Z ∈ hbond list H • f 7 ( BO XH , . . . ) ensure E val → 0 as the covalent bond breaks Θ XHZ • sin 4 ( ) maximized when Θ XHZ = π ensures alignment on a line 2 • crucial for accurately describing water, DNA structure, secondary structures in proteins, etc.

Nonbonded Interaction: Charge Equilibriation (QEq) Minimizing electrostatic energy by redistributing (partial) charges. AQA + 1 ( EA 0 + χ 0 2 J 0 AAQ 2 X X Minimize E ( Q 1 . . . QN ) = A ) + ( JABQAQB ) A A<B N X subject to Qnet = Qi i =1 • Solve the optimization problem using the method of Lagrange multipliers – gives a sparse linear system of equations for finding charges • GMRES with restarts, GMRES (50) – heavy diagonal → diagonal preconditioner – little configurational change between steps → initial guess q t = linear extrapolation ( q t − 1 , q t − 2 ) • Implemented GMRES both with MGS and Householder orthogonalizations – virtually no difference • Implemented CG for comparison – GMRES takes fewer number of matvecs and is faster than CG • Choose tolerance for the norm of the relative residual carefully: – too high → wrong results! – too low → QEq dominates the total computation time! – more on this later . . .