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