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

Merge remote branch 'lammps-ro/master' into lammps-icms

Browse Source 
Resolved Conflicts:
	doc/Manual.html
	doc/Manual.txt
This commit is contained in:
Axel Kohlmeyer
2014-04-12 16:13:47 +02:00
parent 3ec9862f6f 5a26ea8238
commit c6992e8503
19 changed files with 764 additions and 151 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
doc/Manual.html
View File

@ -1,7 +1,7 @@
<HTML>
<HEAD>
<TITLE>LAMMPS-ICMS Users Manual</TITLE>
<META NAME="docnumber" CONTENT="11 Apr 2014 version">
<META NAME="docnumber" CONTENT="12 Apr 2014 version">
<META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
<META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation. This software and manual is distributed under the GNU General Public License.">
</HEAD>
@ -22,7 +22,7 @@
<CENTER><H3>LAMMPS-ICMS Documentation
</H3></CENTER>
<CENTER><H4>11 Apr 2014 version
<CENTER><H4>12 Apr 2014 version
</H4></CENTER>
<H4>Version info:
</H4>

4
doc/Manual.txt
View File

@ -1,7 +1,7 @@
<HEAD>
<TITLE>LAMMPS-ICMS Users Manual</TITLE>
<TITLE>LAMMPS Users Manual</TITLE>
<META NAME="docnumber" CONTENT="11 Apr 2014 version">
<META NAME="docnumber" CONTENT="12 Apr 2014 version">
<META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
<META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation. This software and manual is distributed under the GNU General Public License.">
</HEAD>
@ -19,7 +19,7 @@
<H1></H1>
LAMMPS-ICMS Documentation :c,h3
11 Apr 2014 version :c,h4
12 Apr 2014 version :c,h4
Version info: :h4

54
doc/fix_bond_break.html
View File

@ -47,12 +47,14 @@ specified criteria. This can be used to model the dissolution of a
polymer network due to stretching of the simulation box or other
deformations. In this context, a bond means an interaction between a
pair of atoms computed by the <A HREF = "bond_style.html">bond_style</A> command.
Once the bond is broken it will be permanently deleted. This is
different than a <A HREF = "pair_style.html">pairwise</A> bond-order potential such
as Tersoff or AIREBO which infers bonds and many-body interactions
based on the current geometry of a small cluster of atoms and
effectively creates and destroys bonds from timestep to timestep as
atoms move.
Once the bond is broken it will be permanently deleted, as will all
angle, dihedral, and improper interactions that bond is part of.
</P>
<P>This is different than a <A HREF = "pair_style.html">pairwise</A> bond-order
potential such as Tersoff or AIREBO which infers bonds and many-body
interactions based on the current geometry of a small cluster of atoms
and effectively creates and destroys bonds and higher-order many-body
interactions from timestep to timestep as atoms move.
</P>
<P>A check for possible bond breakage is performed every <I>Nevery</I>
timesteps. If two bonded atoms I,J are further than a distance <I>Rmax</I>
@ -80,17 +82,21 @@ A uniform random number between 0.0 and 1.0 is generated and the
eligible bond is only broken if the random number < fraction.
</P>
<P>When a bond is broken, data structures within LAMMPS that store bond
topology are updated to reflect the breakage. This can also affect
subsequent computation of pairwise interactions involving the atoms in
the bond. See the Restriction section below for additional
information.
topology are updated to reflect the breakage. Likewise, if the bond
is part of a 3-body (angle) or 4-body (dihedral, improper)
interaction, that interaction is removed as well.
</P>
<P>Computationally, each timestep this fix operates, it loops over bond
lists and computes distances between pairs of bonded atoms in the
list. It also communicates between neighboring processors to
coordinate which bonds are broken. Thus it will increase the cost of
a timestep. Thus you should be cautious about invoking this fix too
frequently.
<P>Computationally, each timestep this fix operates, it loops over all
the bonds in the system and computes distances between pairs of bonded
atoms. It also communicates between neighboring processors to
coordinate which bonds are broken. Moreover, if any bonds are broken,
neighbor lists must be immediately updated on the same timestep. This
is to insure that any pairwise interactions that should be turned "on"
due to a bond breaking, because they are no longer excluded by the
presence of the bond and the settings of the
<A HREF = "special_bonds.html">special_bonds</A> command, will be immediately
recognized. All of these operations increase the cost of a timestep.
Thus you should be cautious about invoking this fix too frequently.
</P>
<P>You can dump out snapshots of the current bond topology via the <A HREF = "dump.html">dump
local</A> command.
@ -132,22 +138,6 @@ minimization</A>.
built with that package. See the <A HREF = "Section_start.html#start_3">Making
LAMMPS</A> section for more info.
</P>
<P>Currently, there are 2 restrictions for using this fix. We may relax
these in the future if there are new models that would be enabled by
it.
</P>
<P>When a bond is broken, you might wish to turn off angle and dihedral
interactions that include that bond. However, LAMMPS does not check
for these angles and dihedrals, even if your simulation defines an
<A HREF = "angle_style.html">angle_style</A> or
<A HREF = "dihedral_style.html">dihedral_style</A>.
</P>
<P>This fix requires that the pairwise weightings defined by the
<A HREF = "special_bonds.html">special_bonds</A> command be *,1,1 for 1-3 and 1-4
neighbors within the bond topology (the 1-2 setting is not
restricted). This means that the pairwise interaction of I with J's
other bond partners is unaffected by the breaking of the bond.
</P>
<P><B>Related commands:</B>
</P>
<P><A HREF = "fix_bond_create.html">fix bond/create</A>, <A HREF = "fix_bond_swap.html">fix

54
doc/fix_bond_break.txt
View File

