Brownian dynamics simulations with hard-body interactions: Exact - PowerPoint PPT Presentation

Brownian dynamics simulations with hard-body interactions: Exact numerical treatment PHYSICAL REVIEW E 83 , 065701(R) (2011) Ralf Eichhorn Hard-wall interactions in soft matter systems: Exact numerical treatment Hans Behringer 1,* and Ralf Eichhorn 2, † 1 University of Mainz, Institute of Physics, Staudinger Weg 7, D-55128 Mainz, Germany 2 NORDITA, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden (Received 29 March 2011; published 20 June 2011) THE JOURNAL OF CHEMICAL PHYSICS 137 , 164108 (2012) Brownian dynamics simulations with hard-body interactions: Spherical particles Hans Behringer 1,a) and Ralf Eichhorn 2,b) 1 Johannes Gutenberg-Universität Mainz, Institut für Physik, Staudinger Weg 7, D-55128 Mainz, Germany 2 Nordic Institute for Theoretical Physics (NORDITA), Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 106 91 Stockholm, Sweden

Brownian motion • interacting particles in a suspension (e.g. colloids) • driven by external forces • in a structured environment • solvent : viscous friction and thermal fluctuations microfluidics, biomolecules in the cell, self-assembly, polymers, ... ֒ → numerical simulation of particle interaction with hard walls

Langevin equation √ r ( t ) = 1 2 D � � ˙ F ( � r ( t ) , t ) + ξ ( t ) � η two common idealizations/approximations: • overdamped limit (“ m = 0”) • hard-body interactions (singular!) to represent the extremely short-ranged and strong repulsive contact forces not included in � F ( � r ( t ) , t )

Euler algorithm √ 2 D d t � d � r t = � v t d t + r ( t + d t ) = � r ( t ) + d � G t , � r t v t = 1 η � r ( t ) , t ), and � G t ∈ N (0 , 1) d with � F ( � � r ( t + d t ) v t d t ] 2 � � 1 − [d � r t − � √ p (d � r t ) = 4 πD d t d exp d � r t 4 D d t � r ( t )

Euler algorithm √ 2 D d t � d � r t = � v t d t + G t , � r ( t + d t ) = � r ( t ) + d � r t v t = 1 η � r ( t ) , t ), and � G t ∈ N (0 , 1) d with � F ( � algorithm: � r ( t ) 1) detect unphysical configurations (“collisions”) d � r t r ( t + d t ) � 2) rule to generate physically valid configuration p (d � r t ) = ???

Euler algorithm √ 2 D d t � d � r t = � v t d t + G t , � r ( t + d t ) = � r ( t ) + d � r t v t = 1 η � r ( t ) , t ), and � G t ∈ N (0 , 1) d with � F ( � algorithm: � r ( t ) 1) detect unphysical configurations (“collisions”) d � r t r ( t + d t ) � 2) rule to generate physically valid configuration p (d � r t ) = ???

Heuristic methods (prominent examples) • rejection scheme discard unphysical configurations (advance time or not???) [B. Cichocki and K. Hinsen, Physica A 166 , 473 (1990)] • event-driven scheme propagate fraction of time step ǫ d t until “collision point” use rejection scheme for remaining time step (1 − ǫ )d t [Y.-G. Tao et al., J. Chem. Phys. 124 , 134906 (2006)] lack thorough justification

Euler algorithm √ 2 D d t � d � r t = � v t d t + r ( t + d t ) = � r ( t ) + d � G t , � r t v t = 1 η � r ( t ) , t ), and � G t ∈ N (0 , 1) d with � F ( � � r ( t + d t ) v t d t ] 2 � � 1 − [d � r t − � √ p (d � r t ) = 4 πD d t d exp d � r t 4 D d t � r ( t )

Euler algorithm √ 2 D d t � d � r t = � v t d t + G t , � r ( t + d t ) = � r ( t ) + d � r t v t = 1 η � r ( t ) , t ), and � G t ∈ N (0 , 1) d with � F ( � algorithm: � r ( t ) 1) detect unphysical configurations (“collisions”) d � r t r ( t + d t ) � 2) rule to generate physically valid configuration p (d � r t ) = ???

