diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp
index 27bf4fc70..9947d9a84 100644
--- a/src/GRANULAR/fix_pour.cpp
+++ b/src/GRANULAR/fix_pour.cpp
@@ -1,1038 +1,1038 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "fix_pour.h"
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "update.h"
#include "comm.h"
#include "molecule.h"
#include "modify.h"
#include "fix_gravity.h"
#include "domain.h"
#include "region.h"
#include "region_block.h"
#include "region_cylinder.h"
#include "random_park.h"
#include "math_extra.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum{ATOM,MOLECULE};
enum{ONE,RANGE,POLY};
enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files
#define EPSILON 0.001
/* ---------------------------------------------------------------------- */
FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 6) error->all(FLERR,"Illegal fix pour command");
time_depend = 1;
if (!atom->radius_flag || !atom->rmass_flag)
error->all(FLERR,"Fix pour requires atom attributes radius, rmass");
// required args
ninsert = force->inumeric(FLERR,arg[3]);
ntype = force->inumeric(FLERR,arg[4]);
seed = force->inumeric(FLERR,arg[5]);
if (seed <= 0) error->all(FLERR,"Illegal fix pour command");
// read options from end of input line
options(narg-6,&arg[6]);
// error check on type
if (mode == ATOM && (ntype <= 0 || ntype > atom->ntypes))
error->all(FLERR,"Invalid atom type in fix pour command");
// error checks on region and its extent being inside simulation box
if (iregion == -1) error->all(FLERR,"Must specify a region in fix pour");
if (domain->regions[iregion]->bboxflag == 0)
error->all(FLERR,"Fix pour region does not support a bounding box");
if (domain->regions[iregion]->dynamic_check())
error->all(FLERR,"Fix pour region cannot be dynamic");
if (strcmp(domain->regions[iregion]->style,"block") == 0) {
region_style = 1;
xlo = ((RegBlock *) domain->regions[iregion])->xlo;
xhi = ((RegBlock *) domain->regions[iregion])->xhi;
ylo = ((RegBlock *) domain->regions[iregion])->ylo;
yhi = ((RegBlock *) domain->regions[iregion])->yhi;
zlo = ((RegBlock *) domain->regions[iregion])->zlo;
zhi = ((RegBlock *) domain->regions[iregion])->zhi;
if (xlo < domain->boxlo[0] || xhi > domain->boxhi[0] ||
ylo < domain->boxlo[1] || yhi > domain->boxhi[1] ||
zlo < domain->boxlo[2] || zhi > domain->boxhi[2])
error->all(FLERR,"Insertion region extends outside simulation box");
} else if (strcmp(domain->regions[iregion]->style,"cylinder") == 0) {
region_style = 2;
char axis = ((RegCylinder *) domain->regions[iregion])->axis;
xc = ((RegCylinder *) domain->regions[iregion])->c1;
yc = ((RegCylinder *) domain->regions[iregion])->c2;
rc = ((RegCylinder *) domain->regions[iregion])->radius;
zlo = ((RegCylinder *) domain->regions[iregion])->lo;
zhi = ((RegCylinder *) domain->regions[iregion])->hi;
if (axis != 'z')
error->all(FLERR,"Must use a z-axis cylinder region with fix pour");
if (xc-rc < domain->boxlo[0] || xc+rc > domain->boxhi[0] ||
yc-rc < domain->boxlo[1] || yc+rc > domain->boxhi[1] ||
zlo < domain->boxlo[2] || zhi > domain->boxhi[2])
error->all(FLERR,"Insertion region extends outside simulation box");
} else error->all(FLERR,"Must use a block or cylinder region with fix pour");
if (region_style == 2 && domain->dimension == 2)
error->all(FLERR,
"Must use a block region with fix pour for 2d simulations");
// error check and further setup for mode = MOLECULE
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use fix_pour unless atoms have IDs");
if (mode == MOLECULE) {
for (int i = 0; i < nmol; i++) {
if (onemols[i]->xflag == 0)
error->all(FLERR,"Fix pour molecule must have coordinates");
if (onemols[i]->typeflag == 0)
error->all(FLERR,"Fix pour molecule must have atom types");
if (ntype+onemols[i]->ntypes <= 0 ||
ntype+onemols[i]->ntypes > atom->ntypes)
error->all(FLERR,"Invalid atom type in fix pour mol command");
if (atom->molecular == 2 && onemols != atom->avec->onemols)
error->all(FLERR,"Fix pour molecule template ID must be same "
"as atom style template ID");
onemols[i]->check_attributes(0);
// fix pour uses geoemetric center of molecule for insertion
onemols[i]->compute_center();
}
}
if (rigidflag && mode == ATOM)
error->all(FLERR,"Cannot use fix pour rigid and not molecule");
if (shakeflag && mode == ATOM)
error->all(FLERR,"Cannot use fix pour shake and not molecule");
if (rigidflag && shakeflag)
error->all(FLERR,"Cannot use fix pour rigid and shake");
// setup of coords and imageflags array
if (mode == ATOM) natom_max = 1;
else {
natom_max = 0;
for (int i = 0; i < nmol; i++)
natom_max = MAX(natom_max,onemols[i]->natoms);
}
memory->create(coords,natom_max,4,"pour:coords");
memory->create(imageflags,natom_max,"pour:imageflags");
// find current max atom and molecule IDs if necessary
if (idnext) find_maxid();
// random number generator, same for all procs
random = new RanPark(lmp,seed);
// allgather arrays
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
recvcounts = new int[nprocs];
displs = new int[nprocs];
// grav = gravity in distance/time^2 units
// assume grav = -magnitude at this point, enforce in init()
int ifix;
for (ifix = 0; ifix < modify->nfix; ifix++) {
if (strcmp(modify->fix[ifix]->style,"gravity") == 0) break;
if (strcmp(modify->fix[ifix]->style,"gravity/omp") == 0) break;
}
if (ifix == modify->nfix)
error->all(FLERR,"No fix gravity defined for fix pour");
grav = - ((FixGravity *) modify->fix[ifix])->magnitude * force->ftm2v;
// nfreq = timesteps between insertions
// should be time for a particle to fall from top of insertion region
// to bottom, taking into account that the region may be moving
// set these 2 eqs equal to each other, solve for smallest positive t
// x = zhi + vz*t + 1/2 grav t^2
// x = zlo + rate*t
// gives t = [-(vz-rate) - sqrt((vz-rate)^2 - 2*grav*(zhi-zlo))] / grav
// where zhi-zlo > 0, grav < 0, and vz & rate can be either > 0 or < 0
double v_relative,delta;
if (domain->dimension == 3) {
v_relative = vz - rate;
delta = zhi - zlo;
} else {
v_relative = vy - rate;
delta = yhi - ylo;
}
double t =
(-v_relative - sqrt(v_relative*v_relative - 2.0*grav*delta)) / grav;
nfreq = static_cast<int> (t/update->dt + 0.5);
// 1st insertion on next timestep
force_reneighbor = 1;
next_reneighbor = update->ntimestep + 1;
nfirst = next_reneighbor;
ninserted = 0;
// nper = # to insert each time
// depends on specified volume fraction
// volume = volume of insertion region
// volume_one = volume of inserted particle (with max possible radius)
// in 3d, insure dy >= 1, for quasi-2d simulations
- double volume,volume_one;
+ double volume,volume_one=1.0;
molradius_max = 0.0;
if (mode == MOLECULE) {
for (int i = 0; i < nmol; i++)
molradius_max = MAX(molradius_max,onemols[i]->molradius);
}
if (domain->dimension == 3) {
if (region_style == 1) {
double dy = yhi - ylo;
if (dy < 1.0) dy = 1.0;
volume = (xhi-xlo) * dy * (zhi-zlo);
} else volume = MY_PI*rc*rc * (zhi-zlo);
if (mode == MOLECULE) {
volume_one = 4.0/3.0 * MY_PI * molradius_max*molradius_max*molradius_max;
} else if (dstyle == ONE || dstyle == RANGE) {
volume_one = 4.0/3.0 * MY_PI * radius_max*radius_max*radius_max;
} else if (dstyle == POLY) {
volume_one = 0.0;
for (int i = 0; i < npoly; i++)
volume_one += (4.0/3.0 * MY_PI *
radius_poly[i]*radius_poly[i]*radius_poly[i]) * frac_poly[i];
}
} else {
volume = (xhi-xlo) * (yhi-ylo);
if (mode == MOLECULE) {
volume_one = MY_PI * molradius_max*molradius_max;
} else if (dstyle == ONE || dstyle == RANGE) {
volume_one = MY_PI * radius_max*radius_max;
} else if (dstyle == POLY) {
volume_one = 0.0;
for (int i = 0; i < npoly; i++)
volume_one += (MY_PI * radius_poly[i]*radius_poly[i]) * frac_poly[i];
}
}
nper = static_cast<int> (volfrac*volume/volume_one);
int nfinal = update->ntimestep + 1 + (ninsert-1)/nper * nfreq;
// print stats
if (me == 0) {
if (screen)
fprintf(screen,
"Particle insertion: %d every %d steps, %d by step %d\n",
nper,nfreq,ninsert,nfinal);
if (logfile)
fprintf(logfile,
"Particle insertion: %d every %d steps, %d by step %d\n",
nper,nfreq,ninsert,nfinal);
}
}
/* ---------------------------------------------------------------------- */
FixPour::~FixPour()
{
delete random;
delete [] molfrac;
delete [] idrigid;
delete [] idshake;
delete [] radius_poly;
delete [] frac_poly;
memory->destroy(coords);
memory->destroy(imageflags);
delete [] recvcounts;
delete [] displs;
}
/* ---------------------------------------------------------------------- */
int FixPour::setmask()
{
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixPour::init()
{
if (domain->triclinic)
error->all(FLERR,"Cannot use fix pour with triclinic box");
// insure gravity fix exists
// for 3d must point in -z, for 2d must point in -y
// else insertion cannot work
int ifix;
for (ifix = 0; ifix < modify->nfix; ifix++) {
if (strcmp(modify->fix[ifix]->style,"gravity") == 0) break;
if (strcmp(modify->fix[ifix]->style,"gravity/omp") == 0) break;
}
if (ifix == modify->nfix)
error->all(FLERR,"No fix gravity defined for fix pour");
double xgrav = ((FixGravity *) modify->fix[ifix])->xgrav;
double ygrav = ((FixGravity *) modify->fix[ifix])->ygrav;
double zgrav = ((FixGravity *) modify->fix[ifix])->zgrav;
if (domain->dimension == 3) {
if (fabs(xgrav) > EPSILON || fabs(ygrav) > EPSILON ||
fabs(zgrav+1.0) > EPSILON)
error->all(FLERR,"Gravity must point in -z to use with fix pour in 3d");
} else {
if (fabs(xgrav) > EPSILON || fabs(ygrav+1.0) > EPSILON ||
fabs(zgrav) > EPSILON)
error->all(FLERR,"Gravity must point in -y to use with fix pour in 2d");
}
double gnew = - ((FixGravity *) modify->fix[ifix])->magnitude * force->ftm2v;
if (gnew != grav)
error->all(FLERR,"Gravity changed since fix pour was created");
// if rigidflag defined, check for rigid/small fix
// its molecule template must be same as this one
fixrigid = NULL;
if (rigidflag) {
int ifix = modify->find_fix(idrigid);
if (ifix < 0) error->all(FLERR,"Fix pour rigid fix does not exist");
fixrigid = modify->fix[ifix];
int tmp;
if (onemols != (Molecule **) fixrigid->extract("onemol",tmp))
error->all(FLERR,
"Fix pour and fix rigid/small not using "
"same molecule template ID");
}
// if shakeflag defined, check for SHAKE fix
// its molecule template must be same as this one
fixshake = NULL;
if (shakeflag) {
int ifix = modify->find_fix(idshake);
if (ifix < 0) error->all(FLERR,"Fix pour shake fix does not exist");
fixshake = modify->fix[ifix];
int tmp;
if (onemols != (Molecule **) fixshake->extract("onemol",tmp))
error->all(FLERR,"Fix pour and fix shake not using "
"same molecule template ID");
}
}
/* ----------------------------------------------------------------------
perform particle insertion
------------------------------------------------------------------------- */
void FixPour::pre_exchange()
{
int i,j,m,flag,nlocalprev,imol,natom;
double r[3],rotmat[3][3],quat[4],vnew[3];
double *newcoord;
// just return if should not be called on this timestep
if (next_reneighbor != update->ntimestep) return;
// find current max atom and molecule IDs if necessary
if (!idnext) find_maxid();
// nnew = # of particles (atoms or molecules) to insert this timestep
int nnew = nper;
if (ninserted + nnew > ninsert) nnew = ninsert - ninserted;
// lo/hi current = z (or y) bounds of insertion region this timestep
int dimension = domain->dimension;
if (dimension == 3) {
lo_current = zlo + (update->ntimestep - nfirst) * update->dt * rate;
hi_current = zhi + (update->ntimestep - nfirst) * update->dt * rate;
} else {
lo_current = ylo + (update->ntimestep - nfirst) * update->dt * rate;
hi_current = yhi + (update->ntimestep - nfirst) * update->dt * rate;
}
// ncount = # of my atoms that overlap the insertion region
// nprevious = total of ncount across all procs
int ncount = 0;
for (i = 0; i < atom->nlocal; i++)
if (overlap(i)) ncount++;
int nprevious;
MPI_Allreduce(&ncount,&nprevious,1,MPI_INT,MPI_SUM,world);
// xmine is for my atoms
// xnear is for atoms from all procs + atoms to be inserted
double **xmine,**xnear;
memory->create(xmine,ncount,4,"fix_pour:xmine");
memory->create(xnear,nprevious+nnew*natom_max,4,"fix_pour:xnear");
int nnear = nprevious;
// setup for allgatherv
int n = 4*ncount;
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];
// load up xmine array
double **x = atom->x;
double *radius = atom->radius;
ncount = 0;
for (i = 0; i < atom->nlocal; i++)
if (overlap(i)) {
xmine[ncount][0] = x[i][0];
xmine[ncount][1] = x[i][1];
xmine[ncount][2] = x[i][2];
xmine[ncount][3] = radius[i];
ncount++;
}
// perform allgatherv to acquire list of nearby particles on all procs
double *ptr = NULL;
if (ncount) ptr = xmine[0];
MPI_Allgatherv(ptr,4*ncount,MPI_DOUBLE,
xnear[0],recvcounts,displs,MPI_DOUBLE,world);
// insert new particles into xnear list, one by one
// check against all nearby atoms and previously inserted ones
// if there is an overlap then try again at same z (3d) or y (2d) coord
// else insert by adding to xnear list
// max = maximum # of insertion attempts for all particles
// h = height, biased to give uniform distribution in time of insertion
// for MOLECULE mode:
// coords = coords of all atoms in particle
// perform random rotation around center pt
// apply PBC so final coords are inside box
// store image flag modified due to PBC
int success;
double radtmp,delx,dely,delz,rsq,radsum,rn,h;
double coord[3];
int nfix = modify->nfix;
Fix **fix = modify->fix;
double denstmp;
double *sublo = domain->sublo;
double *subhi = domain->subhi;
int nsuccess = 0;
int attempt = 0;
int maxiter = nnew * maxattempt;
while (nsuccess < nnew) {
rn = random->uniform();
h = hi_current - rn*rn * (hi_current-lo_current);
if (mode == ATOM) radtmp = radius_sample();
success = 0;
while (attempt < maxiter) {
attempt++;
xyz_random(h,coord);
if (mode == ATOM) {
natom = 1;
coords[0][0] = coord[0];
coords[0][1] = coord[1];
coords[0][2] = coord[2];
coords[0][3] = radtmp;
imageflags[0] = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
} else {
double rng = random->uniform();
imol = 0;
while (rng > molfrac[imol]) imol++;
natom = onemols[imol]->natoms;
if (dimension == 3) {
r[0] = random->uniform() - 0.5;
r[1] = random->uniform() - 0.5;
r[2] = random->uniform() - 0.5;
} else {
r[0] = r[1] = 0.0;
r[2] = 1.0;
}
double theta = random->uniform() * MY_2PI;
MathExtra::norm3(r);
MathExtra::axisangle_to_quat(r,theta,quat);
MathExtra::quat_to_mat(quat,rotmat);
for (i = 0; i < natom; i++) {
MathExtra::matvec(rotmat,onemols[imol]->dx[i],coords[i]);
coords[i][0] += coord[0];
coords[i][1] += coord[1];
coords[i][2] += coord[2];
// coords[3] = particle radius
// default to 0.5, if radii not defined in Molecule
// same as atom->avec->create_atom(), invoked below
if (onemols[imol]->radiusflag)
coords[i][3] = onemols[imol]->radius[i];
else coords[i][3] = 0.5;
imageflags[i] = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
domain->remap(coords[i],imageflags[i]);
}
}
// if any pair of atoms overlap, try again
// use minimum_image() to account for PBC
for (m = 0; m < natom; m++) {
for (i = 0; i < nnear; i++) {
delx = coords[m][0] - xnear[i][0];
dely = coords[m][1] - xnear[i][1];
delz = coords[m][2] - xnear[i][2];
domain->minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz;
radsum = coords[m][3] + xnear[i][3];
if (rsq <= radsum*radsum) break;
}
if (i < nnear) break;
}
if (m == natom) {
success = 1;
break;
}
}
if (!success) break;
// proceed with insertion
nsuccess++;
nlocalprev = atom->nlocal;
// add all atoms in particle to xnear
for (m = 0; m < natom; m++) {
xnear[nnear][0] = coords[m][0];
xnear[nnear][1] = coords[m][1];
xnear[nnear][2] = coords[m][2];
xnear[nnear][3] = coords[m][3];
nnear++;
}
// choose random velocity for new particle
// used for every atom in molecule
// z velocity set to what velocity would be if particle
// had fallen from top of insertion region
// this gives continuous stream of atoms
// solution for v from these 2 eqs, after eliminate t:
// v = vz + grav*t
// coord[2] = hi_current + vz*t + 1/2 grav t^2
if (dimension == 3) {
vnew[0] = vxlo + random->uniform() * (vxhi-vxlo);
vnew[1] = vylo + random->uniform() * (vyhi-vylo);
vnew[2] = -sqrt(vz*vz + 2.0*grav*(coord[2]-hi_current));
} else {
vnew[0] = vxlo + random->uniform() * (vxhi-vxlo);
vnew[1] = -sqrt(vy*vy + 2.0*grav*(coord[1]-hi_current));
vnew[2] = 0.0;
}
// check if new atoms are in my sub-box or above it if I am highest proc
// if so, add atom to my list via create_atom()
// initialize additional info about the atoms
// set group mask to "all" plus fix group
for (m = 0; m < natom; m++) {
if (mode == ATOM)
denstmp = density_lo + random->uniform() * (density_hi-density_lo);
newcoord = coords[m];
flag = 0;
if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] &&
newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1;
else if (dimension == 3 && newcoord[2] >= domain->boxhi[2]) {
if (comm->layout != LAYOUT_TILED) {
if (comm->myloc[2] == comm->procgrid[2]-1 &&
newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1;
} else {
if (comm->mysplit[2][1] == 1.0 &&
newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1;
}
} else if (dimension == 2 && newcoord[1] >= domain->boxhi[1]) {
if (comm->layout != LAYOUT_TILED) {
if (comm->myloc[1] == comm->procgrid[1]-1 &&
newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1;
} else {
if (comm->mysplit[1][1] == 1.0 &&
newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1;
}
}
if (flag) {
if (mode == ATOM) atom->avec->create_atom(ntype,coords[m]);
else atom->avec->create_atom(ntype+onemols[imol]->type[m],coords[m]);
int n = atom->nlocal - 1;
atom->tag[n] = maxtag_all + m+1;
if (mode == MOLECULE) {
if (atom->molecule_flag) atom->molecule[n] = maxmol_all+1;
if (atom->molecular == 2) {
atom->molindex[n] = 0;
atom->molatom[n] = m;
}
}
atom->mask[n] = 1 | groupbit;
atom->image[n] = imageflags[m];
atom->v[n][0] = vnew[0];
atom->v[n][1] = vnew[1];
atom->v[n][2] = vnew[2];
if (mode == ATOM) {
radtmp = newcoord[3];
atom->radius[n] = radtmp;
atom->rmass[n] = 4.0*MY_PI/3.0 * radtmp*radtmp*radtmp * denstmp;
} else atom->add_molecule_atom(onemols[imol],m,n,maxtag_all);
for (j = 0; j < nfix; j++)
if (fix[j]->create_attribute) fix[j]->set_arrays(n);
}
}
// FixRigidSmall::set_molecule stores rigid body attributes
// coord is new position of geometric center of mol, not COM
// FixShake::set_molecule stores shake info for molecule
if (rigidflag)
fixrigid->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
else if (shakeflag)
fixshake->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
maxtag_all += natom;
if (mode == MOLECULE && atom->molecule_flag) maxmol_all++;
}
// warn if not successful with all insertions b/c too many attempts
int ninserted_atoms = nnear - nprevious;
int ninserted_mols = ninserted_atoms / natom;
ninserted += ninserted_mols;
if (ninserted_mols < nnew && me == 0)
error->warning(FLERR,"Less insertions than requested",0);
// reset global natoms,nbonds,etc
// increment maxtag_all and maxmol_all if necessary
// if global map exists, reset it now instead of waiting for comm
// since adding atoms messes up ghosts
if (ninserted_atoms) {
atom->natoms += ninserted_atoms;
if (atom->natoms < 0 || atom->natoms > MAXBIGINT)
error->all(FLERR,"Too many total atoms");
if (mode == MOLECULE) {
atom->nbonds += onemols[imol]->nbonds * ninserted_mols;
atom->nangles += onemols[imol]->nangles * ninserted_mols;
atom->ndihedrals += onemols[imol]->ndihedrals * ninserted_mols;
atom->nimpropers += onemols[imol]->nimpropers * ninserted_mols;
}
if (maxtag_all >= MAXTAGINT)
error->all(FLERR,"New atom IDs exceed maximum allowed ID");
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
}
// free local memory
memory->destroy(xmine);
memory->destroy(xnear);
// next timestep to insert
if (ninserted < ninsert) next_reneighbor += nfreq;
else next_reneighbor = 0;
}
/* ----------------------------------------------------------------------
maxtag_all = current max atom ID for all atoms
maxmol_all = current max molecule ID for all atoms
------------------------------------------------------------------------- */
void FixPour::find_maxid()
{
tagint *tag = atom->tag;
tagint *molecule = atom->molecule;
int nlocal = atom->nlocal;
tagint max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]);
MPI_Allreduce(&max,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
if (mode == MOLECULE && molecule) {
max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]);
MPI_Allreduce(&max,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
}
}
/* ----------------------------------------------------------------------
check if particle i could overlap with a particle inserted into region
return 1 if yes, 0 if no
for ATOM mode, use delta with maximum size for inserted atoms
for MOLECULE mode, use delta with max radius of inserted molecules
account for PBC in overlap decision via outside() and minimum_image()
------------------------------------------------------------------------- */
int FixPour::overlap(int i)
{
double delta;
if (mode == ATOM) delta = atom->radius[i] + radius_max;
else delta = atom->radius[i] + molradius_max;
double *x = atom->x[i];
if (domain->dimension == 3) {
if (region_style == 1) {
if (outside(0,x[0],xlo-delta,xhi+delta)) return 0;
if (outside(1,x[1],ylo-delta,yhi+delta)) return 0;
if (outside(2,x[2],lo_current-delta,hi_current+delta)) return 0;
} else {
double delx = x[0] - xc;
double dely = x[1] - yc;
double delz = 0.0;
domain->minimum_image(delx,dely,delz);
double rsq = delx*delx + dely*dely;
double r = rc + delta;
if (rsq > r*r) return 0;
if (outside(2,x[2],lo_current-delta,hi_current+delta)) return 0;
}
} else {
if (outside(0,x[0],xlo-delta,xhi+delta)) return 0;
if (outside(1,x[1],lo_current-delta,hi_current+delta)) return 0;
}
return 1;
}
/* ----------------------------------------------------------------------
check if value is inside/outside lo/hi bounds in dimension
account for PBC if needed
return 1 if value is outside, 0 if inside
------------------------------------------------------------------------- */
int FixPour::outside(int dim, double value, double lo, double hi)
{
double boxlo = domain->boxlo[dim];
double boxhi = domain->boxhi[dim];
if (domain->periodicity[dim]) {
if (lo < boxlo && hi > boxhi) {
return 0;
} else if (lo < boxlo) {
if (value > hi && value < lo + domain->prd[dim]) return 1;
} else if (hi > boxhi) {
if (value > hi - domain->prd[dim] && value < lo) return 1;
} else {
if (value < lo || value > hi) return 1;
}
}
if (value < lo || value > hi) return 1;
return 0;
}
/* ---------------------------------------------------------------------- */
void FixPour::xyz_random(double h, double *coord)
{
if (domain->dimension == 3) {
if (region_style == 1) {
coord[0] = xlo + random->uniform() * (xhi-xlo);
coord[1] = ylo + random->uniform() * (yhi-ylo);
coord[2] = h;
} else {
double r1,r2;
while (1) {
r1 = random->uniform() - 0.5;
r2 = random->uniform() - 0.5;
if (r1*r1 + r2*r2 < 0.25) break;
}
coord[0] = xc + 2.0*r1*rc;
coord[1] = yc + 2.0*r2*rc;
coord[2] = h;
}
} else {
coord[0] = xlo + random->uniform() * (xhi-xlo);
coord[1] = h;
coord[2] = 0.0;
}
}
/* ---------------------------------------------------------------------- */
double FixPour::radius_sample()
{
if (dstyle == ONE) return radius_one;
if (dstyle == RANGE) return radius_lo +
random->uniform()*(radius_hi-radius_lo);
double value = random->uniform();
int i = 0;
double sum = 0.0;
while (sum < value) {
sum += frac_poly[i];
i++;
}
return radius_poly[i-1];
}
/* ----------------------------------------------------------------------
parse optional parameters at end of input line
------------------------------------------------------------------------- */
void FixPour::options(int narg, char **arg)
{
// defaults
iregion = -1;
mode = ATOM;
molfrac = NULL;
rigidflag = 0;
idrigid = NULL;
shakeflag = 0;
idshake = NULL;
idnext = 0;
dstyle = ONE;
radius_max = radius_one = 0.5;
radius_poly = frac_poly = NULL;
density_lo = density_hi = 1.0;
volfrac = 0.25;
maxattempt = 50;
rate = 0.0;
vxlo = vxhi = vylo = vyhi = vy = vz = 0.0;
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1) error->all(FLERR,"Fix pour region ID does not exist");
iarg += 2;
} else if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
int imol = atom->find_molecule(arg[iarg+1]);
if (imol == -1)
error->all(FLERR,"Molecule template ID for fix pour does not exist");
mode = MOLECULE;
onemols = &atom->molecules[imol];
nmol = onemols[0]->nset;
delete [] molfrac;
molfrac = new double[nmol];
molfrac[0] = 1.0/nmol;
for (int i = 1; i < nmol-1; i++) molfrac[i] = molfrac[i-1] + 1.0/nmol;
molfrac[nmol-1] = 1.0;
iarg += 2;
} else if (strcmp(arg[iarg],"molfrac") == 0) {
if (mode != MOLECULE) error->all(FLERR,"Illegal fix deposit command");
if (iarg+nmol+1 > narg) error->all(FLERR,"Illegal fix deposit command");
molfrac[0] = force->numeric(FLERR,arg[iarg+1]);
for (int i = 1; i < nmol; i++)
molfrac[i] = molfrac[i-1] + force->numeric(FLERR,arg[iarg+i+1]);
if (molfrac[nmol-1] < 1.0-EPSILON || molfrac[nmol-1] > 1.0+EPSILON)
error->all(FLERR,"Illegal fix deposit command");
molfrac[nmol-1] = 1.0;
iarg += nmol+1;
} else if (strcmp(arg[iarg],"rigid") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
int n = strlen(arg[iarg+1]) + 1;
delete [] idrigid;
idrigid = new char[n];
strcpy(idrigid,arg[iarg+1]);
rigidflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"shake") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
int n = strlen(arg[iarg+1]) + 1;
delete [] idshake;
idshake = new char[n];
strcpy(idshake,arg[iarg+1]);
shakeflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"id") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
if (strcmp(arg[iarg+1],"max") == 0) idnext = 0;
else if (strcmp(arg[iarg+1],"next") == 0) idnext = 1;
else error->all(FLERR,"Illegal fix pour command");
iarg += 2;
} else if (strcmp(arg[iarg],"diam") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
if (strcmp(arg[iarg+1],"one") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
dstyle = ONE;
radius_one = 0.5 * force->numeric(FLERR,arg[iarg+2]);
radius_max = radius_one;
iarg += 3;
} else if (strcmp(arg[iarg+1],"range") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix pour command");
dstyle = RANGE;
radius_lo = 0.5 * force->numeric(FLERR,arg[iarg+2]);
radius_hi = 0.5 * force->numeric(FLERR,arg[iarg+3]);
radius_max = radius_hi;
iarg += 4;
} else if (strcmp(arg[iarg+1],"poly") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
dstyle = POLY;
npoly = force->inumeric(FLERR,arg[iarg+2]);
if (npoly <= 0) error->all(FLERR,"Illegal fix pour command");
if (iarg+3 + 2*npoly > narg)
error->all(FLERR,"Illegal fix pour command");
radius_poly = new double[npoly];
frac_poly = new double[npoly];
iarg += 3;
radius_max = 0.0;
for (int i = 0; i < npoly; i++) {
radius_poly[i] = 0.5 * force->numeric(FLERR,arg[iarg++]);
frac_poly[i] = force->numeric(FLERR,arg[iarg++]);
if (radius_poly[i] <= 0.0 || frac_poly[i] < 0.0)
error->all(FLERR,"Illegal fix pour command");
radius_max = MAX(radius_max,radius_poly[i]);
}
double sum = 0.0;
for (int i = 0; i < npoly; i++) sum += frac_poly[i];
if (sum != 1.0)
error->all(FLERR,"Fix pour polydisperse fractions do not sum to 1.0");
} else error->all(FLERR,"Illegal fix pour command");
} else if (strcmp(arg[iarg],"dens") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
density_lo = force->numeric(FLERR,arg[iarg+1]);
density_hi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"vol") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
volfrac = force->numeric(FLERR,arg[iarg+1]);
maxattempt = force->inumeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"rate") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
rate = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"vel") == 0) {
if (domain->dimension == 3) {
if (iarg+6 > narg) error->all(FLERR,"Illegal fix pour command");
vxlo = force->numeric(FLERR,arg[iarg+1]);
vxhi = force->numeric(FLERR,arg[iarg+2]);
vylo = force->numeric(FLERR,arg[iarg+3]);
vyhi = force->numeric(FLERR,arg[iarg+4]);
vz = force->numeric(FLERR,arg[iarg+5]);
iarg += 6;
} else {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix pour command");
vxlo = force->numeric(FLERR,arg[iarg+1]);
vxhi = force->numeric(FLERR,arg[iarg+2]);
vy = force->numeric(FLERR,arg[iarg+3]);
vz = 0.0;
iarg += 4;
}
} else error->all(FLERR,"Illegal fix pour command");
}
}
/* ---------------------------------------------------------------------- */
void FixPour::reset_dt()
{
error->all(FLERR,"Cannot change timestep with fix pour");
}
/* ----------------------------------------------------------------------
extract particle radius for atom type = itype
------------------------------------------------------------------------- */
void *FixPour::extract(const char *str, int &itype)
{
if (strcmp(str,"radius") == 0) {
if (mode == ATOM) {
if (itype == ntype) oneradius = radius_max;
else oneradius = 0.0;
} else {
// find a molecule in template with matching type
for (int i = 0; i < nmol; i++) {
if (itype-ntype > onemols[i]->ntypes) continue;
double *radius = onemols[i]->radius;
int *type = onemols[i]->type;
int natoms = onemols[i]->natoms;
// check radii of matching types in Molecule
// default to 0.5, if radii not defined in Molecule
// same as atom->avec->create_atom(), invoked in pre_exchange()
oneradius = 0.0;
for (int i = 0; i < natoms; i++)
if (type[i] == itype-ntype) {
if (radius) oneradius = MAX(oneradius,radius[i]);
else oneradius = 0.5;
}
}
}
itype = 0;
return &oneradius;
}
return NULL;
}
diff --git a/src/MC/fix_bond_create.cpp b/src/MC/fix_bond_create.cpp
index 7c491fdf2..56dc7d1d9 100755
--- a/src/MC/fix_bond_create.cpp
+++ b/src/MC/fix_bond_create.cpp
@@ -1,1441 +1,1444 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "mpi.h"
#include "string.h"
#include "stdlib.h"
#include "fix_bond_create.h"
#include "update.h"
#include "respa.h"
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "pair.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define BIG 1.0e20
#define DELTA 16
/* ---------------------------------------------------------------------- */
FixBondCreate::FixBondCreate(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 8) error->all(FLERR,"Illegal fix bond/create command");
MPI_Comm_rank(world,&me);
nevery = force->inumeric(FLERR,arg[3]);
if (nevery <= 0) error->all(FLERR,"Illegal fix bond/create command");
force_reneighbor = 1;
next_reneighbor = -1;
vector_flag = 1;
size_vector = 2;
global_freq = 1;
extvector = 0;
iatomtype = force->inumeric(FLERR,arg[4]);
jatomtype = force->inumeric(FLERR,arg[5]);
double cutoff = force->numeric(FLERR,arg[6]);
btype = force->inumeric(FLERR,arg[7]);
if (iatomtype < 1 || iatomtype > atom->ntypes ||
jatomtype < 1 || jatomtype > atom->ntypes)
error->all(FLERR,"Invalid atom type in fix bond/create command");
if (cutoff < 0.0) error->all(FLERR,"Illegal fix bond/create command");
if (btype < 1 || btype > atom->nbondtypes)
error->all(FLERR,"Invalid bond type in fix bond/create command");
cutsq = cutoff*cutoff;
// optional keywords
imaxbond = 0;
inewtype = iatomtype;
jmaxbond = 0;
jnewtype = jatomtype;
fraction = 1.0;
int seed = 12345;
atype = dtype = itype = 0;
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg],"iparam") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix bond/create command");
imaxbond = force->inumeric(FLERR,arg[iarg+1]);
inewtype = force->inumeric(FLERR,arg[iarg+2]);
if (imaxbond < 0) error->all(FLERR,"Illegal fix bond/create command");
if (inewtype < 1 || inewtype > atom->ntypes)
error->all(FLERR,"Invalid atom type in fix bond/create command");
iarg += 3;
} else if (strcmp(arg[iarg],"jparam") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix bond/create command");
jmaxbond = force->inumeric(FLERR,arg[iarg+1]);
jnewtype = force->inumeric(FLERR,arg[iarg+2]);
if (jmaxbond < 0) error->all(FLERR,"Illegal fix bond/create command");
if (jnewtype < 1 || jnewtype > atom->ntypes)
error->all(FLERR,"Invalid atom type in fix bond/create command");
iarg += 3;
} else if (strcmp(arg[iarg],"prob") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix bond/create command");
fraction = force->numeric(FLERR,arg[iarg+1]);
seed = force->inumeric(FLERR,arg[iarg+2]);
if (fraction < 0.0 || fraction > 1.0)
error->all(FLERR,"Illegal fix bond/create command");
if (seed <= 0) error->all(FLERR,"Illegal fix bond/create command");
iarg += 3;
} else if (strcmp(arg[iarg],"atype") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/create command");
atype = force->inumeric(FLERR,arg[iarg+1]);
if (atype < 0) error->all(FLERR,"Illegal fix bond/create command");
iarg += 2;
} else if (strcmp(arg[iarg],"dtype") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/create command");
dtype = force->inumeric(FLERR,arg[iarg+1]);
if (dtype < 0) error->all(FLERR,"Illegal fix bond/create command");
iarg += 2;
} else if (strcmp(arg[iarg],"itype") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/create command");
itype = force->inumeric(FLERR,arg[iarg+1]);
if (itype < 0) error->all(FLERR,"Illegal fix bond/create command");
iarg += 2;
} else error->all(FLERR,"Illegal fix bond/create command");
}
// error check
if (atom->molecular != 1)
error->all(FLERR,"Cannot use fix bond/create with non-molecular systems");
if (iatomtype == jatomtype &&
((imaxbond != jmaxbond) || (inewtype != jnewtype)))
error->all(FLERR,
"Inconsistent iparam/jparam values in fix bond/create command");
// initialize Marsaglia RNG with processor-unique seed
random = new RanMars(lmp,seed + me);
// perform initial allocation of atom-based arrays
// register with Atom class
// bondcount values will be initialized in setup()
bondcount = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
countflag = 0;
// set comm sizes needed by this fix
// forward is big due to comm of broken bonds and 1-2 neighbors
comm_forward = MAX(2,2+atom->maxspecial);
comm_reverse = 2;
// allocate arrays local to this fix
nmax = 0;
partner = finalpartner = NULL;
distsq = NULL;
maxcreate = 0;
created = NULL;
// copy = special list for one atom
// size = ms^2 + ms is sufficient
// b/c in rebuild_special() neighs of all 1-2s are added,
// then a dedup(), then neighs of all 1-3s are added, then final dedup()
// this means intermediate size cannot exceed ms^2 + ms
int maxspecial = atom->maxspecial;
copy = new tagint[maxspecial*maxspecial + maxspecial];
// zero out stats
createcount = 0;
createcounttotal = 0;
}
/* ---------------------------------------------------------------------- */
FixBondCreate::~FixBondCreate()
{
// unregister callbacks to this fix from Atom class
atom->delete_callback(id,0);
delete random;
// delete locally stored arrays
memory->destroy(bondcount);
memory->destroy(partner);
memory->destroy(finalpartner);
memory->destroy(distsq);
memory->destroy(created);
delete [] copy;
}
/* ---------------------------------------------------------------------- */
int FixBondCreate::setmask()
{
int mask = 0;
mask |= POST_INTEGRATE;
mask |= POST_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::init()
{
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
// check cutoff for iatomtype,jatomtype
if (force->pair == NULL || cutsq > force->pair->cutsq[iatomtype][jatomtype])
error->all(FLERR,"Fix bond/create cutoff is longer than pairwise cutoff");
// enable angle/dihedral/improper creation if atype/dtype/itype
// option was used and a force field has been specified
if (atype && force->angle) {
angleflag = 1;
if (atype > atom->nangletypes)
error->all(FLERR,"Fix bond/create angle type is invalid");
} else angleflag = 0;
if (dtype && force->dihedral) {
dihedralflag = 1;
if (dtype > atom->ndihedraltypes)
error->all(FLERR,"Fix bond/create dihedral type is invalid");
} else dihedralflag = 0;
if (itype && force->improper) {
improperflag = 1;
if (itype > atom->nimpropertypes)
error->all(FLERR,"Fix bond/create improper type is invalid");
} else improperflag = 0;
if (force->improper) {
if (force->improper_match("class2") || force->improper_match("ring"))
error->all(FLERR,"Cannot yet use fix bond/create with this "
"improper style");
}
// need a half neighbor list, built every Nevery steps
int irequest = neighbor->request((void *) this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->occasional = 1;
lastcheck = -1;
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::init_list(int id, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::setup(int vflag)
{
int i,j,m;
// compute initial bondcount if this is first run
// can't do this earlier, in constructor or init, b/c need ghost info
if (countflag) return;
countflag = 1;
// count bonds stored with each bond I own
// if newton bond is not set, just increment count on atom I
// if newton bond is set, also increment count on atom J even if ghost
// bondcount is long enough to tally ghost atom counts
int *num_bond = atom->num_bond;
int **bond_type = atom->bond_type;
tagint **bond_atom = atom->bond_atom;
int nlocal = atom->nlocal;
int nghost = atom->nghost;
int nall = nlocal + nghost;
int newton_bond = force->newton_bond;
for (i = 0; i < nall; i++) bondcount[i] = 0;
for (i = 0; i < nlocal; i++)
for (j = 0; j < num_bond[i]; j++) {
if (bond_type[i][j] == btype) {
bondcount[i]++;
if (newton_bond) {
m = atom->map(bond_atom[i][j]);
if (m < 0)
error->one(FLERR,"Fix bond/create needs ghost atoms "
"from further away");
bondcount[m]++;
}
}
}
// if newton_bond is set, need to sum bondcount
commflag = 1;
if (newton_bond) comm->reverse_comm_fix(this,1);
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::post_integrate()
{
int i,j,k,m,n,ii,jj,inum,jnum,itype,jtype,n1,n2,n3,possible;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
tagint *slist;
if (update->ntimestep % nevery) return;
// check that all procs have needed ghost atoms within ghost cutoff
// only if neighbor list has changed since last check
// needs to be <= test b/c neighbor list could have been re-built in
// same timestep as last post_integrate() call, but afterwards
// NOTE: no longer think is needed, due to error tests on atom->map()
// NOTE: if delete, can also delete lastcheck and check_ghosts()
//if (lastcheck <= neighbor->lastcall) check_ghosts();
// acquire updated ghost atom positions
// necessary b/c are calling this after integrate, but before Verlet comm
comm->forward_comm();
// forward comm of bondcount, so ghosts have it
commflag = 1;
comm->forward_comm_fix(this,1);
// resize bond partner list and initialize it
// probability array overlays distsq array
// needs to be atom->nmax in length
if (atom->nmax > nmax) {
memory->destroy(partner);
memory->destroy(finalpartner);
memory->destroy(distsq);
nmax = atom->nmax;
memory->create(partner,nmax,"bond/create:partner");
memory->create(finalpartner,nmax,"bond/create:finalpartner");
memory->create(distsq,nmax,"bond/create:distsq");
probability = distsq;
}
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;
for (i = 0; i < nall; i++) {
partner[i] = 0;
finalpartner[i] = 0;
distsq[i] = BIG;
}
// loop over neighbors of my atoms
// each atom sets one closest eligible partner atom ID to bond with
double **x = atom->x;
tagint *tag = atom->tag;
tagint **bond_atom = atom->bond_atom;
int *num_bond = atom->num_bond;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
int *mask = atom->mask;
int *type = atom->type;
neighbor->build_one(list,1);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (!(mask[i] & groupbit)) continue;
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
if (!(mask[j] & groupbit)) continue;
jtype = type[j];
possible = 0;
if (itype == iatomtype && jtype == jatomtype) {
if ((imaxbond == 0 || bondcount[i] < imaxbond) &&
(jmaxbond == 0 || bondcount[j] < jmaxbond))
possible = 1;
} else if (itype == jatomtype && jtype == iatomtype) {
if ((jmaxbond == 0 || bondcount[i] < jmaxbond) &&
(imaxbond == 0 || bondcount[j] < imaxbond))
possible = 1;
}
if (!possible) continue;
// do not allow a duplicate bond to be created
// check 1-2 neighbors of atom I
for (k = 0; k < nspecial[i][0]; k++)
if (special[i][k] == tag[j]) possible = 0;
if (!possible) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq >= cutsq) continue;
if (rsq < distsq[i]) {
partner[i] = tag[j];
distsq[i] = rsq;
}
if (rsq < distsq[j]) {
partner[j] = tag[i];
distsq[j] = rsq;
}
}
}
// reverse comm of distsq and partner
// not needed if newton_pair off since I,J pair was seen by both procs
commflag = 2;
if (force->newton_pair) comm->reverse_comm_fix(this);
// each atom now knows its winning partner
// for prob check, generate random value for each atom with a bond partner
// forward comm of partner and random value, so ghosts have it
if (fraction < 1.0) {
for (i = 0; i < nlocal; i++)
if (partner[i]) probability[i] = random->uniform();
}
commflag = 2;
comm->forward_comm_fix(this,2);
// create bonds for atoms I own
// only if both atoms list each other as winning bond partner
// and probability constraint is satisfied
// if other atom is owned by another proc, it should do same thing
int **bond_type = atom->bond_type;
int newton_bond = force->newton_bond;
ncreate = 0;
for (i = 0; i < nlocal; i++) {
if (partner[i] == 0) continue;
j = atom->map(partner[i]);
if (partner[j] != tag[i]) continue;
// apply probability constraint using RN for atom with smallest ID
if (fraction < 1.0) {
if (tag[i] < tag[j]) {
if (probability[i] >= fraction) continue;
} else {
if (probability[j] >= fraction) continue;
}
}
// if newton_bond is set, only store with I or J
// if not newton_bond, store bond with both I and J
// atom J will also do this consistently, whatever proc it is on
if (!newton_bond || tag[i] < tag[j]) {
if (num_bond[i] == atom->bond_per_atom)
error->one(FLERR,"New bond exceeded bonds per atom in fix bond/create");
bond_type[i][num_bond[i]] = btype;
bond_atom[i][num_bond[i]] = tag[j];
num_bond[i]++;
}
// add a 1-2 neighbor to special bond list for atom I
// atom J will also do this, whatever proc it is on
// need to first remove tag[j] from later in list if it appears
// prevents list from overflowing, will be rebuilt in rebuild_special()
slist = special[i];
n1 = nspecial[i][0];
n2 = nspecial[i][1];
n3 = nspecial[i][2];
for (m = n1; m < n3; m++)
if (slist[m] == tag[j]) break;
if (m < n3) {
for (n = m; n < n3-1; n++) slist[n] = slist[n+1];
n3--;
if (m < n2) n2--;
}
if (n3 == atom->maxspecial)
error->one(FLERR,
"New bond exceeded special list size in fix bond/create");
for (m = n3; m > n1; m--) slist[m] = slist[m-1];
slist[n1] = tag[j];
nspecial[i][0] = n1+1;
nspecial[i][1] = n2+1;
nspecial[i][2] = n3+1;
// increment bondcount, convert atom to new type if limit reached
// atom J will also do this, whatever proc it is on
bondcount[i]++;
if (type[i] == iatomtype) {
if (bondcount[i] == imaxbond) type[i] = inewtype;
} else {
if (bondcount[i] == jmaxbond) type[i] = jnewtype;
}
// store final created bond partners and count the created bond once
finalpartner[i] = tag[j];
finalpartner[j] = tag[i];
if (tag[i] < tag[j]) ncreate++;
}
// tally stats
MPI_Allreduce(&ncreate,&createcount,1,MPI_INT,MPI_SUM,world);
createcounttotal += createcount;
atom->nbonds += createcount;
// trigger reneighboring if any bonds were formed
// this insures neigh lists will immediately reflect the topology changes
// done if any bonds created
if (createcount) next_reneighbor = update->ntimestep;
if (!createcount) return;
// communicate final partner and 1-2 special neighbors
// 1-2 neighs already reflect created bonds
commflag = 3;
comm->forward_comm_fix(this);
// create list of broken bonds that influence my owned atoms
// even if between owned-ghost or ghost-ghost atoms
// finalpartner is now set for owned and ghost atoms so loop over nall
// OK if duplicates in broken list due to ghosts duplicating owned atoms
// check J < 0 to insure a broken bond to unknown atom is included
// i.e. a bond partner outside of cutoff length
ncreate = 0;
for (i = 0; i < nall; i++) {
if (finalpartner[i] == 0) continue;
j = atom->map(finalpartner[i]);
if (j < 0 || tag[i] < tag[j]) {
if (ncreate == maxcreate) {
maxcreate += DELTA;
memory->grow(created,maxcreate,2,"bond/create:created");
}
created[ncreate][0] = tag[i];
created[ncreate][1] = finalpartner[i];
ncreate++;
}
}
// update special neigh lists of all atoms affected by any created bond
// also add angles/dihedrals/impropers induced by created bonds
update_topology();
// DEBUG
//print_bb();
}
/* ----------------------------------------------------------------------
insure all atoms 2 hops away from owned atoms are in ghost list
this allows dihedral 1-2-3-4 to be properly created
and special list of 1 to be properly updated
if I own atom 1, but not 2,3,4, and bond 3-4 is added
then 2,3 will be ghosts and 3 will store 4 as its finalpartner
------------------------------------------------------------------------- */
void FixBondCreate::check_ghosts()
{
int i,j,n;
tagint *slist;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
int nlocal = atom->nlocal;
int flag = 0;
for (i = 0; i < nlocal; i++) {
slist = special[i];
n = nspecial[i][1];
for (j = 0; j < n; j++)
if (atom->map(slist[j]) < 0) flag = 1;
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall)
error->all(FLERR,"Fix3 bond/create needs ghost atoms from further away");
lastcheck = update->ntimestep;
}
/* ----------------------------------------------------------------------
double loop over my atoms and created bonds
influenced = 1 if atom's topology is affected by any created bond
yes if is one of 2 atoms in bond
yes if either atom ID appears in as 1-2 or 1-3 in atom's special list
else no
if influenced by any created bond:
rebuild the atom's special list of 1-2,1-3,1-4 neighs
check for angles/dihedrals/impropers to create due modified special list
------------------------------------------------------------------------- */
void FixBondCreate::update_topology()
{
int i,j,k,n,influence,influenced;
tagint id1,id2;
tagint *slist;
tagint *tag = atom->tag;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
int nlocal = atom->nlocal;
nangles = 0;
ndihedrals = 0;
nimpropers = 0;
overflow = 0;
//printf("NCREATE %d: ",ncreate);
//for (i = 0; i < ncreate; i++)
// printf(" %d %d,",created[i][0],created[i][1]);
//printf("\n");
for (i = 0; i < nlocal; i++) {
influenced = 0;
slist = special[i];
for (j = 0; j < ncreate; j++) {
id1 = created[j][0];
id2 = created[j][1];
influence = 0;
if (tag[i] == id1 || tag[i] == id2) influence = 1;
else {
n = nspecial[i][1];
for (k = 0; k < n; k++)
if (slist[k] == id1 || slist[k] == id2) {
influence = 1;
break;
}
}
if (!influence) continue;
influenced = 1;
}
// rebuild_special first, since used by create_angles, etc
if (influenced) {
rebuild_special(i);
if (angleflag) create_angles(i);
if (dihedralflag) create_dihedrals(i);
if (improperflag) create_impropers(i);
}
}
int overflowall;
MPI_Allreduce(&overflow,&overflowall,1,MPI_INT,MPI_SUM,world);
if (overflowall) error->all(FLERR,"Fix bond/create induced too many "
"angles/dihedrals/impropers per atom");
int newton_bond = force->newton_bond;
int all;
if (angleflag) {
MPI_Allreduce(&nangles,&all,1,MPI_INT,MPI_SUM,world);
if (!newton_bond) all /= 3;
atom->nangles += all;
}
if (dihedralflag) {
MPI_Allreduce(&ndihedrals,&all,1,MPI_INT,MPI_SUM,world);
if (!newton_bond) all /= 4;
atom->ndihedrals += all;
}
if (improperflag) {
MPI_Allreduce(&nimpropers,&all,1,MPI_INT,MPI_SUM,world);
if (!newton_bond) all /= 4;
atom->nimpropers += all;
}
}
/* ----------------------------------------------------------------------
re-build special list of atom M
does not affect 1-2 neighs (already include effects of new bond)
affects 1-3 and 1-4 neighs due to other atom's augmented 1-2 neighs
------------------------------------------------------------------------- */
void FixBondCreate::rebuild_special(int m)
{
int i,j,n,n1,cn1,cn2,cn3;
tagint *slist;
tagint *tag = atom->tag;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
// existing 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]);
if (n < 0)
error->one(FLERR,"Fix bond/create needs ghost atoms from further away");
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);
if (cn2 > atom->maxspecial)
error->one(FLERR,"Special list size exceeded in fix bond/create");
// 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]);
if (n < 0)
error->one(FLERR,"Fix bond/create needs ghost atoms from further away");
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);
if (cn3 > atom->maxspecial)
error->one(FLERR,"Special list size exceeded in fix bond/create");
// 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));
}
/* ----------------------------------------------------------------------
create any angles owned by atom M induced by newly created bonds
walk special list to find all possible angles to create
only add an angle if a new bond is one of its 2 bonds (I-J,J-K)
for newton_bond on, atom M is central atom
for newton_bond off, atom M is any of 3 atoms in angle
------------------------------------------------------------------------- */
void FixBondCreate::create_angles(int m)
{
int i,j,n,i2local,n1,n2;
tagint i1,i2,i3;
tagint *s1list,*s2list;
tagint *tag = atom->tag;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
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];
// atom M is central atom in angle
// double loop over 1-2 neighs
// avoid double counting by 2nd loop as j = i+1,N not j = 1,N
// consider all angles, only add if:
// a new bond is in the angle and atom types match
i2 = tag[m];
n2 = nspecial[m][0];
s2list = special[m];
for (i = 0; i < n2; i++) {
i1 = s2list[i];
for (j = i+1; j < n2; j++) {
i3 = s2list[j];
// angle = i1-i2-i3
for (n = 0; n < ncreate; n++) {
if (created[n][0] == i1 && created[n][1] == i2) break;
if (created[n][0] == i2 && created[n][1] == i1) break;
if (created[n][0] == i2 && created[n][1] == i3) break;
if (created[n][0] == i3 && created[n][1] == i2) break;
}
if (n == ncreate) continue;
// NOTE: this is place to check atom types of i1,i2,i3
if (num_angle < atom->angle_per_atom) {
angle_type[num_angle] = atype;
angle_atom1[num_angle] = i1;
angle_atom2[num_angle] = i2;
angle_atom3[num_angle] = i3;
num_angle++;
nangles++;
} else overflow = 1;
}
}
atom->num_angle[m] = num_angle;
if (force->newton_bond) return;
// for newton_bond off, also consider atom M as atom 1 in angle
i1 = tag[m];
n1 = nspecial[m][0];
s1list = special[m];
for (i = 0; i < n1; i++) {
i2 = s1list[i];
i2local = atom->map(i2);
if (i2local < 0)
error->one(FLERR,"Fix bond/create needs ghost atoms from further away");
s2list = special[i2local];
n2 = nspecial[i2local][0];
for (j = 0; j < n2; j++) {
i3 = s2list[j];
if (i3 == i1) continue;
// angle = i1-i2-i3
for (n = 0; n < ncreate; n++) {
if (created[n][0] == i1 && created[n][1] == i2) break;
if (created[n][0] == i2 && created[n][1] == i1) break;
if (created[n][0] == i2 && created[n][1] == i3) break;
if (created[n][0] == i3 && created[n][1] == i2) break;
}
if (n == ncreate) continue;
// NOTE: this is place to check atom types of i1,i2,i3
if (num_angle < atom->angle_per_atom) {
angle_type[num_angle] = atype;
angle_atom1[num_angle] = i1;
angle_atom2[num_angle] = i2;
angle_atom3[num_angle] = i3;
num_angle++;
nangles++;
} else overflow = 1;
}
}
atom->num_angle[m] = num_angle;
}
/* ----------------------------------------------------------------------
create any dihedrals owned by atom M induced by newly created bonds
walk special list to find all possible dihedrals to create
only add a dihedral if a new bond is one of its 3 bonds (I-J,J-K,K-L)
for newton_bond on, atom M is central atom
for newton_bond off, atom M is any of 4 atoms in dihedral
------------------------------------------------------------------------- */
void FixBondCreate::create_dihedrals(int m)
{
int i,j,k,n,i1local,i2local,i3local,n1,n2,n3;
tagint i1,i2,i3,i4;
tagint *s1list,*s2list,*s3list;
tagint *tag = atom->tag;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
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];
// atom M is 2nd atom in dihedral
// double loop over 1-2 neighs
// two triple loops: one over neighs at each end of triplet
// avoid double counting by 2nd loop as j = i+1,N not j = 1,N
// avoid double counting due to another atom being 2nd atom in same dihedral
// by requiring ID of 2nd atom < ID of 3rd atom
// don't do this if newton bond off since want to double count
// consider all dihedrals, only add if:
// a new bond is in the dihedral and atom types match
i2 = tag[m];
n2 = nspecial[m][0];
s2list = special[m];
for (i = 0; i < n2; i++) {
i1 = s2list[i];
for (j = i+1; j < n2; j++) {
i3 = s2list[j];
if (force->newton_bond && i2 > i3) continue;
i3local = atom->map(i3);
if (i3local < 0)
error->one(FLERR,"Fix bond/create needs ghost atoms from further away");
s3list = special[i3local];
n3 = nspecial[i3local][0];
for (k = 0; k < n3; k++) {
i4 = s3list[k];
if (i4 == i1 || i4 == i2 || i4 == i3) continue;
// dihedral = i1-i2-i3-i4
for (n = 0; n < ncreate; n++) {
if (created[n][0] == i1 && created[n][1] == i2) break;
if (created[n][0] == i2 && created[n][1] == i1) break;
if (created[n][0] == i2 && created[n][1] == i3) break;
if (created[n][0] == i3 && created[n][1] == i2) break;
if (created[n][0] == i3 && created[n][1] == i4) break;
if (created[n][0] == i4 && created[n][1] == i3) break;
}
if (n < ncreate) {
// NOTE: this is place to check atom types of i3,i2,i1,i4
if (num_dihedral < atom->dihedral_per_atom) {
dihedral_type[num_dihedral] = dtype;
dihedral_atom1[num_dihedral] = i1;
dihedral_atom2[num_dihedral] = i2;
dihedral_atom3[num_dihedral] = i3;
dihedral_atom4[num_dihedral] = i4;
num_dihedral++;
ndihedrals++;
} else overflow = 1;
}
}
}
}
for (i = 0; i < n2; i++) {
i1 = s2list[i];
if (force->newton_bond && i2 > i1) continue;
i1local = atom->map(i1);
if (i1local < 0)
error->one(FLERR,"Fix bond/create needs ghost atoms from further away");
s3list = special[i1local];
n3 = nspecial[i1local][0];
for (j = i+1; j < n2; j++) {
i3 = s2list[j];
for (k = 0; k < n3; k++) {
i4 = s3list[k];
if (i4 == i1 || i4 == i2 || i4 == i3) continue;
// dihedral = i3-i2-i1-i4
for (n = 0; n < ncreate; n++) {
if (created[n][0] == i3 && created[n][1] == i2) break;
if (created[n][0] == i2 && created[n][1] == i3) break;
if (created[n][0] == i2 && created[n][1] == i1) break;
if (created[n][0] == i1 && created[n][1] == i2) break;
if (created[n][0] == i1 && created[n][1] == i4) break;
if (created[n][0] == i4 && created[n][1] == i1) break;
}
if (n < ncreate) {
// NOTE: this is place to check atom types of i3,i2,i1,i4
if (num_dihedral < atom->dihedral_per_atom) {
dihedral_type[num_dihedral] = dtype;
dihedral_atom1[num_dihedral] = i3;
dihedral_atom2[num_dihedral] = i2;
dihedral_atom3[num_dihedral] = i1;
dihedral_atom4[num_dihedral] = i4;
num_dihedral++;
ndihedrals++;
} else overflow = 1;
}
}
}
}
atom->num_dihedral[m] = num_dihedral;
if (force->newton_bond) return;
// for newton_bond off, also consider atom M as atom 1 in dihedral
i1 = tag[m];
n1 = nspecial[m][0];
s1list = special[m];
for (i = 0; i < n1; i++) {
i2 = s1list[i];
i2local = atom->map(i2);
if (i2local < 0)
error->one(FLERR,"Fix bond/create needs ghost atoms from further away");
s2list = special[i2local];
n2 = nspecial[i2local][0];
for (j = 0; j < n2; j++) {
i3 = s2list[j];
if (i3 == i1) continue;
i3local = atom->map(i3);
if (i3local < 0)
error->one(FLERR,"Fix bond/create needs ghost atoms from further away");
s3list = special[i3local];
n3 = nspecial[i3local][0];
for (k = 0; k < n3; k++) {
i4 = s3list[k];
if (i4 == i1 || i4 == i2 || i4 == i3) continue;
// dihedral = i1-i2-i3-i4
for (n = 0; n < ncreate; n++) {
if (created[n][0] == i1 && created[n][1] == i2) break;
if (created[n][0] == i2 && created[n][1] == i1) break;
if (created[n][0] == i2 && created[n][1] == i3) break;
if (created[n][0] == i3 && created[n][1] == i2) break;
if (created[n][0] == i3 && created[n][1] == i4) break;
if (created[n][0] == i4 && created[n][1] == i3) break;
}
if (n < ncreate) {
// NOTE: this is place to check atom types of i3,i2,i1,i4
if (num_dihedral < atom->dihedral_per_atom) {
dihedral_type[num_dihedral] = dtype;
dihedral_atom1[num_dihedral] = i1;
dihedral_atom2[num_dihedral] = i2;
dihedral_atom3[num_dihedral] = i3;
dihedral_atom4[num_dihedral] = i4;
num_dihedral++;
ndihedrals++;
} else overflow = 1;
}
}
}
}
}
/* ----------------------------------------------------------------------
create any impropers owned by atom M induced by newly created bonds
walk special list to find all possible impropers to create
only add an improper if a new bond is one of its 3 bonds (I-J,I-K,I-L)
for newton_bond on, atom M is central atom
for newton_bond off, atom M is any of 4 atoms in improper
------------------------------------------------------------------------- */
void FixBondCreate::create_impropers(int m)
{
int i,j,k,n,i1local,n1,n2;
tagint i1,i2,i3,i4;
tagint *s1list,*s2list;
tagint *tag = atom->tag;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
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];
// atom M is central atom in improper
// triple loop over 1-2 neighs
// avoid double counting by 2nd loop as j = i+1,N not j = 1,N
// consider all impropers, only add if:
// a new bond is in the improper and atom types match
i1 = tag[m];
n1 = nspecial[m][0];
s1list = special[m];
for (i = 0; i < n1; i++) {
i2 = s1list[i];
for (j = i+1; j < n1; j++) {
i3 = s1list[j];
for (k = j+1; k < n1; k++) {
i4 = s1list[k];
// improper = i1-i2-i3-i4
for (n = 0; n < ncreate; n++) {
if (created[n][0] == i1 && created[n][1] == i2) break;
if (created[n][0] == i2 && created[n][1] == i1) break;
if (created[n][0] == i1 && created[n][1] == i3) break;
if (created[n][0] == i3 && created[n][1] == i1) break;
if (created[n][0] == i1 && created[n][1] == i4) break;
if (created[n][0] == i4 && created[n][1] == i1) break;
}
if (n == ncreate) continue;
// NOTE: this is place to check atom types of i1,i2,i3,i4
if (num_improper < atom->improper_per_atom) {
improper_type[num_improper] = itype;
improper_atom1[num_improper] = i1;
improper_atom2[num_improper] = i2;
improper_atom3[num_improper] = i3;
improper_atom4[num_improper] = i4;
num_improper++;
nimpropers++;
} else overflow = 1;
}
}
}
atom->num_improper[m] = num_improper;
if (force->newton_bond) return;
// for newton_bond off, also consider atom M as atom 2 in improper
i2 = tag[m];
n2 = nspecial[m][0];
s2list = special[m];
for (i = 0; i < n2; i++) {
i1 = s2list[i];
i1local = atom->map(i1);
if (i1local < 0)
error->one(FLERR,"Fix bond/create needs ghost atoms from further away");
s1list = special[i1local];
n1 = nspecial[i1local][0];
for (j = 0; j < n1; j++) {
i3 = s1list[j];
if (i3 == i1 || i3 == i2) continue;
for (k = j+1; k < n1; k++) {
i4 = s1list[k];
if (i4 == i1 || i4 == i2) continue;
// improper = i1-i2-i3-i4
for (n = 0; n < ncreate; n++) {
if (created[n][0] == i1 && created[n][1] == i2) break;
if (created[n][0] == i2 && created[n][1] == i1) break;
if (created[n][0] == i1 && created[n][1] == i3) break;
if (created[n][0] == i3 && created[n][1] == i1) break;
if (created[n][0] == i1 && created[n][1] == i4) break;
if (created[n][0] == i4 && created[n][1] == i1) break;
}
if (n < ncreate) {
// NOTE: this is place to check atom types of i3,i2,i1,i4
if (num_improper < atom->improper_per_atom) {
improper_type[num_improper] = dtype;
improper_atom1[num_improper] = i1;
improper_atom2[num_improper] = i2;
improper_atom3[num_improper] = i3;
improper_atom4[num_improper] = i4;
num_improper++;
nimpropers++;
} else overflow = 1;
}
}
}
}
}
/* ----------------------------------------------------------------------
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 FixBondCreate::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;
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::post_integrate_respa(int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_integrate();
}
/* ---------------------------------------------------------------------- */
int FixBondCreate::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,k,m,ns;
m = 0;
if (commflag == 1) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = ubuf(bondcount[j]).d;
}
return m;
}
if (commflag == 2) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = ubuf(partner[j]).d;
buf[m++] = probability[j];
}
return m;
}
int **nspecial = atom->nspecial;
tagint **special = atom->special;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = ubuf(finalpartner[j]).d;
ns = nspecial[j][0];
buf[m++] = ubuf(ns).d;
for (k = 0; k < ns; k++)
buf[m++] = ubuf(special[j][k]).d;
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::unpack_forward_comm(int n, int first, double *buf)
{
int i,j,m,ns,last;
m = 0;
last = first + n;
if (commflag == 1) {
for (i = first; i < last; i++)
bondcount[i] = (int) ubuf(buf[m++]).i;
} else if (commflag == 2) {
for (i = first; i < last; i++) {
partner[i] = (tagint) ubuf(buf[m++]).i;
probability[i] = buf[m++];
}
} else {
int **nspecial = atom->nspecial;
tagint **special = atom->special;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
finalpartner[i] = (tagint) ubuf(buf[m++]).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;
}
}
}
/* ---------------------------------------------------------------------- */
int FixBondCreate::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if (commflag == 1) {
for (i = first; i < last; i++)
buf[m++] = ubuf(bondcount[i]).d;
return m;
}
for (i = first; i < last; i++) {
buf[m++] = ubuf(partner[i]).d;
buf[m++] = distsq[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
if (commflag == 1) {
for (i = 0; i < n; i++) {
j = list[i];
bondcount[j] += (int) ubuf(buf[m++]).i;
}
} else {
for (i = 0; i < n; i++) {
j = list[i];
if (buf[m+1] < distsq[j]) {
partner[j] = (tagint) ubuf(buf[m++]).i;
distsq[j] = buf[m++];
} else m += 2;
}
}
}
/* ----------------------------------------------------------------------
allocate local atom-based arrays
------------------------------------------------------------------------- */
void FixBondCreate::grow_arrays(int nmax)
{
memory->grow(bondcount,nmax,"bond/create:bondcount");
}
/* ----------------------------------------------------------------------
copy values within local atom-based arrays
------------------------------------------------------------------------- */
void FixBondCreate::copy_arrays(int i, int j, int delflag)
{
bondcount[j] = bondcount[i];
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for exchange with another proc
------------------------------------------------------------------------- */
int FixBondCreate::pack_exchange(int i, double *buf)
{
buf[0] = bondcount[i];
return 1;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based arrays from exchange with another proc
------------------------------------------------------------------------- */
int FixBondCreate::unpack_exchange(int nlocal, double *buf)
{
bondcount[nlocal] = static_cast<int> (buf[0]);
return 1;
}
/* ---------------------------------------------------------------------- */
double FixBondCreate::compute_vector(int n)
{
if (n == 1) return (double) createcount;
return (double) createcounttotal;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixBondCreate::memory_usage()
{
int nmax = atom->nmax;
double bytes = nmax * sizeof(int);
bytes = 2*nmax * sizeof(tagint);
bytes += nmax * sizeof(double);
return bytes;
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::print_bb()
{
for (int i = 0; i < atom->nlocal; i++) {
printf("TAG " TAGINT_FORMAT ": %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(" " TAGINT_FORMAT,atom->bond_atom[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %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(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT ",",
+ atom->angle_atom1[i][j], atom->angle_atom2[i][j],
+ atom->angle_atom3[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %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],
+ printf(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " "
+ TAGINT_FORMAT ",", atom->dihedral_atom1[i][j],
atom->dihedral_atom2[i][j],atom->dihedral_atom3[i][j],
atom->dihedral_atom4[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %d nimpropers: ",atom->tag[i],atom->num_improper[i]);
for (int j = 0; j < atom->num_improper[i]; j++) {
- printf(" %d %d %d %d,",atom->improper_atom1[i][j],
+ printf(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " "
+ TAGINT_FORMAT ",",atom->improper_atom1[i][j],
atom->improper_atom2[i][j],atom->improper_atom3[i][j],
atom->improper_atom4[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %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(" " TAGINT_FORMAT,atom->special[i][j]);
}
printf("\n");
}
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::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);
+ printf("%s " TAGINT_FORMAT ": %d %d %d nspecial: ",str,m,n1,n2,n3);
for (int j = 0; j < n3; j++) printf(" %d",v[j]);
printf("\n");
}
diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index b014f44ba..9dbf6e826 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -1,3423 +1,3423 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "fix_rigid_small.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "molecule.h"
#include "domain.h"
#include "update.h"
#include "respa.h"
#include "modify.h"
#include "group.h"
#include "comm.h"
#include "force.h"
#include "output.h"
#include "random_mars.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include <map>
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
// allocate space for static class variable
FixRigidSmall *FixRigidSmall::frsptr;
#define MAXLINE 256
#define CHUNK 1024
#define ATTRIBUTE_PERBODY 11
#define TOLERANCE 1.0e-6
#define EPSILON 1.0e-7
#define BIG 1.0e20
#define SINERTIA 0.4 // moment of inertia prefactor for sphere
#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid
#define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment
#define DELTA_BODY 10000
enum{NONE,XYZ,XY,YZ,XZ}; // same as in FixRigid
enum{ISO,ANISO,TRICLINIC}; // same as in FixRigid
enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF};
/* ---------------------------------------------------------------------- */
FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
int i;
scalar_flag = 1;
extscalar = 0;
global_freq = 1;
time_integrate = 1;
rigid_flag = 1;
virial_flag = 1;
create_attribute = 1;
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
// perform initial allocation of atom-based arrays
// register with Atom class
extended = orientflag = dorientflag = 0;
bodyown = NULL;
bodytag = NULL;
atom2body = NULL;
displace = NULL;
eflags = NULL;
orient = NULL;
dorient = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
// parse args for rigid body specification
if (narg < 4) error->all(FLERR,"Illegal fix rigid/small command");
if (strcmp(arg[3],"molecule") != 0)
error->all(FLERR,"Illegal fix rigid/small command");
if (atom->molecule_flag == 0)
error->all(FLERR,"Fix rigid/small requires atom attribute molecule");
if (atom->map_style == 0)
error->all(FLERR,"Fix rigid/small requires an atom map, see atom_modify");
// maxmol = largest molecule #
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int nlocal = atom->nlocal;
maxmol = -1;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) maxmol = MAX(maxmol,molecule[i]);
tagint itmp;
MPI_Allreduce(&maxmol,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
maxmol = itmp;
// parse optional args
int seed;
langflag = 0;
infile = NULL;
onemols = NULL;
tstat_flag = 0;
pstat_flag = 0;
allremap = 1;
id_dilate = NULL;
t_chain = 10;
t_iter = 1;
t_order = 3;
p_chain = 10;
pcouple = NONE;
pstyle = ANISO;
for (int i = 0; i < 3; i++) {
p_start[i] = p_stop[i] = p_period[i] = 0.0;
p_flag[i] = 0;
}
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"langevin") == 0) {
if (iarg+5 > narg) error->all(FLERR,"Illegal fix rigid/small command");
if (strcmp(style,"rigid/small") != 0)
error->all(FLERR,"Illegal fix rigid/small command");
langflag = 1;
t_start = force->numeric(FLERR,arg[iarg+1]);
t_stop = force->numeric(FLERR,arg[iarg+2]);
t_period = force->numeric(FLERR,arg[iarg+3]);
seed = force->inumeric(FLERR,arg[iarg+4]);
if (t_period <= 0.0)
error->all(FLERR,"Fix rigid/small langevin period must be > 0.0");
if (seed <= 0) error->all(FLERR,"Illegal fix rigid/small command");
iarg += 5;
} else if (strcmp(arg[iarg],"infile") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
delete [] infile;
int n = strlen(arg[iarg+1]) + 1;
infile = new char[n];
strcpy(infile,arg[iarg+1]);
restart_file = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
int imol = atom->find_molecule(arg[iarg+1]);
if (imol == -1)
error->all(FLERR,"Molecule template ID for "
"fix rigid/small does not exist");
onemols = &atom->molecules[imol];
nmol = onemols[0]->nset;
iarg += 2;
} else if (strcmp(arg[iarg],"temp") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
if (strcmp(style,"rigid/nvt/small") != 0 &&
strcmp(style,"rigid/npt/small") != 0)
error->all(FLERR,"Illegal fix rigid command");
tstat_flag = 1;
t_start = force->numeric(FLERR,arg[iarg+1]);
t_stop = force->numeric(FLERR,arg[iarg+2]);
t_period = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"iso") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
if (strcmp(style,"rigid/npt/small") != 0 &&
strcmp(style,"rigid/nph/small") != 0)
error->all(FLERR,"Illegal fix rigid/small command");
pcouple = XYZ;
p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
p_period[0] = p_period[1] = p_period[2] =
force->numeric(FLERR,arg[iarg+3]);
p_flag[0] = p_flag[1] = p_flag[2] = 1;
if (domain->dimension == 2) {
p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0;
}
iarg += 4;
} else if (strcmp(arg[iarg],"aniso") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
if (strcmp(style,"rigid/npt/small") != 0 &&
strcmp(style,"rigid/nph/small") != 0)
error->all(FLERR,"Illegal fix rigid/small command");
p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
p_period[0] = p_period[1] = p_period[2] =
force->numeric(FLERR,arg[iarg+3]);
p_flag[0] = p_flag[1] = p_flag[2] = 1;
if (domain->dimension == 2) {
p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0;
}
iarg += 4;
} else if (strcmp(arg[iarg],"x") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
p_start[0] = force->numeric(FLERR,arg[iarg+1]);
p_stop[0] = force->numeric(FLERR,arg[iarg+2]);
p_period[0] = force->numeric(FLERR,arg[iarg+3]);
p_flag[0] = 1;
iarg += 4;
} else if (strcmp(arg[iarg],"y") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
p_start[1] = force->numeric(FLERR,arg[iarg+1]);
p_stop[1] = force->numeric(FLERR,arg[iarg+2]);
p_period[1] = force->numeric(FLERR,arg[iarg+3]);
p_flag[1] = 1;
iarg += 4;
} else if (strcmp(arg[iarg],"z") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
p_start[2] = force->numeric(FLERR,arg[iarg+1]);
p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
p_period[2] = force->numeric(FLERR,arg[iarg+3]);
p_flag[2] = 1;
iarg += 4;
} else if (strcmp(arg[iarg],"couple") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ;
else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY;
else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ;
else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ;
else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE;
else error->all(FLERR,"Illegal fix rigid/small command");
iarg += 2;
} else if (strcmp(arg[iarg],"dilate") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix rigid/small nvt/npt/nph command");
if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
else {
allremap = 0;
delete [] id_dilate;
int n = strlen(arg[iarg+1]) + 1;
id_dilate = new char[n];
strcpy(id_dilate,arg[iarg+1]);
int idilate = group->find(id_dilate);
if (idilate == -1)
error->all(FLERR,"Fix rigid/small nvt/npt/nph dilate group ID "
"does not exist");
}
iarg += 2;
} else if (strcmp(arg[iarg],"tparam") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
if (strcmp(style,"rigid/nvt/small") != 0 &&
strcmp(style,"rigid/npt/small") != 0)
error->all(FLERR,"Illegal fix rigid/small command");
t_chain = force->numeric(FLERR,arg[iarg+1]);
t_iter = force->numeric(FLERR,arg[iarg+2]);
t_order = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"pchain") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
if (strcmp(style,"rigid/npt/small") != 0 &&
strcmp(style,"rigid/nph/small") != 0)
error->all(FLERR,"Illegal fix rigid/small command");
p_chain = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal fix rigid/small command");
}
// error check and further setup for Molecule template
if (onemols) {
for (int i = 0; i < nmol; i++) {
if (onemols[i]->xflag == 0)
error->all(FLERR,"Fix rigid/small molecule must have coordinates");
if (onemols[i]->typeflag == 0)
error->all(FLERR,"Fix rigid/small molecule must have atom types");
// fix rigid/small uses center, masstotal, COM, inertia of molecule
onemols[i]->compute_center();
onemols[i]->compute_mass();
onemols[i]->compute_com();
onemols[i]->compute_inertia();
}
}
// set pstat_flag
pstat_flag = 0;
for (int i = 0; i < 3; i++)
if (p_flag[i]) pstat_flag = 1;
if (pcouple == XYZ || (domain->dimension == 2 && pcouple == XY)) pstyle = ISO;
else pstyle = ANISO;
// create rigid bodies based on molecule ID
// sets bodytag for owned atoms
// body attributes are computed later by setup_bodies()
create_bodies();
// set nlocal_body and allocate bodies I own
tagint *tag = atom->tag;
nlocal_body = nghost_body = 0;
for (i = 0; i < nlocal; i++)
if (bodytag[i] == tag[i]) nlocal_body++;
nmax_body = 0;
while (nmax_body < nlocal_body) nmax_body += DELTA_BODY;
body = (Body *) memory->smalloc(nmax_body*sizeof(Body),
"rigid/small:body");
// set bodyown for owned atoms
nlocal_body = 0;
for (i = 0; i < nlocal; i++)
if (bodytag[i] == tag[i]) {
body[nlocal_body].ilocal = i;
bodyown[i] = nlocal_body++;
} else bodyown[i] = -1;
// bodysize = sizeof(Body) in doubles
bodysize = sizeof(Body)/sizeof(double);
if (bodysize*sizeof(double) != sizeof(Body)) bodysize++;
// set max comm sizes needed by this fix
comm_forward = 1 + bodysize;
comm_reverse = 6;
// bitmasks for properties of extended particles
POINT = 1;
SPHERE = 2;
ELLIPSOID = 4;
LINE = 8;
TRIANGLE = 16;
DIPOLE = 32;
OMEGA = 64;
ANGMOM = 128;
TORQUE = 256;
MINUSPI = -MY_PI;
TWOPI = 2.0*MY_PI;
// atom style pointers to particles that store extra info
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec_line = (AtomVecLine *) atom->style_match("line");
avec_tri = (AtomVecTri *) atom->style_match("tri");
// print statistics
int one = 0;
bigint atomone = 0;
for (int i = 0; i < nlocal; i++) {
if (bodyown[i] >= 0) one++;
if (bodytag[i] > 0) atomone++;
}
MPI_Allreduce(&one,&nbody,1,MPI_INT,MPI_SUM,world);
bigint atomall;
MPI_Allreduce(&atomone,&atomall,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (me == 0) {
if (screen) {
fprintf(screen,"%d rigid bodies with " BIGINT_FORMAT " atoms\n",
nbody,atomall);
fprintf(screen," %g = max distance from body owner to body atom\n",
maxextent);
}
if (logfile) {
fprintf(logfile,"%d rigid bodies with " BIGINT_FORMAT " atoms\n",
nbody,atomall);
fprintf(logfile," %g = max distance from body owner to body atom\n",
maxextent);
}
}
// initialize Marsaglia RNG with processor-unique seed
maxlang = 0;
langextra = NULL;
random = NULL;
if (langflag) random = new RanMars(lmp,seed + comm->me);
// mass vector for granular pair styles
mass_body = NULL;
nmax_mass = 0;
staticflag = 1;
}
/* ---------------------------------------------------------------------- */
FixRigidSmall::~FixRigidSmall()
{
// unregister callbacks to this fix from Atom class
atom->delete_callback(id,0);
// delete locally stored arrays
memory->sfree(body);
memory->destroy(bodyown);
memory->destroy(bodytag);
memory->destroy(atom2body);
memory->destroy(displace);
memory->destroy(eflags);
memory->destroy(orient);
memory->destroy(dorient);
delete random;
delete [] infile;
memory->destroy(langextra);
memory->destroy(mass_body);
}
/* ---------------------------------------------------------------------- */
int FixRigidSmall::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
if (langflag) mask |= POST_FORCE;
mask |= PRE_NEIGHBOR;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::init()
{
int i;
triclinic = domain->triclinic;
// warn if more than one rigid fix
int count = 0;
for (i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"rigid") == 0) count++;
if (count > 1 && me == 0) error->warning(FLERR,"More than one fix rigid");
// error if npt,nph fix comes before rigid fix
for (i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix[i]->style,"npt") == 0) break;
if (strcmp(modify->fix[i]->style,"nph") == 0) break;
}
if (i < modify->nfix) {
for (int j = i; j < modify->nfix; j++)
if (strcmp(modify->fix[j]->style,"rigid") == 0)
error->all(FLERR,"Rigid fix must come before NPT/NPH fix");
}
// timestep info
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dtq = 0.5 * update->dt;
if (strstr(update->integrate_style,"respa"))
step_respa = ((Respa *) update->integrate)->step;
}
/* ----------------------------------------------------------------------
setup static/dynamic properties of rigid bodies, using current atom info
allows resetting of atom properties like mass between runs
only do static initialization once if using an infile
cannot do this until now, b/c requires comm->setup() to have setup stencil
------------------------------------------------------------------------- */
void FixRigidSmall::setup_pre_neighbor()
{
if (!staticflag || !infile) setup_bodies_static();
setup_bodies_dynamic();
}
/* ----------------------------------------------------------------------
compute initial fcm and torque on bodies, also initial virial
reset all particle velocities to be consistent with vcm and omega
------------------------------------------------------------------------- */
void FixRigidSmall::setup(int vflag)
{
int i,n,ibody;
//check(1);
// sum fcm, torque across all rigid bodies
// fcm = force on COM
// torque = torque around COM
double **x = atom->x;
double **f = atom->f;
imageint *image = atom->image;
int nlocal = atom->nlocal;
double *xcm,*fcm,*tcm;
double dx,dy,dz;
double unwrap[3];
for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
fcm = body[ibody].fcm;
fcm[0] = fcm[1] = fcm[2] = 0.0;
tcm = body[ibody].torque;
tcm[0] = tcm[1] = tcm[2] = 0.0;
}
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
fcm = b->fcm;
fcm[0] += f[i][0];
fcm[1] += f[i][1];
fcm[2] += f[i][2];
domain->unmap(x[i],image[i],unwrap);
xcm = b->xcm;
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
tcm = b->torque;
tcm[0] += dy * f[i][2] - dz * f[i][1];
tcm[1] += dz * f[i][0] - dx * f[i][2];
tcm[2] += dx * f[i][1] - dy * f[i][0];
}
// extended particles add their rotation/torque to angmom/torque of body
if (extended) {
double **torque = atom->torque;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
if (eflags[i] & TORQUE) {
tcm = b->torque;
tcm[0] += torque[i][0];
tcm[1] += torque[i][1];
tcm[2] += torque[i][2];
}
}
}
// reverse communicate fcm, torque of all bodies
commflag = FORCE_TORQUE;
comm->reverse_comm_fix(this,6);
// virial setup before call to set_v
if (vflag) v_setup(vflag);
else evflag = 0;
// compute and forward communicate vcm and omega of all bodies
for (ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
b->ez_space,b->inertia,b->omega);
}
commflag = FINAL;
comm->forward_comm_fix(this,10);
// set velocity/rotation of atoms in rigid bodues
set_v();
// guesstimate virial as 2x the set_v contribution
if (vflag_global)
for (n = 0; n < 6; n++) virial[n] *= 2.0;
if (vflag_atom) {
for (i = 0; i < nlocal; i++)
for (n = 0; n < 6; n++)
vatom[i][n] *= 2.0;
}
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::initial_integrate(int vflag)
{
double dtfm;
//check(2);
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
// update vcm by 1/2 step
dtfm = dtf / b->mass;
b->vcm[0] += dtfm * b->fcm[0];
b->vcm[1] += dtfm * b->fcm[1];
b->vcm[2] += dtfm * b->fcm[2];
// update xcm by full step
b->xcm[0] += dtv * b->vcm[0];
b->xcm[1] += dtv * b->vcm[1];
b->xcm[2] += dtv * b->vcm[2];
// update angular momentum by 1/2 step
b->angmom[0] += dtf * b->torque[0];
b->angmom[1] += dtf * b->torque[1];
b->angmom[2] += dtf * b->torque[2];
// compute omega at 1/2 step from angmom at 1/2 step and current q
// update quaternion a full step via Richardson iteration
// returns new normalized quaternion, also updated omega at 1/2 step
// update ex,ey,ez to reflect new quaternion
MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
b->ez_space,b->inertia,b->omega);
MathExtra::richardson(b->quat,b->angmom,b->omega,b->inertia,dtq);
MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space,b->ez_space);
}
// virial setup before call to set_xv
if (vflag) v_setup(vflag);
else evflag = 0;
// forward communicate updated info of all bodies
commflag = INITIAL;
comm->forward_comm_fix(this,26);
// set coords/orient and velocity/rotation of atoms in rigid bodies
set_xv();
}
/* ----------------------------------------------------------------------
apply Langevin thermostat to all 6 DOF of rigid bodies I own
unlike fix langevin, this stores extra force in extra arrays,
which are added in when final_integrate() calculates a new fcm/torque
------------------------------------------------------------------------- */
void FixRigidSmall::post_force(int vflag)
{
double gamma1,gamma2;
// grow langextra if needed
if (nlocal_body > maxlang) {
memory->destroy(langextra);
maxlang = nlocal_body + nghost_body;
memory->create(langextra,maxlang,6,"rigid/small:langextra");
}
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
double t_target = t_start + delta * (t_stop-t_start);
double tsqrt = sqrt(t_target);
double boltz = force->boltz;
double dt = update->dt;
double mvv2e = force->mvv2e;
double ftm2v = force->ftm2v;
double *vcm,*omega,*inertia;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
vcm = body[ibody].vcm;
omega = body[ibody].omega;
inertia = body[ibody].inertia;
gamma1 = -body[ibody].mass / t_period / ftm2v;
gamma2 = sqrt(body[ibody].mass) * tsqrt *
sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
langextra[ibody][0] = gamma1*vcm[0] + gamma2*(random->uniform()-0.5);
langextra[ibody][1] = gamma1*vcm[1] + gamma2*(random->uniform()-0.5);
langextra[ibody][2] = gamma1*vcm[2] + gamma2*(random->uniform()-0.5);
gamma1 = -1.0 / t_period / ftm2v;
gamma2 = tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
langextra[ibody][3] = inertia[0]*gamma1*omega[0] +
sqrt(inertia[0])*gamma2*(random->uniform()-0.5);
langextra[ibody][4] = inertia[1]*gamma1*omega[1] +
sqrt(inertia[1])*gamma2*(random->uniform()-0.5);
langextra[ibody][5] = inertia[2]*gamma1*omega[2] +
sqrt(inertia[2])*gamma2*(random->uniform()-0.5);
}
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::final_integrate()
{
int i,ibody;
double dtfm;
//check(3);
// sum over atoms to get force and torque on rigid body
imageint *image = atom->image;
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
double dx,dy,dz;
double unwrap[3];
double *xcm,*fcm,*tcm;
for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
fcm = body[ibody].fcm;
fcm[0] = fcm[1] = fcm[2] = 0.0;
tcm = body[ibody].torque;
tcm[0] = tcm[1] = tcm[2] = 0.0;
}
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
fcm = b->fcm;
fcm[0] += f[i][0];
fcm[1] += f[i][1];
fcm[2] += f[i][2];
domain->unmap(x[i],image[i],unwrap);
xcm = b->xcm;
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
tcm = b->torque;
tcm[0] += dy*f[i][2] - dz*f[i][1];
tcm[1] += dz*f[i][0] - dx*f[i][2];
tcm[2] += dx*f[i][1] - dy*f[i][0];
}
// extended particles add their torque to torque of body
if (extended) {
double **torque = atom->torque;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
if (eflags[i] & TORQUE) {
tcm = body[atom2body[i]].torque;
tcm[0] += torque[i][0];
tcm[1] += torque[i][1];
tcm[2] += torque[i][2];
}
}
}
// reverse communicate fcm, torque of all bodies
commflag = FORCE_TORQUE;
comm->reverse_comm_fix(this,6);
// include Langevin thermostat forces and torques
if (langflag) {
for (int ibody = 0; ibody < nlocal_body; ibody++) {
fcm = body[ibody].fcm;
fcm[0] += langextra[ibody][0];
fcm[1] += langextra[ibody][1];
fcm[2] += langextra[ibody][2];
tcm = body[ibody].torque;
tcm[0] += langextra[ibody][3];
tcm[1] += langextra[ibody][4];
tcm[2] += langextra[ibody][5];
}
}
// update vcm and angmom, recompute omega
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
// update vcm by 1/2 step
dtfm = dtf / b->mass;
b->vcm[0] += dtfm * b->fcm[0];
b->vcm[1] += dtfm * b->fcm[1];
b->vcm[2] += dtfm * b->fcm[2];
// update angular momentum by 1/2 step
b->angmom[0] += dtf * b->torque[0];
b->angmom[1] += dtf * b->torque[1];
b->angmom[2] += dtf * b->torque[2];
MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
b->ez_space,b->inertia,b->omega);
}
// forward communicate updated info of all bodies
commflag = FINAL;
comm->forward_comm_fix(this,10);
// set velocity/rotation of atoms in rigid bodies
// virial is already setup from initial_integrate
set_v();
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::initial_integrate_respa(int vflag, int ilevel, int iloop)
{
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dtq = 0.5 * step_respa[ilevel];
if (ilevel == 0) initial_integrate(vflag);
else final_integrate();
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::final_integrate_respa(int ilevel, int iloop)
{
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
final_integrate();
}
/* ----------------------------------------------------------------------
atom to processor assignments have changed,
so acquire ghost bodies and setup atom2body
remap xcm of each rigid body back into periodic simulation box
done during pre_neighbor so will be after call to pbc()
and after fix_deform::pre_exchange() may have flipped box
use domain->remap() in case xcm is far away from box
due to 1st definition of rigid body or due to box flip
if don't do this, then atoms of a body which drifts far away
from a triclinic box will be remapped back into box
with huge displacements when the box tilt changes via set_x()
adjust image flag of body and image flags of all atoms in body
------------------------------------------------------------------------- */
void FixRigidSmall::pre_neighbor()
{
// remap xcm and image flags of each body as needed
imageint original,oldimage,newimage;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
original = b->image;
domain->remap(b->xcm,b->image);
if (original == b->image) b->remapflag[3] = 0;
else {
oldimage = original & IMGMASK;
newimage = b->image & IMGMASK;
b->remapflag[0] = newimage - oldimage;
oldimage = (original >> IMGBITS) & IMGMASK;
newimage = (b->image >> IMGBITS) & IMGMASK;
b->remapflag[1] = newimage - oldimage;
oldimage = original >> IMG2BITS;
newimage = b->image >> IMG2BITS;
b->remapflag[2] = newimage - oldimage;
b->remapflag[3] = 1;
}
}
// acquire ghost bodies via forward comm
// also gets new remapflags needed for adjusting atom image flags
// reset atom2body for owned atoms
// forward comm sets atom2body for ghost atoms
nghost_body = 0;
commflag = FULL_BODY;
comm->forward_comm_fix(this);
reset_atom2body();
//check(4);
// adjust image flags of any atom in a rigid body whose xcm was remapped
imageint *image = atom->image;
int nlocal = atom->nlocal;
imageint idim,otherdims;
for (int i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
if (b->remapflag[3] == 0) continue;
if (b->remapflag[0]) {
idim = image[i] & IMGMASK;
otherdims = image[i] ^ idim;
idim -= b->remapflag[0];
idim &= IMGMASK;
image[i] = otherdims | idim;
}
if (b->remapflag[1]) {
idim = (image[i] >> IMGBITS) & IMGMASK;
otherdims = image[i] ^ (idim << IMGBITS);
idim -= b->remapflag[1];
idim &= IMGMASK;
image[i] = otherdims | (idim << IMGBITS);
}
if (b->remapflag[2]) {
idim = image[i] >> IMG2BITS;
otherdims = image[i] ^ (idim << IMG2BITS);
idim -= b->remapflag[2];
idim &= IMGMASK;
image[i] = otherdims | (idim << IMG2BITS);
}
}
}
/* ----------------------------------------------------------------------
count # of DOF removed by rigid bodies for atoms in igroup
return total count of DOF
------------------------------------------------------------------------- */
int FixRigidSmall::dof(int tgroup)
{
int i,j;
// cannot count DOF correctly unless setup_bodies_static() has been called
if (!staticflag) {
if (comm->me == 0)
error->warning(FLERR,"Cannot count rigid body degrees-of-freedom "
"before bodies are fully initialized");
return 0;
}
int tgroupbit = group->bitmask[tgroup];
// counts = 3 values per rigid body I own
// 0 = # of point particles in rigid body and in temperature group
// 1 = # of finite-size particles in rigid body and in temperature group
// 2 = # of particles in rigid body, disregarding temperature group
memory->create(counts,nlocal_body+nghost_body,3,"rigid/small:counts");
for (int i = 0; i < nlocal_body+nghost_body; i++)
counts[i][0] = counts[i][1] = counts[i][2] = 0;
// tally counts from my owned atoms
// 0 = # of point particles in rigid body and in temperature group
// 1 = # of finite-size particles in rigid body and in temperature group
// 2 = # of particles in rigid body, disregarding temperature group
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
j = atom2body[i];
counts[j][2]++;
if (mask[i] & tgroupbit) {
if (extended && eflags[i]) counts[j][1]++;
else counts[j][0]++;
}
}
commflag = DOF;
comm->reverse_comm_fix(this,3);
// nall = count0 = # of point particles in each rigid body
// mall = count1 = # of finite-size particles in each rigid body
// warn if nall+mall != nrigid for any body included in temperature group
int flag = 0;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
if (counts[ibody][0]+counts[ibody][1] > 0 &&
counts[ibody][0]+counts[ibody][1] != counts[ibody][2]) flag = 1;
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
if (flagall && me == 0)
error->warning(FLERR,"Computing temperature of portions of rigid bodies");
// remove appropriate DOFs for each rigid body wholly in temperature group
// N = # of point particles in body
// M = # of finite-size particles in body
// 3d body has 3N + 6M dof to start with
// 2d body has 2N + 3M dof to start with
// 3d point-particle body with all non-zero I should have 6 dof, remove 3N-6
// 3d point-particle body (linear) with a 0 I should have 5 dof, remove 3N-5
// 2d point-particle body should have 3 dof, remove 2N-3
// 3d body with any finite-size M should have 6 dof, remove (3N+6M) - 6
// 2d body with any finite-size M should have 3 dof, remove (2N+3M) - 3
double *inertia;
int n = 0;
if (domain->dimension == 3) {
for (int ibody = 0; ibody < nlocal_body; ibody++) {
if (counts[ibody][0]+counts[ibody][1] == counts[ibody][2]) {
n += 3*counts[ibody][0] + 6*counts[ibody][1] - 6;
inertia = body[ibody].inertia;
if (inertia[0] == 0.0 || inertia[1] == 0.0 || inertia[2] == 0.0) n++;
}
}
} else if (domain->dimension == 2) {
for (int ibody = 0; ibody < nlocal_body; ibody++)
if (counts[ibody][0]+counts[ibody][1] == counts[ibody][2])
n += 2*counts[ibody][0] + 3*counts[ibody][1] - 3;
}
memory->destroy(counts);
int nall;
MPI_Allreduce(&n,&nall,1,MPI_INT,MPI_SUM,world);
return nall;
}
/* ----------------------------------------------------------------------
adjust xcm of each rigid body due to box deformation
called by various fixes that change box size/shape
flag = 0/1 means map from box to lamda coords or vice versa
------------------------------------------------------------------------- */
void FixRigidSmall::deform(int flag)
{
if (flag == 0)
for (int ibody = 0; ibody < nlocal_body; ibody++)
domain->x2lamda(body[ibody].xcm,body[ibody].xcm);
else
for (int ibody = 0; ibody < nlocal_body; ibody++)
domain->lamda2x(body[ibody].xcm,body[ibody].xcm);
}
/* ----------------------------------------------------------------------
set space-frame coords and velocity of each atom in each rigid body
set orientation and rotation of extended particles
x = Q displace + Xcm, mapped back to periodic box
v = Vcm + (W cross (x - Xcm))
------------------------------------------------------------------------- */
void FixRigidSmall::set_xv()
{
int xbox,ybox,zbox;
double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
double ione[3],exone[3],eyone[3],ezone[3],vr[6],p[3][3];
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double xy = domain->xy;
double xz = domain->xz;
double yz = domain->yz;
imageint *image = atom->image;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
// set x and v of each atom
for (int i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
// save old positions and velocities for virial
if (evflag) {
if (triclinic == 0) {
x0 = x[i][0] + xbox*xprd;
x1 = x[i][1] + ybox*yprd;
x2 = x[i][2] + zbox*zprd;
} else {
x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
x1 = x[i][1] + ybox*yprd + zbox*yz;
x2 = x[i][2] + zbox*zprd;
}
v0 = v[i][0];
v1 = v[i][1];
v2 = v[i][2];
}
// x = displacement from center-of-mass, based on body orientation
// v = vcm + omega around center-of-mass
MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,displace[i],x[i]);
v[i][0] = b->omega[1]*x[i][2] - b->omega[2]*x[i][1] + b->vcm[0];
v[i][1] = b->omega[2]*x[i][0] - b->omega[0]*x[i][2] + b->vcm[1];
v[i][2] = b->omega[0]*x[i][1] - b->omega[1]*x[i][0] + b->vcm[2];
// add center of mass to displacement
// map back into periodic box via xbox,ybox,zbox
// for triclinic, add in box tilt factors as well
if (triclinic == 0) {
x[i][0] += b->xcm[0] - xbox*xprd;
x[i][1] += b->xcm[1] - ybox*yprd;
x[i][2] += b->xcm[2] - zbox*zprd;
} else {
x[i][0] += b->xcm[0] - xbox*xprd - ybox*xy - zbox*xz;
x[i][1] += b->xcm[1] - ybox*yprd - zbox*yz;
x[i][2] += b->xcm[2] - zbox*zprd;
}
// virial = unwrapped coords dotted into body constraint force
// body constraint force = implied force due to v change minus f external
// assume f does not include forces internal to body
// 1/2 factor b/c final_integrate contributes other half
// assume per-atom contribution is due to constraint force on that atom
if (evflag) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
vr[0] = 0.5*x0*fc0;
vr[1] = 0.5*x1*fc1;
vr[2] = 0.5*x2*fc2;
vr[3] = 0.5*x0*fc1;
vr[4] = 0.5*x0*fc2;
vr[5] = 0.5*x1*fc2;
v_tally(1,&i,1.0,vr);
}
}
// set orientation, omega, angmom of each extended particle
if (extended) {
double theta_body,theta;
double *shape,*quatatom,*inertiaatom;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
AtomVecLine::Bonus *lbonus;
if (avec_line) lbonus = avec_line->bonus;
AtomVecTri::Bonus *tbonus;
if (avec_tri) tbonus = avec_tri->bonus;
double **omega = atom->omega;
double **angmom = atom->angmom;
double **mu = atom->mu;
int *ellipsoid = atom->ellipsoid;
int *line = atom->line;
int *tri = atom->tri;
for (int i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
if (eflags[i] & SPHERE) {
omega[i][0] = b->omega[0];
omega[i][1] = b->omega[1];
omega[i][2] = b->omega[2];
} else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::quatquat(b->quat,orient[i],quatatom);
MathExtra::qnormalize(quatatom);
ione[0] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,ione,angmom[i]);
} else if (eflags[i] & LINE) {
if (b->quat[3] >= 0.0) theta_body = 2.0*acos(b->quat[0]);
else theta_body = -2.0*acos(b->quat[0]);
theta = orient[i][0] + theta_body;
while (theta <= MINUSPI) theta += TWOPI;
while (theta > MY_PI) theta -= TWOPI;
lbonus[line[i]].theta = theta;
omega[i][0] = b->omega[0];
omega[i][1] = b->omega[1];
omega[i][2] = b->omega[2];
} else if (eflags[i] & TRIANGLE) {
inertiaatom = tbonus[tri[i]].inertia;
quatatom = tbonus[tri[i]].quat;
MathExtra::quatquat(b->quat,orient[i],quatatom);
MathExtra::qnormalize(quatatom);
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,
inertiaatom,angmom[i]);
}
if (eflags[i] & DIPOLE) {
MathExtra::quat_to_mat(b->quat,p);
MathExtra::matvec(p,dorient[i],mu[i]);
MathExtra::snormalize3(mu[i][3],mu[i],mu[i]);
}
}
}
}
/* ----------------------------------------------------------------------
set space-frame velocity of each atom in a rigid body
set omega and angmom of extended particles
v = Vcm + (W cross (x - Xcm))
------------------------------------------------------------------------- */
void FixRigidSmall::set_v()
{
int xbox,ybox,zbox;
double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
double ione[3],exone[3],eyone[3],ezone[3],delta[3],vr[6];
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double xy = domain->xy;
double xz = domain->xz;
double yz = domain->yz;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
imageint *image = atom->image;
int nlocal = atom->nlocal;
// set v of each atom
for (int i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,displace[i],delta);
// save old velocities for virial
if (evflag) {
v0 = v[i][0];
v1 = v[i][1];
v2 = v[i][2];
}
v[i][0] = b->omega[1]*delta[2] - b->omega[2]*delta[1] + b->vcm[0];
v[i][1] = b->omega[2]*delta[0] - b->omega[0]*delta[2] + b->vcm[1];
v[i][2] = b->omega[0]*delta[1] - b->omega[1]*delta[0] + b->vcm[2];
// virial = unwrapped coords dotted into body constraint force
// body constraint force = implied force due to v change minus f external
// assume f does not include forces internal to body
// 1/2 factor b/c initial_integrate contributes other half
// assume per-atom contribution is due to constraint force on that atom
if (evflag) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if (triclinic == 0) {
x0 = x[i][0] + xbox*xprd;
x1 = x[i][1] + ybox*yprd;
x2 = x[i][2] + zbox*zprd;
} else {
x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
x1 = x[i][1] + ybox*yprd + zbox*yz;
x2 = x[i][2] + zbox*zprd;
}
vr[0] = 0.5*x0*fc0;
vr[1] = 0.5*x1*fc1;
vr[2] = 0.5*x2*fc2;
vr[3] = 0.5*x0*fc1;
vr[4] = 0.5*x0*fc2;
vr[5] = 0.5*x1*fc2;
v_tally(1,&i,1.0,vr);
}
}
// set omega, angmom of each extended particle
if (extended) {
double *shape,*quatatom,*inertiaatom;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
AtomVecTri::Bonus *tbonus;
if (avec_tri) tbonus = avec_tri->bonus;
double **omega = atom->omega;
double **angmom = atom->angmom;
int *ellipsoid = atom->ellipsoid;
int *tri = atom->tri;
for (int i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
if (eflags[i] & SPHERE) {
omega[i][0] = b->omega[0];
omega[i][1] = b->omega[1];
omega[i][2] = b->omega[2];
} else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
ione[0] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,ione,
angmom[i]);
} else if (eflags[i] & LINE) {
omega[i][0] = b->omega[0];
omega[i][1] = b->omega[1];
omega[i][2] = b->omega[2];
} else if (eflags[i] & TRIANGLE) {
inertiaatom = tbonus[tri[i]].inertia;
quatatom = tbonus[tri[i]].quat;
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,
inertiaatom,angmom[i]);
}
}
}
}
/* ----------------------------------------------------------------------
one-time identification of which atoms are in which rigid bodies
set bodytag for all owned atoms
------------------------------------------------------------------------- */
void FixRigidSmall::create_bodies()
{
int i,m,n;
double unwrap[3];
// error check on image flags of atoms in rigid bodies
imageint *image = atom->image;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int *periodicity = domain->periodicity;
int xbox,ybox,zbox;
int flag = 0;
for (i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if ((xbox && !periodicity[0]) || (ybox && !periodicity[1]) ||
(zbox && !periodicity[2])) flag = 1;
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall) error->all(FLERR,"Fix rigid/small atom has non-zero image flag "
"in a non-periodic dimension");
// allocate buffer for passing messages around ring of procs
// percount = max number of values to put in buffer for each of ncount
int ncount = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) ncount++;
int percount = 5;
double *buf;
memory->create(buf,ncount*percount,"rigid/small:buf");
// create map hash for storing unique molecule IDs of my atoms
// key = molecule ID
// value = index into per-body data structure
// n = # of entries in hash
hash = new std::map<tagint,int>();
hash->clear();
// setup hash
// key = body ID
// value = index into N-length data structure
// n = count of unique bodies my atoms are part of
tagint *molecule = atom->molecule;
n = 0;
for (i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
if (hash->find(molecule[i]) == hash->end()) (*hash)[molecule[i]] = n++;
}
// bbox = bounding box of each rigid body my atoms are part of
memory->create(bbox,n,6,"rigid/small:bbox");
for (i = 0; i < n; i++) {
bbox[i][0] = bbox[i][2] = bbox[i][4] = BIG;
bbox[i][1] = bbox[i][3] = bbox[i][5] = -BIG;
}
// pack my atoms into buffer as molecule ID, unwrapped coords
double **x = atom->x;
m = 0;
for (i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
domain->unmap(x[i],image[i],unwrap);
buf[m++] = molecule[i];
buf[m++] = unwrap[0];
buf[m++] = unwrap[1];
buf[m++] = unwrap[2];
}
// pass buffer around ring of procs
// func = update bbox with atom coords from every proc
// when done, have full bbox for every rigid body my atoms are part of
frsptr = this;
comm->ring(m,sizeof(double),buf,1,ring_bbox,NULL);
// ctr = center pt of each rigid body my atoms are part of
memory->create(ctr,n,6,"rigid/small:bbox");
for (i = 0; i < n; i++) {
ctr[i][0] = 0.5 * (bbox[i][0] + bbox[i][1]);
ctr[i][1] = 0.5 * (bbox[i][2] + bbox[i][3]);
ctr[i][2] = 0.5 * (bbox[i][4] + bbox[i][5]);
}
// idclose = ID of atom in body closest to center pt (smaller ID if tied)
// rsqclose = distance squared from idclose to center pt
memory->create(idclose,n,"rigid/small:idclose");
memory->create(rsqclose,n,"rigid/small:rsqclose");
for (i = 0; i < n; i++) rsqclose[i] = BIG;
// pack my atoms into buffer as molecule ID, atom ID, unwrapped coords
tagint *tag = atom->tag;
m = 0;
for (i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
domain->unmap(x[i],image[i],unwrap);
buf[m++] = molecule[i];
buf[m++] = ubuf(tag[i]).d;
buf[m++] = unwrap[0];
buf[m++] = unwrap[1];
buf[m++] = unwrap[2];
}
// pass buffer around ring of procs
// func = update idclose,rsqclose with atom IDs from every proc
// when done, have idclose for every rigid body my atoms are part of
frsptr = this;
comm->ring(m,sizeof(double),buf,2,ring_nearest,NULL);
// set bodytag of all owned atoms, based on idclose
// find max value of rsqclose across all procs
double rsqmax = 0.0;
for (i = 0; i < nlocal; i++) {
bodytag[i] = 0;
if (!(mask[i] & groupbit)) continue;
m = hash->find(molecule[i])->second;
bodytag[i] = idclose[m];
rsqmax = MAX(rsqmax,rsqclose[m]);
}
// pack my atoms into buffer as bodytag of owning atom, unwrapped coords
m = 0;
for (i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
domain->unmap(x[i],image[i],unwrap);
buf[m++] = ubuf(bodytag[i]).d;
buf[m++] = unwrap[0];
buf[m++] = unwrap[1];
buf[m++] = unwrap[2];
}
// pass buffer around ring of procs
// func = update rsqfar for atoms belonging to bodies I own
// when done, have rsqfar for all atoms in bodies I own
rsqfar = 0.0;
frsptr = this;
comm->ring(m,sizeof(double),buf,3,ring_farthest,NULL);
// find maxextent of rsqfar across all procs
// if defined, include molecule->maxextent
MPI_Allreduce(&rsqfar,&maxextent,1,MPI_DOUBLE,MPI_MAX,world);
maxextent = sqrt(maxextent);
if (onemols) {
for (int i = 0; i < nmol; i++)
maxextent = MAX(maxextent,onemols[i]->maxextent);
}
// clean up
delete hash;
memory->destroy(buf);
memory->destroy(bbox);
memory->destroy(ctr);
memory->destroy(idclose);
memory->destroy(rsqclose);
}
/* ----------------------------------------------------------------------
process rigid body atoms from another proc
update bounding box for rigid bodies my atoms are part of
------------------------------------------------------------------------- */
void FixRigidSmall::ring_bbox(int n, char *cbuf)
{
std::map<tagint,int> *hash = frsptr->hash;
double **bbox = frsptr->bbox;
double *buf = (double *) cbuf;
int ndatums = n/4;
int j,imol;
double *x;
int m = 0;
for (int i = 0; i < ndatums; i++, m += 4) {
imol = static_cast<int> (buf[m]);
if (hash->find(imol) != hash->end()) {
j = hash->find(imol)->second;
x = &buf[m+1];
bbox[j][0] = MIN(bbox[j][0],x[0]);
bbox[j][1] = MAX(bbox[j][1],x[0]);
bbox[j][2] = MIN(bbox[j][2],x[1]);
bbox[j][3] = MAX(bbox[j][3],x[1]);
bbox[j][4] = MIN(bbox[j][4],x[2]);
bbox[j][5] = MAX(bbox[j][5],x[2]);
}
}
}
/* ----------------------------------------------------------------------
process rigid body atoms from another proc
update nearest atom to body center for rigid bodies my atoms are part of
------------------------------------------------------------------------- */
void FixRigidSmall::ring_nearest(int n, char *cbuf)
{
std::map<tagint,int> *hash = frsptr->hash;
double **ctr = frsptr->ctr;
tagint *idclose = frsptr->idclose;
double *rsqclose = frsptr->rsqclose;
double *buf = (double *) cbuf;
int ndatums = n/5;
int j,imol;
tagint tag;
double delx,dely,delz,rsq;
double *x;
int m = 0;
for (int i = 0; i < ndatums; i++, m += 5) {
imol = static_cast<int> (buf[m]);
if (hash->find(imol) != hash->end()) {
j = hash->find(imol)->second;
tag = (tagint) ubuf(buf[m+1]).i;
x = &buf[m+2];
delx = x[0] - ctr[j][0];
dely = x[1] - ctr[j][1];
delz = x[2] - ctr[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= rsqclose[j]) {
if (rsq == rsqclose[j] && tag > idclose[j]) continue;
idclose[j] = tag;
rsqclose[j] = rsq;
}
}
}
}
/* ----------------------------------------------------------------------
process rigid body atoms from another proc
update rsqfar = distance from owning atom to other atom
------------------------------------------------------------------------- */
void FixRigidSmall::ring_farthest(int n, char *cbuf)
{
double **x = frsptr->atom->x;
imageint *image = frsptr->atom->image;
int nlocal = frsptr->atom->nlocal;
double *buf = (double *) cbuf;
int ndatums = n/4;
int iowner;
tagint tag;
double delx,dely,delz,rsq;
double *xx;
double unwrap[3];
int m = 0;
for (int i = 0; i < ndatums; i++, m += 4) {
tag = (tagint) ubuf(buf[m]).i;
iowner = frsptr->atom->map(tag);
if (iowner < 0 || iowner >= nlocal) continue;
frsptr->domain->unmap(x[iowner],image[iowner],unwrap);
xx = &buf[m+1];
delx = xx[0] - unwrap[0];
dely = xx[1] - unwrap[1];
delz = xx[2] - unwrap[2];
rsq = delx*delx + dely*dely + delz*delz;
frsptr->rsqfar = MAX(frsptr->rsqfar,rsq);
}
}
/* ----------------------------------------------------------------------
initialization of rigid body attributes
called at setup, so body/atom properties can be changed between runs
unless reading from infile, in which case only called once
sets extended flags, masstotal, center-of-mass
sets Cartesian and diagonalized inertia tensor
sets body image flags, but only on first call
------------------------------------------------------------------------- */
void FixRigidSmall::setup_bodies_static()
{
int i,ibody;
// extended = 1 if any particle in a rigid body is finite size
// or has a dipole moment
extended = orientflag = dorientflag = 0;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
AtomVecLine::Bonus *lbonus;
if (avec_line) lbonus = avec_line->bonus;
AtomVecTri::Bonus *tbonus;
if (avec_tri) tbonus = avec_tri->bonus;
double **mu = atom->mu;
double *radius = atom->radius;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *ellipsoid = atom->ellipsoid;
int *line = atom->line;
int *tri = atom->tri;
int *type = atom->type;
int nlocal = atom->nlocal;
if (atom->radius_flag || atom->ellipsoid_flag || atom->line_flag ||
atom->tri_flag || atom->mu_flag) {
int flag = 0;
for (i = 0; i < nlocal; i++) {
if (bodytag[i] == 0) continue;
if (radius && radius[i] > 0.0) flag = 1;
if (ellipsoid && ellipsoid[i] >= 0) flag = 1;
if (line && line[i] >= 0) flag = 1;
if (tri && tri[i] >= 0) flag = 1;
if (mu && mu[i][3] > 0.0) flag = 1;
}
MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world);
}
// extended = 1 if using molecule template with finite-size particles
// require all molecules in template to have consistent radiusflag
if (onemols) {
int radiusflag = onemols[0]->radiusflag;
for (i = 1; i < nmol; i++) {
if (onemols[i]->radiusflag != radiusflag)
error->all(FLERR,"Inconsistent use of finite-size particles "
"by molecule template molecules");
}
if (radiusflag) extended = 1;
}
// grow extended arrays and set extended flags for each particle
// orientflag = 4 if any particle stores ellipsoid or tri orientation
// orientflag = 1 if any particle stores line orientation
// dorientflag = 1 if any particle stores dipole orientation
if (extended) {
- grow_arrays(atom->nmax);
if (atom->ellipsoid_flag) orientflag = 4;
if (atom->line_flag) orientflag = 1;
if (atom->tri_flag) orientflag = 4;
if (atom->mu_flag) dorientflag = 1;
+ grow_arrays(atom->nmax);
for (i = 0; i < nlocal; i++) {
eflags[i] = 0;
if (bodytag[i] == 0) continue;
// set to POINT or SPHERE or ELLIPSOID or LINE
if (radius && radius[i] > 0.0) {
eflags[i] |= SPHERE;
eflags[i] |= OMEGA;
eflags[i] |= TORQUE;
} else if (ellipsoid && ellipsoid[i] >= 0) {
eflags[i] |= ELLIPSOID;
eflags[i] |= ANGMOM;
eflags[i] |= TORQUE;
} else if (line && line[i] >= 0) {
eflags[i] |= LINE;
eflags[i] |= OMEGA;
eflags[i] |= TORQUE;
} else if (tri && tri[i] >= 0) {
eflags[i] |= TRIANGLE;
eflags[i] |= ANGMOM;
eflags[i] |= TORQUE;
} else eflags[i] |= POINT;
// set DIPOLE if atom->mu and mu[3] > 0.0
if (atom->mu_flag && mu[i][3] > 0.0)
eflags[i] |= DIPOLE;
}
}
// acquire ghost bodies via forward comm
// set atom2body for ghost atoms via forward comm
// set atom2body for other owned atoms via reset_atom2body()
nghost_body = 0;
commflag = FULL_BODY;
comm->forward_comm_fix(this);
reset_atom2body();
// compute mass & center-of-mass of each rigid body
double **x = atom->x;
imageint *image = atom->image;
double *xcm;
for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
xcm = body[ibody].xcm;
xcm[0] = xcm[1] = xcm[2] = 0.0;
body[ibody].mass = 0.0;
}
double unwrap[3];
double massone;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
xcm = b->xcm;
xcm[0] += unwrap[0] * massone;
xcm[1] += unwrap[1] * massone;
xcm[2] += unwrap[2] * massone;
b->mass += massone;
}
// reverse communicate xcm, mass of all bodies
commflag = XCM_MASS;
comm->reverse_comm_fix(this,4);
for (ibody = 0; ibody < nlocal_body; ibody++) {
xcm = body[ibody].xcm;
xcm[0] /= body[ibody].mass;
xcm[1] /= body[ibody].mass;
xcm[2] /= body[ibody].mass;
}
// overwrite masstotal and center-of-mass with file values
// inbody[i] = 0/1 if Ith rigid body is initialized by file
int *inbody;
if (infile) {
memory->create(inbody,nlocal_body,"rigid/small:inbody");
for (ibody = 0; ibody < nlocal_body; ibody++) inbody[ibody] = 0;
readfile(0,NULL,inbody);
}
// set image flags for each rigid body to default values
// then remap the xcm of each body back into simulation box if needed
// staticflag check insures this in only done once, not on successive runs
if (!staticflag)
for (ibody = 0; ibody < nlocal_body; ibody++)
body[ibody].image = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
pre_neighbor();
// compute 6 moments of inertia of each body in Cartesian reference frame
// dx,dy,dz = coords relative to center-of-mass
// symmetric 3x3 inertia tensor stored in Voigt notation as 6-vector
memory->create(itensor,nlocal_body+nghost_body,6,"rigid/small:itensor");
for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++)
for (i = 0; i < 6; i++) itensor[ibody][i] = 0.0;
double dx,dy,dz;
double *inertia;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
domain->unmap(x[i],image[i],unwrap);
xcm = b->xcm;
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
inertia = itensor[atom2body[i]];
inertia[0] += massone * (dy*dy + dz*dz);
inertia[1] += massone * (dx*dx + dz*dz);
inertia[2] += massone * (dx*dx + dy*dy);
inertia[3] -= massone * dy*dz;
inertia[4] -= massone * dx*dz;
inertia[5] -= massone * dx*dy;
}
// extended particles may contribute extra terms to moments of inertia
if (extended) {
double ivec[6];
double *shape,*quatatom,*inertiaatom;
double length,theta;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
inertia = itensor[atom2body[i]];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
if (eflags[i] & SPHERE) {
inertia[0] += SINERTIA*massone * radius[i]*radius[i];
inertia[1] += SINERTIA*massone * radius[i]*radius[i];
inertia[2] += SINERTIA*massone * radius[i]*radius[i];
} else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
inertia[0] += ivec[0];
inertia[1] += ivec[1];
inertia[2] += ivec[2];
inertia[3] += ivec[3];
inertia[4] += ivec[4];
inertia[5] += ivec[5];
} else if (eflags[i] & LINE) {
length = lbonus[line[i]].length;
theta = lbonus[line[i]].theta;
MathExtra::inertia_line(length,theta,massone,ivec);
inertia[0] += ivec[0];
inertia[1] += ivec[1];
inertia[2] += ivec[2];
inertia[3] += ivec[3];
inertia[4] += ivec[4];
inertia[5] += ivec[5];
} else if (eflags[i] & TRIANGLE) {
inertiaatom = tbonus[tri[i]].inertia;
quatatom = tbonus[tri[i]].quat;
MathExtra::inertia_triangle(inertiaatom,quatatom,massone,ivec);
inertia[0] += ivec[0];
inertia[1] += ivec[1];
inertia[2] += ivec[2];
inertia[3] += ivec[3];
inertia[4] += ivec[4];
inertia[5] += ivec[5];
}
}
}
// reverse communicate inertia tensor of all bodies
commflag = ITENSOR;
comm->reverse_comm_fix(this,6);
// overwrite Cartesian inertia tensor with file values
if (infile) readfile(1,itensor,inbody);
// diagonalize inertia tensor for each body via Jacobi rotations
// inertia = 3 eigenvalues = principal moments of inertia
// evectors and exzy_space = 3 evectors = principal axes of rigid body
int ierror;
double cross[3];
double tensor[3][3],evectors[3][3];
double *ex,*ey,*ez;
for (ibody = 0; ibody < nlocal_body; ibody++) {
tensor[0][0] = itensor[ibody][0];
tensor[1][1] = itensor[ibody][1];
tensor[2][2] = itensor[ibody][2];
tensor[1][2] = tensor[2][1] = itensor[ibody][3];
tensor[0][2] = tensor[2][0] = itensor[ibody][4];
tensor[0][1] = tensor[1][0] = itensor[ibody][5];
inertia = body[ibody].inertia;
ierror = MathExtra::jacobi(tensor,inertia,evectors);
if (ierror) error->all(FLERR,
"Insufficient Jacobi rotations for rigid body");
ex = body[ibody].ex_space;
ex[0] = evectors[0][0];
ex[1] = evectors[1][0];
ex[2] = evectors[2][0];
ey = body[ibody].ey_space;
ey[0] = evectors[0][1];
ey[1] = evectors[1][1];
ey[2] = evectors[2][1];
ez = body[ibody].ez_space;
ez[0] = evectors[0][2];
ez[1] = evectors[1][2];
ez[2] = evectors[2][2];
// if any principal moment < scaled EPSILON, set to 0.0
double max;
max = MAX(inertia[0],inertia[1]);
max = MAX(max,inertia[2]);
if (inertia[0] < EPSILON*max) inertia[0] = 0.0;
if (inertia[1] < EPSILON*max) inertia[1] = 0.0;
if (inertia[2] < EPSILON*max) inertia[2] = 0.0;
// enforce 3 evectors as a right-handed coordinate system
// flip 3rd vector if needed
MathExtra::cross3(ex,ey,cross);
if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez);
// create initial quaternion
MathExtra::exyz_to_q(ex,ey,ez,body[ibody].quat);
}
// forward communicate updated info of all bodies
commflag = INITIAL;
comm->forward_comm_fix(this,26);
// displace = initial atom coords in basis of principal axes
// set displace = 0.0 for atoms not in any rigid body
// for extended particles, set their orientation wrt to rigid body
double qc[4],delta[3];
double *quatatom;
double theta_body;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) {
displace[i][0] = displace[i][1] = displace[i][2] = 0.0;
continue;
}
Body *b = &body[atom2body[i]];
domain->unmap(x[i],image[i],unwrap);
xcm = b->xcm;
delta[0] = unwrap[0] - xcm[0];
delta[1] = unwrap[1] - xcm[1];
delta[2] = unwrap[2] - xcm[2];
MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space,
delta,displace[i]);
if (extended) {
if (eflags[i] & ELLIPSOID) {
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::qconjugate(b->quat,qc);
MathExtra::quatquat(qc,quatatom,orient[i]);
MathExtra::qnormalize(orient[i]);
} else if (eflags[i] & LINE) {
if (b->quat[3] >= 0.0) theta_body = 2.0*acos(b->quat[0]);
else theta_body = -2.0*acos(b->quat[0]);
orient[i][0] = lbonus[line[i]].theta - theta_body;
while (orient[i][0] <= MINUSPI) orient[i][0] += TWOPI;
while (orient[i][0] > MY_PI) orient[i][0] -= TWOPI;
if (orientflag == 4) orient[i][1] = orient[i][2] = orient[i][3] = 0.0;
} else if (eflags[i] & TRIANGLE) {
quatatom = tbonus[tri[i]].quat;
MathExtra::qconjugate(b->quat,qc);
MathExtra::quatquat(qc,quatatom,orient[i]);
MathExtra::qnormalize(orient[i]);
} else if (orientflag == 4) {
orient[i][0] = orient[i][1] = orient[i][2] = orient[i][3] = 0.0;
} else if (orientflag == 1)
orient[i][0] = 0.0;
if (eflags[i] & DIPOLE) {
MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space,
mu[i],dorient[i]);
MathExtra::snormalize3(mu[i][3],dorient[i],dorient[i]);
} else if (dorientflag)
dorient[i][0] = dorient[i][1] = dorient[i][2] = 0.0;
}
}
// test for valid principal moments & axes
// recompute moments of inertia around new axes
// 3 diagonal moments should equal principal moments
// 3 off-diagonal moments should be 0.0
// extended particles may contribute extra terms to moments of inertia
for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++)
for (i = 0; i < 6; i++) itensor[ibody][i] = 0.0;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
inertia = itensor[atom2body[i]];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
inertia[0] += massone *
(displace[i][1]*displace[i][1] + displace[i][2]*displace[i][2]);
inertia[1] += massone *
(displace[i][0]*displace[i][0] + displace[i][2]*displace[i][2]);
inertia[2] += massone *
(displace[i][0]*displace[i][0] + displace[i][1]*displace[i][1]);
inertia[3] -= massone * displace[i][1]*displace[i][2];
inertia[4] -= massone * displace[i][0]*displace[i][2];
inertia[5] -= massone * displace[i][0]*displace[i][1];
}
if (extended) {
double ivec[6];
double *shape,*inertiaatom;
double length;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
inertia = itensor[atom2body[i]];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
if (eflags[i] & SPHERE) {
inertia[0] += SINERTIA*massone * radius[i]*radius[i];
inertia[1] += SINERTIA*massone * radius[i]*radius[i];
inertia[2] += SINERTIA*massone * radius[i]*radius[i];
} else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
MathExtra::inertia_ellipsoid(shape,orient[i],massone,ivec);
inertia[0] += ivec[0];
inertia[1] += ivec[1];
inertia[2] += ivec[2];
inertia[3] += ivec[3];
inertia[4] += ivec[4];
inertia[5] += ivec[5];
} else if (eflags[i] & LINE) {
length = lbonus[line[i]].length;
MathExtra::inertia_line(length,orient[i][0],massone,ivec);
inertia[0] += ivec[0];
inertia[1] += ivec[1];
inertia[2] += ivec[2];
inertia[3] += ivec[3];
inertia[4] += ivec[4];
inertia[5] += ivec[5];
} else if (eflags[i] & TRIANGLE) {
inertiaatom = tbonus[tri[i]].inertia;
MathExtra::inertia_triangle(inertiaatom,orient[i],massone,ivec);
inertia[0] += ivec[0];
inertia[1] += ivec[1];
inertia[2] += ivec[2];
inertia[3] += ivec[3];
inertia[4] += ivec[4];
inertia[5] += ivec[5];
}
}
}
// reverse communicate inertia tensor of all bodies
commflag = ITENSOR;
comm->reverse_comm_fix(this,6);
// error check that re-computed momemts of inertia match diagonalized ones
// do not do test for bodies with params read from infile
double norm;
for (ibody = 0; ibody < nlocal_body; ibody++) {
if (infile && inbody[ibody]) continue;
inertia = body[ibody].inertia;
if (inertia[0] == 0.0) {
if (fabs(itensor[ibody][0]) > TOLERANCE)
error->all(FLERR,"Fix rigid: Bad principal moments");
} else {
if (fabs((itensor[ibody][0]-inertia[0])/inertia[0]) >
TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
}
if (inertia[1] == 0.0) {
if (fabs(itensor[ibody][1]) > TOLERANCE)
error->all(FLERR,"Fix rigid: Bad principal moments");
} else {
if (fabs((itensor[ibody][1]-inertia[1])/inertia[1]) >
TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
}
if (inertia[2] == 0.0) {
if (fabs(itensor[ibody][2]) > TOLERANCE)
error->all(FLERR,"Fix rigid: Bad principal moments");
} else {
if (fabs((itensor[ibody][2]-inertia[2])/inertia[2]) >
TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
}
norm = (inertia[0] + inertia[1] + inertia[2]) / 3.0;
if (fabs(itensor[ibody][3]/norm) > TOLERANCE ||
fabs(itensor[ibody][4]/norm) > TOLERANCE ||
fabs(itensor[ibody][5]/norm) > TOLERANCE)
error->all(FLERR,"Fix rigid: Bad principal moments");
}
// clean up
memory->destroy(itensor);
if (infile) memory->destroy(inbody);
// static properties have now been initialized once
// used to prevent re-initialization which would re-read infile
staticflag = 1;
}
/* ----------------------------------------------------------------------
one-time initialization of dynamic rigid body attributes
Vcm and angmom, computed explicitly from constituent particles
even if wrong for overlapping particles, is OK,
since is just setting initial time=0 Vcm and angmom of the body
which can be estimated value
------------------------------------------------------------------------- */
void FixRigidSmall::setup_bodies_dynamic()
{
int i,ibody;
double massone,radone;
// sum vcm, angmom across all rigid bodies
// vcm = velocity of COM
// angmom = angular momentum around COM
double **x = atom->x;
double **v = atom->v;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
imageint *image = atom->image;
int nlocal = atom->nlocal;
double *xcm,*vcm,*acm;
double dx,dy,dz;
double unwrap[3];
for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
vcm = body[ibody].vcm;
vcm[0] = vcm[1] = vcm[2] = 0.0;
acm = body[ibody].angmom;
acm[0] = acm[1] = acm[2] = 0.0;
}
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
vcm = b->vcm;
vcm[0] += v[i][0] * massone;
vcm[1] += v[i][1] * massone;
vcm[2] += v[i][2] * massone;
domain->unmap(x[i],image[i],unwrap);
xcm = b->xcm;
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
acm = b->angmom;
acm[0] += dy * massone*v[i][2] - dz * massone*v[i][1];
acm[1] += dz * massone*v[i][0] - dx * massone*v[i][2];
acm[2] += dx * massone*v[i][1] - dy * massone*v[i][0];
}
// extended particles add their rotation to angmom of body
if (extended) {
AtomVecLine::Bonus *lbonus;
if (avec_line) lbonus = avec_line->bonus;
double **omega = atom->omega;
double **angmom = atom->angmom;
double *radius = atom->radius;
int *line = atom->line;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
if (eflags[i] & OMEGA) {
if (eflags[i] & SPHERE) {
radone = radius[i];
acm = b->angmom;
acm[0] += SINERTIA*rmass[i] * radone*radone * omega[i][0];
acm[1] += SINERTIA*rmass[i] * radone*radone * omega[i][1];
acm[2] += SINERTIA*rmass[i] * radone*radone * omega[i][2];
} else if (eflags[i] & LINE) {
radone = lbonus[line[i]].length;
b->angmom[2] += LINERTIA*rmass[i] * radone*radone * omega[i][2];
}
}
if (eflags[i] & ANGMOM) {
acm = b->angmom;
acm[0] += angmom[i][0];
acm[1] += angmom[i][1];
acm[2] += angmom[i][2];
}
}
}
// reverse communicate vcm, angmom of all bodies
commflag = VCM_ANGMOM;
comm->reverse_comm_fix(this,6);
// normalize velocity of COM
for (ibody = 0; ibody < nlocal_body; ibody++) {
vcm = body[ibody].vcm;
vcm[0] /= body[ibody].mass;
vcm[1] /= body[ibody].mass;
vcm[2] /= body[ibody].mass;
}
}
/* ----------------------------------------------------------------------
read per rigid body info from user-provided file
which = 0 to read total mass and center-of-mass
which = 1 to read 6 moments of inertia, store in array
flag inbody = 0 for local bodies whose info is read from file
nlines = # of lines of rigid body info, 0 is OK
one line = rigid-ID mass xcm ycm zcm ixx iyy izz ixy ixz iyz
and rigid-ID = mol-ID for fix rigid/small
------------------------------------------------------------------------- */
void FixRigidSmall::readfile(int which, double **array, int *inbody)
{
int i,j,m,nchunk,eofflag,nlines;
tagint id;
FILE *fp;
char *eof,*start,*next,*buf;
char line[MAXLINE];
// create local hash with key/value pairs
// key = mol ID of bodies my atoms own
// value = index into local body array
int nlocal = atom->nlocal;
hash = new std::map<tagint,int>();
for (i = 0; i < nlocal; i++)
if (bodyown[i] >= 0) (*hash)[atom->molecule[i]] = bodyown[i];
// open file and read header
if (me == 0) {
fp = fopen(infile,"r");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open fix rigid/small infile %s",infile);
error->one(FLERR,str);
}
while (1) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL)
error->one(FLERR,"Unexpected end of fix rigid/small file");
start = &line[strspn(line," \t\n\v\f\r")];
if (*start != '\0' && *start != '#') break;
}
sscanf(line,"%d",&nlines);
}
MPI_Bcast(&nlines,1,MPI_INT,0,world);
char *buffer = new char[CHUNK*MAXLINE];
char **values = new char*[ATTRIBUTE_PERBODY];
int nread = 0;
while (nread < nlines) {
nchunk = MIN(nlines-nread,CHUNK);
eofflag = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eofflag) error->all(FLERR,"Unexpected end of fix rigid/small file");
buf = buffer;
next = strchr(buf,'\n');
*next = '\0';
int nwords = atom->count_words(buf);
*next = '\n';
if (nwords != ATTRIBUTE_PERBODY)
error->all(FLERR,"Incorrect rigid body format in fix rigid/small file");
// loop over lines of rigid body attributes
// tokenize the line into values
// id = rigid body ID = mol-ID
// for which = 0, store mass/com in vec/array
// for which = 1, store interia tensor array, invert 3,4,5 values to Voigt
for (int i = 0; i < nchunk; i++) {
next = strchr(buf,'\n');
values[0] = strtok(buf," \t\n\r\f");
for (j = 1; j < nwords; j++)
values[j] = strtok(NULL," \t\n\r\f");
id = ATOTAGINT(values[0]);
if (id <= 0 || id > maxmol)
error->all(FLERR,"Invalid rigid body ID in fix rigid/small file");
if (hash->find(id) == hash->end()) {
buf = next + 1;
continue;
}
m = (*hash)[id];
inbody[m] = 1;
if (which == 0) {
body[m].mass = atof(values[1]);
body[m].xcm[0] = atof(values[2]);
body[m].xcm[1] = atof(values[3]);
body[m].xcm[2] = atof(values[4]);
} else {
array[m][0] = atof(values[5]);
array[m][1] = atof(values[6]);
array[m][2] = atof(values[7]);
array[m][3] = atof(values[10]);
array[m][4] = atof(values[9]);
array[m][5] = atof(values[8]);
}
buf = next + 1;
}
nread += nchunk;
}
if (me == 0) fclose(fp);
delete [] buffer;
delete [] values;
delete hash;
}
/* ----------------------------------------------------------------------
write out restart info for mass, COM, inertia tensor to file
identical format to infile option, so info can be read in when restarting
each proc contributes info for rigid bodies it owns
------------------------------------------------------------------------- */
void FixRigidSmall::write_restart_file(char *file)
{
FILE *fp;
// do not write file if bodies have not yet been intialized
if (!staticflag) return;
// proc 0 opens file and writes header
if (me == 0) {
char outfile[128];
sprintf(outfile,"%s.rigid",file);
fp = fopen(outfile,"w");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open fix rigid restart file %s",outfile);
error->one(FLERR,str);
}
fprintf(fp,"# fix rigid mass, COM, inertia tensor info for "
"%d bodies on timestep " BIGINT_FORMAT "\n\n",
nbody,update->ntimestep);
fprintf(fp,"%d\n",nbody);
}
// communication buffer for all my rigid body info
// max_size = largest buffer needed by any proc
int ncol = 11;
int sendrow = nlocal_body;
int maxrow;
MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);
double **buf;
if (me == 0) memory->create(buf,MAX(1,maxrow),ncol,"rigid/small:buf");
else memory->create(buf,MAX(1,sendrow),ncol,"rigid/small:buf");
// pack my rigid body info into buf
// compute I tensor against xyz axes from diagonalized I and current quat
// Ispace = P Idiag P_transpose
// P is stored column-wise in exyz_space
double p[3][3],pdiag[3][3],ispace[3][3];
for (int i = 0; i < nlocal_body; i++) {
MathExtra::col2mat(body[i].ex_space,body[i].ey_space,body[i].ez_space,p);
MathExtra::times3_diag(p,body[i].inertia,pdiag);
MathExtra::times3_transpose(pdiag,p,ispace);
buf[i][0] = atom->molecule[body[i].ilocal];
buf[i][1] = body[i].mass;
buf[i][2] = body[i].xcm[0];
buf[i][3] = body[i].xcm[1];
buf[i][4] = body[i].xcm[2];
buf[i][5] = ispace[0][0];
buf[i][6] = ispace[1][1];
buf[i][7] = ispace[2][2];
buf[i][8] = ispace[0][1];
buf[i][9] = ispace[0][2];
buf[i][10] = ispace[1][2];
}
// write one chunk of rigid body info per proc to file
// proc 0 pings each proc, receives its chunk, writes to file
// all other procs wait for ping, send their chunk to proc 0
int tmp,recvrow;
MPI_Status status;
MPI_Request request;
if (me == 0) {
for (int iproc = 0; iproc < nprocs; iproc++) {
if (iproc) {
MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_DOUBLE,iproc,0,world,&request);
MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_DOUBLE,&recvrow);
recvrow /= ncol;
} else recvrow = sendrow;
for (int i = 0; i < recvrow; i++)
fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e "
"%-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e\n",
static_cast<int> (buf[i][0]),buf[i][1],
buf[i][2],buf[i][3],buf[i][4],
buf[i][5],buf[i][6],buf[i][7],buf[i][8],buf[i][9],buf[i][10]);
}
} else {
MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status);
MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_DOUBLE,0,0,world);
}
// clean up and close file
memory->destroy(buf);
if (me == 0) fclose(fp);
}
/* ----------------------------------------------------------------------
allocate local atom-based arrays
------------------------------------------------------------------------- */
void FixRigidSmall::grow_arrays(int nmax)
{
memory->grow(bodyown,nmax,"rigid/small:bodyown");
memory->grow(bodytag,nmax,"rigid/small:bodytag");
memory->grow(atom2body,nmax,"rigid/small:atom2body");
memory->grow(displace,nmax,3,"rigid/small:displace");
if (extended) {
memory->grow(eflags,nmax,"rigid/small:eflags");
if (orientflag) memory->grow(orient,nmax,orientflag,"rigid/small:orient");
if (dorientflag) memory->grow(dorient,nmax,3,"rigid/small:dorient");
}
}
/* ----------------------------------------------------------------------
copy values within local atom-based arrays
------------------------------------------------------------------------- */
void FixRigidSmall::copy_arrays(int i, int j, int delflag)
{
bodytag[j] = bodytag[i];
displace[j][0] = displace[i][0];
displace[j][1] = displace[i][1];
displace[j][2] = displace[i][2];
if (extended) {
eflags[j] = eflags[i];
for (int k = 0; k < orientflag; k++)
orient[j][k] = orient[i][k];
if (dorientflag) {
dorient[j][0] = dorient[i][0];
dorient[j][1] = dorient[i][1];
dorient[j][2] = dorient[i][2];
}
}
// if deleting atom J via delflag and J owns a body, then delete it
if (delflag && bodyown[j] >= 0) {
bodyown[body[nlocal_body-1].ilocal] = bodyown[j];
memcpy(&body[bodyown[j]],&body[nlocal_body-1],sizeof(Body));
nlocal_body--;
}
// if atom I owns a body, reset I's body.ilocal to loc J
// do NOT do this if self-copy (I=J) since I's body is already deleted
if (bodyown[i] >= 0 && i != j) body[bodyown[i]].ilocal = j;
bodyown[j] = bodyown[i];
}
/* ----------------------------------------------------------------------
initialize one atom's array values, called when atom is created
------------------------------------------------------------------------- */
void FixRigidSmall::set_arrays(int i)
{
bodyown[i] = -1;
bodytag[i] = 0;
atom2body[i] = -1;
displace[i][0] = 0.0;
displace[i][1] = 0.0;
displace[i][2] = 0.0;
}
/* ----------------------------------------------------------------------
initialize a molecule inserted by another fix, e.g. deposit or pour
called when molecule is created
nlocalprev = # of atoms on this proc before molecule inserted
tagprev = atom ID previous to new atoms in the molecule
xgeom = geometric center of new molecule
vcm = COM velocity of new molecule
quat = rotation of new molecule (around geometric center)
relative to template in Molecule class
------------------------------------------------------------------------- */
void FixRigidSmall::set_molecule(int nlocalprev, tagint tagprev, int imol,
double *xgeom, double *vcm, double *quat)
{
int m;
double ctr2com[3],ctr2com_rotate[3];
double rotmat[3][3];
int nlocal = atom->nlocal;
if (nlocalprev == nlocal) return;
tagint *tag = atom->tag;
for (int i = nlocalprev; i < nlocal; i++) {
bodytag[i] = tagprev + onemols[imol]->comatom;
if (tag[i]-tagprev == onemols[imol]->comatom) bodyown[i] = nlocal_body;
m = tag[i] - tagprev-1;
displace[i][0] = onemols[imol]->dxbody[m][0];
displace[i][1] = onemols[imol]->dxbody[m][1];
displace[i][2] = onemols[imol]->dxbody[m][2];
eflags[i] = 0;
if (onemols[imol]->radiusflag) {
eflags[i] |= SPHERE;
eflags[i] |= OMEGA;
eflags[i] |= TORQUE;
}
if (bodyown[i] >= 0) {
if (nlocal_body == nmax_body) grow_body();
Body *b = &body[nlocal_body];
b->mass = onemols[imol]->masstotal;
// new COM = Q (onemols[imol]->xcm - onemols[imol]->center) + xgeom
// Q = rotation matrix associated with quat
MathExtra::quat_to_mat(quat,rotmat);
MathExtra::sub3(onemols[imol]->com,onemols[imol]->center,ctr2com);
MathExtra::matvec(rotmat,ctr2com,ctr2com_rotate);
MathExtra::add3(ctr2com_rotate,xgeom,b->xcm);
b->vcm[0] = vcm[0];
b->vcm[1] = vcm[1];
b->vcm[2] = vcm[2];
b->inertia[0] = onemols[imol]->inertia[0];
b->inertia[1] = onemols[imol]->inertia[1];
b->inertia[2] = onemols[imol]->inertia[2];
// final quat is product of insertion quat and original quat
// true even if insertion rotation was not around COM
MathExtra::quatquat(quat,onemols[imol]->quat,b->quat);
MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space,b->ez_space);
b->angmom[0] = b->angmom[1] = b->angmom[2] = 0.0;
b->omega[0] = b->omega[1] = b->omega[2] = 0.0;
b->image = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
b->ilocal = i;
nlocal_body++;
}
}
// increment total # of rigid bodies
nbody++;
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for exchange with another proc
------------------------------------------------------------------------- */
int FixRigidSmall::pack_exchange(int i, double *buf)
{
buf[0] = ubuf(bodytag[i]).d;
buf[1] = displace[i][0];
buf[2] = displace[i][1];
buf[3] = displace[i][2];
// extended attribute info
int m = 4;
if (extended) {
buf[m++] = eflags[i];
for (int j = 0; j < orientflag; j++)
buf[m++] = orient[i][j];
if (dorientflag) {
buf[m++] = dorient[i][0];
buf[m++] = dorient[i][1];
buf[m++] = dorient[i][2];
}
}
// atom not in a rigid body
if (!bodytag[i]) return m;
// atom does not own its rigid body
if (bodyown[i] < 0) {
buf[m++] = 0;
return m;
}
// body info for atom that owns a rigid body
buf[m++] = 1;
memcpy(&buf[m],&body[bodyown[i]],sizeof(Body));
m += bodysize;
return m;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based arrays from exchange with another proc
------------------------------------------------------------------------- */
int FixRigidSmall::unpack_exchange(int nlocal, double *buf)
{
bodytag[nlocal] = (tagint) ubuf(buf[0]).i;
displace[nlocal][0] = buf[1];
displace[nlocal][1] = buf[2];
displace[nlocal][2] = buf[3];
// extended attribute info
int m = 4;
if (extended) {
eflags[nlocal] = static_cast<int> (buf[m++]);
for (int j = 0; j < orientflag; j++)
orient[nlocal][j] = buf[m++];
if (dorientflag) {
dorient[nlocal][0] = buf[m++];
dorient[nlocal][1] = buf[m++];
dorient[nlocal][2] = buf[m++];
}
}
// atom not in a rigid body
if (!bodytag[nlocal]) {
bodyown[nlocal] = -1;
return m;
}
// atom does not own its rigid body
bodyown[nlocal] = static_cast<int> (buf[m++]);
if (bodyown[nlocal] == 0) {
bodyown[nlocal] = -1;
return m;
}
// body info for atom that owns a rigid body
if (nlocal_body == nmax_body) grow_body();
memcpy(&body[nlocal_body],&buf[m],sizeof(Body));
m += bodysize;
body[nlocal_body].ilocal = nlocal;
bodyown[nlocal] = nlocal_body++;
return m;
}
/* ----------------------------------------------------------------------
only pack body info if own or ghost atom owns the body
for FULL_BODY, send 0/1 flag with every atom
------------------------------------------------------------------------- */
int FixRigidSmall::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j;
double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space,*conjqm;
int m = 0;
if (commflag == INITIAL) {
for (i = 0; i < n; i++) {
j = list[i];
if (bodyown[j] < 0) continue;
xcm = body[bodyown[j]].xcm;
buf[m++] = xcm[0];
buf[m++] = xcm[1];
buf[m++] = xcm[2];
vcm = body[bodyown[j]].vcm;
buf[m++] = vcm[0];
buf[m++] = vcm[1];
buf[m++] = vcm[2];
quat = body[bodyown[j]].quat;
buf[m++] = quat[0];
buf[m++] = quat[1];
buf[m++] = quat[2];
buf[m++] = quat[3];
omega = body[bodyown[j]].omega;
buf[m++] = omega[0];
buf[m++] = omega[1];
buf[m++] = omega[2];
ex_space = body[bodyown[j]].ex_space;
buf[m++] = ex_space[0];
buf[m++] = ex_space[1];
buf[m++] = ex_space[2];
ey_space = body[bodyown[j]].ey_space;
buf[m++] = ey_space[0];
buf[m++] = ey_space[1];
buf[m++] = ey_space[2];
ez_space = body[bodyown[j]].ez_space;
buf[m++] = ez_space[0];
buf[m++] = ez_space[1];
buf[m++] = ez_space[2];
conjqm = body[bodyown[j]].conjqm;
buf[m++] = conjqm[0];
buf[m++] = conjqm[1];
buf[m++] = conjqm[2];
buf[m++] = conjqm[3];
}
} else if (commflag == FINAL) {
for (i = 0; i < n; i++) {
j = list[i];
if (bodyown[j] < 0) continue;
vcm = body[bodyown[j]].vcm;
buf[m++] = vcm[0];
buf[m++] = vcm[1];
buf[m++] = vcm[2];
omega = body[bodyown[j]].omega;
buf[m++] = omega[0];
buf[m++] = omega[1];
buf[m++] = omega[2];
conjqm = body[bodyown[j]].conjqm;
buf[m++] = conjqm[0];
buf[m++] = conjqm[1];
buf[m++] = conjqm[2];
buf[m++] = conjqm[3];
}
} else if (commflag == FULL_BODY) {
for (i = 0; i < n; i++) {
j = list[i];
if (bodyown[j] < 0) buf[m++] = 0;
else {
buf[m++] = 1;
memcpy(&buf[m],&body[bodyown[j]],sizeof(Body));
m += bodysize;
}
}
}
return m;
}
/* ----------------------------------------------------------------------
only ghost atoms are looped over
for FULL_BODY, store a new ghost body if this atom owns it
for other commflag values, only unpack body info if atom owns it
------------------------------------------------------------------------- */
void FixRigidSmall::unpack_forward_comm(int n, int first, double *buf)
{
int i,j,last;
double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space,*conjqm;
int m = 0;
last = first + n;
if (commflag == INITIAL) {
for (i = first; i < last; i++) {
if (bodyown[i] < 0) continue;
xcm = body[bodyown[i]].xcm;
xcm[0] = buf[m++];
xcm[1] = buf[m++];
xcm[2] = buf[m++];
vcm = body[bodyown[i]].vcm;
vcm[0] = buf[m++];
vcm[1] = buf[m++];
vcm[2] = buf[m++];
quat = body[bodyown[i]].quat;
quat[0] = buf[m++];
quat[1] = buf[m++];
quat[2] = buf[m++];
quat[3] = buf[m++];
omega = body[bodyown[i]].omega;
omega[0] = buf[m++];
omega[1] = buf[m++];
omega[2] = buf[m++];
ex_space = body[bodyown[i]].ex_space;
ex_space[0] = buf[m++];
ex_space[1] = buf[m++];
ex_space[2] = buf[m++];
ey_space = body[bodyown[i]].ey_space;
ey_space[0] = buf[m++];
ey_space[1] = buf[m++];
ey_space[2] = buf[m++];
ez_space = body[bodyown[i]].ez_space;
ez_space[0] = buf[m++];
ez_space[1] = buf[m++];
ez_space[2] = buf[m++];
conjqm = body[bodyown[i]].conjqm;
conjqm[0] = buf[m++];
conjqm[1] = buf[m++];
conjqm[2] = buf[m++];
conjqm[3] = buf[m++];
}
} else if (commflag == FINAL) {
for (i = first; i < last; i++) {
if (bodyown[i] < 0) continue;
vcm = body[bodyown[i]].vcm;
vcm[0] = buf[m++];
vcm[1] = buf[m++];
vcm[2] = buf[m++];
omega = body[bodyown[i]].omega;
omega[0] = buf[m++];
omega[1] = buf[m++];
omega[2] = buf[m++];
conjqm = body[bodyown[i]].conjqm;
conjqm[0] = buf[m++];
conjqm[1] = buf[m++];
conjqm[2] = buf[m++];
conjqm[3] = buf[m++];
}
} else if (commflag == FULL_BODY) {
for (i = first; i < last; i++) {
bodyown[i] = static_cast<int> (buf[m++]);
if (bodyown[i] == 0) bodyown[i] = -1;
else {
j = nlocal_body + nghost_body;
if (j == nmax_body) grow_body();
memcpy(&body[j],&buf[m],sizeof(Body));
m += bodysize;
body[j].ilocal = i;
bodyown[i] = j;
nghost_body++;
}
}
}
}
/* ----------------------------------------------------------------------
only ghost atoms are looped over
only pack body info if atom owns it
------------------------------------------------------------------------- */
int FixRigidSmall::pack_reverse_comm(int n, int first, double *buf)
{
int i,j,m,last;
double *fcm,*torque,*vcm,*angmom,*xcm;
m = 0;
last = first + n;
if (commflag == FORCE_TORQUE) {
for (i = first; i < last; i++) {
if (bodyown[i] < 0) continue;
fcm = body[bodyown[i]].fcm;
buf[m++] = fcm[0];
buf[m++] = fcm[1];
buf[m++] = fcm[2];
torque = body[bodyown[i]].torque;
buf[m++] = torque[0];
buf[m++] = torque[1];
buf[m++] = torque[2];
}
} else if (commflag == VCM_ANGMOM) {
for (i = first; i < last; i++) {
if (bodyown[i] < 0) continue;
vcm = body[bodyown[i]].vcm;
buf[m++] = vcm[0];
buf[m++] = vcm[1];
buf[m++] = vcm[2];
angmom = body[bodyown[i]].angmom;
buf[m++] = angmom[0];
buf[m++] = angmom[1];
buf[m++] = angmom[2];
}
} else if (commflag == XCM_MASS) {
for (i = first; i < last; i++) {
if (bodyown[i] < 0) continue;
xcm = body[bodyown[i]].xcm;
buf[m++] = xcm[0];
buf[m++] = xcm[1];
buf[m++] = xcm[2];
buf[m++] = body[bodyown[i]].mass;
}
} else if (commflag == ITENSOR) {
for (i = first; i < last; i++) {
if (bodyown[i] < 0) continue;
j = bodyown[i];
buf[m++] = itensor[j][0];
buf[m++] = itensor[j][1];
buf[m++] = itensor[j][2];
buf[m++] = itensor[j][3];
buf[m++] = itensor[j][4];
buf[m++] = itensor[j][5];
}
} else if (commflag == DOF) {
for (i = first; i < last; i++) {
if (bodyown[i] < 0) continue;
j = bodyown[i];
buf[m++] = counts[j][0];
buf[m++] = counts[j][1];
buf[m++] = counts[j][2];
}
}
return m;
}
/* ----------------------------------------------------------------------
only unpack body info if own or ghost atom owns the body
------------------------------------------------------------------------- */
void FixRigidSmall::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,k;
double *fcm,*torque,*vcm,*angmom,*xcm;
int m = 0;
if (commflag == FORCE_TORQUE) {
for (i = 0; i < n; i++) {
j = list[i];
if (bodyown[j] < 0) continue;
fcm = body[bodyown[j]].fcm;
fcm[0] += buf[m++];
fcm[1] += buf[m++];
fcm[2] += buf[m++];
torque = body[bodyown[j]].torque;
torque[0] += buf[m++];
torque[1] += buf[m++];
torque[2] += buf[m++];
}
} else if (commflag == VCM_ANGMOM) {
for (i = 0; i < n; i++) {
j = list[i];
if (bodyown[j] < 0) continue;
vcm = body[bodyown[j]].vcm;
vcm[0] += buf[m++];
vcm[1] += buf[m++];
vcm[2] += buf[m++];
angmom = body[bodyown[j]].angmom;
angmom[0] += buf[m++];
angmom[1] += buf[m++];
angmom[2] += buf[m++];
}
} else if (commflag == XCM_MASS) {
for (i = 0; i < n; i++) {
j = list[i];
if (bodyown[j] < 0) continue;
xcm = body[bodyown[j]].xcm;
xcm[0] += buf[m++];
xcm[1] += buf[m++];
xcm[2] += buf[m++];
body[bodyown[j]].mass += buf[m++];
}
} else if (commflag == ITENSOR) {
for (i = 0; i < n; i++) {
j = list[i];
if (bodyown[j] < 0) continue;
k = bodyown[j];
itensor[k][0] += buf[m++];
itensor[k][1] += buf[m++];
itensor[k][2] += buf[m++];
itensor[k][3] += buf[m++];
itensor[k][4] += buf[m++];
itensor[k][5] += buf[m++];
}
} else if (commflag == DOF) {
for (i = 0; i < n; i++) {
j = list[i];
if (bodyown[j] < 0) continue;
k = bodyown[j];
counts[k][0] += static_cast<int> (buf[m++]);
counts[k][1] += static_cast<int> (buf[m++]);
counts[k][2] += static_cast<int> (buf[m++]);
}
}
}
/* ----------------------------------------------------------------------
grow body data structure
------------------------------------------------------------------------- */
void FixRigidSmall::grow_body()
{
nmax_body += DELTA_BODY;
body = (Body *) memory->srealloc(body,nmax_body*sizeof(Body),
"rigid/small:body");
}
/* ----------------------------------------------------------------------
reset atom2body for all owned atoms
do this via bodyown of atom that owns the body the owned atom is in
atom2body values can point to original body or any image of the body
------------------------------------------------------------------------- */
void FixRigidSmall::reset_atom2body()
{
int iowner;
// iowner = index of atom that owns the body that atom I is in
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
atom2body[i] = -1;
if (bodytag[i]) {
iowner = atom->map(bodytag[i]);
if (iowner == -1) {
char str[128];
sprintf(str,
"Rigid body atoms " TAGINT_FORMAT " " TAGINT_FORMAT
" missing on proc %d at step " BIGINT_FORMAT,
atom->tag[i],bodytag[i],comm->me,update->ntimestep);
error->one(FLERR,str);
}
atom2body[i] = bodyown[iowner];
}
}
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dtq = 0.5 * update->dt;
}
/* ----------------------------------------------------------------------
zero linear momentum of each rigid body
set Vcm to 0.0, then reset velocities of particles via set_v()
------------------------------------------------------------------------- */
void FixRigidSmall::zero_momentum()
{
double *vcm;
for (int ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
vcm = body[ibody].vcm;
vcm[0] = vcm[1] = vcm[2] = 0.0;
}
// forward communicate of vcm to all ghost copies
commflag = FINAL;
comm->forward_comm_fix(this,10);
// set velocity of atoms in rigid bodues
evflag = 0;
set_v();
}
/* ----------------------------------------------------------------------
zero angular momentum of each rigid body
set angmom/omega to 0.0, then reset velocities of particles via set_v()
------------------------------------------------------------------------- */
void FixRigidSmall::zero_rotation()
{
double *angmom,*omega;
for (int ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
angmom = body[ibody].angmom;
angmom[0] = angmom[1] = angmom[2] = 0.0;
omega = body[ibody].omega;
omega[0] = omega[1] = omega[2] = 0.0;
}
// forward communicate of omega to all ghost copies
commflag = FINAL;
comm->forward_comm_fix(this,10);
// set velocity of atoms in rigid bodues
evflag = 0;
set_v();
}
/* ---------------------------------------------------------------------- */
void *FixRigidSmall::extract(const char *str, int &dim)
{
if (strcmp(str,"body") == 0) {
dim = 1;
return atom2body;
}
if (strcmp(str,"onemol") == 0) {
dim = 0;
return onemols;
}
// return vector of rigid body masses, for owned+ghost bodies
// used by granular pair styles, indexed by atom2body
if (strcmp(str,"masstotal") == 0) {
dim = 1;
if (nmax_mass < nmax_body) {
memory->destroy(mass_body);
nmax_mass = nmax_body;
memory->create(mass_body,nmax_mass,"rigid:mass_body");
}
int n = nlocal_body + nghost_body;
for (int i = 0; i < n; i++)
mass_body[i] = body[i].mass;
return mass_body;
}
return NULL;
}
/* ----------------------------------------------------------------------
return translational KE for all rigid bodies
KE = 1/2 M Vcm^2
sum local body results across procs
------------------------------------------------------------------------- */
double FixRigidSmall::extract_ke()
{
double *vcm;
double ke = 0.0;
for (int i = 0; i < nlocal_body; i++) {
vcm = body[i].vcm;
ke += body[i].mass * (vcm[0]*vcm[0] + vcm[1]*vcm[1] + vcm[2]*vcm[2]);
}
double keall;
MPI_Allreduce(&ke,&keall,1,MPI_DOUBLE,MPI_SUM,world);
return 0.5*keall;
}
/* ----------------------------------------------------------------------
return rotational KE for all rigid bodies
Erotational = 1/2 I wbody^2
------------------------------------------------------------------------- */
double FixRigidSmall::extract_erotational()
{
double wbody[3],rot[3][3];
double *inertia;
double erotate = 0.0;
for (int i = 0; i < nlocal_body; i++) {
// for Iw^2 rotational term, need wbody = angular velocity in body frame
// not omega = angular velocity in space frame
inertia = body[i].inertia;
MathExtra::quat_to_mat(body[i].quat,rot);
MathExtra::transpose_matvec(rot,body[i].angmom,wbody);
if (inertia[0] == 0.0) wbody[0] = 0.0;
else wbody[0] /= inertia[0];
if (inertia[1] == 0.0) wbody[1] = 0.0;
else wbody[1] /= inertia[1];
if (inertia[2] == 0.0) wbody[2] = 0.0;
else wbody[2] /= inertia[2];
erotate += inertia[0]*wbody[0]*wbody[0] + inertia[1]*wbody[1]*wbody[1] +
inertia[2]*wbody[2]*wbody[2];
}
double erotateall;
MPI_Allreduce(&erotate,&erotateall,1,MPI_DOUBLE,MPI_SUM,world);
return 0.5*erotateall;
}
/* ----------------------------------------------------------------------
return temperature of collection of rigid bodies
non-active DOF are removed by fflag/tflag and in tfactor
------------------------------------------------------------------------- */
double FixRigidSmall::compute_scalar()
{
double wbody[3],rot[3][3];
double *vcm,*inertia;
double t = 0.0;
for (int i = 0; i < nlocal_body; i++) {
vcm = body[i].vcm;
t += body[i].mass * (vcm[0]*vcm[0] + vcm[1]*vcm[1] + vcm[2]*vcm[2]);
// for Iw^2 rotational term, need wbody = angular velocity in body frame
// not omega = angular velocity in space frame
inertia = body[i].inertia;
MathExtra::quat_to_mat(body[i].quat,rot);
MathExtra::transpose_matvec(rot,body[i].angmom,wbody);
if (inertia[0] == 0.0) wbody[0] = 0.0;
else wbody[0] /= inertia[0];
if (inertia[1] == 0.0) wbody[1] = 0.0;
else wbody[1] /= inertia[1];
if (inertia[2] == 0.0) wbody[2] = 0.0;
else wbody[2] /= inertia[2];
t += inertia[0]*wbody[0]*wbody[0] + inertia[1]*wbody[1]*wbody[1] +
inertia[2]*wbody[2]*wbody[2];
}
double tall;
MPI_Allreduce(&t,&tall,1,MPI_DOUBLE,MPI_SUM,world);
double tfactor = force->mvv2e / (6.0*nbody * force->boltz);
tall *= tfactor;
return tall;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixRigidSmall::memory_usage()
{
int nmax = atom->nmax;
double bytes = 2 * nmax * sizeof(int);
bytes += nmax*3 * sizeof(double);
bytes += maxvatom*6 * sizeof(double); // vatom
if (extended) {
bytes += nmax * sizeof(int);
if (orientflag) bytes = nmax*orientflag * sizeof(double);
if (dorientflag) bytes = nmax*3 * sizeof(double);
}
bytes += nmax_body * sizeof(Body);
return bytes;
}
/* ----------------------------------------------------------------------
debug method for sanity checking of atom/body data pointers
------------------------------------------------------------------------- */
/*
void FixRigidSmall::check(int flag)
{
for (int i = 0; i < atom->nlocal; i++) {
if (bodyown[i] >= 0) {
if (bodytag[i] != atom->tag[i]) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD AAA");
}
if (bodyown[i] < 0 || bodyown[i] >= nlocal_body) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD BBB");
}
if (atom2body[i] != bodyown[i]) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD CCC");
}
if (body[bodyown[i]].ilocal != i) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD DDD");
}
}
}
for (int i = 0; i < atom->nlocal; i++) {
if (bodyown[i] < 0 && bodytag[i] > 0) {
if (atom2body[i] < 0 || atom2body[i] >= nlocal_body+nghost_body) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD EEE");
}
if (bodytag[i] != atom->tag[body[atom2body[i]].ilocal]) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD FFF");
}
}
}
for (int i = atom->nlocal; i < atom->nlocal + atom->nghost; i++) {
if (bodyown[i] >= 0) {
if (bodyown[i] < nlocal_body ||
bodyown[i] >= nlocal_body+nghost_body) {
printf("Values %d %d: %d %d %d\n",
i,atom->tag[i],bodyown[i],nlocal_body,nghost_body);
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD GGG");
}
if (body[bodyown[i]].ilocal != i) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD HHH");
}
}
}
for (int i = 0; i < nlocal_body; i++) {
if (body[i].ilocal < 0 || body[i].ilocal >= atom->nlocal) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD III");
}
if (bodytag[body[i].ilocal] != atom->tag[body[i].ilocal] ||
bodyown[body[i].ilocal] != i) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD JJJ");
}
}
for (int i = nlocal_body; i < nlocal_body + nghost_body; i++) {
if (body[i].ilocal < atom->nlocal ||
body[i].ilocal >= atom->nlocal + atom->nghost) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD KKK");
}
if (bodyown[body[i].ilocal] != i) {
printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
errorx->one(FLERR,"BAD LLL");
}
}
}
*/
diff --git a/src/SNAP/pair_snap.cpp b/src/SNAP/pair_snap.cpp
index 03bee46a8..e008b0c87 100644
--- a/src/SNAP/pair_snap.cpp
+++ b/src/SNAP/pair_snap.cpp
@@ -1,1747 +1,1742 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "pair_snap.h"
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "sna.h"
#include "memory.h"
#include "error.h"
#include "openmp_snap.h"
#include "domain.h"
#include <cmath>
using namespace LAMMPS_NS;
#define MAXLINE 1024
#define MAXWORD 3
/* ---------------------------------------------------------------------- */
PairSNAP::PairSNAP(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
nelements = 0;
elements = NULL;
radelem = NULL;
wjelem = NULL;
coeffelem = NULL;
nmax = 0;
nthreads = 1;
schedule_user = 0;
schedule_time_guided = -1;
schedule_time_dynamic = -1;
ncalls_neigh =-1;
ilistmask_max = 0;
ilistmask = NULL;
ghostinum = 0;
ghostilist_max = 0;
ghostilist = NULL;
ghostnumneigh_max = 0;
ghostnumneigh = NULL;
ghostneighs = NULL;
ghostfirstneigh = NULL;
ghostneighs_total = 0;
ghostneighs_max = 0;
i_max = 0;
i_neighmax = 0;
i_numpairs = 0;
i_rij = NULL;
i_inside = NULL;
i_wj = NULL;
i_rcutij = NULL;
i_ninside = NULL;
i_pairs = NULL;
i_uarraytot_r = NULL;
i_uarraytot_i = NULL;
i_zarray_r = NULL;
i_zarray_i =NULL;
use_shared_arrays = 0;
#ifdef TIMING_INFO
timers[0] = 0;
timers[1] = 0;
timers[2] = 0;
timers[3] = 0;
#endif
// Need to set this because restart not handled by PairHybrid
sna = NULL;
}
/* ---------------------------------------------------------------------- */
PairSNAP::~PairSNAP()
{
if (nelements) {
for (int i = 0; i < nelements; i++)
delete[] elements[i];
delete[] elements;
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(coeffelem);
}
// Need to set this because restart not handled by PairHybrid
if (sna) {
#ifdef TIMING_INFO
double time[5];
double timeave[5];
double timeave_mpi[5];
double timemax_mpi[5];
for (int i = 0; i < 5; i++) {
time[i] = 0;
timeave[i] = 0;
for (int tid = 0; tid<nthreads; tid++) {
if (sna[tid]->timers[i]>time[i])
time[i] = sna[tid]->timers[i];
timeave[i] += sna[tid]->timers[i];
}
timeave[i] /= nthreads;
}
MPI_Reduce(timeave, timeave_mpi, 5, MPI_DOUBLE, MPI_SUM, 0, world);
MPI_Reduce(time, timemax_mpi, 5, MPI_DOUBLE, MPI_MAX, 0, world);
#endif
for (int tid = 0; tid<nthreads; tid++)
delete sna[tid];
delete [] sna;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(map);
}
}
void PairSNAP::compute(int eflag, int vflag)
{
if (use_optimized)
compute_optimized(eflag, vflag);
else
compute_regular(eflag, vflag);
}
/* ----------------------------------------------------------------------
This version is a straightforward implementation
---------------------------------------------------------------------- */
void PairSNAP::compute_regular(int eflag, int vflag)
{
- int i,j,inum,jnum,ninside;
+ int i,j,jnum,ninside;
double delx,dely,delz,evdwl,rsq;
double fij[3];
- int *ilist,*jlist,*numneigh,**firstneigh;
+ int *jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
class SNA* snaptr = sna[0];
- inum = list->inum;
- ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
for (int ii = 0; ii < list->inum; ii++) {
i = list->ilist[ii];
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
const int itype = type[i];
const int ielem = map[itype];
const double radi = radelem[ielem];
jlist = firstneigh[i];
jnum = numneigh[i];
// insure rij, inside, wj, and rcutij are of size jnum
snaptr->grow_rij(jnum);
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
// wj = weights for neighbors of I within cutoff
// rcutij = cutoffs for neighbors of I within cutoff
// note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = x[j][0] - xtmp;
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j];
int jelem = map[jtype];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
snaptr->rij[ninside][0] = delx;
snaptr->rij[ninside][1] = dely;
snaptr->rij[ninside][2] = delz;
snaptr->inside[ninside] = j;
snaptr->wj[ninside] = wjelem[jelem];
snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac;
ninside++;
}
}
// compute Ui, Zi, and Bi for atom I
snaptr->compute_ui(ninside);
snaptr->compute_zi();
if (!gammaoneflag) {
snaptr->compute_bi();
snaptr->copy_bi2bvec();
}
// for neighbors of I within cutoff:
// compute dUi/drj and dBi/drj
// Fij = dEi/dRj = -dEi/dRi => add to Fi, subtract from Fj
double* coeffi = coeffelem[ielem];
for (int jj = 0; jj < ninside; jj++) {
int j = snaptr->inside[jj];
snaptr->compute_duidrj(snaptr->rij[jj],
snaptr->wj[jj],snaptr->rcutij[jj]);
snaptr->compute_dbidrj();
snaptr->copy_dbi2dbvec();
fij[0] = 0.0;
fij[1] = 0.0;
fij[2] = 0.0;
for (int k = 1; k <= ncoeff; k++) {
double bgb;
if (gammaoneflag)
bgb = coeffi[k];
else bgb = coeffi[k]*
gamma*pow(snaptr->bvec[k-1],gamma-1.0);
fij[0] += bgb*snaptr->dbvec[k-1][0];
fij[1] += bgb*snaptr->dbvec[k-1][1];
fij[2] += bgb*snaptr->dbvec[k-1][2];
}
f[i][0] += fij[0];
f[i][1] += fij[1];
f[i][2] += fij[2];
f[j][0] -= fij[0];
f[j][1] -= fij[1];
f[j][2] -= fij[2];
if (evflag)
ev_tally_xyz(i,j,nlocal,newton_pair,0.0,0.0,
fij[0],fij[1],fij[2],
snaptr->rij[jj][0],snaptr->rij[jj][1],
snaptr->rij[jj][2]);
}
if (eflag) {
// evdwl = energy of atom I, sum over coeffs_k * Bi_k
evdwl = coeffi[0];
if (gammaoneflag) {
snaptr->compute_bi();
snaptr->copy_bi2bvec();
for (int k = 1; k <= ncoeff; k++)
evdwl += coeffi[k]*snaptr->bvec[k-1];
} else
for (int k = 1; k <= ncoeff; k++)
evdwl += coeffi[k]*pow(snaptr->bvec[k-1],gamma);
ev_tally_full(i,2.0*evdwl,0.0,0.0,delx,dely,delz);
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
This version is optimized for threading, micro-load balancing
---------------------------------------------------------------------- */
void PairSNAP::compute_optimized(int eflag, int vflag)
{
// if reneighboring took place do load_balance if requested
if (do_load_balance > 0 &&
(neighbor->ncalls != ncalls_neigh)) {
ghostinum = 0;
// reset local ghost neighbor lists
ncalls_neigh = neighbor->ncalls;
if (ilistmask_max < list->inum) {
memory->grow(ilistmask,list->inum,"PairSnap::ilistmask");
ilistmask_max = list->inum;
}
for (int i = 0; i < list->inum; i++)
ilistmask[i] = 1;
//multiple passes for loadbalancing
for (int i = 0; i < do_load_balance; i++)
load_balance();
}
int numpairs = 0;
for (int ii = 0; ii < list->inum; ii++) {
if ((do_load_balance <= 0) || ilistmask[ii]) {
int i = list->ilist[ii];
int jnum = list->numneigh[i];
numpairs += jnum;
}
}
if (do_load_balance)
for (int ii = 0; ii < ghostinum; ii++) {
int i = ghostilist[ii];
int jnum = ghostnumneigh[i];
numpairs += jnum;
}
// optimized schedule setting
int time_dynamic = 0;
int time_guided = 0;
if (schedule_user == 0) schedule_user = 4;
switch (schedule_user) {
case 1:
omp_set_schedule(omp_sched_static,1);
break;
case 2:
omp_set_schedule(omp_sched_dynamic,1);
break;
case 3:
omp_set_schedule(omp_sched_guided,2);
break;
case 4:
omp_set_schedule(omp_sched_auto,0);
break;
case 5:
if (numpairs < 8*nthreads) omp_set_schedule(omp_sched_dynamic,1);
else if (schedule_time_guided < 0.0) {
omp_set_schedule(omp_sched_guided,2);
if (!eflag && !vflag) time_guided = 1;
} else if (schedule_time_dynamic<0.0) {
omp_set_schedule(omp_sched_dynamic,1);
if (!eflag && !vflag) time_dynamic = 1;
} else if (schedule_time_guided<schedule_time_dynamic)
omp_set_schedule(omp_sched_guided,2);
else
omp_set_schedule(omp_sched_dynamic,1);
break;
}
if (use_shared_arrays)
build_per_atom_arrays();
#if defined(_OPENMP)
#pragma omp parallel shared(eflag,vflag,time_dynamic,time_guided) firstprivate(numpairs) default(none)
#endif
{
// begin of pragma omp parallel
int tid = omp_get_thread_num();
int** pairs_tid_unique = NULL;
int** pairs;
if (use_shared_arrays) pairs = i_pairs;
else {
memory->create(pairs_tid_unique,numpairs,4,"numpairs");
pairs = pairs_tid_unique;
}
if (!use_shared_arrays) {
numpairs = 0;
for (int ii = 0; ii < list->inum; ii++) {
if ((do_load_balance <= 0) || ilistmask[ii]) {
int i = list->ilist[ii];
int jnum = list->numneigh[i];
for (int jj = 0; jj<jnum; jj++) {
pairs[numpairs][0] = i;
pairs[numpairs][1] = jj;
pairs[numpairs][2] = -1;
numpairs++;
}
}
}
for (int ii = 0; ii < ghostinum; ii++) {
int i = ghostilist[ii];
int jnum = ghostnumneigh[i];
for (int jj = 0; jj<jnum; jj++) {
pairs[numpairs][0] = i;
pairs[numpairs][1] = jj;
pairs[numpairs][2] = -1;
numpairs++;
}
}
}
int ielem;
- int jj,k,inum,jnum,jtype,ninside;
+ int jj,k,jnum,jtype,ninside;
double delx,dely,delz,evdwl,rsq;
double fij[3];
- int *ilist,*jlist,*numneigh,**firstneigh;
+ int *jlist,*numneigh,**firstneigh;
evdwl = 0.0;
#if defined(_OPENMP)
#pragma omp master
#endif
{
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
}
#if defined(_OPENMP)
#pragma omp barrier
{ ; }
#endif
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
- inum = list->inum;
- ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
#ifdef TIMING_INFO
// only update micro timers after setup
static int count=0;
if (count<2) {
sna[tid]->timers[0] = 0;
sna[tid]->timers[1] = 0;
sna[tid]->timers[2] = 0;
sna[tid]->timers[3] = 0;
sna[tid]->timers[4] = 0;
}
count++;
#endif
// did thread start working on interactions of new atom
int iold = -1;
double starttime, endtime;
if (time_dynamic || time_guided)
starttime = MPI_Wtime();
#if defined(_OPENMP)
#pragma omp for schedule(runtime)
#endif
for (int iijj = 0; iijj < numpairs; iijj++) {
int i = 0;
if (use_shared_arrays) {
i = i_pairs[iijj][0];
if (iold != i) {
set_sna_to_shared(tid,i_pairs[iijj][3]);
ielem = map[type[i]];
}
iold = i;
} else {
i = pairs[iijj][0];
if (iold != i) {
iold = i;
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
const int itype = type[i];
ielem = map[itype];
const double radi = radelem[ielem];
if (i < nlocal) {
jlist = firstneigh[i];
jnum = numneigh[i];
} else {
jlist = ghostneighs+ghostfirstneigh[i];
jnum = ghostnumneigh[i];
}
// insure rij, inside, wj, and rcutij are of size jnum
sna[tid]->grow_rij(jnum);
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
// wj = weights of neighbors of I within cutoff
// rcutij = cutoffs of neighbors of I within cutoff
// note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
ninside = 0;
for (jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
delx = x[j][0] - xtmp; //unitialised
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
int jelem = map[jtype];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) { //unitialised
sna[tid]->rij[ninside][0] = delx;
sna[tid]->rij[ninside][1] = dely;
sna[tid]->rij[ninside][2] = delz;
sna[tid]->inside[ninside] = j;
sna[tid]->wj[ninside] = wjelem[jelem];
sna[tid]->rcutij[ninside] = (radi + radelem[jelem])*rcutfac;
ninside++;
// update index list with inside index
pairs[iijj + (jj - pairs[iijj][1])][2] =
ninside-1; //unitialised
}
}
// compute Ui and Zi for atom I
sna[tid]->compute_ui(ninside); //unitialised
sna[tid]->compute_zi();
}
}
// for neighbors of I within cutoff:
// compute dUi/drj and dBi/drj
// Fij = dEi/dRj = -dEi/dRi => add to Fi, subtract from Fj
// entry into loop if inside index is set
double* coeffi = coeffelem[ielem];
if (pairs[iijj][2] >= 0) {
jj = pairs[iijj][2];
int j = sna[tid]->inside[jj];
sna[tid]->compute_duidrj(sna[tid]->rij[jj],
sna[tid]->wj[jj],sna[tid]->rcutij[jj]);
sna[tid]->compute_dbidrj();
sna[tid]->copy_dbi2dbvec();
if (!gammaoneflag) {
sna[tid]->compute_bi();
sna[tid]->copy_bi2bvec();
}
fij[0] = 0.0;
fij[1] = 0.0;
fij[2] = 0.0;
for (k = 1; k <= ncoeff; k++) {
double bgb;
if (gammaoneflag)
bgb = coeffi[k];
else bgb = coeffi[k]*
gamma*pow(sna[tid]->bvec[k-1],gamma-1.0);
fij[0] += bgb*sna[tid]->dbvec[k-1][0];
fij[1] += bgb*sna[tid]->dbvec[k-1][1];
fij[2] += bgb*sna[tid]->dbvec[k-1][2];
}
#if defined(_OPENMP)
#pragma omp critical
#endif
{
f[i][0] += fij[0];
f[i][1] += fij[1];
f[i][2] += fij[2];
f[j][0] -= fij[0];
f[j][1] -= fij[1];
f[j][2] -= fij[2];
if (evflag)
ev_tally_xyz(i,j,nlocal,newton_pair,0.0,0.0,
fij[0],fij[1],fij[2],
sna[tid]->rij[jj][0],sna[tid]->rij[jj][1],
sna[tid]->rij[jj][2]);
}
}
// evdwl = energy of atom I, sum over coeffs_k * Bi_k
// only call this for first pair of each atom i
// if atom has no pairs, eatom=0, which is wrong
if (eflag&&pairs[iijj][1] == 0) {
evdwl = coeffi[0];
if (gammaoneflag) {
sna[tid]->compute_bi();
sna[tid]->copy_bi2bvec();
for (int k = 1; k <= ncoeff; k++)
evdwl += coeffi[k]*sna[tid]->bvec[k-1];
} else
for (int k = 1; k <= ncoeff; k++)
evdwl += coeffi[k]*pow(sna[tid]->bvec[k-1],gamma);
#if defined(_OPENMP)
#pragma omp critical
#endif
ev_tally_full(i,2.0*evdwl,0.0,0.0,delx,dely,delz);
}
}
if (time_dynamic || time_guided)
endtime = MPI_Wtime();
if (time_dynamic) schedule_time_dynamic = endtime - starttime;
if (time_guided) schedule_time_guided = endtime - starttime;
if (!use_shared_arrays) memory->destroy(pairs);
}// end of pragma omp parallel
if (vflag_fdotr) virial_fdotr_compute();
}
inline int PairSNAP::equal(double* x,double* y)
{
double dist2 =
(x[0]-y[0])*(x[0]-y[0]) +
(x[1]-y[1])*(x[1]-y[1]) +
(x[2]-y[2])*(x[2]-y[2]);
if (dist2 < 1e-20) return 1;
return 0;
}
inline double PairSNAP::dist2(double* x,double* y)
{
return
(x[0]-y[0])*(x[0]-y[0]) +
(x[1]-y[1])*(x[1]-y[1]) +
(x[2]-y[2])*(x[2]-y[2]);
}
// return extra communication cutoff
// extra_cutoff = max(subdomain_length)
double PairSNAP::extra_cutoff()
{
double sublo[3],subhi[3];
if (domain->triclinic == 0) {
for (int dim = 0 ; dim < 3 ; dim++) {
sublo[dim] = domain->sublo[dim];
subhi[dim] = domain->subhi[dim];
}
} else {
domain->lamda2x(domain->sublo_lamda,sublo);
domain->lamda2x(domain->subhi_lamda,subhi);
}
double sub_size[3];
for (int dim = 0; dim < 3; dim++)
sub_size[dim] = subhi[dim] - sublo[dim];
double max_sub_size = 0;
for (int dim = 0; dim < 3; dim++)
max_sub_size = MAX(max_sub_size,sub_size[dim]);
// note: for triclinic, probably need something different
// see Comm::setup()
return max_sub_size;
}
// micro load_balancer: each MPI process will
// check with each of its 26 neighbors,
// whether an imbalance exists in the number
// of atoms to calculate forces for.
// If it does it will set ilistmask of one of
// its local atoms to zero, and send its Tag
// to the neighbor process. The neighboring process
// will check its ghost list for the
// ghost atom with the same Tag which is closest
// to its domain center, and build a
// neighborlist for this ghost atom. For this to work,
// the communication cutoff has to be
// as large as the neighbor cutoff +
// maximum subdomain length.
// Note that at most one atom is exchanged per processor pair.
// Also note that the local atom assignement
// doesn't change. This load balancer will cause
// some ghost atoms to have full neighborlists
// which are unique to PairSNAP.
// They are not part of the generally accessible neighborlist.
// At the same time corresponding local atoms on
// other MPI processes will not be
// included in the force computation since
// their ilistmask is 0. This does not effect
// any other classes which might
// access the same general neighborlist.
// Reverse communication (newton on) of forces is required.
// Currently the load balancer does two passes,
// since its exchanging atoms upstream and downstream.
void PairSNAP::load_balance()
{
double sublo[3],subhi[3];
if (domain->triclinic == 0) {
double* sublotmp = domain->sublo;
double* subhitmp = domain->subhi;
for (int dim = 0 ; dim<3 ; dim++) {
sublo[dim]=sublotmp[dim];
subhi[dim]=subhitmp[dim];
}
} else {
double* sublotmp = domain->sublo_lamda;
double* subhitmp = domain->subhi_lamda;
domain->lamda2x(sublotmp,sublo);
domain->lamda2x(subhitmp,subhi);
}
//if (list->inum==0) list->grow(atom->nmax);
int nlocal = ghostinum;
for (int i=0; i < list->inum; i++)
if (ilistmask[i]) nlocal++;
int ***grid2proc = comm->grid2proc;
int* procgrid = comm->procgrid;
int nlocal_up,nlocal_down;
MPI_Status status;
MPI_Request request;
double sub_mid[3];
for (int dim=0; dim<3; dim++)
sub_mid[dim] = (subhi[dim] + sublo[dim])/2;
if (comm->cutghostuser <
neighbor->cutneighmax+extra_cutoff())
error->all(FLERR,"Communication cutoff is too small "
"for SNAP micro load balancing.\n"
"Typically this can happen, if you change "
"the neighbor skin after your pair_style "
"command or if your box dimensions grow "
"during the run.\n"
"You need to set it via "
"'communicate single cutoff NUMBER' "
"to the needed length.");
int nrecv = ghostinum;
int totalsend = 0;
int nsend = 0;
int depth = 1;
for (int dx = -depth; dx < depth+1; dx++)
for (int dy = -depth; dy < depth+1; dy++)
for (int dz = -depth; dz < depth+1; dz++) {
if (dx == dy && dy == dz && dz == 0) continue;
int sendloc[3] = {comm->myloc[0],
comm->myloc[1], comm->myloc[2]
};
sendloc[0] += dx;
sendloc[1] += dy;
sendloc[2] += dz;
for (int dim = 0; dim < 3; dim++)
if (sendloc[dim] >= procgrid[dim])
sendloc[dim] = sendloc[dim] - procgrid[dim];
for (int dim = 0; dim < 3; dim++)
if (sendloc[dim] < 0)
sendloc[dim] = procgrid[dim] + sendloc[dim];
int recvloc[3] = {comm->myloc[0],
comm->myloc[1], comm->myloc[2]
};
recvloc[0] -= dx;
recvloc[1] -= dy;
recvloc[2] -= dz;
for (int dim = 0; dim < 3; dim++)
if (recvloc[dim] < 0)
recvloc[dim] = procgrid[dim] + recvloc[dim];
for (int dim = 0; dim < 3; dim++)
if (recvloc[dim] >= procgrid[dim])
recvloc[dim] = recvloc[dim] - procgrid[dim];
int sendproc = grid2proc[sendloc[0]][sendloc[1]][sendloc[2]];
int recvproc = grid2proc[recvloc[0]][recvloc[1]][recvloc[2]];
// two stage process, first upstream movement, then downstream
MPI_Sendrecv(&nlocal,1,MPI_INT,sendproc,0,
&nlocal_up,1,MPI_INT,recvproc,0,world,&status);
MPI_Sendrecv(&nlocal,1,MPI_INT,recvproc,0,
&nlocal_down,1,MPI_INT,sendproc,0,world,&status);
nsend = 0;
// send upstream
if (nlocal > nlocal_up+1) {
int i = totalsend++;
while(i < list->inum && ilistmask[i] == 0)
i = totalsend++;
if (i < list->inum)
MPI_Isend(&atom->tag[i],1,MPI_INT,recvproc,0,world,&request);
else {
int j = -1;
MPI_Isend(&j,1,MPI_INT,recvproc,0,world,&request);
}
if (i < list->inum) {
for (int j = 0; j < list->inum; j++)
if (list->ilist[j] == i)
ilistmask[j] = 0;
nsend = 1;
}
}
// recv downstream
if (nlocal < nlocal_down-1) {
nlocal++;
int get_tag = -1;
MPI_Recv(&get_tag,1,MPI_INT,sendproc,0,world,&status);
// if get_tag -1 the other process didnt have local atoms to send
if (get_tag >= 0) {
if (ghostinum >= ghostilist_max) {
memory->grow(ghostilist,ghostinum+10,
"PairSnap::ghostilist");
ghostilist_max = ghostinum+10;
}
if (atom->nlocal + atom->nghost >= ghostnumneigh_max) {
ghostnumneigh_max = atom->nlocal+atom->nghost+100;
memory->grow(ghostnumneigh,ghostnumneigh_max,
"PairSnap::ghostnumneigh");
memory->grow(ghostfirstneigh,ghostnumneigh_max,
"PairSnap::ghostfirstneigh");
}
// find closest ghost image of the transfered particle
double mindist = 1e200;
int closestghost = -1;
for (int j = 0; j < atom->nlocal + atom->nghost; j++)
if (atom->tag[j] == get_tag)
if (dist2(sub_mid, atom->x[j]) < mindist) {
closestghost = j;
mindist = dist2(sub_mid, atom->x[j]);
}
// build neighborlist for this particular
// ghost atom, and add it to list->ilist
if (ghostneighs_max - ghostneighs_total <
neighbor->oneatom) {
memory->grow(ghostneighs,
ghostneighs_total + neighbor->oneatom,
"PairSnap::ghostneighs");
ghostneighs_max = ghostneighs_total + neighbor->oneatom;
}
int j = closestghost;
ghostilist[ghostinum] = j;
ghostnumneigh[j] = 0;
ghostfirstneigh[j] = ghostneighs_total;
ghostinum++;
int* jlist = ghostneighs + ghostfirstneigh[j];
// find all neighbors by looping
// over all local and ghost atoms
for (int k = 0; k < atom->nlocal + atom->nghost; k++)
if (dist2(atom->x[j],atom->x[k]) <
neighbor->cutneighmax*neighbor->cutneighmax) {
jlist[ghostnumneigh[j]] = k;
ghostnumneigh[j]++;
ghostneighs_total++;
}
}
if (get_tag >= 0) nrecv++;
}
// decrease nlocal later, so that it is the
// initial number both for receiving and sending
if (nsend) nlocal--;
// second pass through the grid
MPI_Sendrecv(&nlocal,1,MPI_INT,sendproc,0,
&nlocal_up,1,MPI_INT,recvproc,0,world,&status);
MPI_Sendrecv(&nlocal,1,MPI_INT,recvproc,0,
&nlocal_down,1,MPI_INT,sendproc,0,world,&status);
// send downstream
nsend=0;
if (nlocal > nlocal_down+1) {
int i = totalsend++;
while(i < list->inum && ilistmask[i]==0) i = totalsend++;
if (i < list->inum)
MPI_Isend(&atom->tag[i],1,MPI_INT,sendproc,0,world,&request);
else {
int j =- 1;
MPI_Isend(&j,1,MPI_INT,sendproc,0,world,&request);
}
if (i < list->inum) {
for (int j=0; j<list->inum; j++)
if (list->ilist[j] == i) ilistmask[j] = 0;
nsend = 1;
}
}
// receive upstream
if (nlocal < nlocal_up-1) {
nlocal++;
int get_tag = -1;
MPI_Recv(&get_tag,1,MPI_INT,recvproc,0,world,&status);
if (get_tag >= 0) {
if (ghostinum >= ghostilist_max) {
memory->grow(ghostilist,ghostinum+10,
"PairSnap::ghostilist");
ghostilist_max = ghostinum+10;
}
if (atom->nlocal + atom->nghost >= ghostnumneigh_max) {
ghostnumneigh_max = atom->nlocal + atom->nghost + 100;
memory->grow(ghostnumneigh,ghostnumneigh_max,
"PairSnap::ghostnumneigh");
memory->grow(ghostfirstneigh,ghostnumneigh_max,
"PairSnap::ghostfirstneigh");
}
// find closest ghost image of the transfered particle
double mindist = 1e200;
int closestghost = -1;
for (int j = 0; j < atom->nlocal + atom->nghost; j++)
if (atom->tag[j] == get_tag)
if (dist2(sub_mid,atom->x[j])<mindist) {
closestghost = j;
mindist = dist2(sub_mid,atom->x[j]);
}
// build neighborlist for this particular ghost atom
if (ghostneighs_max-ghostneighs_total < neighbor->oneatom) {
memory->grow(ghostneighs,ghostneighs_total + neighbor->oneatom,
"PairSnap::ghostneighs");
ghostneighs_max = ghostneighs_total + neighbor->oneatom;
}
int j = closestghost;
ghostilist[ghostinum] = j;
ghostnumneigh[j] = 0;
ghostfirstneigh[j] = ghostneighs_total;
ghostinum++;
int* jlist = ghostneighs + ghostfirstneigh[j];
for (int k = 0; k < atom->nlocal + atom->nghost; k++)
if (dist2(atom->x[j],atom->x[k]) <
neighbor->cutneighmax*neighbor->cutneighmax) {
jlist[ghostnumneigh[j]] = k;
ghostnumneigh[j]++;
ghostneighs_total++;
}
}
if (get_tag >= 0) nrecv++;
}
if (nsend) nlocal--;
}
}
void PairSNAP::set_sna_to_shared(int snaid,int i)
{
sna[snaid]->rij = i_rij[i];
sna[snaid]->inside = i_inside[i];
sna[snaid]->wj = i_wj[i];
sna[snaid]->rcutij = i_rcutij[i];
sna[snaid]->zarray_r = i_zarray_r[i];
sna[snaid]->zarray_i = i_zarray_i[i];
sna[snaid]->uarraytot_r = i_uarraytot_r[i];
sna[snaid]->uarraytot_i = i_uarraytot_i[i];
}
void PairSNAP::build_per_atom_arrays()
{
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&starttime);
#endif
- int i_list[list->inum+ghostinum];
int count = 0;
int neighmax = 0;
for (int ii = 0; ii < list->inum; ii++)
if ((do_load_balance <= 0) || ilistmask[ii]) {
neighmax=MAX(neighmax,list->numneigh[list->ilist[ii]]);
- i_list[count++]=list->ilist[ii];
+ ++count;
}
for (int ii = 0; ii < ghostinum; ii++) {
neighmax=MAX(neighmax,ghostnumneigh[ghostilist[ii]]);
- i_list[count++]=ghostilist[ii];
+ ++count;
}
if (i_max < count || i_neighmax < neighmax) {
int i_maxt = MAX(count,i_max);
i_neighmax = MAX(neighmax,i_neighmax);
memory->destroy(i_rij);
memory->destroy(i_inside);
memory->destroy(i_wj);
memory->destroy(i_rcutij);
memory->destroy(i_ninside);
memory->destroy(i_pairs);
memory->create(i_rij,i_maxt,i_neighmax,3,"PairSNAP::i_rij");
memory->create(i_inside,i_maxt,i_neighmax,"PairSNAP::i_inside");
memory->create(i_wj,i_maxt,i_neighmax,"PairSNAP::i_wj");
memory->create(i_rcutij,i_maxt,i_neighmax,"PairSNAP::i_rcutij");
memory->create(i_ninside,i_maxt,"PairSNAP::i_ninside");
memory->create(i_pairs,i_maxt*i_neighmax,4,"PairSNAP::i_pairs");
}
if (i_max < count) {
int jdim = sna[0]->twojmax+1;
memory->destroy(i_uarraytot_r);
memory->destroy(i_uarraytot_i);
memory->create(i_uarraytot_r,count,jdim,jdim,jdim,
"PairSNAP::i_uarraytot_r");
memory->create(i_uarraytot_i,count,jdim,jdim,jdim,
"PairSNAP::i_uarraytot_i");
if (i_zarray_r != NULL)
for (int i = 0; i < i_max; i++) {
memory->destroy(i_zarray_r[i]);
memory->destroy(i_zarray_i[i]);
}
delete [] i_zarray_r;
delete [] i_zarray_i;
i_zarray_r = new double*****[count];
i_zarray_i = new double*****[count];
for (int i = 0; i < count; i++) {
memory->create(i_zarray_r[i],jdim,jdim,jdim,jdim,jdim,
"PairSNAP::i_zarray_r");
memory->create(i_zarray_i[i],jdim,jdim,jdim,jdim,jdim,
"PairSNAP::i_zarray_i");
}
}
if (i_max < count)
i_max = count;
count = 0;
i_numpairs = 0;
for (int ii = 0; ii < list->inum; ii++) {
if ((do_load_balance <= 0) || ilistmask[ii]) {
int i = list->ilist[ii];
int jnum = list->numneigh[i];
int* jlist = list->firstneigh[i];
const double xtmp = atom->x[i][0];
const double ytmp = atom->x[i][1];
const double ztmp = atom->x[i][2];
const int itype = atom->type[i];
const int ielem = map[itype];
const double radi = radelem[ielem];
int ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double delx = atom->x[j][0] - xtmp;
const double dely = atom->x[j][1] - ytmp;
const double delz = atom->x[j][2] - ztmp;
const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = atom->type[j];
int jelem = map[jtype];
i_pairs[i_numpairs][0] = i;
i_pairs[i_numpairs][1] = jj;
i_pairs[i_numpairs][2] = -1;
i_pairs[i_numpairs][3] = count;
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
i_rij[count][ninside][0] = delx;
i_rij[count][ninside][1] = dely;
i_rij[count][ninside][2] = delz;
i_inside[count][ninside] = j;
i_wj[count][ninside] = wjelem[jelem];
i_rcutij[count][ninside] = (radi + radelem[jelem])*rcutfac;
// update index list with inside index
i_pairs[i_numpairs][2] = ninside++;
}
i_numpairs++;
}
i_ninside[count] = ninside;
count++;
}
}
for (int ii = 0; ii < ghostinum; ii++) {
int i = ghostilist[ii];
int jnum = ghostnumneigh[i];
int* jlist = ghostneighs+ghostfirstneigh[i];
const double xtmp = atom->x[i][0];
const double ytmp = atom->x[i][1];
const double ztmp = atom->x[i][2];
const int itype = atom->type[i];
const int ielem = map[itype];
const double radi = radelem[ielem];
int ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double delx = atom->x[j][0] - xtmp;
const double dely = atom->x[j][1] - ytmp;
const double delz = atom->x[j][2] - ztmp;
const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = atom->type[j];
int jelem = map[jtype];
i_pairs[i_numpairs][0] = i;
i_pairs[i_numpairs][1] = jj;
i_pairs[i_numpairs][2] = -1;
i_pairs[i_numpairs][3] = count;
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
i_rij[count][ninside][0] = delx;
i_rij[count][ninside][1] = dely;
i_rij[count][ninside][2] = delz;
i_inside[count][ninside] = j;
i_wj[count][ninside] = wjelem[jelem];
i_rcutij[count][ninside] = (radi + radelem[jelem])*rcutfac;
// update index list with inside index
i_pairs[i_numpairs][2] = ninside++;
}
i_numpairs++;
}
i_ninside[count] = ninside;
count++;
}
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&endtime);
timers[0]+=(endtime.tv_sec-starttime.tv_sec+1.0*
(endtime.tv_nsec-starttime.tv_nsec)/1000000000);
#endif
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&starttime);
#endif
#if defined(_OPENMP)
#pragma omp parallel for shared(count) default(none)
#endif
for (int ii=0; ii < count; ii++) {
int tid = omp_get_thread_num();
set_sna_to_shared(tid,ii);
//sna[tid]->compute_ui(i_ninside[ii]);
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&starttime);
#endif
sna[tid]->compute_ui_omp(i_ninside[ii],MAX(int(nthreads/count),1));
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&endtime);
sna[tid]->timers[0]+=(endtime.tv_sec-starttime.tv_sec+1.0*
(endtime.tv_nsec-starttime.tv_nsec)/1000000000);
#endif
}
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&starttime);
#endif
for (int ii=0; ii < count; ii++) {
int tid = 0;//omp_get_thread_num();
set_sna_to_shared(tid,ii);
sna[tid]->compute_zi_omp(MAX(int(nthreads/count),1));
}
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&endtime);
sna[0]->timers[1]+=(endtime.tv_sec-starttime.tv_sec+1.0*
(endtime.tv_nsec-starttime.tv_nsec)/1000000000);
#endif
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&endtime);
timers[1]+=(endtime.tv_sec-starttime.tv_sec+1.0*
(endtime.tv_nsec-starttime.tv_nsec)/1000000000);
#endif
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSNAP::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(map,n+1,"pair:map");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSNAP::settings(int narg, char **arg)
{
// set default values for optional arguments
nthreads = -1;
use_shared_arrays=-1;
do_load_balance = 0;
use_optimized = 1;
// optional arguments
for (int i=0; i < narg; i++) {
if (i+2>narg) error->all(FLERR,"Illegal pair_style command."
" Too few arguments.");
if (strcmp(arg[i],"nthreads")==0) {
nthreads=force->inumeric(FLERR,arg[++i]);
#if defined(LMP_USER_OMP)
error->all(FLERR,"Please set number of threads via package omp command");
#else
omp_set_num_threads(nthreads);
comm->nthreads=nthreads;
#endif
continue;
}
if (strcmp(arg[i],"optimized")==0) {
use_optimized=force->inumeric(FLERR,arg[++i]);
continue;
}
if (strcmp(arg[i],"shared")==0) {
use_shared_arrays=force->inumeric(FLERR,arg[++i]);
continue;
}
if (strcmp(arg[i],"loadbalance")==0) {
do_load_balance = force->inumeric(FLERR,arg[++i]);
if (do_load_balance) {
double mincutoff = extra_cutoff() +
rcutmax + neighbor->skin;
if (comm->cutghostuser < mincutoff) {
char buffer[255];
//apparently mincutoff is 0 after sprintf command ?????
double tmp = mincutoff + 0.1;
sprintf(buffer, "Communication cutoff is too small "
"for SNAP micro load balancing. "
"It will be increased to: %lf.",mincutoff+0.1);
if (comm->me==0)
error->warning(FLERR,buffer);
comm->cutghostuser = tmp;
}
}
continue;
}
if (strcmp(arg[i],"schedule")==0) {
i++;
if (strcmp(arg[i],"static")==0)
schedule_user = 1;
if (strcmp(arg[i],"dynamic")==0)
schedule_user = 2;
if (strcmp(arg[i],"guided")==0)
schedule_user = 3;
if (strcmp(arg[i],"auto")==0)
schedule_user = 4;
if (strcmp(arg[i],"determine")==0)
schedule_user = 5;
if (schedule_user == 0)
error->all(FLERR,"Illegal pair_style command."
" Illegal schedule argument.");
continue;
}
char buffer[255];
sprintf(buffer, "Illegal pair_style command."
" Unrecognized argument: %s.\n",arg[i]);
error->all(FLERR,buffer);
}
if (nthreads < 0)
nthreads = comm->nthreads;
if (use_shared_arrays < 0) {
if (nthreads > 1 && atom->nlocal <= 2*nthreads)
use_shared_arrays = 1;
else use_shared_arrays = 0;
}
// check if running non-optimized code with
// optimization flags set
if (!use_optimized)
if (nthreads > 1 ||
use_shared_arrays ||
do_load_balance ||
schedule_user)
error->all(FLERR,"Illegal pair_style command."
"Advanced options require setting 'optimized 1'.");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairSNAP::coeff(int narg, char **arg)
{
// read SNAP element names between 2 filenames
// nelements = # of SNAP elements
// elements = list of unique element names
if (narg < 6) error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
if (nelements) {
for (int i = 0; i < nelements; i++)
delete[] elements[i];
delete[] elements;
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(coeffelem);
}
nelements = narg - 4 - atom->ntypes;
if (nelements < 1) error->all(FLERR,"Incorrect args for pair coefficients");
char* type1 = arg[0];
char* type2 = arg[1];
char* coefffilename = arg[2];
char** elemlist = &arg[3];
char* paramfilename = arg[3+nelements];
char** elemtypes = &arg[4+nelements];
// insure I,J args are * *
if (strcmp(type1,"*") != 0 || strcmp(type2,"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
elements = new char*[nelements];
for (int i = 0; i < nelements; i++) {
char* elemname = elemlist[i];
int n = strlen(elemname) + 1;
elements[i] = new char[n];
strcpy(elements[i],elemname);
}
// read snapcoeff and snapparam files
read_files(coefffilename,paramfilename);
// read args that map atom types to SNAP elements
// map[i] = which element the Ith atom type is, -1 if not mapped
// map[0] is not used
for (int i = 1; i <= atom->ntypes; i++) {
char* elemname = elemtypes[i-1];
int jelem;
for (jelem = 0; jelem < nelements; jelem++)
if (strcmp(elemname,elements[jelem]) == 0)
break;
if (jelem < nelements)
map[i] = jelem;
else if (strcmp(elemname,"NULL") == 0) map[i] = -1;
else error->all(FLERR,"Incorrect args for pair coefficients");
}
// clear setflag since coeff() called once with I,J = * *
int n = atom->ntypes;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
sna = new SNA*[nthreads];
// allocate memory for per OpenMP thread data which
// is wrapped into the sna class
#if defined(_OPENMP)
#pragma omp parallel default(none)
#endif
{
int tid = omp_get_thread_num();
sna[tid] = new SNA(lmp,rfac0,twojmax,
diagonalstyle,use_shared_arrays,
rmin0,switchflag);
if (!use_shared_arrays)
sna[tid]->grow_rij(nmax);
}
if (ncoeff != sna[0]->ncoeff) {
printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff);
error->all(FLERR,"Incorrect SNAP parameter file");
}
// Calculate maximum cutoff for all elements
rcutmax = 0.0;
for (int ielem = 0; ielem < nelements; ielem++)
rcutmax = MAX(2.0*radelem[ielem]*rcutfac,rcutmax);
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairSNAP::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style SNAP requires newton pair on");
// need a full neighbor list
int irequest = neighbor->request(this);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
#if defined(_OPENMP)
#pragma omp parallel default(none)
#endif
{
int tid = omp_get_thread_num();
sna[tid]->init();
}
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSNAP::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
return (radelem[map[i]] +
radelem[map[j]])*rcutfac;
}
/* ---------------------------------------------------------------------- */
void PairSNAP::read_files(char *coefffilename, char *paramfilename)
{
// open SNAP ceofficient file on proc 0
FILE *fpcoeff;
if (comm->me == 0) {
fpcoeff = force->open_potential(coefffilename);
if (fpcoeff == NULL) {
char str[128];
sprintf(str,"Cannot open SNAP coefficient file %s",coefffilename);
error->one(FLERR,str);
}
}
char line[MAXLINE],*ptr;
int eof = 0;
int n;
int nwords = 0;
while (nwords == 0) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpcoeff);
if (ptr == NULL) {
eof = 1;
fclose(fpcoeff);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
}
if (nwords != 2)
error->all(FLERR,"Incorrect format in SNAP coefficient file");
// words = ptrs to all words in line
// strip single and double quotes from words
char* words[MAXWORD];
int iword = 0;
words[iword] = strtok(line,"' \t\n\r\f");
iword = 1;
words[iword] = strtok(NULL,"' \t\n\r\f");
int nelemfile = atoi(words[0]);
ncoeff = atoi(words[1])-1;
// Set up element lists
memory->create(radelem,nelements,"pair:radelem");
memory->create(wjelem,nelements,"pair:wjelem");
memory->create(coeffelem,nelements,ncoeff+1,"pair:coeffelem");
int *found = new int[nelements];
for (int ielem = 0; ielem < nelements; ielem++)
found[ielem] = 0;
// Loop over elements in the SNAP coefficient file
for (int ielemfile = 0; ielemfile < nelemfile; ielemfile++) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpcoeff);
if (ptr == NULL) {
eof = 1;
fclose(fpcoeff);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof)
error->all(FLERR,"Incorrect format in SNAP coefficient file");
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
nwords = atom->count_words(line);
if (nwords != 3)
error->all(FLERR,"Incorrect format in SNAP coefficient file");
iword = 0;
words[iword] = strtok(line,"' \t\n\r\f");
iword = 1;
words[iword] = strtok(NULL,"' \t\n\r\f");
iword = 2;
words[iword] = strtok(NULL,"' \t\n\r\f");
char* elemtmp = words[0];
double radtmp = atof(words[1]);
double wjtmp = atof(words[2]);
// skip if element name isn't in element list
int ielem;
for (ielem = 0; ielem < nelements; ielem++)
if (strcmp(elemtmp,elements[ielem]) == 0) break;
if (ielem == nelements) {
if (comm->me == 0)
for (int icoeff = 0; icoeff <= ncoeff; icoeff++)
ptr = fgets(line,MAXLINE,fpcoeff);
continue;
}
// skip if element already appeared
if (found[ielem]) {
if (comm->me == 0)
for (int icoeff = 0; icoeff <= ncoeff; icoeff++)
ptr = fgets(line,MAXLINE,fpcoeff);
continue;
}
found[ielem] = 1;
radelem[ielem] = radtmp;
wjelem[ielem] = wjtmp;
if (comm->me == 0) {
if (screen) fprintf(screen,"SNAP Element = %s, Radius %g, Weight %g \n",
elements[ielem], radelem[ielem], wjelem[ielem]);
if (logfile) fprintf(logfile,"SNAP Element = %s, Radius %g, Weight %g \n",
elements[ielem], radelem[ielem], wjelem[ielem]);
}
for (int icoeff = 0; icoeff <= ncoeff; icoeff++) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpcoeff);
if (ptr == NULL) {
eof = 1;
fclose(fpcoeff);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof)
error->all(FLERR,"Incorrect format in SNAP coefficient file");
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
nwords = atom->count_words(line);
if (nwords != 1)
error->all(FLERR,"Incorrect format in SNAP coefficient file");
iword = 0;
words[iword] = strtok(line,"' \t\n\r\f");
coeffelem[ielem][icoeff] = atof(words[0]);
}
}
// set flags for required keywords
rcutfacflag = 0;
twojmaxflag = 0;
// Set defaults for optional keywords
gamma = 1.0;
gammaoneflag = 1;
rfac0 = 0.99363;
rmin0 = 0.0;
diagonalstyle = 3;
switchflag = 1;
// open SNAP parameter file on proc 0
FILE *fpparam;
if (comm->me == 0) {
fpparam = force->open_potential(paramfilename);
if (fpparam == NULL) {
char str[128];
sprintf(str,"Cannot open SNAP parameter file %s",paramfilename);
error->one(FLERR,str);
}
}
eof = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpparam);
if (ptr == NULL) {
eof = 1;
fclose(fpparam);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
if (nwords != 2)
error->all(FLERR,"Incorrect format in SNAP parameter file");
// words = ptrs to all words in line
// strip single and double quotes from words
char* keywd = strtok(line,"' \t\n\r\f");
char* keyval = strtok(NULL,"' \t\n\r\f");
if (comm->me == 0) {
if (screen) fprintf(screen,"SNAP keyword %s %s \n",keywd,keyval);
if (logfile) fprintf(logfile,"SNAP keyword %s %s \n",keywd,keyval);
}
if (strcmp(keywd,"rcutfac") == 0) {
rcutfac = atof(keyval);
rcutfacflag = 1;
} else if (strcmp(keywd,"twojmax") == 0) {
twojmax = atoi(keyval);
twojmaxflag = 1;
} else if (strcmp(keywd,"gamma") == 0)
gamma = atof(keyval);
else if (strcmp(keywd,"rfac0") == 0)
rfac0 = atof(keyval);
else if (strcmp(keywd,"rmin0") == 0)
rmin0 = atof(keyval);
else if (strcmp(keywd,"diagonalstyle") == 0)
diagonalstyle = atoi(keyval);
else if (strcmp(keywd,"switchflag") == 0)
switchflag = atoi(keyval);
else
error->all(FLERR,"Incorrect SNAP parameter file");
}
if (rcutfacflag == 0 || twojmaxflag == 0)
error->all(FLERR,"Incorrect SNAP parameter file");
if (gamma == 1.0) gammaoneflag = 1;
else gammaoneflag = 0;
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double PairSNAP::memory_usage()
{
double bytes = Pair::memory_usage();
int n = atom->ntypes+1;
bytes += n*n*sizeof(int);
bytes += n*n*sizeof(double);
bytes += 3*nmax*sizeof(double);
bytes += nmax*sizeof(int);
bytes += (2*ncoeff+1)*sizeof(double);
bytes += (ncoeff*3)*sizeof(double);
bytes += sna[0]->memory_usage()*nthreads;
return bytes;
}
diff --git a/src/citeme.cpp b/src/citeme.cpp
index 70b02e83b..66c6eb9dc 100644
--- a/src/citeme.cpp
+++ b/src/citeme.cpp
@@ -1,75 +1,77 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "citeme.h"
#include "version.h"
#include "universe.h"
#include "error.h"
using namespace LAMMPS_NS;
static const char cite_header[] =
"This LAMMPS simulation made specific use of work described in the\n"
"following references. See http://lammps.sandia.gov/cite.html\n"
"for details.\n\n";
static const char cite_nagline[] = "\nPlease see the log.cite file "
"for references relevant to this simulation\n\n";
/* ---------------------------------------------------------------------- */
CiteMe::CiteMe(LAMMPS *lmp) : Pointers(lmp)
{
fp = NULL;
cs = new citeset();
}
/* ----------------------------------------------------------------------
write out nag-line at the end of the regular output and clean up
------------------------------------------------------------------------- */
CiteMe::~CiteMe()
{
if (universe->me || cs->size() == 0) {
delete cs;
return;
}
delete cs;
- if (screen) fprintf(screen,cite_nagline);
- if (logfile) fprintf(logfile,cite_nagline);
+ if (fp) {
+ if (screen) fprintf(screen,cite_nagline);
+ if (logfile) fprintf(logfile,cite_nagline);
- if (fp) fclose(fp);
+ fclose(fp);
+ }
}
/* ----------------------------------------------------------------------
write out and register a citation so it will be written only once
------------------------------------------------------------------------- */
void CiteMe::add(const char *ref)
{
if (universe->me) return;
if (cs->find(ref) != cs->end()) return;
cs->insert(ref);
if (!fp) {
fp = fopen("log.cite","w");
- if (!fp) error->universe_one(FLERR,"Cannot open log.cite file");
+ if (!fp) return;
fputs(cite_header,fp);
fflush(fp);
}
fputs(ref,fp);
fflush(fp);
}
diff --git a/src/fix_recenter.cpp b/src/fix_recenter.cpp
index 2d8ee41f9..5d2e2d5e8 100644
--- a/src/fix_recenter.cpp
+++ b/src/fix_recenter.cpp
@@ -1,216 +1,232 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U)
------------------------------------------------------------------------- */
#include "stdlib.h"
#include "string.h"
#include "fix_recenter.h"
#include "atom.h"
#include "group.h"
+#include "update.h"
#include "domain.h"
#include "lattice.h"
#include "modify.h"
#include "comm.h"
+#include "respa.h"
#include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{BOX,LATTICE,FRACTION};
/* ---------------------------------------------------------------------- */
FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 6) error->all(FLERR,"Illegal fix recenter command");
xcom = ycom = zcom = 0.0;
xflag = yflag = zflag = 1;
xinitflag = yinitflag = zinitflag = 0;
shift[0] = shift[1] = shift[2] = 0.0;
distance = 0.0;
scalar_flag = 1;
vector_flag = 1;
size_vector = 3;
extscalar = 1;
extvector = 1;
global_freq = 1;
/* ---------------------------------------------------------------------- */
if (strcmp(arg[3],"NULL") == 0) xflag = 0;
else if (strcmp(arg[3],"INIT") == 0) xinitflag = 1;
else xcom = force->numeric(FLERR,arg[3]);
if (strcmp(arg[4],"NULL") == 0) yflag = 0;
else if (strcmp(arg[4],"INIT") == 0) yinitflag = 1;
else ycom = force->numeric(FLERR,arg[4]);
if (strcmp(arg[5],"NULL") == 0) zflag = 0;
else if (strcmp(arg[5],"INIT") == 0) zinitflag = 1;
else zcom = force->numeric(FLERR,arg[5]);
// optional args
group2bit = groupbit;
scaleflag = LATTICE;
int iarg = 6;
while (iarg < narg) {
if (strcmp(arg[iarg],"shift") == 0) {
int igroup2 = group->find(arg[iarg+1]);
if (igroup2 < 0) error->all(FLERR,"Could not find fix recenter group ID");
group2bit = group->bitmask[igroup2];
iarg += 2;
} else if (strcmp(arg[iarg],"units") == 0) {
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = BOX;
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = LATTICE;
else if (strcmp(arg[iarg+1],"fraction") == 0) scaleflag = FRACTION;
else error->all(FLERR,"Illegal fix recenter command");
iarg += 2;
} else error->all(FLERR,"Illegal fix recenter command");
}
// scale xcom,ycom,zcom
double xscale,yscale,zscale;
if (scaleflag == LATTICE) {
xscale = domain->lattice->xlattice;
yscale = domain->lattice->ylattice;
zscale = domain->lattice->zlattice;
}
else xscale = yscale = zscale = 1.0;
xcom *= xscale;
ycom *= yscale;
zcom *= zscale;
// cannot have 0 atoms in group
if (group->count(igroup) == 0)
error->all(FLERR,"Fix recenter group has no atoms");
}
/* ---------------------------------------------------------------------- */
int FixRecenter::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
+ mask |= INITIAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixRecenter::init()
{
// warn if any integrate fix comes after this one
int after = 0;
int flag = 0;
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(id,modify->fix[i]->id) == 0) after = 1;
else if ((modify->fmask[i] & INITIAL_INTEGRATE) && after) flag = 1;
}
if (flag && comm->me == 0)
error->warning(FLERR,"Fix recenter should come after all other integration fixes");
masstotal = group->mass(igroup);
// if any components of requested COM were INIT, store initial COM
if (xinitflag || yinitflag || zinitflag) {
double xcm[3];
group->xcm(igroup,masstotal,xcm);
xinit = xcm[0];
yinit = xcm[1];
zinit = xcm[2];
}
+
+ if (strstr(update->integrate_style,"respa"))
+ nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixRecenter::initial_integrate(int vflag)
{
// target COM
// bounding box around domain works for both orthogonal and triclinic
double xtarget,ytarget,ztarget;
double *bboxlo,*bboxhi;
if (scaleflag == FRACTION) {
if (domain->triclinic == 0) {
bboxlo = domain->boxlo;
bboxhi = domain->boxhi;
} else {
bboxlo = domain->boxlo_bound;
bboxhi = domain->boxhi_bound;
}
}
if (xinitflag) xtarget = xinit;
else if (scaleflag == FRACTION)
xtarget = bboxlo[0] + xcom*(bboxhi[0] - bboxlo[0]);
else xtarget = xcom;
if (yinitflag) ytarget = yinit;
else if (scaleflag == FRACTION)
ytarget = bboxlo[1] + ycom*(bboxhi[1] - bboxlo[1]);
else ytarget = ycom;
if (zinitflag) ztarget = zinit;
else if (scaleflag == FRACTION)
ztarget = bboxlo[2] + zcom*(bboxhi[2] - bboxlo[2]);
else ztarget = zcom;
// current COM
double xcm[3];
group->xcm(igroup,masstotal,xcm);
// shift coords by difference between actual COM and requested COM
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
shift[0] = xflag ? (xtarget - xcm[0]) : 0.0;
shift[1] = yflag ? (ytarget - xcm[1]) : 0.0;
shift[2] = zflag ? (ztarget - xcm[2]) : 0.0;
distance = sqrt(shift[0]*shift[0] + shift[1]*shift[1] + shift[2]*shift[2]);
for (int i = 0; i < nlocal; i++)
if (mask[i] & group2bit) {
x[i][0] += shift[0];
x[i][1] += shift[1];
x[i][2] += shift[2];
}
}
/* ---------------------------------------------------------------------- */
double FixRecenter::compute_scalar()
{
return distance;
}
/* ---------------------------------------------------------------------- */
double FixRecenter::compute_vector(int n)
{
return shift[n];
}
+/* ---------------------------------------------------------------------- */
+
+void FixRecenter::initial_integrate_respa(int vflag, int ilevel, int iloop)
+{
+ // outermost level - operate recenter
+ // all other levels - nothing
+
+ if (ilevel == nlevels_respa-1) initial_integrate(vflag);
+}
+