@ -36,12 +36,14 @@ specified criteria. This can be used to model the dissolution of a
polymer network due to stretching of the simulation box or other
deformations. In this context, a bond means an interaction between a
pair of atoms computed by the "bond_style"_bond_style.html command.
Once the bond is broken it will be permanently deleted. This is
different than a "pairwise"_pair_style.html bond-order potential such
as Tersoff or AIREBO which infers bonds and many-body interactions
based on the current geometry of a small cluster of atoms and
effectively creates and destroys bonds from timestep to timestep as
atoms move.
Once the bond is broken it will be permanently deleted, as will all
angle, dihedral, and improper interactions that bond is part of.
This is different than a "pairwise"_pair_style.html bond-order
potential such as Tersoff or AIREBO which infers bonds and many-body
interactions based on the current geometry of a small cluster of atoms
and effectively creates and destroys bonds and higher-order many-body
interactions from timestep to timestep as atoms move.
A check for possible bond breakage is performed every {Nevery}
timesteps. If two bonded atoms I,J are further than a distance {Rmax}
@ -69,17 +71,21 @@ A uniform random number between 0.0 and 1.0 is generated and the
eligible bond is only broken if the random number < fraction.
When a bond is broken, data structures within LAMMPS that store bond
topology are updated to reflect the breakage. This can also affect
subsequent computation of pairwise interactions involving the atoms in
the bond. See the Restriction section below for additional
information.
topology are updated to reflect the breakage. Likewise, if the bond
is part of a 3-body (angle) or 4-body (dihedral, improper)
interaction, that interaction is removed as well.
Computationally, each timestep this fix operates, it loops over bond
lists and computes distances between pairs of bonded atoms in the
list. It also communicates between neighboring processors to
coordinate which bonds are broken. Thus it will increase the cost of
a timestep. Thus you should be cautious about invoking this fix too
frequently.
Computationally, each timestep this fix operates, it loops over all
the bonds in the system and computes distances between pairs of bonded
atoms. It also communicates between neighboring processors to
coordinate which bonds are broken. Moreover, if any bonds are broken,
neighbor lists must be immediately updated on the same timestep. This
is to insure that any pairwise interactions that should be turned "on"
due to a bond breaking, because they are no longer excluded by the
presence of the bond and the settings of the
"special_bonds"_special_bonds.html command, will be immediately
recognized. All of these operations increase the cost of a timestep.
Thus you should be cautious about invoking this fix too frequently.
You can dump out snapshots of the current bond topology via the "dump
local"_dump.html command.
@ -121,22 +127,6 @@ This fix is part of the MC package. It is only enabled if LAMMPS was
built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
Currently, there are 2 restrictions for using this fix. We may relax
these in the future if there are new models that would be enabled by
it.
When a bond is broken, you might wish to turn off angle and dihedral
interactions that include that bond. However, LAMMPS does not check
for these angles and dihedrals, even if your simulation defines an
"angle_style"_angle_style.html or
"dihedral_style"_dihedral_style.html.
This fix requires that the pairwise weightings defined by the
"special_bonds"_special_bonds.html command be *,1,1 for 1-3 and 1-4
neighbors within the bond topology (the 1-2 setting is not
restricted). This means that the pairwise interaction of I with J's
other bond partners is unaffected by the breaking of the bond.
[Related commands:]
"fix bond/create"_fix_bond_create.html, "fix

26
doc/fix_bond_create.html
View File

@ -117,22 +117,20 @@ structures that store this information must have space for it. When
LAMMPS is initialized from a data file, the list of bonds is scanned
and the maximum number of bonds per atom is tallied. If some atom
will acquire more bonds than this limit as this fix operates, then the
"extra bonds per atom" parameter in the data file header must be set
to allow for it. See the <A HREF = "read_data.html">read_data</A> command for more
details. Note that if this parameter needs to be set, it means a data
file must be used to initialize the system, even if it initially has
no bonds. A data file with no atoms can be used if you wish to add
unbonded atoms via the <A HREF = "create_atoms.html">create atoms</A> command,
e.g. for a percolation simulation.
"extra bonds per atom" parameter must be set to allow for it. See the
<A HREF = "read_data.html">read_data</A> or <A HREF = "create_box.html">create_box</A> command for
more details. Note that a data file with no atoms can be used if you
wish to add unbonded atoms via the <A HREF = "create_atoms.html">create atoms</A>
command, e.g. for a percolation simulation.
</P>
<P>IMPORTANT NOTE: LAMMPS also maintains a data structure that stores a
list of 1st, 2nd, and 3rd neighbors of each atom (in the bond topology
of the system) for use in weighting pairwise interactions for bonded
atoms. Adding a bond adds a single entry to this list. The "extra"
keyword of the <A HREF = "special_bonds.html">special_bonds</A> command should be
used to leave space for new bonds if the maximum number of entries for
any atom will be exceeded as this fix operates. See the
<A HREF = "special_bonds.html">special_bonds</A> command for details.
list of 1st, 2nd, and 3rd neighbors of each atom (within the bond
topology of the system) for use in weighting pairwise interactions for
bonded atoms. Adding a bond adds a single entry to this list. The
"extra special per atom" parameter must typically be set to allow for
it. There are 3 ways to do this. See the <A HREF = "read_data.html">read_data</A>
or <A HREF = "create_box.html">create_box</A> or "special_bonds extra" commands for
details.
</P>
<P>Note that even if your simulation starts with no bonds, you must
define a <A HREF = "bond_style.html">bond_style</A> and use the

26
doc/fix_bond_create.txt
View File

