[lammps-users] bond/swap MC attempt limits

Hi everybody! Fix bond/swap imposes a rule that at most one MC swap move per processor can be attempted each time the fix routine is called, i.e. at most one swap move per processor is attempted every N_every timesteps. As I understand things, this was implemented to prevent “reverse moves” (i.e. a swap being executed and reversed in the same timestep). Clearly one doesn’t want to have many reverse moves.

However, the current implementation of fix bond/swap introduces a strong N_proc-dependence into the code. Specifically, N_proc times as many moves per unit time are attempted in a N_proc-processor parallel LAMMPS run as in a serial LAMMPS run, which makes the final state of any such run strongly N_proc-dependent. There are many situations in which this is undesirable.

Looking at the post_integrate() routine in the current (10/29/20) version of fix_bond_swap.cpp, it seems like one could remove the 1-attempt-per-processor restriction by replacing “goto done;” (line 425) with a slightly-modified version of the code that sets “accept” and executes if “accept” is nonzero, i.e.

if (accept) {

naccept++;

next_reneighbor = update->ntimestep;

// change bond partners of affected atoms

// on atom i: bond i-inext changes to i-jnext

// on atom j: bond j-jnext changes to j-inext

// on atom inext: bond inext-i changes to inext-j

// on atom jnext: bond jnext-j changes to jnext-i

for (ibond = 0; ibond < num_bond[i]; ibond++)

if (bond_atom[i][ibond] == tag[inext]) bond_atom[i][ibond] = tag[jnext];

for (jbond = 0; jbond < num_bond[j]; jbond++)

if (bond_atom[j][jbond] == tag[jnext]) bond_atom[j][jbond] = tag[inext];

for (ibond = 0; ibond < num_bond[inext]; ibond++)

if (bond_atom[inext][ibond] == tag[i]) bond_atom[inext][ibond] = tag[j];

for (jbond = 0; jbond < num_bond[jnext]; jbond++)

if (bond_atom[jnext][jbond] == tag[j]) bond_atom[jnext][jbond] = tag[i];

// set global tags of 4 atoms in bonds

itag = tag[i];

inexttag = tag[inext];

jtag = tag[j];

jnexttag = tag[jnext];

// change 1st special neighbors of affected atoms: i,j,inext,jnext

// don’t need to change 2nd/3rd special neighbors for any atom

// since special bonds = 0 1 1 means they are never used

for (m = 0; m < nspecial[i][0]; m++)

if (special[i][m] == inexttag) special[i][m] = jnexttag;

for (m = 0; m < nspecial[j][0]; m++)

if (special[j][m] == jnexttag) special[j][m] = inexttag;

for (m = 0; m < nspecial[inext][0]; m++)

if (special[inext][m] == itag) special[inext][m] = jtag;

for (m = 0; m < nspecial[jnext][0]; m++)

if (special[jnext][m] == jtag) special[jnext][m] = itag;

// done if no angles

if (!angleflag) return;

// set global tags of 4 additional atoms in angles, 0 if no angle

if (iprev >= 0) iprevtag = tag[iprev];

else iprevtag = 0;

if (ilast >= 0) ilasttag = tag[ilast];

else ilasttag = 0;

if (jprev >= 0) jprevtag = tag[jprev];

else jprevtag = 0;

if (jlast >= 0) jlasttag = tag[jlast];

else jlasttag = 0;

// change angle partners of affected atoms

// must check if each angle is stored as a-b-c or c-b-a

// on atom i:

// angle iprev-i-inext changes to iprev-i-jnext

// angle i-inext-ilast changes to i-jnext-jlast

// on atom j:

// angle jprev-j-jnext changes to jprev-j-inext

// angle j-jnext-jlast changes to j-inext-ilast

// on atom inext:

// angle iprev-i-inext changes to jprev-j-inext

// angle i-inext-ilast changes to j-inext-ilast

// on atom jnext:

// angle jprev-j-jnext changes to iprev-i-jnext

// angle j-jnext-jlast changes to i-jnext-jlast

for (iangle = 0; iangle < num_angle[i]; iangle++) {

i1 = angle_atom1[i][iangle];

i2 = angle_atom2[i][iangle];

i3 = angle_atom3[i][iangle];

if (i1 == iprevtag && i2 == itag && i3 == inexttag)

angle_atom3[i][iangle] = jnexttag;

else if (i1 == inexttag && i2 == itag && i3 == iprevtag)

angle_atom1[i][iangle] = jnexttag;

else if (i1 == itag && i2 == inexttag && i3 == ilasttag) {

angle_atom2[i][iangle] = jnexttag;

angle_atom3[i][iangle] = jlasttag;

} else if (i1 == ilasttag && i2 == inexttag && i3 == itag) {

angle_atom1[i][iangle] = jlasttag;

angle_atom2[i][iangle] = jnexttag;

}

}

for (jangle = 0; jangle < num_angle[j]; jangle++) {

j1 = angle_atom1[j][jangle];

j2 = angle_atom2[j][jangle];

j3 = angle_atom3[j][jangle];

if (j1 == jprevtag && j2 == jtag && j3 == jnexttag)

angle_atom3[j][jangle] = inexttag;

else if (j1 == jnexttag && j2 == jtag && j3 == jprevtag)

angle_atom1[j][jangle] = inexttag;

else if (j1 == jtag && j2 == jnexttag && j3 == jlasttag) {

angle_atom2[j][jangle] = inexttag;

angle_atom3[j][jangle] = ilasttag;

} else if (j1 == jlasttag && j2 == jnexttag && j3 == jtag) {

angle_atom1[j][jangle] = ilasttag;

angle_atom2[j][jangle] = inexttag;

}

}

for (iangle = 0; iangle < num_angle[inext]; iangle++) {

i1 = angle_atom1[inext][iangle];

i2 = angle_atom2[inext][iangle];

i3 = angle_atom3[inext][iangle];

if (i1 == iprevtag && i2 == itag && i3 == inexttag) {

angle_atom1[inext][iangle] = jprevtag;

angle_atom2[inext][iangle] = jtag;

} else if (i1 == inexttag && i2 == itag && i3 == iprevtag) {

angle_atom2[inext][iangle] = jtag;

angle_atom3[inext][iangle] = jprevtag;

} else if (i1 == itag && i2 == inexttag && i3 == ilasttag)

angle_atom1[inext][iangle] = jtag;

else if (i1 == ilasttag && i2 == inexttag && i3 == itag)

angle_atom3[inext][iangle] = jtag;

}

for (jangle = 0; jangle < num_angle[jnext]; jangle++) {

j1 = angle_atom1[jnext][jangle];

j2 = angle_atom2[jnext][jangle];

j3 = angle_atom3[jnext][jangle];

if (j1 == jprevtag && j2 == jtag && j3 == jnexttag) {

angle_atom1[jnext][jangle] = iprevtag;

angle_atom2[jnext][jangle] = itag;

} else if (j1 == jnexttag && j2 == jtag && j3 == jprevtag) {

angle_atom2[jnext][jangle] = itag;

angle_atom3[jnext][jangle] = iprevtag;

} else if (j1 == jtag && j2 == jnexttag && j3 == jlasttag)

angle_atom1[jnext][jangle] = itag;

else if (j1 == jlasttag && j2 == jnexttag && j3 == jtag)

angle_atom3[jnext][jangle] = itag;

}

// done if newton bond set

if (newton_bond) return;

// change angles stored by iprev,ilast,jprev,jlast

// on atom iprev: angle iprev-i-inext changes to iprev-i-jnext

// on atom jprev: angle jprev-j-jnext changes to jprev-j-inext

// on atom ilast: angle i-inext-ilast changes to j-inext-ilast

// on atom jlast: angle j-jnext-jlast changes to i-jnext-jlast

for (iangle = 0; iangle < num_angle[iprev]; iangle++) {

i1 = angle_atom1[iprev][iangle];

i2 = angle_atom2[iprev][iangle];

i3 = angle_atom3[iprev][iangle];

if (i1 == iprevtag && i2 == itag && i3 == inexttag)

angle_atom3[iprev][iangle] = jnexttag;

else if (i1 == inexttag && i2 == itag && i3 == iprevtag)

angle_atom1[iprev][iangle] = jnexttag;

}

for (jangle = 0; jangle < num_angle[jprev]; jangle++) {

j1 = angle_atom1[jprev][jangle];

j2 = angle_atom2[jprev][jangle];

j3 = angle_atom3[jprev][jangle];

if (j1 == jprevtag && j2 == jtag && j3 == jnexttag)

angle_atom3[jprev][jangle] = inexttag;

else if (j1 == jnexttag && j2 == jtag && j3 == jprevtag)

angle_atom1[jprev][jangle] = inexttag;

}

for (iangle = 0; iangle < num_angle[ilast]; iangle++) {

i1 = angle_atom1[ilast][iangle];

i2 = angle_atom2[ilast][iangle];

i3 = angle_atom3[ilast][iangle];

if (i1 == itag && i2 == inexttag && i3 == ilasttag)

angle_atom1[ilast][iangle] = jtag;

else if (i1 == ilasttag && i2 == inexttag && i3 == itag)

angle_atom3[ilast][iangle] = jtag;

}

for (jangle = 0; jangle < num_angle[jlast]; jangle++) {

j1 = angle_atom1[jlast][jangle];

j2 = angle_atom2[jlast][jangle];

j3 = angle_atom3[jlast][jangle];

if (j1 == jtag && j2 == jnexttag && j3 == jlasttag)

angle_atom1[jlast][jangle] = itag;

else if (j1 == jlasttag && j2 == jnexttag && j3 == jtag)

angle_atom3[jlast][jangle] = itag;

}

}