Euler algorithm √ 2 D d t � d � r t = � v t d t + G t , � r ( t + d t ) = � r ( t ) + d � r t v t = 1 η � r ( t ) , t ), and � G t ∈ N (0 , 1) d with � F ( � algorithm: � r ( t ) 1) detect unphysical configurations d � r ⊥ (“collisions”) r ( t + d t ) � 2) rule to generate physically valid configuration d � r � p (d � r t ) = ???

Smoluchowski solution [M. V. Smoluchowski, Phys. Z. 17 , 557 (1916)] ∂ 2 ∂ ∂ driven diffusion on a half-line q ∈ [0 , ∞ ) ∂q 2 p − v q ∂tp = D q ∂qp � � ∂ reflecting boundary at q = 0 − ∂qp − v q p = 0 D q q =0 initial position q 0 at time t = 0 solution : p ( q, t ; q 0 ) = p 1 ( q, t ; q 0 ) + p 2 ( q, t ; q 0 ) + p 3 ( q, t ; q 0 ) − ( q − q 0 − v q t ) 2 1 � � with p 1 ( q, t ; q 0 ) = exp � 4 D q t 4 πD q t � � − v q q 0 exp − ( q + q 0 − v q t ) 2 � � D q p 2 ( q, t ; q 0 ) = exp � 4 D q t 4 πD q t � � p 3 ( q, t ; q 0 ) = − v q � v q q � q + q 0 + v q t exp erfc � 2 D q D q 4 D q t

Smoluchowski solution 1.4 1.2 1 p(q,t;q 0 ) 0.8 1 0.6 0.98 0.4 p 0.2 0.96 0 0.1 q 0 0 0.5 1 0 0.5 1 q q v q = − 1 . 0, D q = 1 . 0, q 0 = 0 . 5, t = 0 . 05

Smoluchowski solution 2.5 p p 2 1.05 0.3 p(q,t;q 0 ) q q 0 0.1 0 0.1 1.5 1 0.5 0 1 1 1 0 0.5 0 0.5 0 0.5 q q q q 0 = 0 . 05 q 0 = 0 . 35 q 0 = 0 . 6 v q = 1 . 0, D q = 1 . 0, t = 0 . 05

Smoluchowski solution [M. V. Smoluchowski, Phys. Z. 17 , 557 (1916)] ∂ 2 ∂ ∂ driven diffusion on a half-line q ∈ [0 , ∞ ) ∂q 2 p − v q ∂tp = D q ∂qp � � ∂ reflecting boundary at q = 0 − ∂qp − v q p = 0 D q q =0 initial position q 0 at time t = 0 solution : p ( q, t ; q 0 ) = p 1 ( q, t ; q 0 ) + p 2 ( q, t ; q 0 ) + p 3 ( q, t ; q 0 ) − ( q − q 0 − v q t ) 2 1 � � with p 1 ( q, t ; q 0 ) = exp � 4 D q t 4 πD q t � � − v q q 0 exp − ( q + q 0 − v q t ) 2 � � D q p 2 ( q, t ; q 0 ) = exp � 4 D q t 4 πD q t � � p 3 ( q, t ; q 0 ) = − v q � v q q � q + q 0 + v q t exp erfc � 2 D q D q 4 D q t

The algorithm 1. Perform standard integration as long as the displacements d � r t do not lead to unphysical configurations (“collisions”) 2. If a suggested displacement d � r t results in a “collision”, replace its component along the “collision axis” by a new d q = q − q 0 p 2 ( q, d t ; q 0 ) + p 3 ( q, d t ; q 0 ) drawn from w � ∞ with w = d q [ p 2 ( q, d t ; q 0 ) + p 3 ( q, d t ; q 0 )] 0 � ∞ = 1 − d q p 1 ( q, d t ; q 0 ) 0 � 0 = −∞ d q p 1 ( q, d t ; q 0 )

