From f8212fdb31dc300521870fbe29b987b4314dd0a4 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 6 Dec 2021 13:11:51 -0700 Subject: [PATCH 1/7] add randomization and debug prints --- src/MC/fix_bond_swap.cpp | 71 +++++++++++++++++++++++++++++++++++++--- src/MC/fix_bond_swap.h | 5 +++ 2 files changed, 72 insertions(+), 4 deletions(-) diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp index 5f7ffcbbe5..8332180dc0 100644 --- a/src/MC/fix_bond_swap.cpp +++ b/src/MC/fix_bond_swap.cpp @@ -49,6 +49,8 @@ static const char cite_fix_bond_swap[] = " pages = {12718--12728}\n" "}\n\n"; +#define DELTA_PERMUTE 100 + /* ---------------------------------------------------------------------- */ FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) : @@ -91,11 +93,14 @@ FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) : modify->add_compute(fmt::format("{} all temp",id_temp)); tflag = 1; - // initialize atom list + // initialize two permutation lists nmax = 0; alist = nullptr; + maxpermute = 0; + permute = nullptr; + naccept = foursome = 0; } @@ -111,6 +116,7 @@ FixBondSwap::~FixBondSwap() delete [] id_temp; memory->destroy(alist); + delete [] permute; } /* ---------------------------------------------------------------------- */ @@ -248,8 +254,8 @@ void FixBondSwap::post_integrate() } int tmp; - for (i = 0; i < neligible; i++) { - j = static_cast (random->uniform() * neligible); + for (i = 0; i < neligible-1; i++) { + j = i + static_cast (random->uniform() * (neligible-i)); tmp = alist[i]; alist[i] = alist[j]; alist[j] = tmp; @@ -277,8 +283,35 @@ void FixBondSwap::post_integrate() jlist = firstneigh[i]; jnum = numneigh[i]; + neighbor_permutation(jnum); + + // DEBUG + + int picked[4]; + int count = 0; + printf("BSWAP CANDIDATES nstep %ld itag %d jnum %d jtags", + update->ntimestep,atom->tag[i],jnum); for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; + j = jlist[permute[jj]]; + j &= NEIGHMASK; + printf(" %d",atom->tag[j]); + if (j >= nlocal) continue; + if ((mask[j] & groupbit) == 0) continue; + if (molecule[i] != molecule[j]) continue; + picked[count++] = atom->tag[j]; + } + if (count) { + printf(" picked %d",count); + for (jj = 0; jj < count; jj++) { + printf(" %d",picked[jj]); + } + } + printf("\n"); + + // END DEBUG + + for (jj = 0; jj < jnum; jj++) { + j = jlist[permute[jj]]; j &= NEIGHMASK; if (j >= nlocal) continue; if ((mask[j] & groupbit) == 0) continue; @@ -297,6 +330,8 @@ void FixBondSwap::post_integrate() // already know i != inext, j != jnext // all 4 old and new bonds must have length < cutoff + printf("SWAP candidate jtag\n",atom->tag[j]); + for (ibond = 0; ibond < num_bond[i]; ibond++) { inext = atom->map(bond_atom[i][ibond]); if (inext >= nlocal || inext < 0) continue; @@ -715,6 +750,34 @@ double FixBondSwap::angle_eng(int atype, int i, int j, int k) return force->angle->single(atype,i,j,k); } +/* ---------------------------------------------------------------------- + create a random permutation of one atom's N neighbor list atoms + uses one-pass Fisher-Yates shuffle on an initial identity permutation + output: randomized permute[] vector, used to index neighbors +------------------------------------------------------------------------- */ + +void FixBondSwap::neighbor_permutation(int n) +{ + int i,j,tmp; + + if (n > maxpermute) { + delete [] permute; + maxpermute = n + DELTA_PERMUTE; + permute = new int[maxpermute]; + } + + // Fisher-Yates shuffle + + for (i = 0; i < n; i++) permute[i] = i; + + for (i = 0; i < n-1; i++) { + j = i + static_cast (random->uniform() * (n-i)); + tmp = permute[i]; + permute[i] = permute[j]; + permute[j] = tmp; + } +} + /* ---------------------------------------------------------------------- return bond swapping stats n = 1 is # of swaps diff --git a/src/MC/fix_bond_swap.h b/src/MC/fix_bond_swap.h index e0691233f5..e927616246 100644 --- a/src/MC/fix_bond_swap.h +++ b/src/MC/fix_bond_swap.h @@ -46,6 +46,9 @@ class FixBondSwap : public Fix { int *type; double **x; + int maxpermute; + int *permute; + class NeighList *list; class Compute *temperature; class RanMars *random; @@ -54,6 +57,8 @@ class FixBondSwap : public Fix { double pair_eng(int, int); double bond_eng(int, int, int); double angle_eng(int, int, int, int); + + void neighbor_permutation(int); }; } // namespace LAMMPS_NS From 66c686f73307a34a40d7fd931393df0ff4595afd Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 30 Mar 2022 11:56:07 -0600 Subject: [PATCH 2/7] debug info --- src/MC/fix_bond_swap.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp index 8332180dc0..556c99d93a 100644 --- a/src/MC/fix_bond_swap.cpp +++ b/src/MC/fix_bond_swap.cpp @@ -330,7 +330,7 @@ void FixBondSwap::post_integrate() // already know i != inext, j != jnext // all 4 old and new bonds must have length < cutoff - printf("SWAP candidate jtag\n",atom->tag[j]); + printf(" SWAP candidate jtag %d\n",atom->tag[j]); for (ibond = 0; ibond < num_bond[i]; ibond++) { inext = atom->map(bond_atom[i][ibond]); @@ -351,6 +351,9 @@ void FixBondSwap::post_integrate() if (dist_rsq(i,jnext) >= cutsq) continue; if (dist_rsq(j,inext) >= cutsq) continue; + printf(" BOND jtag %d ibondtag %d jbondtag %d\n", + atom->tag[j],atom->tag[inext],atom->tag[jnext]); + // if angles are enabled: // find other atoms i,inext,j,jnext are in angles with // and angletypes: i/j angletype, i/j nextangletype @@ -468,6 +471,10 @@ void FixBondSwap::post_integrate() factor = exp(-delta/force->boltz/t_current); if (random->uniform() < factor) accept = 1; } + + printf(" MC result jtag %d delta %g accept %d\n", + atom->tag[j],delta,accept); + goto done; } } From 852c5f13fff70da0edf04dac16680173bdf7451b Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 30 Mar 2022 15:35:57 -0600 Subject: [PATCH 3/7] more debugging --- src/MC/fix_bond_swap.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp index 556c99d93a..55f7a112d0 100644 --- a/src/MC/fix_bond_swap.cpp +++ b/src/MC/fix_bond_swap.cpp @@ -330,7 +330,7 @@ void FixBondSwap::post_integrate() // already know i != inext, j != jnext // all 4 old and new bonds must have length < cutoff - printf(" SWAP candidate jtag %d\n",atom->tag[j]); + printf(" SWAP candidate ijtag %d %d\n",atom->tag[i],atom->tag[j]); for (ibond = 0; ibond < num_bond[i]; ibond++) { inext = atom->map(bond_atom[i][ibond]); @@ -351,8 +351,9 @@ void FixBondSwap::post_integrate() if (dist_rsq(i,jnext) >= cutsq) continue; if (dist_rsq(j,inext) >= cutsq) continue; - printf(" BOND jtag %d ibondtag %d jbondtag %d\n", - atom->tag[j],atom->tag[inext],atom->tag[jnext]); + printf(" BOND nstep %ld ijtag %d %d ibondtag %d jbondtag %d\n", + update->ntimestep, + atom->tag[i],atom->tag[j],atom->tag[inext],atom->tag[jnext]); // if angles are enabled: // find other atoms i,inext,j,jnext are in angles with @@ -472,7 +473,7 @@ void FixBondSwap::post_integrate() if (random->uniform() < factor) accept = 1; } - printf(" MC result jtag %d delta %g accept %d\n", + printf(" MC result jtag %d delta %g accept %d\n", atom->tag[j],delta,accept); goto done; From 572e8cc399f1103a4660eff211da0a1302d1209e Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 12 Sep 2022 09:18:06 -0600 Subject: [PATCH 4/7] random permulations of both I and J atoms in double loop --- src/MC/fix_bond_swap.cpp | 44 +++++++++------------------------------- 1 file changed, 10 insertions(+), 34 deletions(-) diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp index 99a5640119..ea31270b92 100644 --- a/src/MC/fix_bond_swap.cpp +++ b/src/MC/fix_bond_swap.cpp @@ -244,6 +244,15 @@ void FixBondSwap::post_integrate() memory->create(alist,nmax,"bondswap:alist"); } + // use randomized permulation of both I and J atoms in double loop below + // this is to avoid any bias in accepted MC swaps based on + // ordering LAMMPS creates on a processor for atoms or their neighbors + + // create a random permutation of list of Neligible atoms + // uses one-pass Fisher-Yates shuffle on an initial identity permutation + // output: randomized alist[] vector, used in outer loop to select an I atom + // similar randomized permulation is created for neighbors of each I atom + int neligible = 0; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; @@ -262,6 +271,7 @@ void FixBondSwap::post_integrate() // examine ntest of my eligible atoms for potential swaps // atom I is randomly selected via atom list // look at all J neighbors of atom I + // done in a randomized permutation, via neighbor_permutation() // J must be on-processor (J < nlocal) // I,J must be in fix group // I,J must have same molecule IDs @@ -283,31 +293,6 @@ void FixBondSwap::post_integrate() neighbor_permutation(jnum); - // DEBUG - - int picked[4]; - int count = 0; - printf("BSWAP CANDIDATES nstep %ld itag %d jnum %d jtags", - update->ntimestep,atom->tag[i],jnum); - for (jj = 0; jj < jnum; jj++) { - j = jlist[permute[jj]]; - j &= NEIGHMASK; - printf(" %d",atom->tag[j]); - if (j >= nlocal) continue; - if ((mask[j] & groupbit) == 0) continue; - if (molecule[i] != molecule[j]) continue; - picked[count++] = atom->tag[j]; - } - if (count) { - printf(" picked %d",count); - for (jj = 0; jj < count; jj++) { - printf(" %d",picked[jj]); - } - } - printf("\n"); - - // END DEBUG - for (jj = 0; jj < jnum; jj++) { j = jlist[permute[jj]]; j &= NEIGHMASK; @@ -328,8 +313,6 @@ void FixBondSwap::post_integrate() // already know i != inext, j != jnext // all 4 old and new bonds must have length < cutoff - printf(" SWAP candidate ijtag %d %d\n",atom->tag[i],atom->tag[j]); - for (ibond = 0; ibond < num_bond[i]; ibond++) { inext = atom->map(bond_atom[i][ibond]); if (inext >= nlocal || inext < 0) continue; @@ -349,10 +332,6 @@ void FixBondSwap::post_integrate() if (dist_rsq(i,jnext) >= cutsq) continue; if (dist_rsq(j,inext) >= cutsq) continue; - printf(" BOND nstep %ld ijtag %d %d ibondtag %d jbondtag %d\n", - update->ntimestep, - atom->tag[i],atom->tag[j],atom->tag[inext],atom->tag[jnext]); - // if angles are enabled: // find other atoms i,inext,j,jnext are in angles with // and angletypes: i/j angletype, i/j nextangletype @@ -471,9 +450,6 @@ void FixBondSwap::post_integrate() if (random->uniform() < factor) accept = 1; } - printf(" MC result jtag %d delta %g accept %d\n", - atom->tag[j],delta,accept); - goto done; } } From 58d6e0d440e5e8b81e7b2220b781a68cf509110f Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 12 Sep 2022 09:41:29 -0600 Subject: [PATCH 5/7] add note to doc page about sticker site model --- doc/src/fix_bond_swap.rst | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/doc/src/fix_bond_swap.rst b/doc/src/fix_bond_swap.rst index 0a237bbf78..2c15664c09 100644 --- a/doc/src/fix_bond_swap.rst +++ b/doc/src/fix_bond_swap.rst @@ -71,6 +71,15 @@ the polymer chains all become interconnected. For this use case, if angles are defined they should not include bonds between sticker sites. +.. note:: + + For the sticker site model, you should set the newton flag for + bonds to "off", via the :doc:`newton on off` command ("on" + is the default for the 2nd argument). This is to ensure + appropriate randomness in bond selection because the I,J bond will + be stored by both atom I and atom J. LAMMPS cannot check for this, + because it is not aware that a sticker site model is being used. + ---------- The bond swapping operation is invoked once every *Nevery* timesteps. @@ -114,11 +123,11 @@ Boltzmann constant, and T is the current temperature of the system. .. note:: - IMPORTANT: Whether the swap is accepted or rejected, no other swaps - are attempted by this processor on this timestep. No other - eligible 4-tuples of atoms are considered. This means that each - processor will perform either a single swap or none on timesteps - this fix is invoked. + Whether the swap is accepted or rejected, no other swaps are + attempted by this processor on this timestep. No other eligible + 4-tuples of atoms are considered. This means that each processor + will perform either a single swap or none on timesteps this fix is + invoked. ---------- From 526cc1566eac8b531e3b3f696804531a44977ccb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 13 Sep 2022 06:13:37 -0400 Subject: [PATCH 6/7] cosmetic --- src/MC/fix_bond_swap.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp index ea31270b92..a1f1c2b307 100644 --- a/src/MC/fix_bond_swap.cpp +++ b/src/MC/fix_bond_swap.cpp @@ -114,10 +114,10 @@ FixBondSwap::~FixBondSwap() // delete temperature if fix created it if (tflag) modify->delete_compute(id_temp); - delete [] id_temp; + delete[] id_temp; memory->destroy(alist); - delete [] permute; + delete[] permute; } /* ---------------------------------------------------------------------- */ @@ -244,14 +244,14 @@ void FixBondSwap::post_integrate() memory->create(alist,nmax,"bondswap:alist"); } - // use randomized permulation of both I and J atoms in double loop below + // use randomized permutation of both I and J atoms in double loop below // this is to avoid any bias in accepted MC swaps based on // ordering LAMMPS creates on a processor for atoms or their neighbors // create a random permutation of list of Neligible atoms // uses one-pass Fisher-Yates shuffle on an initial identity permutation // output: randomized alist[] vector, used in outer loop to select an I atom - // similar randomized permulation is created for neighbors of each I atom + // similar randomized permutation is created for neighbors of each I atom int neligible = 0; for (ii = 0; ii < inum; ii++) { @@ -693,7 +693,7 @@ int FixBondSwap::modify_param(int narg, char **arg) modify->delete_compute(id_temp); tflag = 0; } - delete [] id_temp; + delete[] id_temp; id_temp = utils::strdup(arg[1]); int icompute = modify->find_compute(id_temp); @@ -767,7 +767,7 @@ void FixBondSwap::neighbor_permutation(int n) int i,j,tmp; if (n > maxpermute) { - delete [] permute; + delete[] permute; maxpermute = n + DELTA_PERMUTE; permute = new int[maxpermute]; } From 3770c02752bddfc26ee8c1385c84867f58781a02 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 15 Sep 2022 09:22:20 -0400 Subject: [PATCH 7/7] fix typo --- src/fix_pair.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fix_pair.cpp b/src/fix_pair.cpp index 6f98ad2790..5f9363daaa 100644 --- a/src/fix_pair.cpp +++ b/src/fix_pair.cpp @@ -70,7 +70,7 @@ FixPair::FixPair(LAMMPS *lmp, int narg, char **arg) : for (int ifield = 0; ifield < nfield; ifield++) { if (trigger[ifield] == 0) triggername[ifield] = nullptr; - else triggername[nfield] = utils::strdup(fmt::format("{}_flag", fieldname[ifield])); + else triggername[ifield] = utils::strdup(fmt::format("{}_flag", fieldname[ifield])); } // extract all fields just to get number of per-atom values