@ -105,22 +105,20 @@ structures that store this information must have space for it. When
LAMMPS is initialized from a data file, the list of bonds is scanned
and the maximum number of bonds per atom is tallied. If some atom
will acquire more bonds than this limit as this fix operates, then the
"extra bonds per atom" parameter in the data file header must be set
to allow for it. See the "read_data"_read_data.html command for more
details. Note that if this parameter needs to be set, it means a data
file must be used to initialize the system, even if it initially has
no bonds. A data file with no atoms can be used if you wish to add
unbonded atoms via the "create atoms"_create_atoms.html command,
e.g. for a percolation simulation.
"extra bonds per atom" parameter must be set to allow for it. See the
"read_data"_read_data.html or "create_box"_create_box.html command for
more details. Note that a data file with no atoms can be used if you
wish to add unbonded atoms via the "create atoms"_create_atoms.html
command, e.g. for a percolation simulation.
IMPORTANT NOTE: LAMMPS also maintains a data structure that stores a
list of 1st, 2nd, and 3rd neighbors of each atom (in the bond topology
of the system) for use in weighting pairwise interactions for bonded
atoms. Adding a bond adds a single entry to this list. The "extra"
keyword of the "special_bonds"_special_bonds.html command should be
used to leave space for new bonds if the maximum number of entries for
any atom will be exceeded as this fix operates. See the
"special_bonds"_special_bonds.html command for details.
list of 1st, 2nd, and 3rd neighbors of each atom (within the bond
topology of the system) for use in weighting pairwise interactions for
bonded atoms. Adding a bond adds a single entry to this list. The
"extra special per atom" parameter must typically be set to allow for
it. There are 3 ways to do this. See the "read_data"_read_data.html
or "create_box"_create_box.html or "special_bonds extra" commands for
details.
Note that even if your simulation starts with no bonds, you must
define a "bond_style"_bond_style.html and use the

2
doc/thermo_style.html
View File

@ -29,6 +29,7 @@
emol, elong, etail,
vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi,
xy, xz, yz, xlat, ylat, zlat,
bonds, angles, dihedrals, impropers,
pxx, pyy, pzz, pxy, pxz, pyz,
fmax, fnorm,
cella, cellb, cellc, cellalpha, cellbeta, cellgamma,
@ -68,6 +69,7 @@
xlo,xhi,ylo,yhi,zlo,zhi = box boundaries
xy,xz,yz = box tilt for triclinic (non-orthogonal) simulation boxes
xlat,ylat,zlat = lattice spacings as calculated by <A HREF = "lattice.html">lattice</A> command
bonds,angles,dihedrals,impropers = # of these interactions defined
pxx,pyy,pzz,pxy,pxz,pyz = 6 components of pressure tensor
fmax = max component of force on any atom in any dimension
fnorm = length of force vector for all atoms

2
doc/thermo_style.txt
View File

@ -24,6 +24,7 @@ args = list of arguments for a particular style :l
emol, elong, etail,
vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi,
xy, xz, yz, xlat, ylat, zlat,
bonds, angles, dihedrals, impropers,
pxx, pyy, pzz, pxy, pxz, pyz,
fmax, fnorm,
cella, cellb, cellc, cellalpha, cellbeta, cellgamma,
@ -63,6 +64,7 @@ args = list of arguments for a particular style :l
xlo,xhi,ylo,yhi,zlo,zhi = box boundaries
xy,xz,yz = box tilt for triclinic (non-orthogonal) simulation boxes
xlat,ylat,zlat = lattice spacings as calculated by "lattice"_lattice.html command
bonds,angles,dihedrals,impropers = # of these interactions defined
pxx,pyy,pzz,pxy,pxz,pyz = 6 components of pressure tensor
fmax = max component of force on any atom in any dimension
fnorm = length of force vector for all atoms

557
src/MC/fix_bond_break.cpp
View File