The algorithm “collision axis” ˆ = half-line (with “collision point” at the origin) random number created on the “collision axis”: Q = Q G Θ( Q G ) + Q C Θ( − Q G ) with p G ( q ) = p 1 ( q ) ( q ∈ R ) and p C ( q ) = p 2 ( q ) + p 3 ( q ) ( q ∈ R > 0 ) w � ∞ � ∞ ⇒ p alg ( q ) = −∞ d q 1 d q 2 p G ( q 1 ) p C ( q 2 ) δ ( q − [ q 1 Θ( q 1 ) + q 2 Θ( − q 1 )]) 0 � ∞ � ∞ d q 2 p G ( q 1 ) p C ( q 2 ) [Θ( q 1 ) δ ( q − q 1 ) = −∞ d q 1 0 + Θ( − q 1 ) δ ( q − q 2 )]) � 0 = Θ( q ) p G ( q ) + p C ( q ) −∞ d q 1 p G ( q 1 ) = p 1 ( q ) + p 2 ( q ) + p 3 ( q )

Euler algorithm for hard wall √ r ∗ 2 D d t � d � r t = � v t d t + r ( t + d t ) = � r ( t ) + d � G t , � t v t = 1 η � r ( t ) , t ), and � G t ∈ N (0 , 1) d with � F ( � algorithm: 1) “suggest” d � r t using standard Euler scheme � r ( t ) 2a) no “collision” : d � r ⊥ r ∗ d � t = d � r t 2b) “collision” : � r ( t + d t ) r ∗ d � t = d � r t + ( q − q 0 − � n · d � r t ) � n d � r � n = − d � r ⊥ / | d � r ⊥ | with � note : d � r ⊥ and d � r � are uncorrelated

Euler algorithm for spherical particles � 2 D 1 d t � � 2 D 2 d t � d � r 1 = � v 1 d t + d � r 2 = � v 2 d t + G 1 , G 2 r ∗ r ∗ r 1 ( t + d t ) = � r 1 ( t ) + d � r 2 ( t + d t ) = � r 2 ( t ) + d � � 1 , � 2 algorithm: 1) “suggest” d � r 1 , d � r 2 using standard Euler scheme 2a) no “collision” : r ∗ r ∗ d � 1 = d � r 1 and d � 2 = d � r 2 d � r 2 d � r 1 2b) “collision” : η 2 r ∗ d � 1 = d � r 1 + [(d � r 2 − d � r 1 ) · � e − ( q − q 0 )] � e η 1 + η 2 η 1 � r 2 ( t ) r ∗ d � 2 = d � r 2 − [(d � r 2 − d � r 1 ) · � e − ( q − q 0 )] � e r 1 ( t ) η 1 + η 2 � with � e = ( � r 2 ( t ) − � r 1 ( t )) / | � r 2 ( t ) − � r 1 ( t ) | note : center of friction and relative motion are uncorrelated

Example: Mean first passage time 20 τ comp [a.u.] 2 15 1 0 τ [s] 0 5 10 15 10 v [ µ m/s] � F 5 0 0 5 10 15 v [ µ m/s] � F = ( − f, f ), v = f/η , particle radius 1 µ m, d t = 0 . 01 s

In practice • detection of “collisions” (any integration scheme) • integration time step d t is determined by variations of � F and curvature of structures/particles • generation of random number q according to distribution p C ( q ) = [ p 2 ( q ) + p 3 ( q )] /w � � � � � � q 0 + v q d t v q q q + q 0 + v q d t √ √ erfc − exp erfc � q D q 4 D q d t 4 D q d t 0 d q ′ p C ( q ′ ) = F ( q ) = . � � q 0 + v q d t √ erfc 4 D q d t Then : q = F − 1 ( x ) with x uniformly distributed on [0 , 1] is distributed according to p C ( q ) numerical solution (Brent’s scheme from the GNU scientific library)