Adaptive Brownian Dynamics

Computer simulations have long become an established tool in the investigation of physical phenomena [1, 2, 3]. Complementing experimental results, they build the foundation for the exploration of increasingly complex dynamical systems. From the standpoint of classical statistical mechanics, the simulation of a many-body system consisting of discrete interacting particles can reveal information about its structural correlation as well as its thermodynamic properties. Naturally, this opens up the possibility of tackling many problems in the fields of material science, soft matter and biophysics, such as investigating the dynamics of macromolecules [4], predicting rheological properties of fluids [5, 6, 7], or exploring non-equilibrium processes which occur e.g. in colloidal suspensions under the influence of external forcing [8].

With the ever-increasing capabilities of computer hardware, a variety of different computational methods have emerged since the middle of the last century. Conceptually, at least three distinct classes of particle-based simulation frameworks can be identified: i) Monte-Carlo (MC), which relies on the stochastic exploration of phase space, ii) Molecular Dynamics (MD), in which the set of ordinary differential equations (ODEs) of Hamiltonian dynamics is integrated to obtain particle trajectories, and iii) Langevin Dynamics, where random processes are incorporated into the Newtonian equations of motion so that the evolution of a system is obtained by numerical integration of then stochastic differential equations (SDEs). Brownian Dynamics (BD) can be seen as a special case of iii), since the underlying stochastic Langevin equation is thereby considered in the overdamped limit where particle inertia vanishes and only particle coordinates remain as the sole microscopic degrees of freedom.

Notably, a broad range of refined methods have been developed in all three categories, sometimes even intersecting those. Important examples of such extensions are kinetic Monte-Carlo for the approximation of time evolution from reaction rates [9], the addition of thermostats in MD to model thermodynamic coupling [10, 11], event-driven algorithms which enable both MD and BD in hard-particle systems [12, 13], and the adaptation of molecular algorithms to modern hardware [14]. Improvements in the calculation of observables from the resulting particle configurations have been made as well, e.g. by modifying their generation in MC (umbrella sampling, transition matrix MC [15], Wang-Landau sampling [16]) or by utilizing advanced evaluation schemes in MD and BD, such as force sampling [17, 18, 19] or adaptive resolution MD [20].

The efficiency and accuracy of a certain algorithm are always primary concerns, as these properties are essential for applicability and practicability in real-world problems. One therefore aims to design procedures that are both as fast and as precise as possible – yet it is no surprise that those two goals might often be conflicting. Especially in BD, where stochastic processes complicate the numerical treatment, the development of more sophisticated algorithms apparently lacks behind that of MD, for example, and one often resorts to alternative or combined simulation techniques [21, 22]. If the full dynamical description of BD is indeed considered, the equations of motion are usually integrated with the simple Euler-Maruyama method [23], where stochasticity is accounted for in each equidistant step via normally distributed random numbers. This can lead to inaccuracies and stability problems, making BD seem inferior to other computational methods.

In this work, we propose a novel approach to BD simulations which rectifies the above shortcomings of conventional BD. To achieve this, we employ an adaptive timestepping algorithm that enables the control of the numerical error as follows. The common Euler-Maruyama method is complemented with a higher-order Heun step to obtain an embedded integrator pair for an estimation of the local discretisation error per trial step. By comparison of this error estimate with a prescribed tolerance, the trial step is either accepted or it is rejected and then retried with a smaller stepsize. Particular care is required after rejections so as to not introduce a bias in the random forces which would violate their desired properties. We therefore use Rejection Sampling with Memory (RSwM) [24] to retain a Gaussian random process even in a scenario where already determined random increments may conditionally be discarded. RSwM is a recently developed algorithm for the adaptive generation of random processes in the numerical solution of SDEs, which we improve and specialize to our context of overdamped Brownian motion and thereby formulate a method for adaptive BD simulations.

We demonstrate the practical advantages of adaptive BD over common BD in simulation results for prototypical bulk equilibrium systems and for more involved cases in non-equilibrium. A notable example which we investigate is the drying of colloidal films at planar surfaces. Especially when dealing with non-trivial mixtures, as e.g. present in common paints and coatings, the dynamics of this process can be inherently complex and its quantitative description turns out to be a major challenge [25, 26]. This stands in contrast to the necessity of understanding and predicting stratification processes in those systems. Stratification leads to a dried film that has multiple layers differing in the concentration of constituent particle species, thereby influencing macroscopic properties of the resulting colloidal film. Therefore, controlling this process is an important measure to tailor colloidal suspensions to their field of application. Advances in this area have been made experimentally [27, 28], by utilizing functional many-body frameworks like dynamical density functional theory (DDFT) [29], and with molecular simulations such as conventional BD [30]. By employing the adaptive BD method, we are able to capture the complex dynamical processes occuring in those systems even in the final dense state. Close particle collisions and jammed states are resolved with the required adjustment of the timestep, necessary for the stability and accuracy of the simulation run in those circumstances. This can not be achieved easily with common BD.

