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

– Typeset by FoilTEX –

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 complexity detail comp. demand

approx’s sys size sim. frame

↓

– Post-HF methods ∗ E_{corr} incorporated at great computational expense

∗ semi-empirical methods to the rescue

– 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

• Classical MD methods:

– 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

advantages

reactive

non-reactive

dynamic bonds

static bonds

polarization with QEq fixed charges in general (except for polarizable FF)

challenges

dynamic interaction lists

static lists

complex & costly enegies,forces much simpler energy & force formulas

frequent update of charges (expensive!)

static charges

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

[figure: read system conf, control files; initialize the simulation; generate neighbor lists; compute energy, forces; evolve the system; output energy, trajectroy]

ReaxFF Flowchart

[figure: read geo, control, ffield; initializations; generate neighbors; compute bonds; bonded forces; update charges (qeq); nonbonded forces; evolve the system; compute forces]

ayg: What about valence corrections?

– Typeset by FoilTEX – 5

Implementation: Neighbor Generation

• 3 different neighbor lists:

– near nbrs for bonded forces → bond cut ≈ 4-5 Å, full matrix stored

– hbond list for hydrogen bonds → hbond cut ≈ 6-7.5 Å, only for H

– far nbrs for non-bonded forces → nonb cut ≈ 10 Å, upper-half only

• Bin atoms into 3D grid cells

– grid cell dims ≈ \frac{1}{2}nonb cut

• Verlet lists with delayed re-neighboring not implemented → little benefit

• Compressed adjacency list representation

[figure: …; a_{n−1}; a_{n}; a_{n+1}; …; …; …; nbrs_{n−1}; nbrs_{n}; nbrs_{n+1}]

– Typeset by FoilTEX – 6

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_{system} = E_{bond} + E_{lp} + E_{over} + E_{under} + E_{val} + E_{pen} + E_{3conj}
+ E_{tors} + E_{4conj} + 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{′ \atop ij} = BO{σ′ \atop ij} + BO{π′ \atop ij} + BO{ππ′ \atop ij}
BO{α′ \atop ij}(rij) = exp [aα (\frac{rij}{r0α} )bα]
BOij = BO{′ \atop ij} · f1(Δ′i, Δ′j) · f4(Δ′i, BO{′ \atop ij}) · f5(Δ′j, BO{′ \atop ij}) where Δ′i is the valency of atom i.

Bonded Interaction: Bond Energy

