Mixing
How can I include variable particle number in a Brownian dynamics simulation?
Clash Royale CLAN TAG#URR8PPP
up vote
5
down vote
favorite
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
$$ fracd xdt = fracDk_BT(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.
computational-physics simulations molecular-dynamics brownian-motion
edited 15 mins ago
knzhou
35.3k8100172
asked 7 hours ago
DW147
262
New contributor
DW147 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
add a comment |Â
up vote
5
down vote
favorite
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
$$ fracd xdt = fracDk_BT(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.
computational-physics simulations molecular-dynamics brownian-motion
edited 15 mins ago
knzhou
35.3k8100172
asked 7 hours ago
DW147
262
New contributor
DW147 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
Are you asking a programming question? Or a physics question?
– puppetsock
7 hours ago
1
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.
– DW147
7 hours ago
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.
– Eric Duminil
2 hours ago
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.
– puppetsock
1 hour ago
@EricDuminil I use periodic boundary conditions applied to a rectangular surface, isn't this equivalent to a torus?
– DW147
37 mins ago
add a comment |Â
up vote
5
down vote
favorite
up vote
5
down vote
favorite
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
$$ fracd xdt = fracDk_BT(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.
computational-physics simulations molecular-dynamics brownian-motion
edited 15 mins ago
knzhou
35.3k8100172
asked 7 hours ago
DW147
262
New contributor
DW147 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
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
$$ fracd xdt = fracDk_BT(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.
computational-physics simulations molecular-dynamics brownian-motion
computational-physics simulations molecular-dynamics brownian-motion
edited 15 mins ago
knzhou
35.3k8100172
asked 7 hours ago
DW147
262
New contributor
DW147 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
edited 15 mins ago
knzhou
35.3k8100172
asked 7 hours ago
DW147
262
New contributor
DW147 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
edited 15 mins ago
knzhou
35.3k8100172
edited 15 mins ago
knzhou
35.3k8100172
edited 15 mins ago
knzhou
35.3k8100172
35.3k8100172
asked 7 hours ago
DW147
262
New contributor
DW147 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
asked 7 hours ago
DW147
262
asked 7 hours ago
DW147
262
262
New contributor
DW147 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
New contributor
DW147 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
DW147 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
Are you asking a programming question? Or a physics question?
– puppetsock
7 hours ago
1
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.
– DW147
7 hours ago
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.
– Eric Duminil
2 hours ago
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.
– puppetsock
1 hour ago
@EricDuminil I use periodic boundary conditions applied to a rectangular surface, isn't this equivalent to a torus?
– DW147
37 mins ago
add a comment |Â
Are you asking a programming question? Or a physics question?
– puppetsock
7 hours ago
1
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.
– DW147
7 hours ago
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.
– Eric Duminil
2 hours ago
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.
– puppetsock
1 hour ago
@EricDuminil I use periodic boundary conditions applied to a rectangular surface, isn't this equivalent to a torus?
– DW147
37 mins ago
Are you asking a programming question? Or a physics question?
– puppetsock
7 hours ago
Are you asking a programming question? Or a physics question?
– puppetsock
7 hours ago
1
1
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.
– DW147
7 hours ago
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.
– DW147
7 hours ago
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.
– Eric Duminil
2 hours ago
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.
– Eric Duminil
2 hours ago
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.
– puppetsock
1 hour ago
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.
– puppetsock
1 hour ago
@EricDuminil I use periodic boundary conditions applied to a rectangular surface, isn't this equivalent to a torus?
– DW147
37 mins ago
@EricDuminil I use periodic boundary conditions applied to a rectangular surface, isn't this equivalent to a torus?
– DW147
37 mins ago
add a comment |Â
2 Answers
2
active
oldest
votes
up vote
4
down vote
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
$$
vecr(t+delta t) = vecr + fracDkT vecfdelta t
+ sqrt2Ddelta tvecG
$$
where $D$ is the diffusion coefficient, $T$ the temperature,
$vecf$ the force acting on the particle(s),
and $vecG$ 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, $mathcalO(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
$$
vecr(t+delta t) = vecr + tfrac12(delta t^2/m) vecf
+ sqrtfrackTmdelta tvecG
$$
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,
$mathcalO(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:
- Insert/remove particles at intervals (between time-advancement steps) using the
standard grand canonical Monte Carlo method - 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.
answered 3 hours ago
LonelyProf
1,6981210
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.
– DW147
47 mins ago
add a comment |Â
up vote
2
down vote
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.
answered 6 hours ago
puppetsock
1,40737
add a comment |Â
2 Answers
2
active
oldest
votes
2 Answers
2
active
oldest
votes
active
oldest
votes
active
oldest
votes
up vote
4
down vote
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
$$
vecr(t+delta t) = vecr + fracDkT vecfdelta t
+ sqrt2Ddelta tvecG
$$
where $D$ is the diffusion coefficient, $T$ the temperature,
$vecf$ the force acting on the particle(s),
and $vecG$ 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, $mathcalO(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
$$
vecr(t+delta t) = vecr + tfrac12(delta t^2/m) vecf
+ sqrtfrackTmdelta tvecG
$$
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,
$mathcalO(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:
- Insert/remove particles at intervals (between time-advancement steps) using the
standard grand canonical Monte Carlo method - 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.
answered 3 hours ago
LonelyProf
1,6981210
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.
– DW147
47 mins ago
add a comment |Â
up vote
4
down vote
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
$$
vecr(t+delta t) = vecr + fracDkT vecfdelta t
+ sqrt2Ddelta tvecG
$$
where $D$ is the diffusion coefficient, $T$ the temperature,
$vecf$ the force acting on the particle(s),
and $vecG$ 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, $mathcalO(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
$$
vecr(t+delta t) = vecr + tfrac12(delta t^2/m) vecf
+ sqrtfrackTmdelta tvecG
$$
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,
$mathcalO(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:
- Insert/remove particles at intervals (between time-advancement steps) using the
standard grand canonical Monte Carlo method - 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.
answered 3 hours ago
LonelyProf
1,6981210
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.
– DW147
47 mins ago
add a comment |Â
up vote
4
down vote
up vote
4
down vote
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
$$
vecr(t+delta t) = vecr + fracDkT vecfdelta t
+ sqrt2Ddelta tvecG
$$
where $D$ is the diffusion coefficient, $T$ the temperature,
$vecf$ the force acting on the particle(s),
and $vecG$ 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, $mathcalO(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
$$
vecr(t+delta t) = vecr + tfrac12(delta t^2/m) vecf
+ sqrtfrackTmdelta tvecG
$$
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,
$mathcalO(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:
- Insert/remove particles at intervals (between time-advancement steps) using the
standard grand canonical Monte Carlo method - 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.
answered 3 hours ago
LonelyProf
1,6981210
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
$$
vecr(t+delta t) = vecr + fracDkT vecfdelta t
+ sqrt2Ddelta tvecG
$$
where $D$ is the diffusion coefficient, $T$ the temperature,
$vecf$ the force acting on the particle(s),
and $vecG$ 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, $mathcalO(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
$$
vecr(t+delta t) = vecr + tfrac12(delta t^2/m) vecf
+ sqrtfrackTmdelta tvecG
$$
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,
$mathcalO(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:
- Insert/remove particles at intervals (between time-advancement steps) using the
standard grand canonical Monte Carlo method - 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.
answered 3 hours ago
LonelyProf
1,6981210
edited 1 hour ago
answered 3 hours ago
LonelyProf
1,6981210
answered 3 hours ago
LonelyProf
1,6981210
answered 3 hours ago
LonelyProf
1,6981210
1,6981210
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.
– DW147
47 mins ago
add a comment |Â
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.
– DW147
47 mins ago
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.
– DW147
47 mins ago
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.
– DW147
47 mins ago
add a comment |Â
up vote
2
down vote
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.
answered 6 hours ago
puppetsock
1,40737
add a comment |Â
up vote
2
down vote
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.
answered 6 hours ago
puppetsock
1,40737
add a comment |Â
up vote
2
down vote
up vote
2
down vote
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.
answered 6 hours ago
puppetsock
1,40737
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.
answered 6 hours ago
puppetsock
1,40737
answered 6 hours ago
puppetsock
1,40737
answered 6 hours ago
puppetsock
1,40737
answered 6 hours ago
puppetsock
1,40737
1,40737
add a comment |Â
add a comment |Â
DW147 is a new contributor. Be nice, and check out our Code of Conduct.
DW147 is a new contributor. Be nice, and check out our Code of Conduct.
DW147 is a new contributor. Be nice, and check out our Code of Conduct.
DW147 is a new contributor. Be nice, and check out our Code of Conduct.
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fphysics.stackexchange.com%2fquestions%2f430761%2fhow-can-i-include-variable-particle-number-in-a-brownian-dynamics-simulation%23new-answer', 'question_page');
);
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Are you asking a programming question? Or a physics question?
– puppetsock
7 hours ago
1
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.
– DW147
7 hours ago
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.
– Eric Duminil
2 hours ago
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.
– puppetsock
1 hour ago
@EricDuminil I use periodic boundary conditions applied to a rectangular surface, isn't this equivalent to a torus?
– DW147
37 mins ago