Hi, I am new to LAMMPS and I’m exploring the use of the bond/create and bond/break fixes. I was testing them out using a simple system of two atoms in a hard-walled box, and I noticed something interesting: the probability of a bond forming does not equal the prescribed probability from the “fix bond/create” call. To illustrate this, here’s the relevant parts of my input file:
fix createfix all bond/create 1 1 1 3.0 1 prob ${nominal_prob} 123
fix breakfix all bond/break 1 1 0.0 prob 1 123456
Running the code for different nominal probabilities produces the following plot of Measured Probability of Forming a Bond vs. Nominal Probability of formation: http://i.imgur.com/mjJ2IMK.png
Does anyone have any info on this?
Thanks much for your time,
Scott Grindy
Always a good idea to check. Your results demonstrate an approximate
version of the central limit theorem. In this case, the average of two
uniform variables is a non-uniform distribution with depleted
probability on both ends and enhanced probability in the middle. In
fact, this probability density has a very simple form. It's a
triangle, instead of a rectangle (Thanks to Toby and Jonathan for
figuring that out).
The relevant piece of code in fix_bond_create.cpp reads:
// apply probability constraint
// MIN,MAX insures values are added in same order on different procs
if (fraction < 1.0) {
min = MIN(probability[i],probability[j]);
max = MAX(probability[i],probability[j]);
if (0.5*(min+max) >= fraction) continue;
}
I can see how this might be problematic for values of the keyword prob
that are very small or very large. I see a couple of solutions, but
Steve may ultimately have a better one. In the meantime, try this
change in fix_bond_create.cpp and let us know that it satisfies your
test:
// if (0.5*(min+max) >= fraction) continue;
if (fraction < 0.5) {
if (0.5*(min+max) >= sqrt(fraction/2.0)) continue;
} else {
if (0.5*(min+max) >= 1.0 - sqrt((1.0-fraction)/2.0)) continue;
}
These formula come from inverting the joint cumulative probability
function (your blue points):
P(x) = 2 x^2, 0 <= x <= 1/2
P(x) = 1 - 2 (1-x)^2, 1/2 <= x <= 1
Aidan
Ah yes that works perfectly! (http://i.imgur.com/dd3bqEX.png) Thanks much for your help.
-Scott
Posted a 26May patch for this. This version
does not apply the correction to the triangular
probability distribution - it just makes a simple
choice about which of the 2 random numbers
to use for the check. So you should verify
that it still gives the right answer.
Steve
Yes, this fixes the issue. Thanks much for your help!
-Scott
