Event Driven Langevin simulations of Hard Spheres

Event Driven Langevin simulations of Hard Spheres

A. Scala ISC-CNR Dipartimento di Fisica, Sapienza Università di Roma Piazzale Moro 5, 00185 Roma, Italy London Institute of Mathematical Sciences, 22 South Audley St Mayfair London W1K 2NY, UK
July 15, 2019

The blossoming of interest in colloids and nano-particles has given renewed impulse to the study of hard-body systems. In particular, hard spheres have become a real test system for theories and experiments. It is therefore necessary to study the complex dynamics of such systems in presence of a solvent; disregarding hydrodynamic interactions, the simplest model is the Langevin equation. Unfortunately, standard algorithms for the numerical integration of the Langevin equation require that interactions are slowly varying during an integration time-step. This in not the case for hard-body systems, where there is no clearcut between the correlation time of the noise and the time-scale of the interactions. Starting first from a splitting of the Fokker-Plank operator associated with the Langevin dynamics, and then from an approximation of the two-body Green’s function, we introduce and test two new algorithms for the simulation of the Langevin dynamics of hard-spheres.


I Introduction

Hard spheres (HS) are a reference system for structural and dynamical theories of fluids Andersen71 (); HansenMcDonaldBook (), but idealized: the infinitely steep potential is essentially a way of capturing the effects of steric interactions. On the atomic or the molecular scale two body interactions are mostly modelled by Lennard-Jones or Coulumb potentials; experiments on colloids shift the length scales of interest up to roughly to where objects can behave as hard bodies and are still small enough to exhibit thermal or Brownian motion in a solvent. Dynamical light scattering BerneBook (); DeGiorgioBook () has already provided a rich collection of data for such systems, encouraging a considerable effort in understanding the dynamics; the possibility of following single particle trajectories via confocal microscopy of latex particle vanBlaaderen95 () has allowed a direct view on an experimental realization of HS systems and their dynamics Kegel00 (); Weeks00 ().

The simplest model of a suspension of neutral particles is to consider a system of HS in an ideal solvent with no hydrodynamic interactions; real suspensions are often described in terms of their deviations from such ideal system. This is the most interesting model for theoreticians and many results have been derived: the two body case (and hence the low density case) has been solved exactly Hanna82 (); Ackerson82 (), while at moderate and high packing fractions various Enskog-like Felderhof83A (); Cichocki90theo () or Mode Coupling theories Hess83 (); FuchsPRL2002 () have been applied to understand the dynamics. While hydrodynamic interactions (HI) are well understood at low particle densities, much less is known at high densities, and theories often proceed by disregarding them Felderhof83B (). As an example, theories regarding glass transition often disregard HI effect, like in the case of the Mode Coupling theory for Brownian hard spheresFuchsPRL2002 (); BraderPRL2007 (); BraderPRL2008 () or Brownian hard discs in shear flowHenrich2009 ().

Non-HI simulations therefore have their place in testing such theories, and in circumventing the huge effort needed to simulate HI Brady88 (); Rjoy05 (); Groot97 (); Tanaka00 ().

In order to validate non-HI theories for HSs it is necessary to use computer simulations, as only a qualitative agreement is to be expected among non-HI theories and data for real suspension. Standard simulation methods for Brownian dynamic like the well-known Ermak-McCammon Ermak78 () require continuous potentials; to circumvent such problem several algorithms have been introduced with various degrees of justification Cichocki90sim (); Heyes93 (); Sillescu94 (); Foss2000 () for the over-damped dynamics; only recently it has been recognized that in the case of hard interactions such simulations are better performed by event-driven (ED) codes Strating99 (); Tao06 (); BD4HS (). We introduce two new ED algorithms that go beyond the over-damped approximation and allow for the simulation of the full Brownian dynamics of HSs.

Ii Methods

We consider a system of HSs governed by the Langevin equation


for the positions and the velocities ; here is the friction constant, the acceleration, is the mass of the HSs, is the potential energy and are the zero-mean random forces due to the solvent. We assume that such random forces are delta correlated and satisfy the fluctuation-dissipation theorem


In the case of continuous interactions, it is possible to define stochastic Taylor expansions KannanBook (); correspondingly, integration schemes of the -th order with errors of order in the time-step can be introduced KloedenPlatenBook (). In the case of hard-body interactions, all the standard machinery of stochastic calculus breaks down due to the singular nature of the interaction potential and new methods must be developed.

We consider the Fokker-Plank equation associated to the SDE (1) (Kramers’ equation Kramers1940 ())


where is the probability distribution function (PDF) for the positions and the velocities of the particles, relates to the temperature and


is the Kramer operator. Integrating the SDE (1) for a finite time-step corresponds to extracting a configuration according to the probability .

