# How can I include variable particle number in a Brownian dynamics simulation?

I programmed a Brownian dynamics simulation in two dimensions. (Coarse-grained proteins on surfaces with interaction potentials i.e. patchy particles.) Now I want to allow particles to leave or enter the system, in other words particle number doesn't have to be conserved.

Does anybody know a good way to do that? What are some common pitfalls?

The basis of my simulation is the overdamped Langevin equation $$\frac{d x}{dt} = \frac{D}{k_{B}T}(F^{int} + F^{s}).$$ In my case I have anisotropic diffusion, so I basically have a langevin equation for each degree of freedom of each particle; 2 directional and 1 rotational. The random movement of particles is included in the stochastic force generated from a standard Wiener process. Point patches attached to the particles mediate attractive interactions whereas repulsion is described by a standard soft-sphere potential. The system is held at constant temperature (included in the moments of the random motion). I'm interested in the self-assembled structures at simulation end i.e. clusters of different sizes.

So far my simulation passed important tests and agrees with analytical predictions from statistical calculations (canonical ensemble). In real experiments it is not so easy to fix the particle number since particles diffuse away and onto the surface throughout the experiment. I want to include this observation into my computer simulation which means that I have to remove particles and add particles to my system.

I had the idea to define a rates at which randomly selected particles are removed and added to the system. Now I fear that this approach is to naive and wanted to know if anybody has a suggestion for a proper solution to this or at least some hints. Maybe there is something I'm missing.

I would prefer to have a theoretical backup for my idea such as for example grand canonical Monte Carlo simulations have.

• Are you asking a programming question? Or a physics question?
– user93146
Sep 25, 2018 at 15:01
• Actually a combination of both. I know that there are algorithms for grand canonical ensembles for Monte-Carlo approaches. I'm searching for something similar for BD simulations. Sep 25, 2018 at 15:06
• I'm just thinking out loud. If you have been working with a square or rectangle area until now, you could maybe try it with a torus : a molecule going through the right edge would reappear on the left edge. Same goes for top and bottom edges. You could then crop a smaller rectangle out of the torus surface, and you'd get a variable number of particles. Sep 25, 2018 at 19:43
• I want to thank you for the term "grand canonical Monet-Carlo." Searching for that has caused me to come across some interesting things I should learn about. You taught me something with your question.
– user93146
Sep 25, 2018 at 20:50
• @EricDuminil I use periodic boundary conditions applied to a rectangular surface, isn't this equivalent to a torus? Sep 25, 2018 at 21:31

Introducing particle creation and annihilation to Brownian dynamics will involve similar issues to doing it for molecular dynamics. You can bolt on the standard grand canonical Monte Carlo (GCMC) moves to your dynamical scheme. The practical danger is that, having accepted a move, the consequences for the dynamics will sometimes be dramatic: very large forces between the particles, and hence very large displacements. To avoid this, schemes have been concocted to insert particles gradually, or even introduce a continuous parameter in the Lagrangian which controls the appearance of the extra particle. In Agarwal et al, New J Phys, 17, 083042 (2015), which is open access, they review a few of these methods. As they point out, though, such approaches have not been widely used and are somewhat fiddly. I think that the same is true of the approach proposed by Agarwal. I wouldn't recommend going down that route, but at least you can consider these alternatives, with a suitable adaptation from molecular dynamics to Brownian dynamics.

Here's another possibility. Use Monte Carlo instead of Brownian dynamics. The timescale is somewhat fictitious, but the particles will still diffuse around realistically, and arguably you are abandoning completely realistic dynamics anyway by adding GCMC moves which allow particles to appear and disappear.

There is an intermediate solution. Brownian dynamics without inertia, uses an algorithm $$\vec{r}(t+\delta t) = \vec{r} + \frac{D}{kT} \vec{f}\delta t + \sqrt{2D\delta t}\vec{G}$$ where $$D$$ is the diffusion coefficient, $$T$$ the temperature, $$\vec{f}$$ the force acting on the particle(s), and $$\vec{G}$$ a set of independent normalized Gaussian random numbers. This can be shown to be (almost) equivalent to "Smart Monte Carlo" (SMC), going back to an old paper by Rossky et al, J Chem Phys, 69, 4628 (1978). The essential difference is that SMC applies an acceptance/rejection stage to the above step, based on a slightly complicated formula involving the forces at the start and the end. This guarantees sampling the correct ensemble, and may be the way to save your simulation from the consequences of large forces if you have just added/removed particles using the usual GCMC procedure. It depends whether you are prepared to reject a (small, $$\mathcal{O}(\delta t^2)$$) fraction of advancement steps.

This approach can also be related to "Hybrid Monte Carlo" (HMC), Duane et al, Phys Lett B, 195, 216 (1987). Simply substitute $$D=(kT/2m)\delta t$$ in the above equation to give $$\vec{r}(t+\delta t) = \vec{r} + \tfrac{1}{2}(\delta t^2/m) \vec{f} + \sqrt{\frac{kT}{m}}\delta t\vec{G}$$ This has the form of a standard algorithm (velocity Verlet) for molecular dynamics, where the last term corresponds to choosing random velocities at the start of each step, from the Maxwell-Boltzmann distribution. In HMC, you also write down the rest of the velocity Verlet algorithm, to advance the velocities, compute the change in kinetic energy over the step, add it to the change in potential energy, and use the resultant total energy in a Metropolis acceptance/rejection criterion. The velocities are then thrown away, and generated afresh at the start of the next step. This turns out to be exactly the same as SMC, but simpler to write down. The reason for the small fraction of rejected moves, $$\mathcal{O}(\delta t^2)$$, is more obvious: the total energy is conserved to this order.

If your Brownian dynamics is of the non-inertial kind, this might be the approach I'd recommend:

1. Insert/remove particles at intervals (between time-advancement steps) using the standard grand canonical Monte Carlo method
2. Include a Metropolis acceptance/rejection decision on the time-advancement steps, as described above, meaning that a small fraction of moves will be rejected, but guaranteeing that the right ensemble will be sampled.

I have a simpler suggestion, which might be preferable or not depending on your circumstances. Stick with your existing Brownian dynamics, but make your simulation three dimensional, rather than two dimensional, with a flat surface that attracts your coarse-grained particles. The rest of the system would be a low-density gas acting as a reservoir of particles. A lot would depend on whether you could adjust the parameters to make the thermodynamically stable state of your system a monolayer adsorbed on the surface. If this is possible, you should see particles arriving and departing from the surface in a physically reasonable way.

• Thank you very much for taking the time to answer my question. I decided to do a BD simulation to have easy access to each particle position. With BD it is also possible to account for flexibility of the binding interface naturally. Your last suggestion might work but I think the implementation (changes in my source code) will be very time consuming. Nevertheless, you mentioned a few interesting things and I check them one by one. Sep 25, 2018 at 21:22