diff --git a/doc/Manual.html b/doc/Manual.html index 1849ffdc32..4726983eb1 100644 --- a/doc/Manual.html +++ b/doc/Manual.html @@ -1,7 +1,7 @@ LAMMPS-ICMS Users Manual - + @@ -22,7 +22,7 @@

LAMMPS-ICMS Documentation

-

11 Apr 2014 version +

12 Apr 2014 version

Version info:

diff --git a/doc/Manual.txt b/doc/Manual.txt index 087549652d..1d75783e1d 100644 --- a/doc/Manual.txt +++ b/doc/Manual.txt @@ -1,7 +1,7 @@ LAMMPS-ICMS Users Manual LAMMPS Users Manual - + @@ -19,7 +19,7 @@

LAMMPS-ICMS Documentation :c,h3 -11 Apr 2014 version :c,h4 +12 Apr 2014 version :c,h4 Version info: :h4 diff --git a/doc/fix_bond_break.html b/doc/fix_bond_break.html index 4cae77cb34..dd9d4d9983 100644 --- a/doc/fix_bond_break.html +++ b/doc/fix_bond_break.html @@ -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 bond_style command. -Once the bond is broken it will be permanently deleted. This is -different than a pairwise 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 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 @@ -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.

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 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 command. @@ -132,22 +138,6 @@ minimization. built with that package. See the Making LAMMPS 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 or -dihedral_style. -

-

This fix requires that the pairwise weightings defined by the -special_bonds 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 diff --git a/doc/fix_bond_break.txt b/doc/fix_bond_break.txt index c3b1a502c2..556bfc1dcb 100644 --- a/doc/fix_bond_break.txt +++ b/doc/fix_bond_break.txt @@ -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 diff --git a/doc/fix_bond_create.html b/doc/fix_bond_create.html index 079bc359e8..3098617f90 100644 --- a/doc/fix_bond_create.html +++ b/doc/fix_bond_create.html @@ -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 read_data 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 command, -e.g. for a percolation simulation. +"extra bonds per atom" parameter must be set to allow for it. See the +read_data or create_box 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 +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 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 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 +or create_box or "special_bonds extra" commands for +details.

Note that even if your simulation starts with no bonds, you must define a bond_style and use the diff --git a/doc/fix_bond_create.txt b/doc/fix_bond_create.txt index 89fb3a0c88..c8bbdab829 100644 --- a/doc/fix_bond_create.txt +++ b/doc/fix_bond_create.txt @@ -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 diff --git a/doc/thermo_style.html b/doc/thermo_style.html index 2c43fde256..f8c8030406 100644 --- a/doc/thermo_style.html +++ b/doc/thermo_style.html @@ -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 lattice 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 diff --git a/doc/thermo_style.txt b/doc/thermo_style.txt index 0fbcbce8bb..f0adef0ab2 100644 --- a/doc/thermo_style.txt +++ b/doc/thermo_style.txt @@ -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 diff --git a/src/MC/fix_bond_break.cpp b/src/MC/fix_bond_break.cpp index 7f00bf37d5..99f78db824 100644 --- a/src/MC/fix_bond_break.cpp +++ b/src/MC/fix_bond_break.cpp @@ -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 (outbuf[m++]); + brokenall[i][1] = static_cast (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 (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 (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; } diff --git a/src/MC/fix_bond_break.h b/src/MC/fix_bond_break.h index 2d1d202284..e7b0b25065 100644 --- a/src/MC/fix_bond_break.h +++ b/src/MC/fix_bond_break.h @@ -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 *); }; } diff --git a/src/comm.cpp b/src/comm.cpp index 1bb82bd8c4..bf5213af53 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -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 diff --git a/src/comm.h b/src/comm.h index 83b247984b..e9c1bc68dd 100644 --- a/src/comm.h +++ b/src/comm.h @@ -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 diff --git a/src/fix_property_atom.cpp b/src/fix_property_atom.cpp index 0c573fba10..f74e135125 100644 --- a/src/fix_property_atom.cpp +++ b/src/fix_property_atom.cpp @@ -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]]; diff --git a/src/force.cpp b/src/force.cpp index 04dea56633..7e8e189e84 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/force.h b/src/force.h index bd72567812..a9efb08fbf 100644 --- a/src/force.h +++ b/src/force.h @@ -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 &); diff --git a/src/neighbor.cpp b/src/neighbor.cpp index c38d8db7aa..37305effca 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -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) { diff --git a/src/thermo.cpp b/src/thermo.cpp index 003fc0881b..6708758b82 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -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]; diff --git a/src/thermo.h b/src/thermo.h index 8c4e7435cb..7bfe9e394d 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -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(); diff --git a/src/version.h b/src/version.h index 7e739f872f..ee9cf1af8e 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "11 Apr 2014" +#define LAMMPS_VERSION "12 Apr 2014"