Fix bond/create and fix bond/break probabilities

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:

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

      // 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


Ah yes that works perfectly! ( Thanks much for your help.


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.


Yes, this fixes the issue. Thanks much for your help!