This paper is structured as follows. In Sec. II, a brief and non-rigorous mathematical introduction to the numerical solution of SDEs is given. Particularly, we illustrate the prerequisites for adaptive and higher-order algorithms in the case of general SDEs and emphasize certain pitfalls. In Sec. III, these considerations are applied to the case of Brownian motion. We construct the embedded Heun-Euler integration scheme in Sec. III.1 and incorporate RSwM in Sec. III.2, which yields the adaptive BD method. Observables can then be sampled from the resulting particle trajectories with the means illustrated in Sec. III.3. In Sec. IV, simulation results of the above mentioned Lennard-Jones systems are shown and the practical use of adaptive BD is confirmed. In Sec. V, we conclude with a summary of the shown concepts, propose possible improvements for the adaptation of timesteps, and present further ideas for use cases.

Brownian Dynamics of a classical many-body system of $N$ particles in $d$ spatial dimensions with positions $\mathbf{r}^{N}(t)=(\mathbf{r}^{(1)}(t),\dots,\mathbf{r}^{(N)}(t))$ at time $t$ and temperature $T$ is described by the overdamped Langevin equation. The trajectory of particle $i$ satisfies

where $\mathbf{F}^{(i)}(\mathbf{r}^{N}(t))$ is the total force (composed of external and interparticle contributions) acting on particle $i$, $\gamma^{(i)}$ is the friction coefficient of particle $i$ and $k_{B}$ is Boltzmann’s constant; the dot denotes a time derivative. In eq. (1), the right-hand side consists of a deterministic (first summand) and a random contribution (second summand). The random forces are modeled via multivariate Gaussian white noise processes $\mathbf{R}^{(i)}(t)$ which satisfy

where $\langle\cdot\rangle$ denotes an average over realizations of the random process, $\mathbf{I}$ is the $d\times d$ unit matrix, $\delta_{ij}$ denotes the Kronecker delta, and $\delta(\cdot)$ is the Dirac delta function.

One can recognize that eq. (1) has the typical form of an SDE

if the dependent random variable $X$ is identified with the particle positions $\mathbf{r}^{N}$ and $W$ is a Wiener process corresponding to the integral of the Gaussian processes $\mathbf{R}^{N}=(\mathbf{R}^{(1)},\dots,\mathbf{R}^{(N)})$. As we do not consider hydrodynamic interactions, the random forces in eq. (1) are obtained by a mere scaling of $\mathbf{R}^{(i)}(t)$ with the constant prefactors $\sqrt{2k_{B}T/\gamma^{(i)}}$. Therefore, the noise in BD is additive, since $g(X(t),t)=\mathrm{const.}$ in the sense of eq. (4). This is a crucial property for the construction of a simple higher-order integrator below in Sec. III.

In computer simulations, particle trajectories are obtained from eq. (1) by numerical integration. Contrary to the numerics of ODEs, where higher-order schemes and adaptivity are textbook material, the derivation of corresponding methods for SDEs poses several difficulties which we address below. Due to the complications, SDEs of type (4) are often integrated via the Euler-Maruyama method instead, which follows the notion of the explicit Euler scheme for ODEs. Thereby, the true solution of eq. (4) with initial value $X(0)=X_{0}$ is approximated in $t\in[0,T]$ by partitioning the time interval into $n$ equidistant subintervals of length $\Delta t=T/n$. Then, for $0\leq k<n$, a timestep is defined by

with Wiener increments $\Delta W_{k}$. An Euler-Maruyama step is also incorporated in the adaptive BD method which we construct below, applying eq. (5) to the overdamped Langevin equation (1).

Crucially, the random increments $\Delta W_{k}$ in each Euler-Maruyama step (5) have to be constructed from independent and identically distributed normal random variables with expectation value $\mathop{}\!\mathrm{E}(\Delta W_{k})=0$ and variance $\mathop{}\!\mathrm{Var}(\Delta W_{k})=\Delta t$. In practice, this is realized by drawing a new random number $R$ (or vector thereof) from a pseudo-random number generator obeying the normal distribution $\mathcal{N}(0,\Delta t)=\sqrt{\Delta t}\mathcal{N}(0,1)$ in each step $k$ and setting $\Delta W_{k}=R$ in eq. (5). The process of obtaining such a scalar (or vectorial) random increment $R$ will be denoted in the following by