@ -19,6 +19,7 @@
#include "update.h"
#include "respa.h"
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
@ -30,6 +31,8 @@
using namespace LAMMPS_NS;
using namespace FixConst;
#define DELTA 16
/* ---------------------------------------------------------------------- */
FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) :
@ -38,6 +41,7 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) :
if (narg < 6) error->all(FLERR,"Illegal fix bond/break command");
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
nevery = force->inumeric(FLERR,arg[3]);
if (nevery <= 0) error->all(FLERR,"Illegal fix bond/break command");
@ -50,7 +54,7 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) :
extvector = 0;
btype = force->inumeric(FLERR,arg[4]);
double cutoff = force->numeric(FLERR,arg[5]);
cutoff = force->numeric(FLERR,arg[5]);
if (btype < 1 || btype > atom->nbondtypes)
error->all(FLERR,"Invalid bond type in fix bond/break command");
@ -86,8 +90,9 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) :
random = new RanMars(lmp,seed + me);
// set comm sizes needed by this fix
// forward is big due to comm of 1-2 neighbors
comm_forward = 2;
comm_forward = MAX(2,1+atom->maxspecial);
comm_reverse = 2;
// allocate arrays local to this fix
@ -96,6 +101,22 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) :
partner = NULL;
distsq = NULL;
maxbreak = maxbreakall = 0;
broken = brokenall = NULL;
inbuf = NULL;
maxinfluenced = 0;
influenced = NULL;
recvcounts = NULL;
displs = NULL;
// copy = special list for one atom
// may contain 1-2 neighs of all 1-3 neighs before dedup() shrinks it
int maxspecial = atom->maxspecial;
copy = new tagint[maxspecial*maxspecial + maxspecial];
// zero out stats
breakcount = 0;
@ -112,6 +133,13 @@ FixBondBreak::~FixBondBreak()
memory->destroy(partner);
memory->destroy(distsq);
memory->destroy(broken);
memory->destroy(brokenall);
memory->destroy(inbuf);
memory->destroy(influenced);
memory->destroy(recvcounts);
memory->destroy(displs);
delete [] copy;
}
/* ---------------------------------------------------------------------- */
@ -128,33 +156,36 @@ int FixBondBreak::setmask()
void FixBondBreak::init()
{
// require special bonds = *,1,1
// [0] can be anything b/c I,J are removed from each other's special list
// [1],[2] must be 1.0 b/c only special lists of I,J are updated when
// bond I-J is broken, not special lists of neighbors of I,J,etc
int flag = 0;
if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) flag = 1;
if (force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) flag = 1;
if (flag) error->all(FLERR,"Fix bond/break requires special_bonds = *,1,1");
// warn if angles, dihedrals, impropers are being used
if (force->angle || force->dihedral || force->improper) {
if (me == 0)
error->warning(FLERR,"Broken bonds will not alter angles, "
"dihedrals, or impropers");
}
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
// commextent = 3*bondcutoff
// use 3 b/c 4 atom in 1-2-3-4 needs to know 1-2 bond has broken
// and that info could be known by atom 1 with
// atom 1 on one proc and atoms 2,3,4 on the other proc
// using "cutoff" as bond length is a guesstimate of what's OK
// see NOTE below for possible issues with this
commextent = 3.0*cutoff;
// improper class2 and ring styles not allowed for now
// due to different ordering of improper topology (not atom I centric)
if (force->improper) {
if (force->improper_match("class2") || force->improper_match("ring"))
error->all(FLERR,"Cannot yet use fix bond/break with this "
"improper style");
}
// DEBUG
// print_bb();
}
/* ---------------------------------------------------------------------- */
void FixBondBreak::post_integrate()
{
int i,j,k,m,n,i1,i2,n1,n3,type;
int i,j,k,m,n,i1,i2,n1,n2,n3,type;
double delx,dely,delz,rsq;
tagint *slist;
@ -231,6 +262,7 @@ void FixBondBreak::post_integrate()
if (partner[i]) probability[i] = random->uniform();
}
commflag = 1;
comm->forward_comm_fix(this);
// break bonds
@ -243,7 +275,7 @@ void FixBondBreak::post_integrate()
int **nspecial = atom->nspecial;
tagint **special = atom->special;
int nbreak = 0;
nbreak = 0;
for (i = 0; i < nlocal; i++) {
if (partner[i] == 0) continue;
j = atom->map(partner[i]);
@ -278,17 +310,26 @@ void FixBondBreak::post_integrate()
slist = special[i];
n1 = nspecial[i][0];
n3 = nspecial[i][2];
for (m = 0; m < n1; m++)
if (slist[m] == partner[i]) break;
n3 = nspecial[i][2];
for (; m < n3-1; m++) slist[m] = slist[m+1];
nspecial[i][0]--;
nspecial[i][1]--;
nspecial[i][2]--;
// count the broken bond once
// count the broken bond once and store in broken list
if (tag[i] < tag[j]) nbreak++;
if (tag[i] < tag[j]) {
if (nbreak == maxbreak) {
maxbreak += DELTA;
memory->grow(broken,maxbreak,2,"bond/break:broken");
memory->grow(inbuf,2*maxbreak,"bond/break:inbuf");
}
broken[nbreak][0] = tag[i];
broken[nbreak][1] = tag[j];
nbreak++;
}
}
// tally stats
@ -297,9 +338,370 @@ void FixBondBreak::post_integrate()
breakcounttotal += breakcount;
atom->nbonds -= breakcount;
// trigger reneighboring if any bonds were formed
// trigger reneighboring if any bonds were broken
// this insures neigh lists will immediately reflect the topology changes
// done if no bonds broken
if (breakcount) next_reneighbor = update->ntimestep;
if (!breakcount) return;
// communicate broken bonds to procs that need them
// local comm via exchange_variable() if commextent < prob sub-domain
// else global comm via MPI_Allgatherv()
// NOTE: not fully happy with this test,
// but want to avoid Allgather if at all possible
// issue is there is no simple way to guarantee that local comm is OK
// what is really needed is every 4 atom in 1-2-3-4 knows about 1-2 bond
// this is hard to predict dynamically b/c current lengths of 2-3,3-4
// bonds are not known
// could use 3*(domain->maxbondall+BONDSTRETCH) for better estimate?
// am not doing the local comm via ghost atoms, so ghost cutoff
// is irrelevant
// this test is also for orthogonal boxes and equi-partitioning
double *outbuf;
int local = 1;
if (domain->xprd/comm->procgrid[0] < commextent) local = 0;
if (domain->yprd/comm->procgrid[1] < commextent) local = 0;
if (domain->dimension == 3 && domain->zprd/comm->procgrid[2] < commextent)
local = 0;
if (local) {
m = 0;
for (i = 0; i < nbreak; i++) {
inbuf[m++] = broken[i][0];
inbuf[m++] = broken[i][1];
}
int ntotal = comm->exchange_variable(2*nbreak,inbuf,outbuf);
nbreakall = ntotal/2;
} else MPI_Allreduce(&nbreak,&nbreakall,1,MPI_INT,MPI_SUM,world);
if (nbreakall > maxbreakall) {
while (maxbreakall < nbreakall) maxbreakall += DELTA;
memory->destroy(brokenall);
memory->grow(brokenall,maxbreakall,2,"bond/break:brokenall");
}
if (local) {
m = 0;
for (i = 0; i < nbreakall; i++) {
brokenall[i][0] = static_cast<int> (outbuf[m++]);
brokenall[i][1] = static_cast<int> (outbuf[m++]);
}
} else {
if (!recvcounts) {
memory->create(recvcounts,nprocs,"bond/break:recvcounts");
memory->create(displs,nprocs,"bond/break:displs");
}
int n = 2*nbreak;
MPI_Allgather(&n,1,MPI_INT,recvcounts,1,MPI_INT,world);
displs[0] = 0;
for (int iproc = 1; iproc < nprocs; iproc++)
displs[iproc] = displs[iproc-1] + recvcounts[iproc-1];
int *ptr = NULL;
if (nbreak) ptr = &broken[0][0];
MPI_Allgatherv(ptr,2*nbreak,MPI_INT,
&brokenall[0][0],recvcounts,displs,MPI_INT,world);
}
// communicate 1-2 special neighs of ghost atoms
// 1-2 neighs already reflect broken bonds
commflag = 2;
comm->forward_comm_variable_fix(this);
// update special neigh lists of all atoms affected by any broken bond
// also remove angles/dihedrals/impropers broken by broken bonds
update_topology();
// DEBUG
// print_bb();
}
/* ----------------------------------------------------------------------
double loop over my atoms and broken bonds
influenced = 1 if atom's topology is affected by any broken bond
yes if is one of 2 atoms in bond
yes if both atom IDs appear in atom's special list
else no
if influenced:
check for angles/dihedrals/impropers to break due to broken bond
rebuild the atom's special list of 1-2,1-3,1-4 neighs
------------------------------------------------------------------------- */
void FixBondBreak::update_topology()
{
int i,j,k,n,influence,found;
tagint id1,id2;
tagint *slist;
int angles_allow = atom->avec->angles_allow;
int dihedrals_allow = atom->avec->dihedrals_allow;
int impropers_allow = atom->avec->impropers_allow;
tagint *tag = atom->tag;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
int nlocal = atom->nlocal;
if (nlocal > maxinfluenced) {
maxinfluenced = atom->nmax;
memory->destroy(influenced);
memory->create(influenced,maxinfluenced,"bond/break:influenced");
}
nangles = 0;
ndihedrals = 0;
nimpropers = 0;
for (i = 0; i < nlocal; i++) {
influenced[i] = 0;
slist = special[i];
for (j = 0; j < nbreakall; j++) {
id1 = brokenall[j][0];
id2 = brokenall[j][1];
influence = 0;
if (tag[i] == id1 || tag[i] == id2) influence = 1;
else {
n = nspecial[i][2];
found = 0;
for (k = 0; k < n; k++)
if (slist[k] == id1 || slist[k] == id2) found++;
if (found == 2) influence = 1;
}
if (!influence) continue;
influenced[i] = 1;
if (angles_allow) break_angles(i,id1,id2);
if (dihedrals_allow) break_dihedrals(i,id1,id2);
if (impropers_allow) break_impropers(i,id1,id2);
}
}
int newton_bond = force->newton_bond;
int all;
if (angles_allow) {
MPI_Allreduce(&nangles,&all,1,MPI_INT,MPI_SUM,world);
if (!newton_bond) all /= 3;
atom->nangles -= all;
}
if (dihedrals_allow) {
MPI_Allreduce(&ndihedrals,&all,1,MPI_INT,MPI_SUM,world);
if (!newton_bond) all /= 4;
atom->ndihedrals -= all;
}
if (impropers_allow) {
MPI_Allreduce(&nimpropers,&all,1,MPI_INT,MPI_SUM,world);
if (!newton_bond) all /= 4;
atom->nimpropers -= all;
}
for (i = 0; i < nlocal; i++)
if (influenced[i]) rebuild_special(i,id1,id2);
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void FixBondBreak::rebuild_special(int m, tagint id1, tagint id2)
{
int i,j,n,n1,n2,n3,cn1,cn2,cn3;
tagint *slist;
tagint *tag = atom->tag;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
// new 1-2 neighs of atom M
slist = special[m];
n1 = nspecial[m][0];
cn1 = 0;
for (i = 0; i < n1; i++)
copy[cn1++] = slist[i];
// new 1-3 neighs of atom M, based on 1-2 neighs of 1-2 neighs
// exclude self
// remove duplicates after adding all possible 1-3 neighs
cn2 = cn1;
for (i = 0; i < cn1; i++) {
n = atom->map(copy[i]);
slist = special[n];
n1 = nspecial[n][0];
for (j = 0; j < n1; j++)
if (slist[j] != tag[m]) copy[cn2++] = slist[j];
}
cn2 = dedup(cn1,cn2,copy);
// new 1-4 neighs of atom M, based on 1-2 neighs of 1-3 neighs
// exclude self
// remove duplicates after adding all possible 1-4 neighs
cn3 = cn2;
for (i = cn1; i < cn2; i++) {
n = atom->map(copy[i]);
slist = special[n];
n1 = nspecial[n][0];
for (j = 0; j < n1; j++)
if (slist[j] != tag[m]) copy[cn3++] = slist[j];
}
cn3 = dedup(cn2,cn3,copy);
// store new special list with atom M
nspecial[m][0] = cn1;
nspecial[m][1] = cn2;
nspecial[m][2] = cn3;
memcpy(special[m],copy,cn3*sizeof(int));
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void FixBondBreak::break_angles(int m, tagint id1, tagint id2)
{
int j,found;
int num_angle = atom->num_angle[m];
int *angle_type = atom->angle_type[m];
tagint *angle_atom1 = atom->angle_atom1[m];
tagint *angle_atom2 = atom->angle_atom2[m];
tagint *angle_atom3 = atom->angle_atom3[m];
int i = 0;
while (i < num_angle) {
found = 0;
if (angle_atom1[i] == id1 && angle_atom2[i] == id2) found = 1;
else if (angle_atom2[i] == id1 && angle_atom3[i] == id2) found = 1;
else if (angle_atom1[i] == id2 && angle_atom2[i] == id1) found = 1;
else if (angle_atom2[i] == id2 && angle_atom3[i] == id1) found = 1;
if (!found) i++;
else {
for (j = i; j < num_angle-1; j++) {
angle_type[j] = angle_type[j+1];
angle_atom1[j] = angle_atom1[j+1];
angle_atom2[j] = angle_atom2[j+1];
angle_atom3[j] = angle_atom3[j+1];
}
num_angle--;
nangles++;
}
}
atom->num_angle[m] = num_angle;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void FixBondBreak::break_dihedrals(int m, tagint id1, tagint id2)
{
int j,found;
int num_dihedral = atom->num_dihedral[m];
int *dihedral_type = atom->dihedral_type[m];
tagint *dihedral_atom1 = atom->dihedral_atom1[m];
tagint *dihedral_atom2 = atom->dihedral_atom2[m];
tagint *dihedral_atom3 = atom->dihedral_atom3[m];
tagint *dihedral_atom4 = atom->dihedral_atom4[m];
int i = 0;
while (i < num_dihedral) {
found = 0;
if (dihedral_atom1[i] == id1 && dihedral_atom2[i] == id2) found = 1;
else if (dihedral_atom2[i] == id1 && dihedral_atom3[i] == id2) found = 1;
else if (dihedral_atom3[i] == id1 && dihedral_atom4[i] == id2) found = 1;
else if (dihedral_atom1[i] == id2 && dihedral_atom2[i] == id1) found = 1;
else if (dihedral_atom2[i] == id2 && dihedral_atom3[i] == id1) found = 1;
else if (dihedral_atom3[i] == id2 && dihedral_atom4[i] == id1) found = 1;
if (!found) i++;
else {
for (j = i; j < num_dihedral-1; j++) {
dihedral_type[j] = dihedral_type[j+1];
dihedral_atom1[j] = dihedral_atom1[j+1];
dihedral_atom2[j] = dihedral_atom2[j+1];
dihedral_atom3[j] = dihedral_atom3[j+1];
dihedral_atom4[j] = dihedral_atom4[j+1];
}
num_dihedral--;
ndihedrals++;
}
}
atom->num_dihedral[m] = num_dihedral;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void FixBondBreak::break_impropers(int m, tagint id1, tagint id2)
{
int j,found;
int num_improper = atom->num_improper[m];
int *improper_type = atom->improper_type[m];
tagint *improper_atom1 = atom->improper_atom1[m];
tagint *improper_atom2 = atom->improper_atom2[m];
tagint *improper_atom3 = atom->improper_atom3[m];
tagint *improper_atom4 = atom->improper_atom4[m];
int i = 0;
while (i < num_improper) {
found = 0;
if (improper_atom1[i] == id1 && improper_atom2[i] == id2) found = 1;
else if (improper_atom1[i] == id1 && improper_atom3[i] == id2) found = 1;
else if (improper_atom1[i] == id1 && improper_atom4[i] == id2) found = 1;
else if (improper_atom1[i] == id2 && improper_atom2[i] == id1) found = 1;
else if (improper_atom1[i] == id2 && improper_atom3[i] == id1) found = 1;
else if (improper_atom1[i] == id2 && improper_atom4[i] == id1) found = 1;
if (!found) i++;
else {
for (j = i; j < num_improper-1; j++) {
improper_type[j] = improper_type[j+1];
improper_atom1[j] = improper_atom1[j+1];
improper_atom2[j] = improper_atom2[j+1];
improper_atom3[j] = improper_atom3[j+1];
improper_atom4[j] = improper_atom4[j+1];
}
num_improper--;
nimpropers++;
}
}
atom->num_improper[m] = num_improper;
}
/* ----------------------------------------------------------------------
remove all ID duplicates in copy from Nstart:Nstop-1
compare to all previous values in copy
return N decremented by any discarded duplicates
------------------------------------------------------------------------- */
int FixBondBreak::dedup(int nstart, int nstop, tagint *copy)
{
int i;
int m = nstart;
while (m < nstop) {
for (i = 0; i < m; i++)
if (copy[i] == copy[m]) {
copy[m] = copy[nstop-1];
nstop--;
break;
}
if (i == m) m++;
}
return nstop;
}
/* ---------------------------------------------------------------------- */
@ -312,30 +714,61 @@ void FixBondBreak::post_integrate_respa(int ilevel, int iloop)
/* ---------------------------------------------------------------------- */
int FixBondBreak::pack_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int pbc_flag, int *pbc)
{
int i,j,m;
int i,j,k,m,ns;
if (commflag == 1) {
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = ubuf(partner[j]).d;
buf[m++] = probability[j];
}
return 2;
}
int **nspecial = atom->nspecial;
int **special = atom->special;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = partner[j];
buf[m++] = probability[j];
ns = nspecial[j][0];
buf[m++] = ubuf(ns).d;
for (k = 0; k < ns; k++)
buf[m++] = ubuf(special[j][k]).d;
}
return 2;
return m;
}
/* ---------------------------------------------------------------------- */
void FixBondBreak::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
int i,j,m,ns,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
partner[i] = static_cast<tagint> (buf[m++]);
probability[i] = buf[m++];
if (commflag == 1) {
m = 0;
last = first + n;
for (i = first; i < last; i++) {
partner[i] = (tagint) ubuf(buf[m++]).i;
probability[i] = buf[m++];
}
} else {
int **nspecial = atom->nspecial;
int **special = atom->special;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
ns = (int) ubuf(buf[m++]).i;
nspecial[i][0] = ns;
for (j = 0; j < ns; j++)
special[i][j] = (tagint) ubuf(buf[m++]).i;
}
}
}
@ -348,7 +781,7 @@ int FixBondBreak::pack_reverse_comm(int n, int first, double *buf)
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = partner[i];
buf[m++] = ubuf(partner[i]).d;
buf[m++] = distsq[i];
}
return 2;
@ -364,12 +797,55 @@ void FixBondBreak::unpack_reverse_comm(int n, int *list, double *buf)
for (i = 0; i < n; i++) {
j = list[i];
if (buf[m+1] > distsq[j]) {
partner[j] = static_cast<tagint> (buf[m++]);
partner[j] = (tagint) ubuf(buf[m++]).i;
distsq[j] = buf[m++];
} else m += 2;
}
}
/* ---------------------------------------------------------------------- */
void FixBondBreak::print_bb()
{
for (int i = 0; i < atom->nlocal; i++) {
printf("TAG %i: %d nbonds: ",atom->tag[i],atom->num_bond[i]);
for (int j = 0; j < atom->num_bond[i]; j++) {
printf(" %d",atom->bond_atom[i][j]);
}
printf("\n");
printf("TAG %i: %d nangles: ",atom->tag[i],atom->num_angle[i]);
for (int j = 0; j < atom->num_angle[i]; j++) {
printf(" %d %d %d,",atom->angle_atom1[i][j],
atom->angle_atom2[i][j],atom->angle_atom3[i][j]);
}
printf("\n");
printf("TAG %i: %d ndihedrals: ",atom->tag[i],atom->num_dihedral[i]);
for (int j = 0; j < atom->num_dihedral[i]; j++) {
printf(" %d %d %d %d,",atom->dihedral_atom1[i][j],
atom->dihedral_atom2[i][j],atom->dihedral_atom3[i][j],
atom->dihedral_atom4[i][j]);
}
printf("\n");
printf("TAG %i: %d %d %d nspecial: ",atom->tag[i],
atom->nspecial[i][0],atom->nspecial[i][1],atom->nspecial[i][2]);
for (int j = 0; j < atom->nspecial[i][2]; j++) {
printf(" %d",atom->special[i][j]);
}
printf("\n");
}
}
/* ---------------------------------------------------------------------- */
void FixBondBreak::print_copy(const char *str, tagint m,
int n1, int n2, int n3, int *v)
{
printf("%s %i: %d %d %d nspecial: ",str,m,n1,n2,n3);
for (int j = 0; j < n3; j++) printf(" %d",v[j]);
printf("\n");
}
/* ---------------------------------------------------------------------- */
double FixBondBreak::compute_vector(int n)
@ -385,7 +861,8 @@ double FixBondBreak::compute_vector(int n)
double FixBondBreak::memory_usage()
{
int nmax = atom->nmax;
double bytes = nmax * sizeof(int);
double bytes = nmax * sizeof(tagint);
bytes += nmax * sizeof(double);
bytes += maxinfluenced * sizeof(int);
return bytes;
}

