I saw on the 8/22 patch that Monte Carlo bond swapping will be added soon. Will this be the algorithm for equilibrating long-chain polymers and brushes that was in the FORTRAN LAMMPS, or something similar? That would be great!
I want to write a Monte Carlo routines for forming and breaking bonds between the ends of bead-spring chains. Looking at bond_quartic.cpp, it seems that when a bond is broken, it is not deleted, but rather set to type 0. Now, suppose you wanted to reform the bond - would you just reset it to type 1? (assuming the normal bonds are type 1)
What I have in mind, as a way to avoid messing with the bond list that is created during startup, is to list all pairs of ends as bonded in my initial configuration file, but have all these bonds initially set to type 0. Then, when a pair of chains ends passes close to one another, "activate" the interactions for that bond (if it is favorable in a Monte Carlo sense) by switching it to type 1, which will be a FENE bond, offset so that U_bind = -1 (or something of order kT). Then, keep the bond intact until it is favorable in a Monte Carlo sense to break it, at which point switch it back to type zero. This could be repeated arbirarily many times.
Then after each breaking or formation I would have to either do this or the reverse of this, respectively:
if (rsq > rc[type]*rc[type]) {
// IF criterion would have to be suitably modified
bondlist[n][2] = 0;
for (m = 0; m < atom->num_bond[i1]; m++)
if (atom->bond_atom[i1][m] == atom->tag[i2])
atom->bond_type[i1][m] = 0;
// change "0" to "1" if forming bond?
if (i2 < atom->nlocal)
for (m = 0; m < atom->num_bond[i2]; m++)
if (atom->bond_atom[i2][m] == atom->tag[i1])
atom->bond_type[i2][m] = 0;
continue;
}
Am I on the right track? Any advice or pitfalls, in particular for multiprocessor runs? I know interior bonds on chains would have to be type 2 or something so they would not be involved.
Thanks,
Rob