ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2888 from lammps/sticker-bonds

Browse Source 
Sticker bonds
This commit is contained in:
Axel Kohlmeyer
2021-08-20 14:28:57 -04:00
committed by GitHub
parent 3404920e98 fcc6cc628e
commit 1eb506002b
2 changed files with 138 additions and 76 deletions
Download Patch File Download Diff File Expand all files Collapse all files

171
doc/src/fix_bond_swap.rst
View File

@ -27,21 +27,26 @@ Examples
Description
"""""""""""
In a simulation of polymer chains, this command attempts to swap bonds
between two different chains, effectively grafting the end of one
chain onto another chain and vice versa. This is done via Monte Carlo
rules using the Boltzmann acceptance criterion. The purpose is to
equilibrate the polymer chain conformations more rapidly than dynamics
alone would do it, by enabling instantaneous large conformational
changes in a dense polymer melt. The polymer chains should thus more
rapidly converge to the proper end-to-end distances and radii of
gyration. It is designed for use with systems of
:doc:`FENE <bond_fene>` or :doc:`harmonic <bond_harmonic>` bead-spring
polymer chains where each polymer is a linear chain of monomers, but
LAMMPS does not enforce this requirement, i.e. any
:doc:`bond_style <bond_style>` can be used.
In a simulation of polymer chains this command attempts to swap a pair
of bonds, as illustrated below. This is done via Monte Carlo rules
using the Boltzmann acceptance criterion, typically with the goal of
equilibrating the polymer system more quickly. This fix is designed
for use with idealized bead-spring polymer chains where each polymer
is a linear chain of monomers, but LAMMPS does not check that is the
case for your system.
A schematic of the kinds of bond swaps that can occur is shown here:
Here are two use cases for this fix.
The first use case is for swapping bonds on two different chains,
effectively grafting the end of one chain onto the other chain and
vice versa. The purpose is to equilibrate the polymer chain
conformations more rapidly than dynamics alone would do it, by
enabling instantaneous large conformational changes in a dense polymer
melt. The polymer chains should thus more rapidly converge to the
proper end-to-end distances and radii of gyration.
A schematic of the kinds of bond swaps that can occur in this use case
is shown here:
.. image:: JPG/bondswap.jpg
:align: center
@ -53,38 +58,76 @@ attempt to delete the A1-A2 and B1-B2 bonds and replace them with
A1-B2 and B1-A2 bonds. If the swap is energetically favorable, the
two chains on the right are the result and each polymer chain has
undergone a dramatic conformational change. This reference,
:ref:`(Sides) <Sides>` provides more details on how the algorithm works and
its application:
:ref:`(Sides) <Sides>` provides more details on the algorithm's
effectiveness for this use case.
The bond swapping operation is invoked every *Nevery* timesteps. If
any bond is swapped, a re-build of the neighbor lists is triggered,
since a swap alters the list of which neighbors are considered for
pairwise interaction. At each invocation, each processor considers a
random specified *fraction* of its atoms as potential swapping
monomers for this timestep. Choosing a small *fraction* value can
reduce the likelihood of a reverse swap occurring soon after an
initial swap.
The second use case is a collection of polymer chains with some
fraction of their sites identified as "sticker" sites. Initially each
polymer chain is isolated from the others in a topological sense, and
there is an intra-chain bond between every pair of sticker sites on
the same chain. Over time, bonds swap so that inter-moleculer sticker
bonds are created. This models a vitrification-style process whereby
the polymer chains all become interconnected. For this use case, if
angles are defined they should not include bonds between sticker
sites.
For each monomer A1, its neighbors are examined to find a possible B1
monomer. Both A1 and B1 must be in the fix group, their separation
must be less than the specified *cutoff*, and the molecule IDs of A1
and B1 must be the same (see below). If a suitable partner is found,
the energy change due to swapping the 2 bonds is computed. This
includes changes in pairwise, bond, and angle energies due to the
altered connectivity of the 2 chains. Dihedral and improper
interactions are not allowed to be defined when this fix is used.
----------
The bond swapping operation is invoked once every *Nevery* timesteps.
If any bond in the entire system is swapped, a re-build of the
neighbor lists is triggered, since a swap alters the list of which
neighbors are considered for pairwise interaction. At each
invocation, each processor considers a random specified *fraction* of
its atoms as potential swapping monomers for this timestep. Choosing
a small *fraction* value can reduce the likelihood of a reverse swap
occurring soon after an initial swap.
For each monomer A1, its neighbors are looped over as B1 monomers.
For each A1,B1 an additional double loop of bond partners A2 of A1,
and bond partners B2 of B1 a is performed. For each pair of A1-A2 and
B1-B2 bonds to be eligible for swapping, the following 4 criteria must
be met:
(1) All 4 monomers must be in the fix group.
(2) All 4 monomers must be owned by the processor (not ghost atoms).
This insures that another processor does not attempt to swap bonds
involving the same atoms on the same timestep. Note that this also
means that bond pairs which straddle processor boundaries are not
eligible for swapping on this step.
(3) The distances between 4 pairs of atoms -- (A1,A2), (B1,B2),
(A1,B2), (B1,A2) -- must all be less thann the specified *cutoff*\ .
(4) The molecule IDs of A1 and B1 must be the same (see below).
If an eligible B1 partner is found, the energy change due to swapping
the 2 bonds is computed. This includes changes in pairwise, bond, and
angle energies due to the altered connectivity of the 2 chains.
Dihedral and improper interactions are not allowed to be defined when
this fix is used.
If the energy decreases due to the swap operation, the bond swap is
accepted. If the energy increases it is accepted with probability
exp(-delta/kT) where delta is the increase in energy, k is the
Boltzmann constant, and T is the current temperature of the system.
Whether the swap is accepted or rejected, no other swaps are attempted
by this processor on this timestep.
The criterion for matching molecule IDs is how bond swaps performed by
this fix conserve chain length. To use this features you must setup
the molecule IDs for your polymer chains in a certain way, typically
in the data file, read by the :doc:`read_data <read_data>` command.
.. note::
IMPORTANT: Whether the swap is accepted or rejected, no other swaps
are attempted by this processor on this timestep. No other
eliglble 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.
----------
The criterion for matching molecule IDs is how the first use case
described above can be simulated while conserving chain lengths. This
is done by setting up the molecule IDs for the polymer chains in a
specific way, typically in the data file, read by the :doc:`read_data
<read_data>` command.
Consider a system of 6-mer chains. You have 2 choices. If the
molecule IDs for monomers on each chain are set to 1,2,3,4,5,6 then
swaps will conserve chain length. For a particular monomer there will
@ -113,6 +156,9 @@ ends of a chain swap with each other.
running dynamics, but can affect calculation of some diagnostic
quantities or the printing of unwrapped coordinates to a dump file.
For the second use case described above, the molecule IDs for all
sticker sites should be the same.
----------
This fix computes a temperature each time it is invoked for use by the
@ -129,27 +175,28 @@ appended and the group for the new compute is "all", so that the
temperature of the entire system is used.
Note that this is NOT the compute used by thermodynamic output (see
the :doc:`thermo_style <thermo_style>` command) with ID = *thermo_temp*.
This means you can change the attributes of this fix's temperature
(e.g. its degrees-of-freedom) via the
:doc:`compute_modify <compute_modify>` command or print this temperature
during thermodynamic output via the :doc:`thermo_style custom <thermo_style>` command using the appropriate compute-ID.
It also means that changing attributes of *thermo_temp* will have no
effect on this fix.
the :doc:`thermo_style <thermo_style>` command) with ID =
*thermo_temp*. This means you can change the attributes of this fix's
temperature (e.g. its degrees-of-freedom) via the :doc:`compute_modify
<compute_modify>` command or print this temperature during
thermodynamic output via the :doc:`thermo_style custom <thermo_style>`
command using the appropriate compute-ID. It also means that changing
attributes of *thermo_temp* will have no effect on this fix.
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`. Because the state of the random number generator
is not saved in restart files, this means you cannot do "exact"
restarts with this fix, where the simulation continues on the same as
if no restart had taken place. However, in a statistical sense, a
restarted simulation should produce the same behavior. Also note that
each processor generates possible swaps independently of other
processors. Thus if you repeat the same simulation on a different number
of processors, the specific swaps performed will be different.
No information about this fix is written to :doc:`binary restart files
<restart>`. Because the state of the random number generator is not
saved in restart files, this means you cannot do "exact" restarts with
this fix, where the simulation continues on the same as if no restart
had taken place. However, in a statistical sense, a restarted
simulation should produce the same behavior. Also note that each
processor generates possible swaps independently of other processors.
Thus if you repeat the same simulation on a different number of
processors, the specific swaps performed will be different.
The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
fix. You can use it to assign a :doc:`compute <compute>` you have
@ -157,16 +204,18 @@ defined to this fix which will be used to compute the temperature for
the Boltzmann criterion.
This fix computes two statistical quantities as a global 2-vector of
output, which can be accessed by various :doc:`output commands <Howto_output>`. The first component of the vector is the
cumulative number of swaps performed by all processors. The second
component of the vector is the cumulative number of swaps attempted
(whether accepted or rejected). Note that a swap "attempt" only
occurs when swap partners meeting the criteria described above are
found on a particular timestep. The vector values calculated by this
fix are "intensive".
output, which can be accessed by various :doc:`output commands
<Howto_output>`. The first component of the vector is the cumulative
number of swaps performed by all processors. The second component of
the vector is the cumulative number of swaps attempted (whether
accepted or rejected). Note that a swap "attempt" only occurs when
swap partners meeting the criteria described above are found on a
particular timestep. The vector values calculated by this fix are
"intensive".
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minimization <minimize>`.
the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""

43
src/MC/fix_bond_swap.cpp
View File

@ -148,7 +148,8 @@ void FixBondSwap::init()
error->all(FLERR,"Pair style does not support fix bond/swap");
if (force->angle == nullptr && atom->nangles > 0 && comm->me == 0)
error->warning(FLERR,"Fix bond/swap will ignore defined angles");
error->warning(FLERR,"Fix bond/swap will not preserve correct angle "
"topology because no angle_style is defined");
if (force->dihedral || force->improper)
error->all(FLERR,"Fix bond/swap cannot use dihedral or improper styles");
@ -255,12 +256,18 @@ 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
// atom j must be on-processor (j < nlocal)
// atom j must be in fix group
// i and j must be same distance from chain end (mol[i] = mol[j])
// NOTE: must use extra parens in if test on mask[j] & groupbit
// atom I is randomly selected via atom list
// look at all J neighbors of atom I
// J must be on-processor (J < nlocal)
// I,J must be in fix group
// I,J must have same molecule IDs
// use case 1 (see doc page):
// if user defines mol IDs appropriately for linear chains,
// this will mean they are same distance from (either) chain end
// use case 2 (see doc page):
// if user defines a unique mol ID for desired bond sites (on any chain)
// and defines the fix group as these sites,
// this will mean they are eligible bond sites
int ntest = static_cast<int> (fraction * neligible);
int accept = 0;
@ -277,23 +284,29 @@ void FixBondSwap::post_integrate()
if ((mask[j] & groupbit) == 0) continue;
if (molecule[i] != molecule[j]) continue;
// look at all bond partners of atoms i and j
// use num_bond for this, not special list, so also find bondtypes
// inext,jnext = bonded atoms
// loop over all bond partners of atoms I and J
// use num_bond for this, not special list, so also have bondtypes
// inext,jnext = atoms bonded to I,J
// inext,jnext must be on-processor (inext,jnext < nlocal)
// inext,jnext must be same dist from chain end (mol[inext] = mol[jnext])
// since swaps may occur between two ends of a single chain, insure
// the 4 atoms are unique (no duplicates): inext != jnext, inext != j
// inext,jnext must be in fix group
// inext,jnext must have same molecule IDs
// in use cases above ...
// for case 1: this insures chain length is preserved
// for case 2: always satisfied b/c fix group = bond-able atoms
// 4 atoms must be unique (no duplicates): inext != jnext, inext != j
// already know i != inext, j != jnext
// all 4 old and new bonds must have length < cutoff
for (ibond = 0; ibond < num_bond[i]; ibond++) {
inext = atom->map(bond_atom[i][ibond]);
if (inext >= nlocal || inext < 0) continue;
if ((mask[inext] & groupbit) == 0) continue;
ibondtype = bond_type[i][ibond];
for (jbond = 0; jbond < num_bond[j]; jbond++) {
jnext = atom->map(bond_atom[j][jbond]);
if (jnext >= nlocal || jnext < 0) continue;
if ((mask[jnext] & groupbit) == 0) continue;
jbondtype = bond_type[j][jbond];
if (molecule[inext] != molecule[jnext]) continue;
@ -306,7 +319,7 @@ void FixBondSwap::post_integrate()
// if angles are enabled:
// find other atoms i,inext,j,jnext are in angles with
// and angletypes: i/j angletype, i/j nextangletype
// use num_angle for this, not special list, so also find angletypes
// use num_angle for this, not special list, so also have angletypes
// 4 atoms consecutively along 1st chain: iprev,i,inext,ilast
// 4 atoms consecutively along 2nd chain: jprev,j,jnext,jlast
// prev or last atom can be non-existent at end of chain
@ -428,7 +441,7 @@ void FixBondSwap::post_integrate()
done:
// trigger immediate reneighboring if any swaps occurred
// trigger immediate reneighboring if swaps occurred on one or more procs
int accept_any;
MPI_Allreduce(&accept,&accept_any,1,MPI_INT,MPI_SUM,world);