where $\mathcal{N}(\mu,\eta)$ is a scalar (or multivariate) normal distribution with expectation value $\mu$ and variance $\eta$.

As in the case of ODEs, an important measure for the quality of an integration method is its order of convergence. However, what convergence exactly means in the case of SDEs must be carefully reconsidered due to their stochasticity. We refer the reader to the pertinent literature (see e.g. Ref. 23) and only summarize the key concepts and main results in the following.

Since both the approximation $X_{k}$ and the true solution $X(t_{k})$ are random variables, one can define two distinct convergence criteria. For a certain method with discretisation $\Delta t\rightarrow 0$, weak convergence captures the error of average values whereas strong convergence assesses the error of following a specific realization of a random process. One can show that the Euler-Maruyama method has a strong convergence order of 0.5, i.e. when increasing or decreasing the stepsize $\Delta t$, the error of the numerical solution only scales with $\Delta t^{0.5}$. For general $g(X(t),t)$ in eq. (4), the construction of schemes with higher strong order is complicated due to the occurence of higher-order stochastic integrals. Practically, this means that the careful evaluation of additional random variables is necessary in each iteration step. These random variables then enable the approximation of the stochastic integrals. There exist schemes of Runge-Kutta type with strong orders up to 2, although only strong order 1 and 1.5 Runge-Kutta methods are mostly used due to practical concerns [31, 32].

In order to incorporate adaptivity, one needs a means of comparison of two integration schemes of different strong order to formulate a local error estimate for the proposed step. If the error that occurs in a specific timestep is too large (we define below precisely what we mean by that), the step is rejected and a retrial with a smaller value of $\Delta t$ is performed. Otherwise, the step is accepted and, based on the previously calculated error, a new optimized stepsize is chosen for the next step, which hence can be larger than the previous one. This protocol makes an optimal and automatic control of the stepsize possible, meaning that $\Delta t$ can both be reduced when the propagation of the SDE is difficult and it can be increased when the error estimate allows to do so. Similar to the case of ODEs, it is computationally advantageous to construct so-called embedded integration schemes, analogous to e.g. Runge-Kutta-Fehlberg integrators [33], which minimize the number of costly evaluations of the right-hand side of eq. (4). Developments have been made in this direction e.g. with embedded stochastic Runge-Kutta methods [34].

There is still one caveat to consider when rejecting a step, in that one has to be careful to preserve the properties of the Wiener process. In the naive approach of simply redrawing new uncorrelated random increments, rejections would alter the random process implicitly. The reason lies in the introduction of an undesired bias. Since large random increments (generally causing larger errors) get rejected more often, the variance of the Wiener process would be systematically decreased, ultimately violating its desired properties. To avoid this effect, it must be guaranteed that once a random increment is chosen, it will not be discarded until the time interval it originally applied to has passed. Consequently, when rejecting a trial step and retrying the numerical propagation of the SDE with smaller time intervals, new random increments cannot be drawn independently for those substeps anymore. The new random increments must instead be created based on the rejected original timestep such that an unbiased Brownian path is still followed.

The formal framework to the above procedure is given by the so-called Brownian bridge theorem [35], which interpolates a Wiener process between known values at two timepoints. If $W(0)=0$ and $W(\Delta t)=R$ are given (e.g. due to the previous rejection of a timestep of length $\Delta t$ where the random increment $R$ has been drawn), then a new intermediate random value must be constructed by

such that the statistical properties of the to-be-simulated Wiener process are left intact and a substep $q\Delta t$ can be tried. The value of $q$ thereby sets the fraction of the original time interval to which the Wiener process shall be interpolated. Equation (7) extends naturally (i.e. component-wise) to the multivariate case and it hence enables the construction of interpolating random vectors in a straightforward manner.

With this idea in mind, several adaptive timestepping algorithms for SDEs have been designed [36, 37, 38, 39, 40, 41, 42]. Still, most of these approaches are quite restrictive in the choice of timesteps (e.g. only allowing discrete variations such as halving or doubling) [36, 42], involve the deterministic or random part only separately into an a-priori error estimation [39, 41, 43], or store rejected timesteps in a costly manner, e.g. in the form of Brownian trees [36]. In particular, the above methods rely on precomputed Brownian paths and do not illustrate an ad hoc generation, which is desirable from a performance and memory consumption standpoint in a high-dimensional setting like BD.

