Improving Particle Methods Robert Krasny (Michigan) collaborators - PowerPoint PPT Presentation

March 12-16, 2018 ICERM Workshop on “Fast Algorithms for Generating Static and Dynamically Changing Point Configurations” Improving Particle Methods Robert Krasny (Michigan) collaborators Weihua Geng (Southern Methodist University) Peter Bosler (Sandia National Laboratories) Ling Xu (Michigan) support NSF DMS-1418966 ONR N00014-14-1-0075 MICDE (Michigan Institute for Computational Discovery and Engineering)

outline 1. protein-solvent electrostatics 2. incompressible fluid flow

1. protein-solvent electrostatics atomic charges Van der Waals radius water molecules molecular surface dissolved salt ions goal φ ( x ) : electrostatic potential E sol : electrostatic solvation energy

Poisson-Boltzmann implicit solvent model y k Ω 1 , ✏ 1 Ω 2 , ✏ 2 Γ N c X Ω 1 : protein domain : ✏ 1 r 2 � ( x ) = � q k � ( x � y k ) k =1 Ω 2 : solvent domain : ✏ 2 r 2 � ( x ) � ¯  2 � ( x ) = 0 : PB equation @� 1 @� 2 Γ : molecular surface : � 1 = � 2 , ✏ 1 @ n = ✏ 2 @ n | x | →∞ � ( x ) = 0 lim far-field boundary condition :

boundary integral formulation Ju ff er et al. (1991) 1 G 0 ( x, y ) = 4 ⇡ | x − y | : Coulomb potential G  ( x, y ) = e −  | x − y | 4 ⇡ | x − y | : screened Coulomb potential ,  2 = ¯  2 / ✏ 2 Z ⇣ ⌘ h i 1 1 + ✏ 2 � ( x ) = K 1 ( x, y ) @ n � ( y ) + K 2 ( x, y ) � ( y ) dS y + S 1 ( x ) 2 ✏ 1 Γ Z ⇣ ⌘ h i 1 1 + ✏ 1 @ n � ( x ) = K 3 ( x, y ) @ n � ( y ) + K 4 ( x, y ) � ( y ) dS y + S 2 ( x ) 2 ✏ 2 Γ K 2 = ✏ 2 K 1 = G 0 − G  , @ n y G  − @ n y G 0 ✏ 1 K 4 = @ 2 n x n y ( G  − G 0 ) , K 3 = @ n x G 0 − ✏ 1 ✏ 2 @ n x G  N c N c S 1 ( x ) = 1 q k G 0 ( x, y k ) , S 2 ( x ) = 1 X X q k @ n x G 0 ( x, y k ) ✏ 1 ✏ 1 k =1 k =1

discretization triangulation + centroid collocation x i : centroid , A i : area , i = 1 : N N ⇣ ⌘ h i X 1 1 + ✏ 2 φ i = K 1 ( x i , x j ) ∂ n φ j + K 2 ( x i , x j ) φ j A j + S 1 ( x i ) 2 ✏ 1 j =1 j 6 = i N ⇣ ⌘ h i X 1 1 + ✏ 1 ∂ n φ i = K 3 ( x i , x j ) ∂ n φ j + K 4 ( x i , x j ) φ j A j + S 2 ( x i ) 2 ✏ 2 j =1 j 6 = i ⇒ linear system : Ax = b - GMRES - matrix-vector product at each step - fast summation → treecode

treecode : higher-order version of Barnes-Hut (1986) N q j G ( x i , x j ) , i = 1 , . . . , N : particle-particle , O ( N 2 ) X φ i = j =1 j 6 = i p X X a k ( x i , x c ) M k ( c ) : particle-cluster , O ( N log N ) ≈ c k =0 - Cartesian coordinates - recurrence relations for a k ( x i , x c ) = 1 k ! D k y G ( x, x c ) - geometrically adaptive - low memory usage, good parallel e ffi ciency Li-Johnston-K (2009) JCP 228

example : protein 1A63 RNA binding domain of E. coli rho factor 2069 atoms , triangulation by MSMS with N = 132196

