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

Question

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.

Answer 1

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.

Answer 2

Your question is quite broad. So my answer is going to be quite broad.

Some of the usual things to worry about in Monte-Carlo calculations include the following. But there are more.

• Have you carefully recorded the assumptions you will be using to form your calculations? Maybe things like constant temperature, no fluid flow, gravity neglected, etc. You should be looking for choices you made that affect how you did the calculations as opposed to trivial details that don't make any difference.
• Have you paid adequate attention to the usual mundane issues of software quality and correctness? Yada yada, useful variable names, structured program, basic testing, at least minimal documenting of program behavior, ease of input, sensible arrangement of output, etc.
• Have you correctly normalized your calculation? You may have 1 million particles in your calculation, but there are many more molecules involved in a real experiment. Have you got a reasonable way to scale that? This applies to each aspect of your calculation. How do you set the temperature of the working material? How set the density? How to scale tallies? And so on.
• Have you conserved the quantities that should be conserved? You want to change particle number. Does it match the physical process in your experiment? For example, particles passing the edge of your experiment. Or you might have a place you inject particles. Maybe you need to worry about energy or total mass or ratio of particle type.
• Have you got a handle on how the things that should change are actually changing? So you inject new particles and watch them go through Brownian motion. Does that change local conditions in any important fashion? Density, viscosity, temperature, etc.
• Have you got tallies that make sense? Maybe you are looking for particle number density at various locations. Are you counting the right things, and counting them reasonably?
• Have you got the system conditions at reasonable values? So the particles move under Brownian motion. Have you got reasonable values for the speed they move? The average distance they move before colliding? The angle they move after? The energy change in a collision?
• Have you got some kind of authoritative data source for validation? Some frequently used sources are experiment, calculations from some other program that is known to be accurate, or hand calculations that are of known accuracy.