Iii Splitted Brownian Dynamics

To obtain a numerical approximation, a powerful approach is to split the evolution operator in a product of exactly-integrable operators Forbert00 () ensuring that that the decomposition is positive (i.e. all ) Chin05 (). Therefore, to each splitting corresponds an algorithm in which in a single time-step the operators are applied in sequence.

We first choose to split into the reversible (or streaming) operator and the irreversible (or collision) operator RiskenBook (); we indicate the corresponding algorithm as Splitted Brownian Dynamics (SBD).

The operator is the Liouvillian associated to the Hamiltonian . In the case of step potentials, the associated reversible equation of motion can be integrated via event-driven molecular dynamics (EDMD) RapaBook () with a precision limited only by the numerical round-off errors; therefore the propagator can be implemented with extreme accuracy.

The operator corresponds to the interaction with the bath; the associated SDE is can be exactly integrated giving an explicit formula for the evolution :


where is a unitary Gaussian random variable and .

The algorithm for the single SBD time-step consists therefore in an EDMD simulation RapaBook () of length followed by a thermalization of the velocities according to eq.(5). We notice that the error is at most quadratic (as can be checked via Taylor expansion ) and regards only in the dynamics; in fact, SBD is equivalent (upon identifying the angle mixing reversible and irreversible evolution with ) to the Generalized Hybrid Monte Carlo Kennedy99 () and therefore explores the canonical ensemble as long as the propagation steps , can be exactly implemented (as in our case).

It is therefore of interest to give some physical bounds on the magnitude of the feasible time-step . First, we notice that for the dynamics reduces to and MD simulations where velocities are extracted each from a Maxwellian; therefore if the time-step is much bigger than the average inter-particle collision time, results of classical MD are to be expected. Accordingly, we find that for big the algorithm overestimates the diffusion coefficient (fig. 1 ); this is to be expected as the mean free path (in absence of collisions) of a particle is of order instead of . Second, the magnitude of is naturally bounded the damping time ; therefore the SBD is not well indicated for simulations in the over-damped limit . Accordingly, we find that SBD overestimates diffusion coefficients for (fig.1); it is therefore necessary to develop an alternative approach for the simulation of systems with high damping.

Iv Approximate Green’s function dynamics

It has been shown in BD4HS () that the over-damped limit of eq.(1) can be simulated efficiently using ED codesBD4HS (). The algorithm relies on considering time steps small enough so that mostly binary collisions are relevant, i.e. the average displacement should be less than the average inter-particle separation. Moreover, average displacement should be smaller than the HSs’ radii in order to map the interaction of two nearby HSs in the problem of a random walk near a reflective wall. Under such approximations, the true two-body stochastic dynamics for over-damped Brownian HSs can be implemented by algorithm of BD4HS () in which each step consists in predicting the displacements of the HSs via the free propagator, introducing fictive velocities , and performing an EDMD with such fictive velocities during and . We extend such approach to the general Brownian case.