example : protein 1A63 compare 2 codes TABI : N = 20K → 500K APBS : grid = 65 × 33 2 → 513 × 321 2 , h max = 1 . 63 ˚ A → 0 . 13 ˚ A CPU time / error memory / error TABI − PB 4 TABI − PB 10 3 TABI − P TABI − P 10 3% APBS − PB APBS − PB APBS − P APBS − P memory (MB) CPU time (s) 3 10 2 10 2 10 6% 1 10 1 10 0 1 0 1 10 10 10 10 error in E sol (%) error in E sol (%) Geng-K (2013) JCP 247

example : protein 1A63 TABI parallel performance , N = 132196 , error in E sol ≈ 1 . 3% Poisson-Boltzmann # processors run time (s) speedup parallel e ffi ciency (%) 1 799.3 1.00 100.0 2 410.0 1.95 97.5 4 223.8 3.57 89.3 8 123.7 6.46 80.9 Poisson # processors run time (s) speedup parallel e ffi ciency (%) 1 324.4 1.00 100.0 2 166.7 1.95 97.3 4 92.6 3.50 87.6 8 51.0 6.36 79.5

PROTEIN SCIENCE 2018 VOL 27:112-128 TOOLS FOR PROTEIN SCIENCE Improvements to the APBS biomolecular solvation software suite Elizabeth Jurrus, 1 Dave Engel, 1 Keith Star, 1 Kyle Monson, 1 Juan Brandi, 1 Lisa E. Felberg, 2 David H. Brookes, 2 Leighton Wilson, 3 Jiahui Chen, 4 Karina Liles, 1 Minju Chun, 1 Peter Li, 1 David W. Gohara, 5 Todd Dolinsky, 6 Robert Konecny, 7 David R. Koes , 8 Jens Erik Nielsen, 9 Teresa Head-Gordon, 2 Weihua Geng, 4 Robert Krasny, 3 Guo-Wei Wei, 10 Michael J. Holst, 7 J. Andrew McCammon, 7 and Nathan A. Baker 1,11 * 1 Pacific Northwest National Laboratory, Richland, Washington 2 University of California, Berkeley, California 3 University of Michigan, Ann Arbor, Michigan 4 Southern Methodist University, Dallas, Texas 5 St. Louis University, St. Louis, Missouri 6 FoodLogiQ, Durham, North Carolina 7 University of California San Diego, San Diego, California 8 University of Pittsburgh, Pittsburgh, Pennsylvania 9 Protein Engineering, Novozymes A/S, Copenhagen, Denmark 10 Michigan State University, East Lansing, Michigan 11 Brown University, Providence, Rhode Island

comparison of so>ware for molecular surface NanoShaper MSMS N = 29224, iter = 9 N = 29867, iter = 62 E_sol = -10443 kJ/mol E_sol = -10675 kJ/mol 93.9 Surface Potential (kJ/mol/e) 80 60 40 20 0 -20 -40 PDBID: 3dfg -53.8

outline 1. protein-solvent electrostatics 2. incompressible fluid flow - tracer transport on a sphere - vortex dynamics on a sphere - vortex dynamics in 2D

motivation : atmosphere

tracer transport on a sphere tracer transport given : q 0 ( x ) , u ( x, t ) , x 2 S problem : determine q ( x, t ) for t > 0 Eulerian form ∂ q ∂ t ( x, t ) + u · r q ( x, t ) = 0 Lagrangian form flow map : a ! x ( a, t ) , a 2 S ∂ ∂ tx ( a, t ) = u ( x ( a, t ) , t ) , x ( a, 0) = a q ( x ( a, t ) , t ) = q 0 ( a )

Lagrangian particle method particles : x j ( t ) ≈ x ( a j , t ) , x j (0) = a j , j = 1 : M panels : P k ( t ) = x ( P k (0) , t ) , S = ∪ N k =1 P k ( t ) , k = 1 : N icosahedral triangulation N 20 80 320 1280 5120 20480 81920 327680 M 32 122 482 1922 7682 30722 122882 491522 ∆ λ 63 . 43 � 33 . 87 � 17 . 22 � 8 . 64 � 4 . 33 � 2 . 16 � 1 . 08 � 0 . 54 � Bosler-Kent-K-Jablonowski (2017) JCP 340

particle advection d dtx j ( t ) = u ( x j , t ) , x j (0) = a j q j = q 0 ( a j ) problem : particles become disordered for t > 0 solution : remeshing : { x old j } → { x new } , { q new } = ? j j

test case 1 : moving vortices flow - no remeshing movie