31
src/MC/fix_bond_break.h
View File

@ -41,17 +41,44 @@ class FixBondBreak : public Fix {
double memory_usage();
private:
int me;
int me,nprocs;
int btype,seed;
double cutsq,fraction;
double cutoff,cutsq,fraction;
double commextent;
int breakcount,breakcounttotal;
int nmax;
tagint *partner;
double *distsq,*probability;
int nbreak,nbreakall;
int maxbreak,maxbreakall;
tagint **broken,**brokenall;
double *inbuf;
int maxinfluenced;
int *influenced;
int *recvcounts,*displs;
tagint *copy;
class RanMars *random;
int nlevels_respa;
int commflag;
int nbroken;
int nangles,ndihedrals,nimpropers;
void update_topology();
void rebuild_special(int, tagint, tagint);
void break_angles(int, tagint, tagint);
void break_dihedrals(int, tagint, tagint);
void break_impropers(int, tagint, tagint);
int dedup(int, int, tagint *);
// DEBUG
void print_bb();
void print_copy(const char *, tagint, int, int, int, int *);
};
}

55
src/comm.cpp
View File

@ -1518,6 +1518,61 @@ void Comm::forward_comm_array(int n, double **array)
}
}
/* ----------------------------------------------------------------------
exchange info provided with all 6 stencil neighbors
------------------------------------------------------------------------- */
int Comm::exchange_variable(int n, double *inbuf, double *&outbuf)
{
int nsend,nrecv,nrecv1,nrecv2;
MPI_Request request;
MPI_Status status;
nrecv = n;
if (nrecv > maxrecv) grow_recv(nrecv);
memcpy(buf_recv,inbuf,nrecv*sizeof(double));
// loop over dimensions
for (int dim = 0; dim < 3; dim++) {
// no exchange if only one proc in a dimension
if (procgrid[dim] == 1) continue;
// send/recv info in both directions using same buf_recv
// if 2 procs in dimension, single send/recv
// if more than 2 procs in dimension, send/recv to both neighbors
nsend = nrecv;
MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][0],0,
&nrecv1,1,MPI_INT,procneigh[dim][1],0,world,&status);
nrecv += nrecv1;
if (procgrid[dim] > 2) {
MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][1],0,
&nrecv2,1,MPI_INT,procneigh[dim][0],0,world,&status);
nrecv += nrecv2;
} else nrecv2 = 0;
if (nrecv > maxrecv) grow_recv(nrecv);
MPI_Irecv(&buf_recv[nsend],nrecv1,MPI_DOUBLE,procneigh[dim][1],0,
world,&request);
MPI_Send(buf_recv,nsend,MPI_DOUBLE,procneigh[dim][0],0,world);
MPI_Wait(&request,&status);
if (procgrid[dim] > 2) {
MPI_Irecv(&buf_recv[nsend+nrecv1],nrecv2,MPI_DOUBLE,procneigh[dim][0],0,
world,&request);
MPI_Send(buf_recv,nsend,MPI_DOUBLE,procneigh[dim][1],0,world);
MPI_Wait(&request,&status);
}
}
outbuf = buf_recv;
return nrecv;
}
/* ----------------------------------------------------------------------
communicate inbuf around full ring of processors with messtag
nbytes = size of inbuf = n datums * nper bytes

2
src/comm.h
View File

@ -63,6 +63,8 @@ class Comm : protected Pointers {
virtual void reverse_comm_dump(class Dump *); // reverse comm from a Dump
void forward_comm_array(int, double **); // forward comm of array
int exchange_variable(int, double *, // exchange of info on neigh stencil
double *&);
void ring(int, int, void *, int, void (*)(int, char *), // ring comm
void *, int self = 1);
int read_lines_from_file(FILE *, int, int, char *); // read/bcast file lines

2
src/fix_property_atom.cpp
View File

@ -429,7 +429,7 @@ int FixPropertyAtom::pack_border(int n, int *list, double *buf)
int *ivector = atom->ivector[index[k]];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = ubuf(ivector[j]).i;
buf[m++] = ubuf(ivector[j]).d;
}
} else if (style[k] == DOUBLE) {
double *dvector = atom->dvector[index[k]];

15
src/force.cpp
View File

@ -485,6 +485,21 @@ Improper *Force::new_improper(const char *style, const char *suffix, int &sflag)
return NULL;
}
/* ----------------------------------------------------------------------
return ptr to current improper class or hybrid sub-class if matches style
------------------------------------------------------------------------- */
Improper *Force::improper_match(const char *style)
{
if (strcmp(improper_style,style) == 0) return improper;
else if (strcmp(improper_style,"hybrid") == 0) {
ImproperHybrid *hybrid = (ImproperHybrid *) bond;
for (int i = 0; i < hybrid->nstyles; i++)
if (strcmp(hybrid->keywords[i],style) == 0) return hybrid->styles[i];
}
return NULL;
}
/* ----------------------------------------------------------------------
new kspace style
------------------------------------------------------------------------- */