In contrast, a very flexible and performant class of adaptive timestepping algorithms called Rejection Sampling with Memory (RSwM) has recently been developed by Rackauckas and Nie [24]. Their work provides the arguably optimal means for the adaptive numerical solution of SDEs while still being computationally efficient in the generation of random increments as well as in the handling of rejections. We therefore use RSwM in the construction of an adaptive algorithm and specialize the method to Brownian motion in the following.

Based on the remarks of Sec. II, we next proceed to apply the general framework to the case of BD with the overdamped Langevin equation (1) forming the underlying SDE. An embedded integration scheme is constructed which allows the derivation of an error estimate and an acceptance criterion in each step. Furthermore, the application of RSwM for handling rejected timesteps in BD is shown and discussed. We also illustrate how the calculation of observables from sampling of phase space functions has to be altered in a variable timestep scenario.

To test and illustrate the adaptive BD algorithm, we investigate the truncated and shifted Lennard-Jones fluid with interaction potential

where

and $r$ is the interparticle distance. We set a cutoff radius of $r_{c}=2.5\sigma$ throughout the next sections and use reduced Lennard-Jones units which yield the reduced timescale $\tau=\sigma^{2}\gamma/\epsilon$.

A common problem in conventional BD simulations is the choice of an appropriate timestep. Obviously, a too small value of $\Delta t$ – while leading to accurate trajectories – has a strong performance impact, hindering runs which would reveal long time behavior and prohibiting extensive sampling periods, which are desirable from the viewpoint of the time average (20). Still, $\Delta t$ must be kept below a certain threshold above which results might be wrong or the simulation becomes unstable. Unfortunately, due to the absence of any intrinsic error estimates, judgement of a chosen $\Delta t$ is generally not straightforward. For instance, one can merely observe the stability of a single simulation run and accept $\Delta t$ if sensible output is produced and certain properties of the system (such as its energy) are well-behaved. Another possibility is the costly conduction of several simulation runs with differing timesteps, thereby cross-validating gathered results. Consequently, a true a priori choice of the timestep is not possible in general and test runs cannot always be avoided.

With adaptive timestepping, this problem is entirely prevented as one does not need to make a conscious choice for $\Delta t$ at all. Instead, the maximum local error of a step is restricted by the user-defined tolerance (11), ensuring correctness of results up to a controllable discretisation error. This does come at the moderate cost of overhead due to the additional operations per step necessary in the embedded Heun-Euler method and the RSwM algorithm. However, as we demonstrate in the following, the benefits of this method far outweigh the cost even in simple situations.

In this work, we have constructed a novel method for BD simulations by employing recently developed algorithms for the adaptive numerical solution of SDEs to the case of Brownian motion as described on the level of the overdamped Langevin equation (1). For the evaluation of a local error estimate in each trial step, we have complemented the simple Euler-Maruyama scheme (8) found in common BD with a higher-order Heun step (9). By comparison of their discrepancy with a user-defined tolerance (11) composed of an absolute and a relative contribution, we were able to impose a criterion (16) for the acceptance or rejection of the trial step and for the adaptation of $\Delta t$. Special care was thereby required in the reduction from particle-wise errors (10) to a global scalar error estimate (12).

Due to the stochastic nature of Brownian motion, the rejection and subsequent retrial of a timestep could not be done naively by redrawing uncorrelated random vectors in this case. Instead, a sophisticated method had to be employed to ensure the validity of the statistical properties of the random process in eq. (1) and hence to avoid biased random forces from a physical point of view. We have illustrated that the Brownian bridge theorem (7) resolves this problem formally, and that RSwM [24] provides a feasible implementation method based on this theorem. Hence, we specialized RSwM to the case of Brownian motion in Sec. III.2 and constructed the fully adaptive BD scheme outlined in alg. 2. A correction to the original algorithm of Ref. 24 is given in Appendix B and a generalization of adaptive BD to non-overdamped Langevin dynamics is outlined in Appendix A.

To test our framework, we applied the adaptive BD method to both equilibrium and non-equilibrium Lennard-Jones systems focusing on the analysis of individual trajectories. Even in the standard case of a phase-separating bulk fluid, we could verify that the use of RSwM induced no significant computational overhead and that a performance gain could be achieved compared to the fixed-timestep Euler-Maruyama method. This is complemented by practical convenience, since $\Delta t$ needs not be chosen a priori.