Is this correct? If not, could someone please specify how best to remove the 1-attempt-per-processor restriction?

Also, if we implement this, could we reduce the # of reverse moves by reducing the value of “fraction”? Or better yet, is there a straightforward way to explicitly forbid reverse moves by adding some additional code?

Thanks,
Rob

i cannot comment on the subject itself, since I currently don’t have time to look into these kinds of issues.
but this kind of comment/suggestion/question is probably better submitted as a feature request issue on the LAMMPS github project.
mailing list emails have the disadvantage to “fade out of focus” as more emails are sent to the list while project issues require to be actively closed to be removed from view.

axel.

Hi Axel. That’s a good suggestion, I’ve never used GitHub but apparently it’s time to start :slight_smile:

So far I’ve figured out that at the very least the “return” statements within the if(accept) {} block need to be replaced with goto’s…

Anyway, since I did post the question to this list, anyone else have any suggestions?

Thanks,
Rob

Can you effectively turn off the dependence on the number of procs
by adjusting the frequency of invocation or other params. And thus
get similar behavior for any number of procs?

Steve

Hi Steve! So far it seems that we can do this by adjusting N_every in such a way that N_every*N_proc remains constant; this seems to give a nearly-constant # of successful swaps, as one would expect. I think it would still be useful at some point to allow multiple swap attempts per processor per timestep, but I guess we don’t NEED to tackle that right now.

Thanks!
Rob