1
src/force.h
View File

@ -93,6 +93,7 @@ class Force : protected Pointers {
void create_improper(const char *, const char *suffix = NULL);
class Improper *new_improper(const char *, const char *, int &);
class Improper *improper_match(const char *);
void create_kspace(int, char **, const char *suffix = NULL);
class KSpace *new_kspace(int, char **, const char *, int &);

19
src/neighbor.cpp
View File

@ -565,11 +565,16 @@ void Neighbor::init()
if (processed) continue;
// pair and half: if there is a full non-occasional non-skip list
// pair and half and newton != 2:
// if there is a full non-occasional non-skip list
// change this list to half_from_full and point at the full list
// parent could be copy list or pair or fix
// could remove newton != 2 check if added half_from_full_no_newton_ghost
// option in neigh_derive.cpp and below in choose_build()
// this would require full list had ghost info
// would be useful when reax/c used in hybrid mode, e.g. with airebo
if (requests[i]->pair && requests[i]->half) {
if (requests[i]->pair && requests[i]->half && requests[i]->newton != 2) {
for (j = 0; j < nrequest; j++) {
if (!lists[j]) continue;
if (requests[j]->full && requests[j]->occasional == 0 &&
@ -924,8 +929,14 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
else pb = &Neighbor::skip_from;
} else if (rq->half_from_full) {
if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton;
else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton;
if (rq->newton == 0) {
if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton;
else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton;
} else if (rq->newton == 1) {
pb = &Neighbor::half_from_full_newton;
} else if (rq->newton == 2) {
pb = &Neighbor::half_from_full_no_newton;
}
} else if (rq->half) {
if (style == NSQ) {

54
src/thermo.cpp
View File

@ -55,6 +55,7 @@ using namespace MathConst;
// evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail
// vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz,
// xlat, ylat, zlat
// bonds, angles, dihedrals, impropers,
// pxx, pyy, pzz, pxy, pxz, pyz
// fmax, fnorm
// cella, cellb, cellc, cellalpha, cellbeta, cellgamma
@ -767,6 +768,15 @@ void Thermo::parse_fields(char *str)
} else if (strcmp(word,"zlat") == 0) {
addfield("Zlat",&Thermo::compute_zlat,FLOAT);
} else if (strcmp(word,"bonds") == 0) {
addfield("Bonds",&Thermo::compute_bonds,BIGINT);
} else if (strcmp(word,"angles") == 0) {
addfield("Angles",&Thermo::compute_angles,BIGINT);
} else if (strcmp(word,"dihedrals") == 0) {
addfield("Diheds",&Thermo::compute_dihedrals,BIGINT);
} else if (strcmp(word,"impropers") == 0) {
addfield("Impros",&Thermo::compute_impropers,BIGINT);
} else if (strcmp(word,"pxx") == 0) {
addfield("Pxx",&Thermo::compute_pxx,FLOAT);
index_press_vector = add_compute(id_press,VECTOR);
@ -1269,14 +1279,16 @@ int Thermo::evaluate_keyword(char *word, double *answer)
else if (strcmp(word,"xz") == 0) compute_xz();
else if (strcmp(word,"yz") == 0) compute_yz();
else if (strcmp(word,"xlat") == 0) {
compute_xlat();
} else if (strcmp(word,"ylat") == 0) {
compute_ylat();
} else if (strcmp(word,"zlat") == 0) {
compute_zlat();
else if (strcmp(word,"xlat") == 0) compute_xlat();
else if (strcmp(word,"ylat") == 0) compute_ylat();
else if (strcmp(word,"zlat") == 0) compute_zlat();
} else if (strcmp(word,"pxx") == 0) {
else if (strcmp(word,"bonds") == 0) compute_bonds();
else if (strcmp(word,"angles") == 0) compute_angles();
else if (strcmp(word,"dihedrals") == 0) compute_dihedrals();
else if (strcmp(word,"impropers") == 0) compute_impropers();
else if (strcmp(word,"pxx") == 0) {
if (!pressure)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init press");
@ -1860,6 +1872,34 @@ void Thermo::compute_zlat()
/* ---------------------------------------------------------------------- */
void Thermo::compute_bonds()
{
bivalue = atom->nbonds;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_angles()
{
bivalue = atom->nangles;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_dihedrals()
{
bivalue = atom->ndihedrals;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_impropers()
{
bivalue = atom->nimpropers;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_pxx()
{
dvalue = pressure->vector[0];

5
src/thermo.h
View File

@ -169,6 +169,11 @@ class Thermo : protected Pointers {
void compute_ylat();
void compute_zlat();
void compute_bonds();
void compute_angles();
void compute_dihedrals();
void compute_impropers();
void compute_pxx();
void compute_pyy();
void compute_pzz();

2
src/version.h
View File

@ -1 +1 @@
#define LAMMPS_VERSION "11 Apr 2014"
#define LAMMPS_VERSION "12 Apr 2014"