particle advection d dtx j ( t ) = u ( x j , t ) , x j (0) = a j q j = q 0 ( a j ) problem : particles become disordered for t > 0 solution : remeshing : { x old j } → { x new } , { q new } = ? j j 1. direct remeshing { q new } = I ( { x new } ; { x old j } , { q old j } ) j j 2. indirect remeshing { a new } = I ( { x new } ; { x old j } , { a old j } ) j j { q new } = q 0 ( { a new } ) j j I : STRIPACK/SSRFPACK (Renka 1997)

test case 1 : moving vorZces flow movie

tracer error ∆ λ = 8 . 64 � test case 1 test case 2 ∞ 2 2 ∞ 1 0 10 10 0 10 − 1 10 − 1 10 − 2 − 2 10 10 LR, l ∞ LR, l ∞ − 3 LR, l 2 LR, l 2 10 − 3 10 − 4 LPM − d, l ∞ LPM − d, l ∞ 10 ∆ λ = 4 . 33 � − 4 10 LPM − d, l 2 LPM − d, l 2 − 5 10 LPM − i, l ∞ LPM − i, l ∞ − 6 − 5 10 10 O( ∆λ 4 ) O( ∆λ 4 ) LPM − i, l 2 LPM − i, l 2 − 7 10 − 6 10 0.5 1 2 4 8 16 32 0.5 1 2 4 8 16 32 ∆λ (degrees) ∆λ (degrees) ∆ λ = 2 . 16 �

test case 1 : moving vorZces flow - adapZve panel refinement ⇣ ( x k ) A k < ✏ Γ max movie

test case 2 : reversing-deformational flow - Gaussian hills tracer movie

test case 2 : reversing-deformational flow Gaussian-hills tracer , grid spacing ∆ λ = 8 . 64 � direct remeshing indirect remeshing

tracer error test case 1 test case 2 ∞ 2 2 ∞ 1 0 10 10 0 10 − 1 10 − 1 10 − 2 − 2 10 10 LR, l ∞ LR, l ∞ − 3 LR, l 2 LR, l 2 10 − 3 10 − 4 LPM − d, l ∞ LPM − d, l ∞ 10 − 4 10 LPM − d, l 2 LPM − d, l 2 − 5 10 LPM − i, l ∞ LPM − i, l ∞ − 6 − 5 10 10 O( ∆λ 4 ) O( ∆λ 4 ) LPM − i, l 2 LPM − i, l 2 − 7 10 − 6 10 0.5 1 2 4 8 16 32 0.5 1 2 4 8 16 32 ∆λ (degrees) ∆λ (degrees)

vortex dynamics on a sphere (partial list!) point vortices . . . … Kimura-Okamoto (1987), Newton-Sakajo (2007) vortex patches . . . Dritschel (1989, 2004), Dritschel-Polvani (1992) … Crowdy-Cloke (2003), Surana-Crowdy (2008) vortex sheets . . . … Sakajo (2009) present work : smooth vorticity

barotropic vorticity equation : Eulerian form velocity : u = r ψ ⇥ x vorticity : ζ = r ⇥ u Poisson equation : r 2 ψ = � ζ ∂ζ vortex dynamics : ∂ t + u · r ( ζ + f ) = 0 Coriolis parameter : f = 2 Ω sin θ Bosler-Wang-Jablonowski-K (2014) FDR 46

barotropic vorticity equation : Lagrangian form flow map : x = x ( ~ ↵ , t ) vorticity : ζ = ζ ( x , t ) ψ ( x ) = − 1 Z stream function : log(1 − x · ˜ x ) ζ (˜ x ) dA (˜ x ) 4 π S dt = − 1 x × ˜ d x Z Biot-Savart law : x x ζ (˜ x , t ) dA (˜ x ) 4 π 1 − x · ˜ S d ζ dt = − 2 Ω dz vortex dynamics : dt

discretization S = ∪ N k =1 A k panels : cubed sphere particles : x j ( t ) = x ( ~ ↵ j , t ) d x j dt = − 1 x j × x k X ζ k A k 4 π 1 − x j · x k k =1 k 6 = j d ζ j dt = − 2 Ω dz j vorticity : dt

Rossby-Haurwitz wave - no remeshing Rossby-Haurwitz wave – no remeshing vorticity panels movie

RH4 wave , remesh + restart Rossby-Haurwitz wave – remeshing vorticity panels movie