First, we need to predict the positions of the HSs after a time-step according to their free propagation, i.e. the solution of eq.(1) with no interaction (:


The particle displacements contain both systematic parts , and stochastic displacements. The stochastic displacements , are zero-mean correlated gaussian variables with variances , and cross-correlation ATBook ().

If we consider a time-step such that the average displacement is less than the average inter-particle separation, we can consider only the corrections due to two-body interactions. In the limit of small , a couple of HSs will interact only when they start from nearby positions. In particular, if , i.e. the average free displacement is much smaller than the diameter of the HSs, the dynamics of two particles and can be approximated as the Langevin dynamics of a point particle at a distance from a flat wall. It is possible to solve such problem with a straightforward generalization of the image method applied in BD4HS (). In fact, the solution given by the free particle Green’s function plus an image particle with a reflected velocity beyond the reflective wall (fig. 2) correctly satisfies the zero-current boundary condition , where is the normal to the wall and is the probability current for the position.

Such solution can be implemented exactly by predicting the new positions and velocities , according to eq.(6), defining fictive velocities and performing an EDMD simulation with such fictive velocities during ; if a collision happens, the component of the relative velocity normal to the contact point must be reflected for both the fictive and the predicted velocities . We indicate such algorithm as the approximate Green’s function dynamics (AGD). In the over-damped limit, the prediction of the velocities and positions decorrelates and the algorithm correctly reduces to the over-damped case of BD4HS ().

As for the SBD algorithm, it can be proven that the AGD scheme respects detailed balance and ergodicity and therefore explores the correct ensemble for HSs; hence, errors are again only in dynamic quantities. At difference with SBD, we have no analytic estimate for the error; nevertheless, we expect that the the mean-free path in absence of collisions must be smaller than the radius of the HSs in order to satisfy the flat-wall approximation, and must be smaller than the average inter-particle distance in order to avoid multiple collisions (hence higher than two-body effects) during .

In order to check that the behaviour of AGS is driven just by geometrical considerations, we have simulated HS systems at different and varying the time-step in the range (reduced units). At difference with SBD where diffusion can vary even by a order of magnitude in such a range, the values of measured from AGD vary a few percent over the range and long simulations are been necessary to have enough statistics to detect the behaviour of that would otherwise look flat. In fig. 3 we show that the measured diffusion coefficient versus the AGS simulation time-step displays a plateau (i.e. fluctuations become much smaller than ) already for regardless of and .

V Conclusions

Hard spheres, and in general hard body systems in suspension, have become a realistic model due to the developments of experimental techniques for the investigation of colloidal systems and nano-particles; yet the dynamics of such systems is hard to simulate via the standard Brownian dynamics algorithms. In fact, classical continuous-time algorithms fail due to instantaneous character of the interactions; we have shown instead how it is possible to simulate the full Langevin dynamics of Hard Spheres.

First, we have shown how the simplest splitting of the stochastic evolution operator (a technique often referred to as ”Trotterization” from Trotter’s seminal workTrotterPAMS59 ()) allows to write an algorithm (the Splitted Brownian Dynamics - SBD). The SBD algorithm becomes inefficient of high viscosities but via the operator-splitting technique could easily take account for the interaction with external fields or with the presence of fluxes (like shear) in the surrounding fluid.

Second, we have shown how by considering the two body dynamics of Brownian Hard Spheres it is possible do develop an algorithm (the Approximate Green’s function Dynamics AGD) that overcomes such problem and works equally well for a wide range of packing fractions and viscosities. To develop the AGD algorithm, we have solved the problem of the Langevin dynamics of a point particle in presence of a reflective wall by extending the classical Image Method solution for the over-damped Brownian dynamics of a point particle in presence of a reflective wall (here , are noises). The AGD algorithm is Event Driven and considers fictive collisions between Hard Spheres. While it should possible to take into account the polydispersity of a system by considering also effective masses in the fictive collisions as in BD4HS (), including shear or external fields in the AGD algorithm looks more complicated as it would require the solution of the particle - reflective wall problem with external fields/shear.

Both SBD and AGD simulations explore the canonical ensemble for Hard Spheres and therefore reproduce the correct equilibrium thermodynamics. They belong to the class of Asynchronous Event-Driven Particle AlgorithmsDonevSIM2009 () and can be easily implemented by adapting existing codes for ED dynamics RapaBook () or Brownian Dynamics HardBrown () of Hard Spheres.

The author thanks Th. Voigtmann for long and useful discussions.


  • (1) H. C. Andersen, J. D. Weeks, and D. Chandler, Physical Review A 4, 1597 (Oct. 1971)
  • (2) J. P. Hansen and I. R. McDonald, Theory of Simple Liquid, 2nd ed. (Academic Press, New York, 1989)
  • (3) B. J. Berne and R. Pecora, Dynamic Light Scattering: with Applications to Chemistry Biology, and Physics (John Wiley & Sons, New York, 1976)
  • (4) V. De Giorgio, M. Corti, and M. Giglio, Light Scattering in Liquids and Macromolecular Solutions (Plenum, New York, 1980)
  • (5) A. van Blaaderen and P. Wiltzius, Science 270, 1177 (1995)
  • (6) W. K. Kegel and A. van Blaaderen, Science 287, 290 (2000)
  • (7) E. R. Weeks, J. C. Crocker, A. Levitt, A.C.and Schofield, and D. A. Weitz, Science 287, 627 (2000)
  • (8) S. Hanna, W. Hess, and R. Klein, Physica A 111, 181 (Mar. 1982)
  • (9) B. J. Ackerson and L. Fleishman, Journal of Chemical Physics 76, 2675 (Mar. 1982)
  • (10) B. Felderhof and R. Jones, Physica A Statistical Mechanics and its Applications 121, 329 (Aug. 1983)
  • (11) B. Cichocki and B. U. Felderhof, Physical Review A 42, 6024 (Nov. 1990)
  • (12) W. Hess and R. Klein, Adv. Phys. 32, 173 (1983)
  • (13) M. Fuchs and M. E. Cates, Phys. Rev. Lett. 89, 248304 (Nov 2002)
  • (14) B. Felderhof and R. Jones, Faraday Discuss. Chem. Soc. 76, 179 (1983)
  • (15) J. M. Brader, T. Voigtmann, M. E. Cates, and M. Fuchs, Phys. Rev. Lett. 98, 058301 (Jan 2007)
  • (16) J. M. Brader, M. E. Cates, and M. Fuchs, Phys. Rev. Lett. 101, 138301 (Sep 2008)
  • (17) O. Henrich, F. Weysser, M. E. Cates, and M. Fuchs, Physical and Engineering Sciences 367, 5033 (Dec. 2009)
  • (18) J. F. Brady and G. Bossis, Annual Review of Fluid Mechanics 20, 111 (1988)
  • (19) R. Adhikari, K. Stratford, M. E. Cates, and A. J. Wagner, Europhysics Letters 71, 473 (2005)
  • (20) R. D. Groot and P. B. Warren, Journal of Chemical Physics 107, 4423 (1997)
  • (21) H. Tanaka and T. Araki, Phys. Rev. Lett. 85, 1338 (2000)
  • (22) D. L. Ermak and J. A. McCammon, Journal of Chemical Physics 69, 1352 (Aug. 1978)
  • (23) B. Cichocki and H. K., Physica A 166, 473 (Jul. 1990)
  • (24) D. M. Heyes and J. R. Melrose, Journal of Non-Newtonian Fluid Mechanics 46, 1 (Jan. 1993)
  • (25) W. Schaertl and H. Sillescu, Journal of Statistical Physics 74, 687 (Feb. 1994)
  • (26) D. R. Foss and J. F. Brady, Journal of Fluid Mechanics 407, 167 (Mar. 2000)
  • (27) P. Strating, Physical Review E 59, 2175 (Feb. 1999)
  • (28) Y.-G. Tao, W. K. den Otter, J. K. G. Dhont, and W. J. Briels, Journal of Chemical Physics 124, 134906 (2006)
  • (29) A. Scala, C. De Michele, and T. Voigtmann, Journal of Chemical Physics 126, 134109 (Apr. 2007)
  • (30) D. Kannan and V. Lakshmikantham, Handbook of stochastic analysis and applications (Marcel Dekker, 2002) Chap. 5
  • (31) P. E. Kloeden and E. Platen, Numerical solution of stochastic differential equations, 3rd ed., Applications of Mathematics, Vol. 23 (Springer, 1999)
  • (32) H. Kramers, Physica 7, 284 (1940), ISSN 0031-8914
  • (33) H. A. Forbert and S. A. Chin, Physical Review E 63, 016703 (Dec. 2000)
  • (34) S. A. Chin, Physical Review E 71, 016703 (Jan. 2005)
  • (35) H. Risken, The Fokker-Planck equation. Methods of solution and applications (Springer Series in Synergetics, Berlin, New York: Springer, —c1989, 2nd ed., 1989)
  • (36) D. C. Rapaport, The Art of Molecular Dynamics Simulation (Cambridge University Press, 2004)
  • (37) A. D. Kennedy, Parallel Computing 7, 284 (Apr. 1999)
  • (38) M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, 2nd ed. (Clarendon Press, Oxford, 1987)
  • (39) H. F. Trotter, Proc. Amer. Math. Soc. 10, 545 (1959)
  • (40) A. Donev, SIMULATION 85, 229 (2009)
  • (41) A. Scala, “Simulation of hard brownian and granular particles,” (2008), http://gna.org/projects/hardbrown
Figure 1: Effect of the damping coefficient on the size of the simulation step (all quantities in reduced units). The diffusion coefficient from simulations is plotted versus the time-step size for various ’s. As expected, the system approaches the value for diffusion regardless of for . The “true” value of is obtained for . We observe at small ’s a plateau in the  vs  plot for , signalling that the “true” value of is approached. Results are presented for packing fraction ; a completely analogous behaviour is found at a low packing fraction and an high packing fraction .
Figure 2: A two body problem for hard spheres can be mapped into the problem of a point particle interacting with a larger sphere. When particles are very near, the problem further simplifies to the interaction of a Langevin particle with a reflective flat wall, whose solution can be derived by applying the Image Method to the Langevin equation. In fact, the Green function must zero inside the wall and must satisfy the no-flux boundary conditions at the wall. Combining the free Green function of the particle in its initial position and the free Green function of its image (with the normal-to-the-wall component of the initial velocity reflected) satisfies both Kramer’ equations and reflective boundary conditions giving the correct solution.
Figure 3: Effects of packing fraction (upper panel) and of the damping coefficient (lower panel) on the time-step for the AGF algorithm. All quantities in reduced units; thick lines are just a guide for the eye. Diffusion is calculated averaging over independent trajectories for particle systems; simulations are long at least times the structural correlation time. In the upper panel, results are shown for at fixed damping . In the lower panel, results are shown for at fixed packing fraction . Notice that the estimated diffusion coefficient has a small relative variation in the wide range of dampings s and packing fractions s analysed. As a rule of thumb, to estimate with an accuracy much smaller than time-step of order are already enough.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description