E_{bond} = −D_{e}^{σ} · BO_{ij}^{σ} · exp {p_{be1} (1 − (BO_{ij}^{σ})^{pbe2})}
−D_{e}^{π} · BO_{ij}^{π} − D_{e}^{ππ} · BO_{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

– Δ_{i}^{lp} = n_{opt}^{lp} − n_{i}^{lp}

– 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

– Δ_{i} = ∑_{j∈nbrs(i)} bo (i, j) − V al_{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 \atop val} = f7(BOij, pval3, pval4) · f7(BOjk, pval3, pval4) · f8(Δj, pval5, pval6, pval7) ·

(pval1 − pval1 · exp {−pval2 · (Θ0 − Θijk)2})

• Θ_{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

[figure: …; ai; ai+1; …; …; …; bonds i+1; …; bo i,jn; …; bo i,jlast; atom list; bond list; 3body list; …; …; < j0, i, x; …; < jn, i, x; …; < jlast, i, x; < x, i + 1, y; x ∈ bo listi; x, y ∈ bo listi+1; bo i,j0]

Bonded Interaction: Dihedral Energy

E_{tors} = \frac{1}{2} · f_{10}(BO_{ij}, BO_{jk}, BO_{kl}, p_{tor2}, 1) · sinΘ_{ijk} · sinΘ_{jkl} · V_{123}(ω_{ijkl})

• 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

EXHZ = phb1 · f7(BOXH, phb2, 1) · sin4 (\frac{ΘXHZ}{2} ) · exp \{−phb3 · (\frac{r{0 \atop hb}}{rHZ} + \frac{rHZ}{r{0 \atop hb}} − 2)\}

• 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

• sin^{4}(\frac{ΘXHZ}{2} ) maximized when Θ_{XHZ} = π ensures alignment on a line

• crucial for accurately describing water, DNA structure, secondary structures in proteins, etc.

Nonbonded Interaction: Charge Equilibriation (QEq)

Minimizing electrostatic energy by redistributing (partial) charges.

Minimize E(Q1 … QN) = ∑ ; A(EA0 + χ0 AQA + \frac{1}{2}J{0 \atop AA}Q2 A) + ∑ ; A<B(JABQAQB)

subject to Qnet = N ; ∑ ; i=1 Qi

• 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 …

Nonbonded Interactions: Coulomb & van der Waals Energy

ECoulomb = C · T ap(rij) · \frac{qi · qj}{[r{3 \atop ij} + γ{−3 \atop ij}]\frac{1}{3}}
EvdW aals = T ap(rij) · Dij ·

[exp \{αij · (1 − \frac{f13(rij)}{rvdW})\} − 2 · exp \{\frac{1}{2} · αij · (1 − \frac{f13(rij)}{rvdW})\}]

• Shielding prevents energies from increasing drastically at close distances

• Long range interactions with cutoffs

– Taper term ensures smooth vanishing of energies after the cutoff

• no 1 − 2, 1 − 3 or 1 − 4 exclusions → smooth bond forming/breaking

• Takes up a large portion of the total computation time

– tabulate long range energy & forces

– linear interpolation

– large table → good approximations, reasonable memory usage

– more on our gains later …

Summing Alltogether: Net Force

• Let E_{i} be the sum of energies from all interactions involving atom i

• Let r_{i} denote the position of atom i

• F_{i} = \frac{∂Ei}{∂ri}

• Problem:

– bonded energy expressions include BO_{ij}(BO′_{ij}, Δ′_{i}, Δ′_{j}) terms

– \frac{∂BOij}{∂rk} arise in every bonded interaction

– \frac{∂BOij}{∂rk} = c_{1} · \frac{∂BO{′ \atop ij}}{∂rk} + c_{2} · \frac{∂Δ′i}{∂rk} + c_{3} · \frac{∂Δ′j}{∂rk}

– \frac{∂BOij}{∂rk} ≠ 0, ∀k ∈ bonds _{i} ∪ bonds _{j}, huge memory overhead!

– even if we choose to store them, very time consuming to compute each single bonded interaction!

Summing Alltogether: Net Force - Solution

• idea: distribution law of multiplication over summation

• let C_{0}, …, C_{n} be the coefficients of \frac{∂BOij}{∂rk} arising in different interactions

• re-write ∑_{t} C_{t} × \frac{∂BO_{ij}}{∂r_{k}} as \frac{∂BOij}{∂rk} × ∑_{t} C_{t}

• while computing interactions, accumulate C_{t}’s in C_{ijk}

• delay computation of \frac{∂BOij}{∂rk} × ∑_{t} C_{t} and forces due to them until C_{ijk}’s are determined

• no additional storage, important savings in CPU time

Implementation: Additional Features

• Modular implementation

– a different force field can be adopted by plugging-in new interaction routines

• NVE, NVT and NPT ensembles ayg: explain these terms first and what it takes to implement them

• Compressible custom trajectory format

• Tools for performing common analysis

– detection of reactions (on-the-fly)

– property calculations such as drift coefficient, dipole moment (on-the-fly)

– distributions of bond lengths, strengths, valence angles, charges, etc. (over the trajectory file)

Applications

• Hexane Simulations: Validation and Performance Analysis

– Preparation of systems

– Effects of QEq tolerance on accuracy and performance

– Effects of tabulating long range interactions on accuracy and performance

– Hexane structure verification

– Scalabilty of ReaxFF compared to ab-initio and classical MD

• Corrosion of silica surface in water (in collaboration with Dr Pandit’s group)

• Measuring the strain tensor of Si/Ge nanobar (in collaboration with Dr Strachan’s group)

Hexane Simulations: Preparation od the Systems

• Hexane: C_{6}H_{14}, hydrocarbon, constituent of gasoline

• Initial configuration setup:

– Very large box compared to the ideal volume → reduces overlaps!

– Randomly spread copies of a model hexane molecule

– Rotations of the model molecule around x,y and z axis to increase randomness – Various system sizes for scalability analysis (343,512,1000,1728, 3375 molecules)

• Energy minimization and NPT simulations using Gromacs

– brings the systems to the ideal volume quickly

– output to be used for ReaxFF studies and scalability analysis

• Energy minimization and NVT equilibration using ReaxFF

– added H to Gromacs output conf using the Avogadro program

– energy minimization for 2.5 ps

– NVT equilibration at 200 K for 2.5 ps

– QEq tolerance set to 1e − 8 to be safe

Hexane Simulations: Effect of QEq Tolerance on Accuracy

• Chosen: hexane system with 343 molecules = 6860 atoms

• Restart from the system equilibrated at 200 K

• How to determine the “right” tolerance:

– observe how the same system evolves over time at different QEq tolerances

– pick the highest one with reasonable accuracy

– tol_{1} = 1e − 3, tol_{2} = 1e − 4, tol_{3} = 1e − 8 → control run

Hexane Simulations: Effect of Tabulation on Accuracy

tol_{2} = 1e − 4 looks good enough, now turn on tabulation of long range interactions, too!

More in depth comparison shows they are almost identical:

property  tol3 = 1e − 8  tol2 = 1e − 4  tol2 = 1e − 4 with opt.
C-H bond  1.09 ± 0.01  1.09 ± 0.01  1.09 ± 0.01
C-C bond  1.57 ± 0.01  1.57 ± 0.01  1.57 ± 0.01
<C-C-C  108.0 ± 2.9  107.9 ± 2.9  108.0 ± 2.9
<C-C-H  111.0 ± 0.0  111.0 ± 0.0  111.0 ± 0.0
<H-C-H  106.6 ± 0.0  106.6 ± 0.0  106.6 ± 0.0
qC-tip  −0.171  −0.171  −0.171
qC-mid  −0.080  −0.080  −0.080
qH-tip  0.040  0.040  0.040
qH-mid  0.040  0.040  0.040

Hexane Simulations: Profiling Analysis & Scalability

• Different qeq tolerances, with/without tabulation

• Used the head.cs cluster except for the first case

ayg: What are the numbers here.. are they times in seconds?

total neighbors  bonded  nonb  QEq matvecs  QEq%
tol=1e − 4 w/opt§  0.93  0.22  0.14  0.22 0.34  9.9  37%
tol=1e − 4 w/opt  3.26  0.47  0.32  0.65 1.78  9.9  55%
tol=1e − 4  4.41  0.45  0.31  1.92 1.67  9.9  38%
tol=1e − 6 w/opt  3.35  0.45  0.31  0.59 1.97  13.6  59%
tol=1e − 6  4.65  0.45  0.31  1.92 1.95  13.5  42%
tol=1e − 8  7.56  0.44  0.30  1.91 4.89  46.8  65%

ayg: Dont say Dell Studio.. instead, say what the processor is, etc. §Architecture can make a huge difference. Same system with same parameters on a Dell Studio XPS with 2.67GHz quad-core i7 processor and 1066MHz memory.

• Lessons learnt:

– QEq tolerance is crucial for accuracy

– arbitrarily large QEq tolerance might cause QEq domination ayg: What does the above bullet mean?

– QEq is just a precursor to electrostatics yet takes up at least one third of total time!

– QEq must be improved to make ReaxFF scalable.

Hexane Simulations: Validation

Compare the structure of our hexane molecules to those of experimental results in the literature and ab-initio simulation:

property  ReaxFF  experimental†  ab-initio‡
C-H bond  1.09 ± 0.01 1.118 ± 0.006  1.100
C-C bond  1.57 ± 0.01 1.533 ± 0.003  1.533
<C-C-C  108.0 ± 2.9  111.9 ± 0.4  114.2
<C-C-H  111.0 ± 0.0  109.5 ± 0.5  109.5
<H-C-H  106.6 ± 0.0  NA  106.5
q_{C-tip}  −0.171  NA  −0.205
q_{C-mid}  −0.080  NA  0.033
q_{H-tip}  0.040  NA  0.047
q_{H-mid}  0.040  NA  −0.10 ∼ 0.10

†R. A. Bonham, L. S. Bartell, and D. A. Kohl. “The Molecular Structures of n-Pentane, n-Hexane and n-Heptane” J. Am. Chem. Soc., 1959, 81 (18), 4765 − 4769 ‡geometry optimization of an isolated hexane using CPMD v3.13.2 with P BE Troullier-Martins pseudopotentials

Hexane Simulations: Scalability Analysis

– Typeset by FoilTEX – 25

Silica Surface Corrosion in Water

– Typeset by FoilTEX – 26

Si/Ge Nanobar

Ongoing & Future Work

• ParallelReax

– complete & verify the implementation

– a large scale application (maybe with the PRISM device)

• Integration into LAMMPS (a well-known, widely used MD package from SNL)

– QEq integrated independently ∗ opens the door for polarizable-ff inside LAMMPS

∗ compatibility issues to be sorted out!

• Better solvers for QEq

– block Jacobi type pre-conditioner

– inner-outer schemes (use a lower cutoff inner solve to precondition outer solve)

– use a fast multipole-type preconditioner.

• Tight relations among items on the agenda

– better QEq solvers necessary for scalable ParallelReax

– better QEq solvers necessary for QEq in LAMMPS

– completion of ParallelReax necessary for Reax in LAMMPS

Parallelization of ReaxFF

• A draft version of ParallelReax

– domain decomposition technique

∗ repeat: get my share of atoms → communicate boundaries → compute forces →

move my atoms

– not fully verified

– inefficient handling of processor boundaries

Two big challenges:

• Parallelization of QEq

– even CG needs at least 4 communications per iteration!

– parallel GMRES (with Householder orth.) is even worse

– QEq will dominate even more

– definitely need: better solvers for QEq

• Processor boundaries

– avoid double computation at boundaries!

– avoid thick boundaries, retain accuracy for any system!

Parallelization: Our Solution Approaches

• Parallelization of QEq: better solvers for QEq

– block Jacobi type pre-conditioner

– inner-outer schemes

• Avoid double computation

– coordination through the mid-point rule ∗ bond(i, j): owner(\frac{1}{2}(r_{i} + r_{j}))

∗ < (i, j, k): owner(r_{j})

∗ dihedral(i, j, k, l): owner(\frac{1}{2}(r_{j} + r_{k}))

∗ hbond(X, H, Z): owner(r_{H})

∗ nonbonded(i, j): owner(\frac{1}{2}(r_{i} + r_{j}))

dihedral bond nonbonded hbond angle P0  P1 P2  P3 P4

Parallelization: Our Solution Approaches

• Why avoid thick boundaries?

– a 20 Å cubic box ≈ 1000 atoms

– assume nonb cut long boundaries → a few thousand atoms/processor → upto a couple of seconds per iteration!

– definitely need: thin boundaries

• How to avoid thick boundaries?

– Mid-point rule to the rescue!

– boundary thickness: max(\frac{3}{2}bond cut, hbond cut, \frac{1}{2}nonb cut)

– gets better if no hydrogen bonds present: max(\frac{3}{2}bond cut, \frac{1}{2}nonb cut)

worst case scenarios P1  P2  bond→ bond cut ∼ 4A angle→ bond cut ∼ 4A dihedral → \frac{3}{2}bond cut ∼ 6A hbond → hbond cut ∼ 6 − 7.5A long range → \frac{1}{2}nonb cut ∼ 5A

Conclusions

– Typeset by FoilTEX – 32

Current Userbase

• Dr Pandit’s group at USF

– silica-water systems

• Dr Strachan’s group at Purdue

– Si/Ge nanobar

• Dr Buehler’s group at MIT

– Silica cracking with strain

• Dr van Duin at PennState

• Dr Goddard’s group at Caltech

• Dr Aluru’s group at UIUC