The real advantages of adaptive BD become clear in more demanding situations where a fixed timestep would ultimately lead to inaccuracy and instability of the simulation without manual intervention. We have shown this with non-equilibrium systems of drying colloidal suspensions and have modeled the rapid evaporation of the implicit solvent by a moving interface as in Ref. 29. With our method, we achieved efficient timestepping and unconditional stability of the simulation even in the final stages where the packing fraction increases rapidly and the final structural configuration of the dried film develops. Particularly, in a single-species Lennard-Jones system, the freezing transition could be captured effortlessly even if crystal defects remained as a result of strong external driving. In the more elaborate case of a binary mixture of different-size particles, the intricate structure of the colloidal film could be resolved accurately. Here, we reported the development of a dual layer of the big particle species both at the substrate as well as at the top interface, while a concentration gradient of the small particle species could be observed in-between. This shows that the method can be used to predict the arrangement of stratified films efficiently and with great detail, which could be helpful in the determination of their macroscopic properties.

We expect similar benefits in all sorts of situations where fundamental changes within the system call for more accurate or more relaxed numerical timestepping. For instance, this could be the case in sheared systems [59, 60], for sedimenting [61, 62] or cluster-forming [63] colloidal suspensions, in microrheological simulations [64], for phase transitions and especially the glass transition [65], in the formation of crystals and quasicrystals [66, 67], in active matter [68, 69, 70, 71], possibly including active crystallization [72, 73], as well as for responsive colloids [74]. In any respect, adaptive BD can help in a more accurate measurement of observables, since the evaluation of the average (20) is not hampered by erroneous particle configurations, which might occur in conventional BD. Especially for quantities with a one-sided bias such as absolute or square values of particle forces and velocities, the abundance of outliers in the set of samples used in eq. (20) is essential to yield valid results.

We finally illustrate possible improvements of the adaptive BD method itself, where several aspects come to mind. Firstly, the choices of the parameters $\alpha$ in eq. (16) as well as $q_{\mathrm{min}}$ and $q_{\mathrm{max}}$ in eq. (17) were made mostly heuristically. We note that in Ref. 24, the influence of those parameters on the performance of the algorithm has been investigated, but emphasize that results could vary for our high-dimensional setting of Brownian motion.

Secondly, recently developed adaptation methods for the stepsize which employ control theory could be helpful. Thereby, $\Delta t$ is not merely scaled after each trial step by a locally determined factor $q$ as evaluated in eq. (16). Instead, one constructs a proportional-integral-derivative (PID) controller for the selection of new timesteps. This means that the adaptation is now not only influenced proportionally by the error estimate of the current step, but rather it is also determined by integral and differential contributions, i.e. the memory of previous stepsizes and the local change of the stepsize. Adaptive timestepping with control theory has already been applied to ODEs [75, 76, 77] and SDEs [40]. This method could help in our case to reduce the number of unneccessarily small steps even further. In our present implementation, a conservatively chosen $q_{\mathrm{max}}$ in eq. (17) keeps the number of rejections low because it restricts the growth of $\Delta t$ after moves with coincidentally low error. But it also has the effect of preventing a fast relaxation of $\Delta t$ after a sudden drop, e.g. due to an unfortunate random event. Then, many steps are needed for $\Delta t$ to grow back to its nominal value, because the maximum gain in each step is limited by $q_{\mathrm{max}}$. With a control theory approach, a mechanism could be established which permits sudden drops to low values of $\Delta t$ but still ensures a rapid relaxation afterwards.

Finally, the adaptive BD method could be augmented to include hydrodynamic interactions. This requires a careful treatment of the random forces, as they now incorporate the particle configuration $\mathbf{r}^{N}$. Therefore, the noise is no longer additive, which complicates the numerical scheme necessary for a correct discretisation of the overdamped Langevin equation. Still, instead of eqs. (8) and (9), existing methods for SDEs with non-additive noise could be employed within the presented framework to yield an adaptive timestepping procedure for BD with hydrodynamic interactions.

In summary, to capture the dynamics of systems that undergo structural changes such as the ones shown above, and to gain further understanding of the implications regarding e.g. the resulting macroscopic properties, sophisticated simulation methods are needed. Adaptive BD provides the means to treat these problems with great accuracy while still being computationally feasible and robust, which is a crucial trait from a practical standpoint.

The data that support the findings of this study are available from the corresponding author upon reasonable request.