diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp deleted file mode 100644 index e5e59acc53..0000000000 --- a/src/fix_deposit.cpp +++ /dev/null @@ -1,704 +0,0 @@ -/* ---------------------------------------------------------------------- - 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_deposit.h" -#include "atom.h" -#include "atom_vec.h" -#include "molecule.h" -#include "force.h" -#include "update.h" -#include "modify.h" -#include "fix.h" -#include "comm.h" -#include "domain.h" -#include "lattice.h" -#include "region.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}; - -/* ---------------------------------------------------------------------- */ - -FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg < 7) error->all(FLERR,"Illegal fix deposit command"); - - restart_global = 1; - time_depend = 1; - - // required args - - ninsert = force->inumeric(FLERR,arg[3]); - ntype = force->inumeric(FLERR,arg[4]); - nfreq = force->inumeric(FLERR,arg[5]); - seed = force->inumeric(FLERR,arg[6]); - - if (seed <= 0) error->all(FLERR,"Illegal fix deposit command"); - - // read options from end of input line - - options(narg-7,&arg[7]); - - // error check on type - - if (mode == ATOM && (ntype <= 0 || ntype > atom->ntypes)) - error->all(FLERR,"Invalid atom type in fix deposit command"); - - // error checks on region and its extent being inside simulation box - - if (iregion == -1) error->all(FLERR,"Must specify a region in fix deposit"); - if (domain->regions[iregion]->bboxflag == 0) - error->all(FLERR,"Fix deposit region does not support a bounding box"); - if (domain->regions[iregion]->dynamic_check()) - error->all(FLERR,"Fix deposit region cannot be dynamic"); - - xlo = domain->regions[iregion]->extent_xlo; - xhi = domain->regions[iregion]->extent_xhi; - ylo = domain->regions[iregion]->extent_ylo; - yhi = domain->regions[iregion]->extent_yhi; - zlo = domain->regions[iregion]->extent_zlo; - zhi = domain->regions[iregion]->extent_zhi; - - if (domain->triclinic == 0) { - 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,"Deposition region extends outside simulation box"); - } else { - if (xlo < domain->boxlo_bound[0] || xhi > domain->boxhi_bound[0] || - ylo < domain->boxlo_bound[1] || yhi > domain->boxhi_bound[1] || - zlo < domain->boxlo_bound[2] || zhi > domain->boxhi_bound[2]) - error->all(FLERR,"Deposition region extends outside simulation box"); - } - - // error check and further setup for mode = MOLECULE - - if (atom->tag_enable == 0) - error->all(FLERR,"Cannot use fix_deposit unless atoms have IDs"); - - if (mode == MOLECULE) { - if (atom->molecule_flag == 0) - error->all(FLERR,"Fix deposit requires atom attribute molecule"); - if (onemol->xflag == 0) - error->all(FLERR,"Fix deposit molecule must have coordinates"); - if (onemol->typeflag == 0) - error->all(FLERR,"Fix deposit molecule must have atom types"); - if (ntype+onemol->maxtype <= 0 || ntype+onemol->maxtype > atom->ntypes) - error->all(FLERR,"Invalid atom type in fix deposit mol command"); - - // fix deposit uses geoemetric center of molecule for insertion - - onemol->compute_center(); - } - - if (rigidflag && mode == ATOM) - error->all(FLERR,"Cannot use fix deposit rigid and not molecule"); - if (shakeflag && mode == ATOM) - error->all(FLERR,"Cannot use fix deposit shake and not molecule"); - if (rigidflag && shakeflag) - error->all(FLERR,"Cannot use fix deposit rigid and shake"); - - // setup of coords and imageflags array - - if (mode == ATOM) natom = 1; - else natom = onemol->natoms; - memory->create(coords,natom,3,"deposit:coords"); - memory->create(imageflags,natom,"deposit:imageflags"); - - // setup scaling - - double xscale,yscale,zscale; - if (scaleflag) { - xscale = domain->lattice->xlattice; - yscale = domain->lattice->ylattice; - zscale = domain->lattice->zlattice; - } - else xscale = yscale = zscale = 1.0; - - // apply scaling to all input parameters with dist/vel units - - if (domain->dimension == 2) { - lo *= yscale; - hi *= yscale; - rate *= yscale; - } else { - lo *= zscale; - hi *= zscale; - rate *= zscale; - } - deltasq *= xscale*xscale; - nearsq *= xscale*xscale; - vxlo *= xscale; - vxhi *= xscale; - vylo *= yscale; - vyhi *= yscale; - vzlo *= zscale; - vzhi *= zscale; - tx *= xscale; - ty *= yscale; - tz *= zscale; - - // 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); - - // set up reneighboring - - force_reneighbor = 1; - next_reneighbor = update->ntimestep + 1; - nfirst = next_reneighbor; - ninserted = 0; -} - -/* ---------------------------------------------------------------------- */ - -FixDeposit::~FixDeposit() -{ - delete random; - delete [] idrigid; - delete [] idshake; - delete [] idregion; - memory->destroy(coords); - memory->destroy(imageflags); -} - -/* ---------------------------------------------------------------------- */ - -int FixDeposit::setmask() -{ - int mask = 0; - mask |= PRE_EXCHANGE; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixDeposit::init() -{ - // set index and check validity of region - - iregion = domain->find_region(idregion); - if (iregion == -1) - error->all(FLERR,"Region ID for fix deposit does not exist"); - - // 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 (onemol != (Molecule *) fixrigid->extract("onemol",tmp)) - error->all(FLERR, - "Fix deposit and fix rigid/small not using same molecule 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 deposit shake fix does not exist"); - fixshake = modify->fix[ifix]; - int tmp; - if (onemol != (Molecule *) fixshake->extract("onemol",tmp)) - error->all(FLERR,"Fix deposit and fix shake not using same molecule ID"); - } -} - -/* ---------------------------------------------------------------------- - perform particle insertion -------------------------------------------------------------------------- */ - -void FixDeposit::pre_exchange() -{ - int i,j,m,n,nlocalprev,flag,flagall; - double coord[3],lamda[3],delx,dely,delz,rsq; - double alpha,beta,gamma; - double r[3],vnew[3],rotmat[3][3],quat[4]; - double *newcoord; - - // just return if should not be called on this timestep - - if (next_reneighbor != update->ntimestep) return; - - // compute current offset = bottom of insertion volume - - double offset = 0.0; - if (rateflag) offset = (update->ntimestep - nfirst) * update->dt * rate; - - double *sublo,*subhi; - if (domain->triclinic == 0) { - sublo = domain->sublo; - subhi = domain->subhi; - } else { - sublo = domain->sublo_lamda; - subhi = domain->subhi_lamda; - } - - // find current max atom and molecule IDs if necessary - - if (!idnext) find_maxid(); - - // attempt an insertion until successful - - int dimension = domain->dimension; - int nfix = modify->nfix; - Fix **fix = modify->fix; - - int success = 0; - int attempt = 0; - while (attempt < maxattempt) { - attempt++; - - // choose random position for new particle within region - - coord[0] = xlo + random->uniform() * (xhi-xlo); - coord[1] = ylo + random->uniform() * (yhi-ylo); - coord[2] = zlo + random->uniform() * (zhi-zlo); - while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { - coord[0] = xlo + random->uniform() * (xhi-xlo); - coord[1] = ylo + random->uniform() * (yhi-ylo); - coord[2] = zlo + random->uniform() * (zhi-zlo); - } - - // adjust vertical coord by offset - - if (dimension == 2) coord[1] += offset; - else coord[2] += offset; - - // if global, reset vertical coord to be lo-hi above highest atom - // if local, reset vertical coord to be lo-hi above highest "nearby" atom - // local computation computes lateral distance between 2 particles w/ PBC - // when done, have final coord of atom or center pt of molecule - - if (globalflag || localflag) { - int dim; - double max,maxall,delx,dely,delz,rsq; - - if (dimension == 2) { - dim = 1; - max = domain->boxlo[1]; - } else { - dim = 2; - max = domain->boxlo[2]; - } - - double **x = atom->x; - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) { - if (localflag) { - delx = coord[0] - x[i][0]; - dely = coord[1] - x[i][1]; - delz = 0.0; - domain->minimum_image(delx,dely,delz); - if (dimension == 2) rsq = delx*delx; - else rsq = delx*delx + dely*dely; - if (rsq > deltasq) continue; - } - if (x[i][dim] > max) max = x[i][dim]; - } - - MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world); - if (dimension == 2) - coord[1] = maxall + lo + random->uniform()*(hi-lo); - else - coord[2] = maxall + lo + random->uniform()*(hi-lo); - } - - // coords = coords of all atoms - // for molecule, perform random rotation around center pt - // apply PBC so final coords are inside box - // also store image flag modified due to PBC - - if (mode == ATOM) { - coords[0][0] = coord[0]; - coords[0][1] = coord[1]; - coords[0][2] = coord[2]; - } else { - 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,onemol->dx[i],coords[i]); - coords[i][0] += coord[0]; - coords[i][1] += coord[1]; - coords[i][2] += coord[2]; - - imageflags[i] = ((tagint) IMGMAX << IMG2BITS) | - ((tagint) IMGMAX << IMGBITS) | IMGMAX; - domain->remap(coords[i],imageflags[i]); - } - } - - // if distance to any inserted atom is less than near, try again - // use minimum_image() to account for PBC - - double **x = atom->x; - int nlocal = atom->nlocal; - - flag = 0; - for (m = 0; m < natom; m++) { - for (i = 0; i < nlocal; i++) { - delx = coords[m][0] - x[i][0]; - dely = coords[m][1] - x[i][1]; - delz = coords[m][2] - x[i][2]; - domain->minimum_image(delx,dely,delz); - rsq = delx*delx + dely*dely + delz*delz; - if (rsq < nearsq) flag = 1; - } - } - MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world); - if (flagall) continue; - - // proceed with insertion - - nlocalprev = atom->nlocal; - - // choose random velocity for new particle - // used for every atom in molecule - - vnew[0] = vxlo + random->uniform() * (vxhi-vxlo); - vnew[1] = vylo + random->uniform() * (vyhi-vylo); - vnew[2] = vzlo + random->uniform() * (vzhi-vzlo); - - // if target specified, change velocity vector accordingly - - if (targetflag) { - double vel = sqrt(vnew[0]*vnew[0] + vnew[1]*vnew[1] + vnew[2]*vnew[2]); - delx = tx - coord[0]; - dely = ty - coord[1]; - delz = tz - coord[2]; - double rsq = delx*delx + dely*dely + delz*delz; - if (rsq > 0.0) { - double rinv = sqrt(1.0/rsq); - vnew[0] = delx*rinv*vel; - vnew[1] = dely*rinv*vel; - vnew[2] = delz*rinv*vel; - } - } - - // 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 (domain->triclinic) { - domain->x2lamda(coords[m],lamda); - newcoord = lamda; - } else 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] && - 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 (dimension == 2 && newcoord[1] >= domain->boxhi[1] && - comm->myloc[1] == comm->procgrid[1]-1 && - 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+onemol->type[m],coords[m]); - n = atom->nlocal - 1; - atom->tag[n] = maxtag_all + m+1; - if (mode == MOLECULE) atom->molecule[n] = maxmol_all; - 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 == MOLECULE) atom->add_molecule_atom(onemol,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,coord,vnew,quat); - else if (shakeflag) - fixshake->set_molecule(nlocalprev,maxtag_all,coord,vnew,quat); - } - - // old code: unsuccessful if no proc performed insertion of an atom - // don't think that check is necessary - // if get this far, should always be succesful - // would be hard to undo partial insertion for a molecule - // better to check how many atoms could be inserted (w/out inserting) - // then sum to insure all are inserted, before doing actual insertion - // MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world); - - success = 1; - break; - } - - // warn if not successful b/c too many attempts - - if (!success && comm->me == 0) - error->warning(FLERR,"Particle deposition was unsuccessful",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 (success) { - atom->natoms += natom; - if (mode == MOLECULE) { - atom->nbonds += onemol->nbonds; - atom->nangles += onemol->nangles; - atom->ndihedrals += onemol->ndihedrals; - atom->nimpropers += onemol->nimpropers; - } - if (idnext) { - maxtag_all += natom; - if (mode == MOLECULE) maxmol_all++; - } - if (atom->map_style) { - atom->nghost = 0; - atom->map_init(); - atom->map_set(); - } - } - - // next timestep to insert - // next_reneighbor = 0 if done - - if (success) ninserted++; - 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 FixDeposit::find_maxid() -{ - int *tag = atom->tag; - int *molecule = atom->molecule; - int nlocal = atom->nlocal; - - int max = 0; - for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]); - MPI_Allreduce(&max,&maxtag_all,1,MPI_INT,MPI_MAX,world); - - if (mode == MOLECULE) { - max = 0; - for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]); - MPI_Allreduce(&max,&maxmol_all,1,MPI_INT,MPI_MAX,world); - } -} - -/* ---------------------------------------------------------------------- - parse optional parameters at end of input line -------------------------------------------------------------------------- */ - -void FixDeposit::options(int narg, char **arg) -{ - // defaults - - iregion = -1; - idregion = NULL; - mode = ATOM; - rigidflag = 0; - idrigid = NULL; - shakeflag = 0; - idshake = NULL; - idnext = 0; - globalflag = localflag = 0; - lo = hi = deltasq = 0.0; - nearsq = 0.0; - maxattempt = 10; - rateflag = 0; - vxlo = vxhi = vylo = vyhi = vzlo = vzhi = 0.0; - scaleflag = 1; - targetflag = 0; - - int iarg = 0; - while (iarg < narg) { - if (strcmp(arg[iarg],"region") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - iregion = domain->find_region(arg[iarg+1]); - if (iregion == -1) - error->all(FLERR,"Region ID for fix deposit does not exist"); - int n = strlen(arg[iarg+1]) + 1; - idregion = new char[n]; - strcpy(idregion,arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"mol") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - int imol = atom->find_molecule(arg[iarg+1]); - if (imol == -1) - error->all(FLERR,"Molecule ID for fix deposit does not exist"); - mode = MOLECULE; - onemol = atom->molecules[imol]; - iarg += 2; - } else if (strcmp(arg[iarg],"rigid") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit 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 deposit 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 deposit 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 deposit command"); - iarg += 2; - } else if (strcmp(arg[iarg],"global") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command"); - globalflag = 1; - localflag = 0; - lo = force->numeric(FLERR,arg[iarg+1]); - hi = force->numeric(FLERR,arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg],"local") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command"); - localflag = 1; - globalflag = 0; - lo = force->numeric(FLERR,arg[iarg+1]); - hi = force->numeric(FLERR,arg[iarg+2]); - deltasq = force->numeric(FLERR,arg[iarg+3]) * - force->numeric(FLERR,arg[iarg+3]); - iarg += 4; - - } else if (strcmp(arg[iarg],"near") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - nearsq = force->numeric(FLERR,arg[iarg+1]) * - force->numeric(FLERR,arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"attempt") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - maxattempt = force->inumeric(FLERR,arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"rate") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - rateflag = 1; - rate = force->numeric(FLERR,arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"vx") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command"); - vxlo = force->numeric(FLERR,arg[iarg+1]); - vxhi = force->numeric(FLERR,arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg],"vy") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command"); - vylo = force->numeric(FLERR,arg[iarg+1]); - vyhi = force->numeric(FLERR,arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg],"vz") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command"); - vzlo = force->numeric(FLERR,arg[iarg+1]); - vzhi = force->numeric(FLERR,arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg],"units") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; - else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; - else error->all(FLERR,"Illegal fix deposit command"); - iarg += 2; - } else if (strcmp(arg[iarg],"target") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command"); - tx = force->numeric(FLERR,arg[iarg+1]); - ty = force->numeric(FLERR,arg[iarg+2]); - tz = force->numeric(FLERR,arg[iarg+3]); - targetflag = 1; - iarg += 4; - } else error->all(FLERR,"Illegal fix deposit command"); - } -} - -/* ---------------------------------------------------------------------- - pack entire state of Fix into one write -------------------------------------------------------------------------- */ - -void FixDeposit::write_restart(FILE *fp) -{ - int n = 0; - double list[4]; - list[n++] = random->state(); - list[n++] = ninserted; - list[n++] = nfirst; - list[n++] = next_reneighbor; - - if (comm->me == 0) { - int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(list,sizeof(double),n,fp); - } -} - -/* ---------------------------------------------------------------------- - use state info from restart file to restart the Fix -------------------------------------------------------------------------- */ - -void FixDeposit::restart(char *buf) -{ - int n = 0; - double *list = (double *) buf; - - seed = static_cast (list[n++]); - ninserted = static_cast (list[n++]); - nfirst = static_cast (list[n++]); - next_reneighbor = static_cast (list[n++]); - - random->reset(seed); -} diff --git a/src/fix_deposit.h b/src/fix_deposit.h deleted file mode 100644 index dfac1bcbec..0000000000 --- a/src/fix_deposit.h +++ /dev/null @@ -1,110 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#ifdef FIX_CLASS - -FixStyle(deposit,FixDeposit) - -#else - -#ifndef LMP_FIX_DEPOSIT_H -#define LMP_FIX_DEPOSIT_H - -#include "stdio.h" -#include "fix.h" - -namespace LAMMPS_NS { - -class FixDeposit : public Fix { - public: - int ntype; // type of deposited atom, visible to PairGran - - FixDeposit(class LAMMPS *, int, char **); - ~FixDeposit(); - int setmask(); - void init(); - void pre_exchange(); - void write_restart(FILE *); - void restart(char *); - - private: - int ninsert,nfreq,seed; - int iregion,globalflag,localflag,maxattempt,rateflag,scaleflag,targetflag; - int mode,rigidflag,shakeflag,idnext; - double lo,hi,deltasq,nearsq,rate; - double vxlo,vxhi,vylo,vyhi,vzlo,vzhi; - double xlo,xhi,ylo,yhi,zlo,zhi; - double tx,ty,tz; - char *idregion; - char *idrigid,*idshake; - - class Molecule *onemol; - int natom; - double **coords; - int *imageflags; - class Fix *fixrigid,*fixshake; - - int nfirst,ninserted; - int maxtag_all,maxmol_all; - class RanPark *random; - - void find_maxid(); - void options(int, char **); -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: Must specify a region in fix deposit - -The region keyword must be specified with this fix. - -E: Fix deposit region does not support a bounding box - -Not all regions represent bounded volumes. You cannot use -such a region with the fix deposit command. - -E: Fix deposit region cannot be dynamic - -Only static regions can be used with fix deposit. - -E: Deposition region extends outside simulation box - -Self-explanatory. - -E: Region ID for fix deposit does not exist - -Self-explanatory. - -W: Particle deposition was unsuccessful - -The fix deposit command was not able to insert as many atoms as -needed. The requested volume fraction may be too high, or other atoms -may be in the insertion region. - -U: Use of fix deposit with undefined lattice - -Must use lattice command with compute fix deposit command if units -option is set to lattice. - -*/ diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp deleted file mode 100644 index 301cfd815a..0000000000 --- a/src/fix_shake.cpp +++ /dev/null @@ -1,2578 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 "lmptype.h" -#include "mpi.h" -#include "math.h" -#include "stdlib.h" -#include "string.h" -#include "stdio.h" -#include "fix_shake.h" -#include "atom.h" -#include "atom_vec.h" -#include "molecule.h" -#include "update.h" -#include "respa.h" -#include "modify.h" -#include "domain.h" -#include "force.h" -#include "bond.h" -#include "angle.h" -#include "comm.h" -#include "group.h" -#include "fix_respa.h" -#include "math_const.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace FixConst; -using namespace MathConst; - -// allocate space for static class variable - -FixShake *FixShake::fsptr; - -#define BIG 1.0e20 -#define MASSDELTA 0.1 - -/* ---------------------------------------------------------------------- */ - -FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - MPI_Comm_rank(world,&me); - MPI_Comm_size(world,&nprocs); - - virial_flag = 1; - create_attribute = 1; - - // error check - - if (atom->molecular == 0) - error->all(FLERR,"Cannot use fix shake with non-molecular system"); - - // perform initial allocation of atom-based arrays - // register with Atom class - - shake_flag = NULL; - shake_atom = shake_type = NULL; - xshake = NULL; - - grow_arrays(atom->nmax); - atom->add_callback(0); - - // set comm size needed by this fix - - comm_forward = 3; - - // parse SHAKE args - - if (narg < 8) error->all(FLERR,"Illegal fix shake command"); - - tolerance = force->numeric(FLERR,arg[3]); - max_iter = force->inumeric(FLERR,arg[4]); - output_every = force->inumeric(FLERR,arg[5]); - - // parse SHAKE args for bond and angle types - // will be used by find_clusters - // store args for "b" "a" "t" as flags in (1:n) list for fast access - // store args for "m" in list of length nmass for looping over - // for "m" verify that atom masses have been set - - bond_flag = new int[atom->nbondtypes+1]; - for (int i = 1; i <= atom->nbondtypes; i++) bond_flag[i] = 0; - angle_flag = new int[atom->nangletypes+1]; - for (int i = 1; i <= atom->nangletypes; i++) angle_flag[i] = 0; - type_flag = new int[atom->ntypes+1]; - for (int i = 1; i <= atom->ntypes; i++) type_flag[i] = 0; - mass_list = new double[atom->ntypes]; - nmass = 0; - - char mode = '\0'; - int next = 6; - while (next < narg) { - if (strcmp(arg[next],"b") == 0) mode = 'b'; - else if (strcmp(arg[next],"a") == 0) mode = 'a'; - else if (strcmp(arg[next],"t") == 0) mode = 't'; - else if (strcmp(arg[next],"m") == 0) { - mode = 'm'; - atom->check_mass(); - - // break if keyword that is not b,a,t,m - - } else if (isalpha(arg[next][0])) break; - - // read numeric args of b,a,t,m - - else if (mode == 'b') { - int i = force->inumeric(FLERR,arg[next]); - if (i < 1 || i > atom->nbondtypes) - error->all(FLERR,"Invalid bond type index for fix shake"); - bond_flag[i] = 1; - - } else if (mode == 'a') { - int i = force->inumeric(FLERR,arg[next]); - if (i < 1 || i > atom->nangletypes) - error->all(FLERR,"Invalid angle type index for fix shake"); - angle_flag[i] = 1; - - } else if (mode == 't') { - int i = force->inumeric(FLERR,arg[next]); - if (i < 1 || i > atom->ntypes) - error->all(FLERR,"Invalid atom type index for fix shake"); - type_flag[i] = 1; - - } else if (mode == 'm') { - double massone = force->numeric(FLERR,arg[next]); - if (massone == 0.0) error->all(FLERR,"Invalid atom mass for fix shake"); - if (nmass == atom->ntypes) - error->all(FLERR,"Too many masses for fix shake"); - mass_list[nmass++] = massone; - - } else error->all(FLERR,"Illegal fix shake command"); - next++; - } - - // parse optional args - - onemol = NULL; - - int iarg = next; - while (iarg < narg) { - if (strcmp(arg[next],"mol") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix shake command"); - int imol = atom->find_molecule(arg[iarg+1]); - if (imol == -1) - error->all(FLERR,"Molecule ID for fix shake does not exist"); - onemol = atom->molecules[imol]; - iarg += 2; - } else error->all(FLERR,"Illegal fix shake command"); - } - - // error check for Molecule template - - if (onemol && onemol->shakeflag == 0) - error->all(FLERR,"Fix shake molecule must have shake info"); - - // allocate bond and angle distance arrays, indexed from 1 to n - - bond_distance = new double[atom->nbondtypes+1]; - angle_distance = new double[atom->nangletypes+1]; - - // allocate statistics arrays - - if (output_every) { - int nb = atom->nbondtypes + 1; - b_count = new int[nb]; - b_count_all = new int[nb]; - b_ave = new double[nb]; - b_ave_all = new double[nb]; - b_max = new double[nb]; - b_max_all = new double[nb]; - b_min = new double[nb]; - b_min_all = new double[nb]; - - int na = atom->nangletypes + 1; - a_count = new int[na]; - a_count_all = new int[na]; - a_ave = new double[na]; - a_ave_all = new double[na]; - a_max = new double[na]; - a_max_all = new double[na]; - a_min = new double[na]; - a_min_all = new double[na]; - } - - // identify all SHAKE clusters - - find_clusters(); - - // initialize list of SHAKE clusters to constrain - - maxlist = 0; - list = NULL; -} - -/* ---------------------------------------------------------------------- */ - -FixShake::~FixShake() -{ - // unregister callbacks to this fix from Atom class - - atom->delete_callback(id,0); - - // set bond_type and angle_type back to positive for SHAKE clusters - // must set for all SHAKE bonds and angles stored by each atom - - int **bond_type = atom->bond_type; - int **angle_type = atom->angle_type; - int nlocal = atom->nlocal; - - int n; - for (int i = 0; i < nlocal; i++) { - if (shake_flag[i] == 0) continue; - else if (shake_flag[i] == 1) { - n = bondfind(i,shake_atom[i][0],shake_atom[i][1]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - n = bondfind(i,shake_atom[i][0],shake_atom[i][2]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - n = anglefind(i,shake_atom[i][1],shake_atom[i][2]); - if (n >= 0) angle_type[i][n] = -angle_type[i][n]; - } else if (shake_flag[i] == 2) { - n = bondfind(i,shake_atom[i][0],shake_atom[i][1]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - } else if (shake_flag[i] == 3) { - n = bondfind(i,shake_atom[i][0],shake_atom[i][1]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - n = bondfind(i,shake_atom[i][0],shake_atom[i][2]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - } else if (shake_flag[i] == 4) { - n = bondfind(i,shake_atom[i][0],shake_atom[i][1]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - n = bondfind(i,shake_atom[i][0],shake_atom[i][2]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - n = bondfind(i,shake_atom[i][0],shake_atom[i][3]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - } - } - - // delete locally stored arrays - - memory->destroy(shake_flag); - memory->destroy(shake_atom); - memory->destroy(shake_type); - memory->destroy(xshake); - - delete [] bond_flag; - delete [] angle_flag; - delete [] type_flag; - delete [] mass_list; - - delete [] bond_distance; - delete [] angle_distance; - - if (output_every) { - delete [] b_count; - delete [] b_count_all; - delete [] b_ave; - delete [] b_ave_all; - delete [] b_max; - delete [] b_max_all; - delete [] b_min; - delete [] b_min_all; - - delete [] a_count; - delete [] a_count_all; - delete [] a_ave; - delete [] a_ave_all; - delete [] a_max; - delete [] a_max_all; - delete [] a_min; - delete [] a_min_all; - } - - memory->destroy(list); -} - -/* ---------------------------------------------------------------------- */ - -int FixShake::setmask() -{ - int mask = 0; - mask |= PRE_NEIGHBOR; - mask |= POST_FORCE; - mask |= POST_FORCE_RESPA; - return mask; -} - -/* ---------------------------------------------------------------------- - set bond and angle distances - this init must happen after force->bond and force->angle inits -------------------------------------------------------------------------- */ - -void FixShake::init() -{ - int i,m,flag,flag_all,type1,type2,bond1_type,bond2_type; - double rsq,angle; - - // error if more than one shake fix - - int count = 0; - for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"shake") == 0) count++; - if (count > 1) error->all(FLERR,"More than one fix shake"); - - // cannot use with minimization since SHAKE turns off bonds - // that should contribute to potential energy - - if (update->whichflag == 2) - error->all(FLERR,"Fix shake cannot be used with minimization"); - - // error if npt,nph fix comes before shake 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,"shake") == 0) - error->all(FLERR,"Shake fix must come before NPT/NPH fix"); - } - - // if rRESPA, find associated fix that must exist - // could have changed locations in fix list since created - // set ptrs to rRESPA variables - - if (strstr(update->integrate_style,"respa")) { - for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"RESPA") == 0) ifix_respa = i; - nlevels_respa = ((Respa *) update->integrate)->nlevels; - loop_respa = ((Respa *) update->integrate)->loop; - step_respa = ((Respa *) update->integrate)->step; - } - - // set equilibrium bond distances - - if (force->bond == NULL) - error->all(FLERR,"Bond potential must be defined for SHAKE"); - for (i = 1; i <= atom->nbondtypes; i++) - bond_distance[i] = force->bond->equilibrium_distance(i); - - // set equilibrium angle distances - - int nlocal = atom->nlocal; - - for (i = 1; i <= atom->nangletypes; i++) { - if (angle_flag[i] == 0) continue; - if (force->angle == NULL) - error->all(FLERR,"Angle potential must be defined for SHAKE"); - - // scan all atoms for a SHAKE angle cluster - // extract bond types for the 2 bonds in the cluster - // bond types must be same in all clusters of this angle type, - // else set error flag - - flag = 0; - bond1_type = bond2_type = 0; - for (m = 0; m < nlocal; m++) { - if (shake_flag[m] != 1) continue; - if (shake_type[m][2] != i) continue; - type1 = MIN(shake_type[m][0],shake_type[m][1]); - type2 = MAX(shake_type[m][0],shake_type[m][1]); - if (bond1_type > 0) { - if (type1 != bond1_type || type2 != bond2_type) { - flag = 1; - break; - } - } - bond1_type = type1; - bond2_type = type2; - } - - // error check for any bond types that are not the same - - MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_MAX,world); - if (flag_all) error->all(FLERR,"Shake angles have different bond types"); - - // insure all procs have bond types - - MPI_Allreduce(&bond1_type,&flag_all,1,MPI_INT,MPI_MAX,world); - bond1_type = flag_all; - MPI_Allreduce(&bond2_type,&flag_all,1,MPI_INT,MPI_MAX,world); - bond2_type = flag_all; - - // if bond types are 0, no SHAKE angles of this type exist - // just skip this angle - - if (bond1_type == 0) { - angle_distance[i] = 0.0; - continue; - } - - // compute the angle distance as a function of 2 bond distances - - angle = force->angle->equilibrium_angle(i); - rsq = 2.0*bond_distance[bond1_type]*bond_distance[bond2_type] * - (1.0-cos(angle)); - angle_distance[i] = sqrt(rsq); - } -} - -/* ---------------------------------------------------------------------- - SHAKE as pre-integrator constraint -------------------------------------------------------------------------- */ - -void FixShake::setup(int vflag) -{ - pre_neighbor(); - - if (output_every) stats(); - - // setup SHAKE output - - bigint ntimestep = update->ntimestep; - if (output_every) { - next_output = ntimestep + output_every; - if (ntimestep % output_every != 0) - next_output = (ntimestep/output_every)*output_every + output_every; - } else next_output = -1; - - // half timestep constraint on pre-step, full timestep thereafter - - if (strstr(update->integrate_style,"verlet")) { - dtv = update->dt; - dtfsq = 0.5 * update->dt * update->dt * force->ftm2v; - post_force(vflag); - dtfsq = update->dt * update->dt * force->ftm2v; - } else { - dtv = step_respa[0]; - dtf_innerhalf = 0.5 * step_respa[0] * force->ftm2v; - dtf_inner = dtf_innerhalf; - - // apply correction to all rRESPA levels - - for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) { - ((Respa *) update->integrate)->copy_flevel_f(ilevel); - post_force_respa(vflag,ilevel,loop_respa[ilevel]-1); - ((Respa *) update->integrate)->copy_f_flevel(ilevel); - } - - dtf_inner = step_respa[0] * force->ftm2v; - } -} - -/* ---------------------------------------------------------------------- - build list of SHAKE clusters to constrain - if one or more atoms in cluster are on this proc, - this proc lists the cluster exactly once -------------------------------------------------------------------------- */ - -void FixShake::pre_neighbor() -{ - int atom1,atom2,atom3,atom4; - - // local copies of atom quantities - // used by SHAKE until next re-neighboring - - x = atom->x; - v = atom->v; - f = atom->f; - mass = atom->mass; - rmass = atom->rmass; - type = atom->type; - nlocal = atom->nlocal; - - // extend size of SHAKE list if necessary - - if (nlocal > maxlist) { - maxlist = nlocal; - memory->destroy(list); - memory->create(list,maxlist,"shake:list"); - } - - // build list of SHAKE clusters I compute - - nlist = 0; - - for (int i = 0; i < nlocal; i++) - if (shake_flag[i]) { - if (shake_flag[i] == 2) { - atom1 = atom->map(shake_atom[i][0]); - atom2 = atom->map(shake_atom[i][1]); - if (atom1 == -1 || atom2 == -1) { - char str[128]; - sprintf(str, - "Shake atoms %d %d missing on proc %d at step " BIGINT_FORMAT, - shake_atom[i][0],shake_atom[i][1],me,update->ntimestep); - error->one(FLERR,str); - } - if (i <= atom1 && i <= atom2) list[nlist++] = i; - } else if (shake_flag[i] % 2 == 1) { - atom1 = atom->map(shake_atom[i][0]); - atom2 = atom->map(shake_atom[i][1]); - atom3 = atom->map(shake_atom[i][2]); - if (atom1 == -1 || atom2 == -1 || atom3 == -1) { - char str[128]; - sprintf(str, - "Shake atoms %d %d %d missing on proc %d at step " - BIGINT_FORMAT, - shake_atom[i][0],shake_atom[i][1],shake_atom[i][2], - me,update->ntimestep); - error->one(FLERR,str); - } - if (i <= atom1 && i <= atom2 && i <= atom3) list[nlist++] = i; - } else { - atom1 = atom->map(shake_atom[i][0]); - atom2 = atom->map(shake_atom[i][1]); - atom3 = atom->map(shake_atom[i][2]); - atom4 = atom->map(shake_atom[i][3]); - if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1) { - char str[128]; - sprintf(str, - "Shake atoms %d %d %d %d missing on proc %d at step " - BIGINT_FORMAT, - shake_atom[i][0],shake_atom[i][1], - shake_atom[i][2],shake_atom[i][3], - me,update->ntimestep); - error->one(FLERR,str); - } - if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4) - list[nlist++] = i; - } - } -} - -/* ---------------------------------------------------------------------- - compute the force adjustment for SHAKE constraint -------------------------------------------------------------------------- */ - -void FixShake::post_force(int vflag) -{ - if (update->ntimestep == next_output) stats(); - - // xshake = unconstrained move with current v,f - // communicate results if necessary - - unconstrained_update(); - if (nprocs > 1) comm->forward_comm_fix(this); - - // virial setup - - if (vflag) v_setup(vflag); - else evflag = 0; - - // loop over clusters to add constraint forces - - int m; - for (int i = 0; i < nlist; i++) { - m = list[i]; - if (shake_flag[m] == 2) shake(m); - else if (shake_flag[m] == 3) shake3(m); - else if (shake_flag[m] == 4) shake4(m); - else shake3angle(m); - } -} - -/* ---------------------------------------------------------------------- - enforce SHAKE constraints from rRESPA - xshake prediction portion is different than Verlet -------------------------------------------------------------------------- */ - -void FixShake::post_force_respa(int vflag, int ilevel, int iloop) -{ - // call stats only on outermost level - - if (ilevel == nlevels_respa-1 && update->ntimestep == next_output) stats(); - - // might be OK to skip enforcing SHAKE constraings - // on last iteration of inner levels if pressure not requested - // however, leads to slightly different trajectories - - //if (ilevel < nlevels_respa-1 && iloop == loop_respa[ilevel]-1 && !vflag) - // return; - - // xshake = unconstrained move with current v,f as function of level - // communicate results if necessary - - unconstrained_update_respa(ilevel); - if (nprocs > 1) comm->forward_comm_fix(this); - - // virial setup only needed on last iteration of innermost level - // and if pressure is requested - // virial accumulation happens via evflag at last iteration of each level - - if (ilevel == 0 && iloop == loop_respa[ilevel]-1 && vflag) v_setup(vflag); - if (iloop == loop_respa[ilevel]-1) evflag = 1; - else evflag = 0; - - // loop over clusters to add constraint forces - - int m; - for (int i = 0; i < nlist; i++) { - m = list[i]; - if (shake_flag[m] == 2) shake(m); - else if (shake_flag[m] == 3) shake3(m); - else if (shake_flag[m] == 4) shake4(m); - else shake3angle(m); - } -} - -/* ---------------------------------------------------------------------- - count # of degrees-of-freedom removed by SHAKE for atoms in igroup -------------------------------------------------------------------------- */ - -int FixShake::dof(int igroup) -{ - int groupbit = group->bitmask[igroup]; - - int *mask = atom->mask; - int *tag = atom->tag; - int nlocal = atom->nlocal; - - // count dof in a cluster if and only if - // the central atom is in group and atom i is the central atom - - int n = 0; - for (int i = 0; i < nlocal; i++) { - if (!(mask[i] & groupbit)) continue; - if (shake_flag[i] == 0) continue; - if (shake_atom[i][0] != tag[i]) continue; - if (shake_flag[i] == 1) n += 3; - else if (shake_flag[i] == 2) n += 1; - else if (shake_flag[i] == 3) n += 2; - else if (shake_flag[i] == 4) n += 3; - } - - int nall; - MPI_Allreduce(&n,&nall,1,MPI_INT,MPI_SUM,world); - return nall; -} - -/* ---------------------------------------------------------------------- - identify whether each atom is in a SHAKE cluster - only include atoms in fix group and those bonds/angles specified in input - test whether all clusters are valid - set shake_flag, shake_atom, shake_type values - set bond,angle types negative so will be ignored in neighbor lists -------------------------------------------------------------------------- */ - -void FixShake::find_clusters() -{ - int i,j,m,n; - int flag,flag_all,messtag,loop,nbuf,nbufmax,size; - double massone; - int *buf,*bufcopy; - MPI_Request request; - MPI_Status status; - - if (me == 0 && screen) fprintf(screen,"Finding SHAKE clusters ...\n"); - - // local copies of atom ptrs - - int *tag = atom->tag; - int *type = atom->type; - int *mask = atom->mask; - double *mass = atom->mass; - double *rmass = atom->rmass; - int **bond_type = atom->bond_type; - int **angle_type = atom->angle_type; - int **nspecial = atom->nspecial; - int **special = atom->special; - int nlocal = atom->nlocal; - - int angles_allow = atom->avec->angles_allow; - - // setup ring of procs - - int next = me + 1; - int prev = me -1; - if (next == nprocs) next = 0; - if (prev < 0) prev = nprocs - 1; - - // ----------------------------------------------------- - // allocate arrays for self (1d) and bond partners (2d) - // max = max # of bond partners for owned atoms = 2nd dim of partner arrays - // npartner[i] = # of bonds attached to atom i - // nshake[i] = # of SHAKE bonds attached to atom i - // partner_tag[i][] = global IDs of each partner - // partner_mask[i][] = mask of each partner - // partner_type[i][] = type of each partner - // partner_massflag[i][] = 1 if partner meets mass criterion, 0 if not - // partner_bondtype[i][] = type of bond attached to each partner - // partner_shake[i][] = 1 if SHAKE bonded to partner, 0 if not - // partner_nshake[i][] = nshake value for each partner - // ----------------------------------------------------- - - int max = 0; - for (i = 0; i < nlocal; i++) max = MAX(max,nspecial[i][0]); - - int *npartner; - memory->create(npartner,nlocal,"shake:npartner"); - memory->create(nshake,nlocal,"shake:nshake"); - - int **partner_tag,**partner_mask,**partner_type,**partner_massflag; - int ** partner_bondtype,**partner_shake,**partner_nshake; - memory->create(partner_tag,nlocal,max,"shake:partner_tag"); - memory->create(partner_mask,nlocal,max,"shake:partner_mask"); - memory->create(partner_type,nlocal,max,"shake:partner_type"); - memory->create(partner_massflag,nlocal,max,"shake:partner_massflag"); - memory->create(partner_bondtype,nlocal,max,"shake:partner_bondtype"); - memory->create(partner_shake,nlocal,max,"shake:partner_shake"); - memory->create(partner_nshake,nlocal,max,"shake:partner_nshake"); - - // ----------------------------------------------------- - // set npartner and partner_tag from special arrays - // ----------------------------------------------------- - - for (i = 0; i < nlocal; i++) { - npartner[i] = nspecial[i][0]; - for (j = 0; j < npartner[i]; j++) partner_tag[i][j] = special[i][j]; - } - - // ----------------------------------------------------- - // set partner_mask, partner_type, partner_massflag, partner_bondtype - // for bonded partners - // requires communication for off-proc partners - // ----------------------------------------------------- - - // fill in mask, type, massflag, bondtype if own bond partner - // info to store in buf for each off-proc bond = nper = 6 - // 2 atoms IDs in bond, space for mask, type, massflag, bondtype - // nbufmax = largest buffer needed to hold info from any proc - - int nper = 6; - - nbuf = 0; - for (i = 0; i < nlocal; i++) { - for (j = 0; j < npartner[i]; j++) { - partner_mask[i][j] = 0; - partner_type[i][j] = 0; - partner_massflag[i][j] = 0; - partner_bondtype[i][j] = 0; - - m = atom->map(partner_tag[i][j]); - if (m >= 0 && m < nlocal) { - partner_mask[i][j] = mask[m]; - partner_type[i][j] = type[m]; - if (nmass) { - if (rmass) massone = rmass[m]; - else massone = mass[type[m]]; - partner_massflag[i][j] = masscheck(massone); - } - n = bondfind(i,tag[i],partner_tag[i][j]); - if (n >= 0) partner_bondtype[i][j] = bond_type[i][n]; - else { - n = bondfind(m,tag[i],partner_tag[i][j]); - if (n >= 0) partner_bondtype[i][j] = bond_type[m][n]; - } - } else nbuf += nper; - } - } - - memory->create(buf,nbuf,"shake:buf"); - - // fill buffer with info - - size = 0; - for (i = 0; i < nlocal; i++) { - for (j = 0; j < npartner[i]; j++) { - m = atom->map(partner_tag[i][j]); - if (m < 0 || m >= nlocal) { - buf[size] = tag[i]; - buf[size+1] = partner_tag[i][j]; - buf[size+2] = 0; - buf[size+3] = 0; - buf[size+4] = 0; - n = bondfind(i,tag[i],partner_tag[i][j]); - if (n >= 0) buf[size+5] = bond_type[i][n]; - else buf[size+5] = 0; - size += nper; - } - } - } - - // cycle buffer around ring of procs back to self - - fsptr = this; - comm->ring(size,sizeof(int),buf,1,ring_bonds,buf); - - // store partner info returned to me - - m = 0; - while (m < size) { - i = atom->map(buf[m]); - for (j = 0; j < npartner[i]; j++) - if (buf[m+1] == partner_tag[i][j]) break; - partner_mask[i][j] = buf[m+2]; - partner_type[i][j] = buf[m+3]; - partner_massflag[i][j] = buf[m+4]; - partner_bondtype[i][j] = buf[m+5]; - m += nper; - } - - memory->destroy(buf); - - // error check for unfilled partner info - // if partner_type not set, is an error - // partner_bondtype may not be set if special list is not consistent - // with bondatom (e.g. due to delete_bonds command) - // this is OK if one or both atoms are not in fix group, since - // bond won't be SHAKEn anyway - // else it's an error - - flag = 0; - for (i = 0; i < nlocal; i++) - for (j = 0; j < npartner[i]; j++) { - if (partner_type[i][j] == 0) flag = 1; - if (!(mask[i] & groupbit)) continue; - if (!(partner_mask[i][j] & groupbit)) continue; - if (partner_bondtype[i][j] == 0) flag = 1; - } - - MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); - if (flag_all) error->all(FLERR,"Did not find fix shake partner info"); - - // ----------------------------------------------------- - // identify SHAKEable bonds - // set nshake[i] = # of SHAKE bonds attached to atom i - // set partner_shake[i][] = 1 if SHAKE bonded to partner, 0 if not - // both atoms must be in group, bondtype must be > 0 - // check if bondtype is in input bond_flag - // check if type of either atom is in input type_flag - // check if mass of either atom is in input mass_list - // ----------------------------------------------------- - - int np; - - for (i = 0; i < nlocal; i++) { - nshake[i] = 0; - np = npartner[i]; - for (j = 0; j < np; j++) { - partner_shake[i][j] = 0; - - if (!(mask[i] & groupbit)) continue; - if (!(partner_mask[i][j] & groupbit)) continue; - if (partner_bondtype[i][j] <= 0) continue; - - if (bond_flag[partner_bondtype[i][j]]) { - partner_shake[i][j] = 1; - nshake[i]++; - continue; - } - if (type_flag[type[i]] || type_flag[partner_type[i][j]]) { - partner_shake[i][j] = 1; - nshake[i]++; - continue; - } - if (nmass) { - if (partner_massflag[i][j]) { - partner_shake[i][j] = 1; - nshake[i]++; - continue; - } else { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - if (masscheck(massone)) { - partner_shake[i][j] = 1; - nshake[i]++; - continue; - } - } - } - } - } - - // ----------------------------------------------------- - // set partner_nshake for bonded partners - // requires communication for off-proc partners - // ----------------------------------------------------- - - // fill in partner_nshake if own bond partner - // info to store in buf for each off-proc bond = - // 2 atoms IDs in bond, space for nshake value - // nbufmax = largest buffer needed to hold info from any proc - - nbuf = 0; - for (i = 0; i < nlocal; i++) { - for (j = 0; j < npartner[i]; j++) { - m = atom->map(partner_tag[i][j]); - if (m >= 0 && m < nlocal) partner_nshake[i][j] = nshake[m]; - else nbuf += 3; - } - } - - memory->create(buf,nbuf,"shake:buf"); - - // fill buffer with info - - size = 0; - for (i = 0; i < nlocal; i++) { - for (j = 0; j < npartner[i]; j++) { - m = atom->map(partner_tag[i][j]); - if (m < 0 || m >= nlocal) { - buf[size] = tag[i]; - buf[size+1] = partner_tag[i][j]; - size += 3; - } - } - } - - // cycle buffer around ring of procs back to self - - fsptr = this; - comm->ring(size,sizeof(int),buf,2,ring_nshake,buf); - - // store partner info returned to me - - m = 0; - while (m < size) { - i = atom->map(buf[m]); - for (j = 0; j < npartner[i]; j++) - if (buf[m+1] == partner_tag[i][j]) break; - partner_nshake[i][j] = buf[m+2]; - m += 3; - } - - memory->destroy(buf); - - // ----------------------------------------------------- - // error checks - // no atom with nshake > 3 - // no connected atoms which both have nshake > 1 - // ----------------------------------------------------- - - flag = 0; - for (i = 0; i < nlocal; i++) if (nshake[i] > 3) flag = 1; - MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); - if (flag_all) error->all(FLERR,"Shake cluster of more than 4 atoms"); - - flag = 0; - for (i = 0; i < nlocal; i++) { - if (nshake[i] <= 1) continue; - for (j = 0; j < npartner[i]; j++) - if (partner_shake[i][j] && partner_nshake[i][j] > 1) flag = 1; - } - MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); - if (flag_all) error->all(FLERR,"Shake clusters are connected"); - - // ----------------------------------------------------- - // set SHAKE arrays that are stored with atoms & add angle constraints - // zero shake arrays for all owned atoms - // if I am central atom set shake_flag & shake_atom & shake_type - // for 2-atom clusters, I am central atom if my atom ID < partner ID - // for 3-atom clusters, test for angle constraint - // angle will be stored by this atom if it exists - // if angle type matches angle_flag, then it is angle-constrained - // shake_flag[] = 0 if atom not in SHAKE cluster - // 2,3,4 = size of bond-only cluster - // 1 = 3-atom angle cluster - // shake_atom[][] = global IDs of 2,3,4 atoms in cluster - // central atom is 1st - // for 2-atom cluster, lowest ID is 1st - // shake_type[][] = bondtype of each bond in cluster - // for 3-atom angle cluster, 3rd value is angletype - // ----------------------------------------------------- - - for (i = 0; i < nlocal; i++) { - shake_flag[i] = 0; - shake_atom[i][0] = 0; - shake_atom[i][1] = 0; - shake_atom[i][2] = 0; - shake_atom[i][3] = 0; - shake_type[i][0] = 0; - shake_type[i][1] = 0; - shake_type[i][2] = 0; - - if (nshake[i] == 1) { - for (j = 0; j < npartner[i]; j++) - if (partner_shake[i][j]) break; - if (partner_nshake[i][j] == 1 && tag[i] < partner_tag[i][j]) { - shake_flag[i] = 2; - shake_atom[i][0] = tag[i]; - shake_atom[i][1] = partner_tag[i][j]; - shake_type[i][0] = partner_bondtype[i][j]; - } - } - - if (nshake[i] > 1) { - shake_flag[i] = 1; - shake_atom[i][0] = tag[i]; - for (j = 0; j < npartner[i]; j++) - if (partner_shake[i][j]) { - m = shake_flag[i]; - shake_atom[i][m] = partner_tag[i][j]; - shake_type[i][m-1] = partner_bondtype[i][j]; - shake_flag[i]++; - } - } - - if (nshake[i] == 2 && angles_allow) { - n = anglefind(i,shake_atom[i][1],shake_atom[i][2]); - if (n < 0) continue; - if (angle_type[i][n] < 0) continue; - if (angle_flag[angle_type[i][n]]) { - shake_flag[i] = 1; - shake_type[i][2] = angle_type[i][n]; - } - } - } - - // ----------------------------------------------------- - // set shake_flag,shake_atom,shake_type for non-central atoms - // requires communication for off-proc atoms - // ----------------------------------------------------- - - // fill in shake arrays for each bond partner I own - // info to store in buf for each off-proc bond = - // all values from shake_flag, shake_atom, shake_type - // nbufmax = largest buffer needed to hold info from any proc - - nbuf = 0; - for (i = 0; i < nlocal; i++) { - if (shake_flag[i] == 0) continue; - for (j = 0; j < npartner[i]; j++) { - if (partner_shake[i][j] == 0) continue; - m = atom->map(partner_tag[i][j]); - if (m >= 0 && m < nlocal) { - shake_flag[m] = shake_flag[i]; - shake_atom[m][0] = shake_atom[i][0]; - shake_atom[m][1] = shake_atom[i][1]; - shake_atom[m][2] = shake_atom[i][2]; - shake_atom[m][3] = shake_atom[i][3]; - shake_type[m][0] = shake_type[i][0]; - shake_type[m][1] = shake_type[i][1]; - shake_type[m][2] = shake_type[i][2]; - } else nbuf += 9; - } - } - - memory->create(buf,nbuf,"shake:buf"); - - // fill buffer with info - - size = 0; - for (i = 0; i < nlocal; i++) { - if (shake_flag[i] == 0) continue; - for (j = 0; j < npartner[i]; j++) { - if (partner_shake[i][j] == 0) continue; - m = atom->map(partner_tag[i][j]); - if (m < 0 || m >= nlocal) { - buf[size] = partner_tag[i][j]; - buf[size+1] = shake_flag[i]; - buf[size+2] = shake_atom[i][0]; - buf[size+3] = shake_atom[i][1]; - buf[size+4] = shake_atom[i][2]; - buf[size+5] = shake_atom[i][3]; - buf[size+6] = shake_type[i][0]; - buf[size+7] = shake_type[i][1]; - buf[size+8] = shake_type[i][2]; - size += 9; - } - } - } - - // cycle buffer around ring of procs back to self - - fsptr = this; - comm->ring(size,sizeof(int),buf,3,ring_shake,NULL); - - memory->destroy(buf); - - // ----------------------------------------------------- - // free local memory - // ----------------------------------------------------- - - memory->destroy(npartner); - memory->destroy(nshake); - memory->destroy(partner_tag); - memory->destroy(partner_mask); - memory->destroy(partner_type); - memory->destroy(partner_massflag); - memory->destroy(partner_bondtype); - memory->destroy(partner_shake); - memory->destroy(partner_nshake); - - // ----------------------------------------------------- - // set bond_type and angle_type negative for SHAKE clusters - // must set for all SHAKE bonds and angles stored by each atom - // ----------------------------------------------------- - - for (i = 0; i < nlocal; i++) { - if (shake_flag[i] == 0) continue; - else if (shake_flag[i] == 1) { - n = bondfind(i,shake_atom[i][0],shake_atom[i][1]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - n = bondfind(i,shake_atom[i][0],shake_atom[i][2]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - n = anglefind(i,shake_atom[i][1],shake_atom[i][2]); - if (n >= 0) angle_type[i][n] = -angle_type[i][n]; - } else if (shake_flag[i] == 2) { - n = bondfind(i,shake_atom[i][0],shake_atom[i][1]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - } else if (shake_flag[i] == 3) { - n = bondfind(i,shake_atom[i][0],shake_atom[i][1]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - n = bondfind(i,shake_atom[i][0],shake_atom[i][2]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - } else if (shake_flag[i] == 4) { - n = bondfind(i,shake_atom[i][0],shake_atom[i][1]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - n = bondfind(i,shake_atom[i][0],shake_atom[i][2]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - n = bondfind(i,shake_atom[i][0],shake_atom[i][3]); - if (n >= 0) bond_type[i][n] = -bond_type[i][n]; - } - } - - // ----------------------------------------------------- - // print info on SHAKE clusters - // ----------------------------------------------------- - - int count1,count2,count3,count4; - count1 = count2 = count3 = count4 = 0; - for (i = 0; i < nlocal; i++) { - if (shake_flag[i] == 1) count1++; - else if (shake_flag[i] == 2) count2++; - else if (shake_flag[i] == 3) count3++; - else if (shake_flag[i] == 4) count4++; - } - - int tmp; - tmp = count1; - MPI_Allreduce(&tmp,&count1,1,MPI_INT,MPI_SUM,world); - tmp = count2; - MPI_Allreduce(&tmp,&count2,1,MPI_INT,MPI_SUM,world); - tmp = count3; - MPI_Allreduce(&tmp,&count3,1,MPI_INT,MPI_SUM,world); - tmp = count4; - MPI_Allreduce(&tmp,&count4,1,MPI_INT,MPI_SUM,world); - - if (me == 0) { - if (screen) { - fprintf(screen," %d = # of size 2 clusters\n",count2/2); - fprintf(screen," %d = # of size 3 clusters\n",count3/3); - fprintf(screen," %d = # of size 4 clusters\n",count4/4); - fprintf(screen," %d = # of frozen angles\n",count1/3); - } - if (logfile) { - fprintf(logfile," %d = # of size 2 clusters\n",count2/2); - fprintf(logfile," %d = # of size 3 clusters\n",count3/3); - fprintf(logfile," %d = # of size 4 clusters\n",count4/4); - fprintf(logfile," %d = # of frozen angles\n",count1/3); - } - } -} - -/* ---------------------------------------------------------------------- - when receive buffer, scan bond partner IDs for atoms I own - if I own partner: - fill in mask and type and massflag - search for bond with 1st atom and fill in bondtype -------------------------------------------------------------------------- */ - -void FixShake::ring_bonds(int ndatum, char *cbuf) -{ - Atom *atom = fsptr->atom; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *mask = atom->mask; - int **bond_type = atom->bond_type; - int *type = atom->type; - int nlocal = atom->nlocal; - int nmass = fsptr->nmass; - - int *buf = (int *) cbuf; - int m,n; - double massone; - - for (int i = 0; i < ndatum; i += 6) { - m = atom->map(buf[i+1]); - if (m >= 0 && m < nlocal) { - buf[i+2] = mask[m]; - buf[i+3] = type[m]; - if (nmass) { - if (rmass) massone = rmass[m]; - else massone = mass[type[m]]; - buf[i+4] = fsptr->masscheck(massone); - } - if (buf[i+5] == 0) { - n = fsptr->bondfind(m,buf[i],buf[i+1]); - if (n >= 0) buf[i+5] = bond_type[m][n]; - } - } - } -} - -/* ---------------------------------------------------------------------- - when receive buffer, scan bond partner IDs for atoms I own - if I own partner, fill in nshake value -------------------------------------------------------------------------- */ - -void FixShake::ring_nshake(int ndatum, char *cbuf) -{ - Atom *atom = fsptr->atom; - int nlocal = atom->nlocal; - - int *nshake = fsptr->nshake; - - int *buf = (int *) cbuf; - int m; - - for (int i = 0; i < ndatum; i += 3) { - m = atom->map(buf[i+1]); - if (m >= 0 && m < nlocal) buf[i+2] = nshake[m]; - } -} - -/* ---------------------------------------------------------------------- - when receive buffer, scan bond partner IDs for atoms I own - if I own partner, fill in nshake value -------------------------------------------------------------------------- */ - -void FixShake::ring_shake(int ndatum, char *cbuf) -{ - Atom *atom = fsptr->atom; - int nlocal = atom->nlocal; - - int *shake_flag = fsptr->shake_flag; - int **shake_atom = fsptr->shake_atom; - int **shake_type = fsptr->shake_type; - - int *buf = (int *) cbuf; - int m; - - for (int i = 0; i < ndatum; i += 9) { - m = atom->map(buf[i]); - if (m >= 0 && m < nlocal) { - shake_flag[m] = buf[i+1]; - shake_atom[m][0] = buf[i+2]; - shake_atom[m][1] = buf[i+3]; - shake_atom[m][2] = buf[i+4]; - shake_atom[m][3] = buf[i+5]; - shake_type[m][0] = buf[i+6]; - shake_type[m][1] = buf[i+7]; - shake_type[m][2] = buf[i+8]; - } - } -} - -/* ---------------------------------------------------------------------- - check if massone is within MASSDELTA of any mass in mass_list - return 1 if yes, 0 if not -------------------------------------------------------------------------- */ - -int FixShake::masscheck(double massone) -{ - for (int i = 0; i < nmass; i++) - if (fabs(mass_list[i]-massone) <= MASSDELTA) return 1; - return 0; -} - -/* ---------------------------------------------------------------------- - update the unconstrained position of each atom - only for SHAKE clusters, else set to 0.0 - assumes NVE update, seems to be accurate enough for NVT,NPT,NPH as well -------------------------------------------------------------------------- */ - -void FixShake::unconstrained_update() -{ - double dtfmsq; - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (shake_flag[i]) { - dtfmsq = dtfsq / rmass[i]; - xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; - xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; - xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; - } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; - } - } else { - for (int i = 0; i < nlocal; i++) { - if (shake_flag[i]) { - dtfmsq = dtfsq / mass[type[i]]; - xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; - xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; - xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; - } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; - } - } -} - -/* ---------------------------------------------------------------------- - update the unconstrained position of each atom in a rRESPA step - only for SHAKE clusters, else set to 0.0 - assumes NVE update, seems to be accurate enough for NVT,NPT,NPH as well -------------------------------------------------------------------------- */ - -void FixShake::unconstrained_update_respa(int ilevel) -{ - // xshake = atom coords after next x update in innermost loop - // depends on rRESPA level - // for levels > 0 this includes more than one velocity update - // xshake = predicted position from call to this routine at level N = - // x + dt0 (v + dtN/m fN + 1/2 dt(N-1)/m f(N-1) + ... + 1/2 dt0/m f0) - // also set dtfsq = dt0*dtN so that shake,shake3,etc can use it - - double ***f_level = ((FixRespa *) modify->fix[ifix_respa])->f_level; - dtfsq = dtf_inner * step_respa[ilevel]; - - double invmass,dtfmsq; - int jlevel; - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (shake_flag[i]) { - invmass = 1.0 / rmass[i]; - dtfmsq = dtfsq * invmass; - xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; - xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; - xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; - for (jlevel = 0; jlevel < ilevel; jlevel++) { - dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; - xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; - xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; - xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; - } - } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; - } - - } else { - for (int i = 0; i < nlocal; i++) { - if (shake_flag[i]) { - invmass = 1.0 / mass[type[i]]; - dtfmsq = dtfsq * invmass; - xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; - xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; - xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; - for (jlevel = 0; jlevel < ilevel; jlevel++) { - dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; - xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; - xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; - xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; - } - } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixShake::shake(int m) -{ - int nlist,list[2]; - double v[6]; - double invmass0,invmass1; - - // local atom IDs and constraint distances - - int i0 = atom->map(shake_atom[m][0]); - int i1 = atom->map(shake_atom[m][1]); - double bond1 = bond_distance[shake_type[m][0]]; - - // r01 = distance vec between atoms, with PBC - - double r01[3]; - r01[0] = x[i0][0] - x[i1][0]; - r01[1] = x[i0][1] - x[i1][1]; - r01[2] = x[i0][2] - x[i1][2]; - domain->minimum_image(r01); - - // s01 = distance vec after unconstrained update, with PBC - - double s01[3]; - s01[0] = xshake[i0][0] - xshake[i1][0]; - s01[1] = xshake[i0][1] - xshake[i1][1]; - s01[2] = xshake[i0][2] - xshake[i1][2]; - domain->minimum_image(s01); - - // scalar distances between atoms - - double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2]; - double s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2]; - - // a,b,c = coeffs in quadratic equation for lamda - - if (rmass) { - invmass0 = 1.0/rmass[i0]; - invmass1 = 1.0/rmass[i1]; - } else { - invmass0 = 1.0/mass[type[i0]]; - invmass1 = 1.0/mass[type[i1]]; - } - - double a = (invmass0+invmass1)*(invmass0+invmass1) * r01sq; - double b = 2.0 * (invmass0+invmass1) * - (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); - double c = s01sq - bond1*bond1; - - // error check - - double determ = b*b - 4.0*a*c; - if (determ < 0.0) { - error->warning(FLERR,"Shake determinant < 0.0",0); - determ = 0.0; - } - - // exact quadratic solution for lamda - - double lamda,lamda1,lamda2; - lamda1 = (-b+sqrt(determ)) / (2.0*a); - lamda2 = (-b-sqrt(determ)) / (2.0*a); - - if (fabs(lamda1) <= fabs(lamda2)) lamda = lamda1; - else lamda = lamda2; - - // update forces if atom is owned by this processor - - lamda /= dtfsq; - - if (i0 < nlocal) { - f[i0][0] += lamda*r01[0]; - f[i0][1] += lamda*r01[1]; - f[i0][2] += lamda*r01[2]; - } - - if (i1 < nlocal) { - f[i1][0] -= lamda*r01[0]; - f[i1][1] -= lamda*r01[1]; - f[i1][2] -= lamda*r01[2]; - } - - if (evflag) { - nlist = 0; - if (i0 < nlocal) list[nlist++] = i0; - if (i1 < nlocal) list[nlist++] = i1; - - v[0] = lamda*r01[0]*r01[0]; - v[1] = lamda*r01[1]*r01[1]; - v[2] = lamda*r01[2]*r01[2]; - v[3] = lamda*r01[0]*r01[1]; - v[4] = lamda*r01[0]*r01[2]; - v[5] = lamda*r01[1]*r01[2]; - - v_tally(nlist,list,2.0,v); - } -} - -/* ---------------------------------------------------------------------- */ - -void FixShake::shake3(int m) -{ - int nlist,list[3]; - double v[6]; - double invmass0,invmass1,invmass2; - - // local atom IDs and constraint distances - - int i0 = atom->map(shake_atom[m][0]); - int i1 = atom->map(shake_atom[m][1]); - int i2 = atom->map(shake_atom[m][2]); - double bond1 = bond_distance[shake_type[m][0]]; - double bond2 = bond_distance[shake_type[m][1]]; - - // r01,r02 = distance vec between atoms, with PBC - - double r01[3]; - r01[0] = x[i0][0] - x[i1][0]; - r01[1] = x[i0][1] - x[i1][1]; - r01[2] = x[i0][2] - x[i1][2]; - domain->minimum_image(r01); - - double r02[3]; - r02[0] = x[i0][0] - x[i2][0]; - r02[1] = x[i0][1] - x[i2][1]; - r02[2] = x[i0][2] - x[i2][2]; - domain->minimum_image(r02); - - // s01,s02 = distance vec after unconstrained update, with PBC - - double s01[3]; - s01[0] = xshake[i0][0] - xshake[i1][0]; - s01[1] = xshake[i0][1] - xshake[i1][1]; - s01[2] = xshake[i0][2] - xshake[i1][2]; - domain->minimum_image(s01); - - double s02[3]; - s02[0] = xshake[i0][0] - xshake[i2][0]; - s02[1] = xshake[i0][1] - xshake[i2][1]; - s02[2] = xshake[i0][2] - xshake[i2][2]; - domain->minimum_image(s02); - - // scalar distances between atoms - - double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2]; - double r02sq = r02[0]*r02[0] + r02[1]*r02[1] + r02[2]*r02[2]; - double s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2]; - double s02sq = s02[0]*s02[0] + s02[1]*s02[1] + s02[2]*s02[2]; - - // matrix coeffs and rhs for lamda equations - - if (rmass) { - invmass0 = 1.0/rmass[i0]; - invmass1 = 1.0/rmass[i1]; - invmass2 = 1.0/rmass[i2]; - } else { - invmass0 = 1.0/mass[type[i0]]; - invmass1 = 1.0/mass[type[i1]]; - invmass2 = 1.0/mass[type[i2]]; - } - - double a11 = 2.0 * (invmass0+invmass1) * - (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); - double a12 = 2.0 * invmass0 * - (s01[0]*r02[0] + s01[1]*r02[1] + s01[2]*r02[2]); - double a21 = 2.0 * invmass0 * - (s02[0]*r01[0] + s02[1]*r01[1] + s02[2]*r01[2]); - double a22 = 2.0 * (invmass0+invmass2) * - (s02[0]*r02[0] + s02[1]*r02[1] + s02[2]*r02[2]); - - // inverse of matrix - - double determ = a11*a22 - a12*a21; - if (determ == 0.0) error->one(FLERR,"Shake determinant = 0.0"); - double determinv = 1.0/determ; - - double a11inv = a22*determinv; - double a12inv = -a12*determinv; - double a21inv = -a21*determinv; - double a22inv = a11*determinv; - - // quadratic correction coeffs - - double r0102 = (r01[0]*r02[0] + r01[1]*r02[1] + r01[2]*r02[2]); - - double quad1_0101 = (invmass0+invmass1)*(invmass0+invmass1) * r01sq; - double quad1_0202 = invmass0*invmass0 * r02sq; - double quad1_0102 = 2.0 * (invmass0+invmass1)*invmass0 * r0102; - - double quad2_0202 = (invmass0+invmass2)*(invmass0+invmass2) * r02sq; - double quad2_0101 = invmass0*invmass0 * r01sq; - double quad2_0102 = 2.0 * (invmass0+invmass2)*invmass0 * r0102; - - // iterate until converged - - double lamda01 = 0.0; - double lamda02 = 0.0; - int niter = 0; - int done = 0; - - double quad1,quad2,b1,b2,lamda01_new,lamda02_new; - - while (!done && niter < max_iter) { - quad1 = quad1_0101 * lamda01*lamda01 + quad1_0202 * lamda02*lamda02 + - quad1_0102 * lamda01*lamda02; - quad2 = quad2_0101 * lamda01*lamda01 + quad2_0202 * lamda02*lamda02 + - quad2_0102 * lamda01*lamda02; - - b1 = bond1*bond1 - s01sq - quad1; - b2 = bond2*bond2 - s02sq - quad2; - - lamda01_new = a11inv*b1 + a12inv*b2; - lamda02_new = a21inv*b1 + a22inv*b2; - - done = 1; - if (fabs(lamda01_new-lamda01) > tolerance) done = 0; - if (fabs(lamda02_new-lamda02) > tolerance) done = 0; - - lamda01 = lamda01_new; - lamda02 = lamda02_new; - niter++; - } - - // update forces if atom is owned by this processor - - lamda01 = lamda01/dtfsq; - lamda02 = lamda02/dtfsq; - - if (i0 < nlocal) { - f[i0][0] += lamda01*r01[0] + lamda02*r02[0]; - f[i0][1] += lamda01*r01[1] + lamda02*r02[1]; - f[i0][2] += lamda01*r01[2] + lamda02*r02[2]; - } - - if (i1 < nlocal) { - f[i1][0] -= lamda01*r01[0]; - f[i1][1] -= lamda01*r01[1]; - f[i1][2] -= lamda01*r01[2]; - } - - if (i2 < nlocal) { - f[i2][0] -= lamda02*r02[0]; - f[i2][1] -= lamda02*r02[1]; - f[i2][2] -= lamda02*r02[2]; - } - - if (evflag) { - nlist = 0; - if (i0 < nlocal) list[nlist++] = i0; - if (i1 < nlocal) list[nlist++] = i1; - if (i2 < nlocal) list[nlist++] = i2; - - v[0] = lamda01*r01[0]*r01[0] + lamda02*r02[0]*r02[0]; - v[1] = lamda01*r01[1]*r01[1] + lamda02*r02[1]*r02[1]; - v[2] = lamda01*r01[2]*r01[2] + lamda02*r02[2]*r02[2]; - v[3] = lamda01*r01[0]*r01[1] + lamda02*r02[0]*r02[1]; - v[4] = lamda01*r01[0]*r01[2] + lamda02*r02[0]*r02[2]; - v[5] = lamda01*r01[1]*r01[2] + lamda02*r02[1]*r02[2]; - - v_tally(nlist,list,3.0,v); - } -} - -/* ---------------------------------------------------------------------- */ - -void FixShake::shake4(int m) -{ - int nlist,list[4]; - double v[6]; - double invmass0,invmass1,invmass2,invmass3; - - // local atom IDs and constraint distances - - int i0 = atom->map(shake_atom[m][0]); - int i1 = atom->map(shake_atom[m][1]); - int i2 = atom->map(shake_atom[m][2]); - int i3 = atom->map(shake_atom[m][3]); - double bond1 = bond_distance[shake_type[m][0]]; - double bond2 = bond_distance[shake_type[m][1]]; - double bond3 = bond_distance[shake_type[m][2]]; - - // r01,r02,r03 = distance vec between atoms, with PBC - - double r01[3]; - r01[0] = x[i0][0] - x[i1][0]; - r01[1] = x[i0][1] - x[i1][1]; - r01[2] = x[i0][2] - x[i1][2]; - domain->minimum_image(r01); - - double r02[3]; - r02[0] = x[i0][0] - x[i2][0]; - r02[1] = x[i0][1] - x[i2][1]; - r02[2] = x[i0][2] - x[i2][2]; - domain->minimum_image(r02); - - double r03[3]; - r03[0] = x[i0][0] - x[i3][0]; - r03[1] = x[i0][1] - x[i3][1]; - r03[2] = x[i0][2] - x[i3][2]; - domain->minimum_image(r03); - - // s01,s02,s03 = distance vec after unconstrained update, with PBC - - double s01[3]; - s01[0] = xshake[i0][0] - xshake[i1][0]; - s01[1] = xshake[i0][1] - xshake[i1][1]; - s01[2] = xshake[i0][2] - xshake[i1][2]; - domain->minimum_image(s01); - - double s02[3]; - s02[0] = xshake[i0][0] - xshake[i2][0]; - s02[1] = xshake[i0][1] - xshake[i2][1]; - s02[2] = xshake[i0][2] - xshake[i2][2]; - domain->minimum_image(s02); - - double s03[3]; - s03[0] = xshake[i0][0] - xshake[i3][0]; - s03[1] = xshake[i0][1] - xshake[i3][1]; - s03[2] = xshake[i0][2] - xshake[i3][2]; - domain->minimum_image(s03); - - // scalar distances between atoms - - double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2]; - double r02sq = r02[0]*r02[0] + r02[1]*r02[1] + r02[2]*r02[2]; - double r03sq = r03[0]*r03[0] + r03[1]*r03[1] + r03[2]*r03[2]; - double s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2]; - double s02sq = s02[0]*s02[0] + s02[1]*s02[1] + s02[2]*s02[2]; - double s03sq = s03[0]*s03[0] + s03[1]*s03[1] + s03[2]*s03[2]; - - // matrix coeffs and rhs for lamda equations - - if (rmass) { - invmass0 = 1.0/rmass[i0]; - invmass1 = 1.0/rmass[i1]; - invmass2 = 1.0/rmass[i2]; - invmass3 = 1.0/rmass[i3]; - } else { - invmass0 = 1.0/mass[type[i0]]; - invmass1 = 1.0/mass[type[i1]]; - invmass2 = 1.0/mass[type[i2]]; - invmass3 = 1.0/mass[type[i3]]; - } - - double a11 = 2.0 * (invmass0+invmass1) * - (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); - double a12 = 2.0 * invmass0 * - (s01[0]*r02[0] + s01[1]*r02[1] + s01[2]*r02[2]); - double a13 = 2.0 * invmass0 * - (s01[0]*r03[0] + s01[1]*r03[1] + s01[2]*r03[2]); - double a21 = 2.0 * invmass0 * - (s02[0]*r01[0] + s02[1]*r01[1] + s02[2]*r01[2]); - double a22 = 2.0 * (invmass0+invmass2) * - (s02[0]*r02[0] + s02[1]*r02[1] + s02[2]*r02[2]); - double a23 = 2.0 * invmass0 * - (s02[0]*r03[0] + s02[1]*r03[1] + s02[2]*r03[2]); - double a31 = 2.0 * invmass0 * - (s03[0]*r01[0] + s03[1]*r01[1] + s03[2]*r01[2]); - double a32 = 2.0 * invmass0 * - (s03[0]*r02[0] + s03[1]*r02[1] + s03[2]*r02[2]); - double a33 = 2.0 * (invmass0+invmass3) * - (s03[0]*r03[0] + s03[1]*r03[1] + s03[2]*r03[2]); - - // inverse of matrix; - - double determ = a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - - a11*a23*a32 - a12*a21*a33 - a13*a22*a31; - if (determ == 0.0) error->one(FLERR,"Shake determinant = 0.0"); - double determinv = 1.0/determ; - - double a11inv = determinv * (a22*a33 - a23*a32); - double a12inv = -determinv * (a12*a33 - a13*a32); - double a13inv = determinv * (a12*a23 - a13*a22); - double a21inv = -determinv * (a21*a33 - a23*a31); - double a22inv = determinv * (a11*a33 - a13*a31); - double a23inv = -determinv * (a11*a23 - a13*a21); - double a31inv = determinv * (a21*a32 - a22*a31); - double a32inv = -determinv * (a11*a32 - a12*a31); - double a33inv = determinv * (a11*a22 - a12*a21); - - // quadratic correction coeffs - - double r0102 = (r01[0]*r02[0] + r01[1]*r02[1] + r01[2]*r02[2]); - double r0103 = (r01[0]*r03[0] + r01[1]*r03[1] + r01[2]*r03[2]); - double r0203 = (r02[0]*r03[0] + r02[1]*r03[1] + r02[2]*r03[2]); - - double quad1_0101 = (invmass0+invmass1)*(invmass0+invmass1) * r01sq; - double quad1_0202 = invmass0*invmass0 * r02sq; - double quad1_0303 = invmass0*invmass0 * r03sq; - double quad1_0102 = 2.0 * (invmass0+invmass1)*invmass0 * r0102; - double quad1_0103 = 2.0 * (invmass0+invmass1)*invmass0 * r0103; - double quad1_0203 = 2.0 * invmass0*invmass0 * r0203; - - double quad2_0101 = invmass0*invmass0 * r01sq; - double quad2_0202 = (invmass0+invmass2)*(invmass0+invmass2) * r02sq; - double quad2_0303 = invmass0*invmass0 * r03sq; - double quad2_0102 = 2.0 * (invmass0+invmass2)*invmass0 * r0102; - double quad2_0103 = 2.0 * invmass0*invmass0 * r0103; - double quad2_0203 = 2.0 * (invmass0+invmass2)*invmass0 * r0203; - - double quad3_0101 = invmass0*invmass0 * r01sq; - double quad3_0202 = invmass0*invmass0 * r02sq; - double quad3_0303 = (invmass0+invmass3)*(invmass0+invmass3) * r03sq; - double quad3_0102 = 2.0 * invmass0*invmass0 * r0102; - double quad3_0103 = 2.0 * (invmass0+invmass3)*invmass0 * r0103; - double quad3_0203 = 2.0 * (invmass0+invmass3)*invmass0 * r0203; - - // iterate until converged - - double lamda01 = 0.0; - double lamda02 = 0.0; - double lamda03 = 0.0; - int niter = 0; - int done = 0; - - double quad1,quad2,quad3,b1,b2,b3,lamda01_new,lamda02_new,lamda03_new; - - while (!done && niter < max_iter) { - quad1 = quad1_0101 * lamda01*lamda01 + - quad1_0202 * lamda02*lamda02 + - quad1_0303 * lamda03*lamda03 + - quad1_0102 * lamda01*lamda02 + - quad1_0103 * lamda01*lamda03 + - quad1_0203 * lamda02*lamda03; - - quad2 = quad2_0101 * lamda01*lamda01 + - quad2_0202 * lamda02*lamda02 + - quad2_0303 * lamda03*lamda03 + - quad2_0102 * lamda01*lamda02 + - quad2_0103 * lamda01*lamda03 + - quad2_0203 * lamda02*lamda03; - - quad3 = quad3_0101 * lamda01*lamda01 + - quad3_0202 * lamda02*lamda02 + - quad3_0303 * lamda03*lamda03 + - quad3_0102 * lamda01*lamda02 + - quad3_0103 * lamda01*lamda03 + - quad3_0203 * lamda02*lamda03; - - b1 = bond1*bond1 - s01sq - quad1; - b2 = bond2*bond2 - s02sq - quad2; - b3 = bond3*bond3 - s03sq - quad3; - - lamda01_new = a11inv*b1 + a12inv*b2 + a13inv*b3; - lamda02_new = a21inv*b1 + a22inv*b2 + a23inv*b3; - lamda03_new = a31inv*b1 + a32inv*b2 + a33inv*b3; - - done = 1; - if (fabs(lamda01_new-lamda01) > tolerance) done = 0; - if (fabs(lamda02_new-lamda02) > tolerance) done = 0; - if (fabs(lamda03_new-lamda03) > tolerance) done = 0; - - lamda01 = lamda01_new; - lamda02 = lamda02_new; - lamda03 = lamda03_new; - niter++; - } - - // update forces if atom is owned by this processor - - lamda01 = lamda01/dtfsq; - lamda02 = lamda02/dtfsq; - lamda03 = lamda03/dtfsq; - - if (i0 < nlocal) { - f[i0][0] += lamda01*r01[0] + lamda02*r02[0] + lamda03*r03[0]; - f[i0][1] += lamda01*r01[1] + lamda02*r02[1] + lamda03*r03[1]; - f[i0][2] += lamda01*r01[2] + lamda02*r02[2] + lamda03*r03[2]; - } - - if (i1 < nlocal) { - f[i1][0] -= lamda01*r01[0]; - f[i1][1] -= lamda01*r01[1]; - f[i1][2] -= lamda01*r01[2]; - } - - if (i2 < nlocal) { - f[i2][0] -= lamda02*r02[0]; - f[i2][1] -= lamda02*r02[1]; - f[i2][2] -= lamda02*r02[2]; - } - - if (i3 < nlocal) { - f[i3][0] -= lamda03*r03[0]; - f[i3][1] -= lamda03*r03[1]; - f[i3][2] -= lamda03*r03[2]; - } - - if (evflag) { - nlist = 0; - if (i0 < nlocal) list[nlist++] = i0; - if (i1 < nlocal) list[nlist++] = i1; - if (i2 < nlocal) list[nlist++] = i2; - if (i3 < nlocal) list[nlist++] = i3; - - v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda03*r03[0]*r03[0]; - v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda03*r03[1]*r03[1]; - v[2] = lamda01*r01[2]*r01[2]+lamda02*r02[2]*r02[2]+lamda03*r03[2]*r03[2]; - v[3] = lamda01*r01[0]*r01[1]+lamda02*r02[0]*r02[1]+lamda03*r03[0]*r03[1]; - v[4] = lamda01*r01[0]*r01[2]+lamda02*r02[0]*r02[2]+lamda03*r03[0]*r03[2]; - v[5] = lamda01*r01[1]*r01[2]+lamda02*r02[1]*r02[2]+lamda03*r03[1]*r03[2]; - - v_tally(nlist,list,4.0,v); - } -} - -/* ---------------------------------------------------------------------- */ - -void FixShake::shake3angle(int m) -{ - int nlist,list[3]; - double v[6]; - double invmass0,invmass1,invmass2; - - // local atom IDs and constraint distances - - int i0 = atom->map(shake_atom[m][0]); - int i1 = atom->map(shake_atom[m][1]); - int i2 = atom->map(shake_atom[m][2]); - double bond1 = bond_distance[shake_type[m][0]]; - double bond2 = bond_distance[shake_type[m][1]]; - double bond12 = angle_distance[shake_type[m][2]]; - - // r01,r02,r12 = distance vec between atoms, with PBC - - double r01[3]; - r01[0] = x[i0][0] - x[i1][0]; - r01[1] = x[i0][1] - x[i1][1]; - r01[2] = x[i0][2] - x[i1][2]; - domain->minimum_image(r01); - - double r02[3]; - r02[0] = x[i0][0] - x[i2][0]; - r02[1] = x[i0][1] - x[i2][1]; - r02[2] = x[i0][2] - x[i2][2]; - domain->minimum_image(r02); - - double r12[3]; - r12[0] = x[i1][0] - x[i2][0]; - r12[1] = x[i1][1] - x[i2][1]; - r12[2] = x[i1][2] - x[i2][2]; - domain->minimum_image(r12); - - // s01,s02,s12 = distance vec after unconstrained update, with PBC - - double s01[3]; - s01[0] = xshake[i0][0] - xshake[i1][0]; - s01[1] = xshake[i0][1] - xshake[i1][1]; - s01[2] = xshake[i0][2] - xshake[i1][2]; - domain->minimum_image(s01); - - double s02[3]; - s02[0] = xshake[i0][0] - xshake[i2][0]; - s02[1] = xshake[i0][1] - xshake[i2][1]; - s02[2] = xshake[i0][2] - xshake[i2][2]; - domain->minimum_image(s02); - - double s12[3]; - s12[0] = xshake[i1][0] - xshake[i2][0]; - s12[1] = xshake[i1][1] - xshake[i2][1]; - s12[2] = xshake[i1][2] - xshake[i2][2]; - domain->minimum_image(s12); - - // scalar distances between atoms - - double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2]; - double r02sq = r02[0]*r02[0] + r02[1]*r02[1] + r02[2]*r02[2]; - double r12sq = r12[0]*r12[0] + r12[1]*r12[1] + r12[2]*r12[2]; - double s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2]; - double s02sq = s02[0]*s02[0] + s02[1]*s02[1] + s02[2]*s02[2]; - double s12sq = s12[0]*s12[0] + s12[1]*s12[1] + s12[2]*s12[2]; - - // matrix coeffs and rhs for lamda equations - - if (rmass) { - invmass0 = 1.0/rmass[i0]; - invmass1 = 1.0/rmass[i1]; - invmass2 = 1.0/rmass[i2]; - } else { - invmass0 = 1.0/mass[type[i0]]; - invmass1 = 1.0/mass[type[i1]]; - invmass2 = 1.0/mass[type[i2]]; - } - - double a11 = 2.0 * (invmass0+invmass1) * - (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); - double a12 = 2.0 * invmass0 * - (s01[0]*r02[0] + s01[1]*r02[1] + s01[2]*r02[2]); - double a13 = - 2.0 * invmass1 * - (s01[0]*r12[0] + s01[1]*r12[1] + s01[2]*r12[2]); - double a21 = 2.0 * invmass0 * - (s02[0]*r01[0] + s02[1]*r01[1] + s02[2]*r01[2]); - double a22 = 2.0 * (invmass0+invmass2) * - (s02[0]*r02[0] + s02[1]*r02[1] + s02[2]*r02[2]); - double a23 = 2.0 * invmass2 * - (s02[0]*r12[0] + s02[1]*r12[1] + s02[2]*r12[2]); - double a31 = - 2.0 * invmass1 * - (s12[0]*r01[0] + s12[1]*r01[1] + s12[2]*r01[2]); - double a32 = 2.0 * invmass2 * - (s12[0]*r02[0] + s12[1]*r02[1] + s12[2]*r02[2]); - double a33 = 2.0 * (invmass1+invmass2) * - (s12[0]*r12[0] + s12[1]*r12[1] + s12[2]*r12[2]); - - // inverse of matrix - - double determ = a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - - a11*a23*a32 - a12*a21*a33 - a13*a22*a31; - if (determ == 0.0) error->one(FLERR,"Shake determinant = 0.0"); - double determinv = 1.0/determ; - - double a11inv = determinv * (a22*a33 - a23*a32); - double a12inv = -determinv * (a12*a33 - a13*a32); - double a13inv = determinv * (a12*a23 - a13*a22); - double a21inv = -determinv * (a21*a33 - a23*a31); - double a22inv = determinv * (a11*a33 - a13*a31); - double a23inv = -determinv * (a11*a23 - a13*a21); - double a31inv = determinv * (a21*a32 - a22*a31); - double a32inv = -determinv * (a11*a32 - a12*a31); - double a33inv = determinv * (a11*a22 - a12*a21); - - // quadratic correction coeffs - - double r0102 = (r01[0]*r02[0] + r01[1]*r02[1] + r01[2]*r02[2]); - double r0112 = (r01[0]*r12[0] + r01[1]*r12[1] + r01[2]*r12[2]); - double r0212 = (r02[0]*r12[0] + r02[1]*r12[1] + r02[2]*r12[2]); - - double quad1_0101 = (invmass0+invmass1)*(invmass0+invmass1) * r01sq; - double quad1_0202 = invmass0*invmass0 * r02sq; - double quad1_1212 = invmass1*invmass1 * r12sq; - double quad1_0102 = 2.0 * (invmass0+invmass1)*invmass0 * r0102; - double quad1_0112 = - 2.0 * (invmass0+invmass1)*invmass1 * r0112; - double quad1_0212 = - 2.0 * invmass0*invmass1 * r0212; - - double quad2_0101 = invmass0*invmass0 * r01sq; - double quad2_0202 = (invmass0+invmass2)*(invmass0+invmass2) * r02sq; - double quad2_1212 = invmass2*invmass2 * r12sq; - double quad2_0102 = 2.0 * (invmass0+invmass2)*invmass0 * r0102; - double quad2_0112 = 2.0 * invmass0*invmass2 * r0112; - double quad2_0212 = 2.0 * (invmass0+invmass2)*invmass2 * r0212; - - double quad3_0101 = invmass1*invmass1 * r01sq; - double quad3_0202 = invmass2*invmass2 * r02sq; - double quad3_1212 = (invmass1+invmass2)*(invmass1+invmass2) * r12sq; - double quad3_0102 = - 2.0 * invmass1*invmass2 * r0102; - double quad3_0112 = - 2.0 * (invmass1+invmass2)*invmass1 * r0112; - double quad3_0212 = 2.0 * (invmass1+invmass2)*invmass2 * r0212; - - // iterate until converged - - double lamda01 = 0.0; - double lamda02 = 0.0; - double lamda12 = 0.0; - int niter = 0; - int done = 0; - - double quad1,quad2,quad3,b1,b2,b3,lamda01_new,lamda02_new,lamda12_new; - - while (!done && niter < max_iter) { - quad1 = quad1_0101 * lamda01*lamda01 + - quad1_0202 * lamda02*lamda02 + - quad1_1212 * lamda12*lamda12 + - quad1_0102 * lamda01*lamda02 + - quad1_0112 * lamda01*lamda12 + - quad1_0212 * lamda02*lamda12; - - quad2 = quad2_0101 * lamda01*lamda01 + - quad2_0202 * lamda02*lamda02 + - quad2_1212 * lamda12*lamda12 + - quad2_0102 * lamda01*lamda02 + - quad2_0112 * lamda01*lamda12 + - quad2_0212 * lamda02*lamda12; - - quad3 = quad3_0101 * lamda01*lamda01 + - quad3_0202 * lamda02*lamda02 + - quad3_1212 * lamda12*lamda12 + - quad3_0102 * lamda01*lamda02 + - quad3_0112 * lamda01*lamda12 + - quad3_0212 * lamda02*lamda12; - - b1 = bond1*bond1 - s01sq - quad1; - b2 = bond2*bond2 - s02sq - quad2; - b3 = bond12*bond12 - s12sq - quad3; - - lamda01_new = a11inv*b1 + a12inv*b2 + a13inv*b3; - lamda02_new = a21inv*b1 + a22inv*b2 + a23inv*b3; - lamda12_new = a31inv*b1 + a32inv*b2 + a33inv*b3; - - done = 1; - if (fabs(lamda01_new-lamda01) > tolerance) done = 0; - if (fabs(lamda02_new-lamda02) > tolerance) done = 0; - if (fabs(lamda12_new-lamda12) > tolerance) done = 0; - - lamda01 = lamda01_new; - lamda02 = lamda02_new; - lamda12 = lamda12_new; - niter++; - } - - // update forces if atom is owned by this processor - - lamda01 = lamda01/dtfsq; - lamda02 = lamda02/dtfsq; - lamda12 = lamda12/dtfsq; - - if (i0 < nlocal) { - f[i0][0] += lamda01*r01[0] + lamda02*r02[0]; - f[i0][1] += lamda01*r01[1] + lamda02*r02[1]; - f[i0][2] += lamda01*r01[2] + lamda02*r02[2]; - } - - if (i1 < nlocal) { - f[i1][0] -= lamda01*r01[0] - lamda12*r12[0]; - f[i1][1] -= lamda01*r01[1] - lamda12*r12[1]; - f[i1][2] -= lamda01*r01[2] - lamda12*r12[2]; - } - - if (i2 < nlocal) { - f[i2][0] -= lamda02*r02[0] + lamda12*r12[0]; - f[i2][1] -= lamda02*r02[1] + lamda12*r12[1]; - f[i2][2] -= lamda02*r02[2] + lamda12*r12[2]; - } - - if (evflag) { - nlist = 0; - if (i0 < nlocal) list[nlist++] = i0; - if (i1 < nlocal) list[nlist++] = i1; - if (i2 < nlocal) list[nlist++] = i2; - - v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda12*r12[0]*r12[0]; - v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda12*r12[1]*r12[1]; - v[2] = lamda01*r01[2]*r01[2]+lamda02*r02[2]*r02[2]+lamda12*r12[2]*r12[2]; - v[3] = lamda01*r01[0]*r01[1]+lamda02*r02[0]*r02[1]+lamda12*r12[0]*r12[1]; - v[4] = lamda01*r01[0]*r01[2]+lamda02*r02[0]*r02[2]+lamda12*r12[0]*r12[2]; - v[5] = lamda01*r01[1]*r01[2]+lamda02*r02[1]*r02[2]+lamda12*r12[1]*r12[2]; - - v_tally(nlist,list,3.0,v); - } -} - -/* ---------------------------------------------------------------------- - print-out bond & angle statistics -------------------------------------------------------------------------- */ - -void FixShake::stats() -{ - int i,j,m,n,iatom,jatom,katom; - double delx,dely,delz; - double r,r1,r2,r3,angle; - - // zero out accumulators - - int nb = atom->nbondtypes + 1; - int na = atom->nangletypes + 1; - - for (i = 0; i < nb; i++) { - b_count[i] = 0; - b_ave[i] = b_max[i] = 0.0; - b_min[i] = BIG; - } - for (i = 0; i < na; i++) { - a_count[i] = 0; - a_ave[i] = a_max[i] = 0.0; - a_min[i] = BIG; - } - - // log stats for each bond & angle - // OK to double count since are just averaging - - double **x = atom->x; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) { - if (shake_flag[i] == 0) continue; - - // bond stats - - n = shake_flag[i]; - if (n == 1) n = 3; - iatom = atom->map(shake_atom[i][0]); - for (j = 1; j < n; j++) { - jatom = atom->map(shake_atom[i][j]); - delx = x[iatom][0] - x[jatom][0]; - dely = x[iatom][1] - x[jatom][1]; - delz = x[iatom][2] - x[jatom][2]; - domain->minimum_image(delx,dely,delz); - r = sqrt(delx*delx + dely*dely + delz*delz); - - m = shake_type[i][j-1]; - b_count[m]++; - b_ave[m] += r; - b_max[m] = MAX(b_max[m],r); - b_min[m] = MIN(b_min[m],r); - } - - // angle stats - - if (shake_flag[i] == 1) { - iatom = atom->map(shake_atom[i][0]); - jatom = atom->map(shake_atom[i][1]); - katom = atom->map(shake_atom[i][2]); - - delx = x[iatom][0] - x[jatom][0]; - dely = x[iatom][1] - x[jatom][1]; - delz = x[iatom][2] - x[jatom][2]; - domain->minimum_image(delx,dely,delz); - r1 = sqrt(delx*delx + dely*dely + delz*delz); - - delx = x[iatom][0] - x[katom][0]; - dely = x[iatom][1] - x[katom][1]; - delz = x[iatom][2] - x[katom][2]; - domain->minimum_image(delx,dely,delz); - r2 = sqrt(delx*delx + dely*dely + delz*delz); - - delx = x[jatom][0] - x[katom][0]; - dely = x[jatom][1] - x[katom][1]; - delz = x[jatom][2] - x[katom][2]; - domain->minimum_image(delx,dely,delz); - r3 = sqrt(delx*delx + dely*dely + delz*delz); - - angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2)); - angle *= 180.0/MY_PI; - m = shake_type[i][2]; - a_count[m]++; - a_ave[m] += angle; - a_max[m] = MAX(a_max[m],angle); - a_min[m] = MIN(a_min[m],angle); - } - } - - // sum across all procs - - MPI_Allreduce(b_count,b_count_all,nb,MPI_INT,MPI_SUM,world); - MPI_Allreduce(b_ave,b_ave_all,nb,MPI_DOUBLE,MPI_SUM,world); - MPI_Allreduce(b_max,b_max_all,nb,MPI_DOUBLE,MPI_MAX,world); - MPI_Allreduce(b_min,b_min_all,nb,MPI_DOUBLE,MPI_MIN,world); - - MPI_Allreduce(a_count,a_count_all,na,MPI_INT,MPI_SUM,world); - MPI_Allreduce(a_ave,a_ave_all,na,MPI_DOUBLE,MPI_SUM,world); - MPI_Allreduce(a_max,a_max_all,na,MPI_DOUBLE,MPI_MAX,world); - MPI_Allreduce(a_min,a_min_all,na,MPI_DOUBLE,MPI_MIN,world); - - // print stats only for non-zero counts - - if (me == 0) { - - if (screen) { - fprintf(screen, - "SHAKE stats (type/ave/delta) on step " BIGINT_FORMAT "\n", - update->ntimestep); - for (i = 1; i < nb; i++) - if (b_count_all[i]) - fprintf(screen," %d %g %g %d\n",i, - b_ave_all[i]/b_count_all[i],b_max_all[i]-b_min_all[i], - b_count_all[i]); - for (i = 1; i < na; i++) - if (a_count_all[i]) - fprintf(screen," %d %g %g\n",i, - a_ave_all[i]/a_count_all[i],a_max_all[i]-a_min_all[i]); - } - if (logfile) { - fprintf(logfile, - "SHAKE stats (type/ave/delta) on step " BIGINT_FORMAT "\n", - update->ntimestep); - for (i = 0; i < nb; i++) - if (b_count_all[i]) - fprintf(logfile," %d %g %g\n",i, - b_ave_all[i]/b_count_all[i],b_max_all[i]-b_min_all[i]); - for (i = 0; i < na; i++) - if (a_count_all[i]) - fprintf(logfile," %d %g %g\n",i, - a_ave_all[i]/a_count_all[i],a_max_all[i]-a_min_all[i]); - } - } - - // next timestep for stats - - next_output += output_every; -} - -/* ---------------------------------------------------------------------- - find a bond between global tags n1 and n2 stored with local atom i - return -1 if don't find it - return bond index if do find it -------------------------------------------------------------------------- */ - -int FixShake::bondfind(int i, int n1, int n2) -{ - int *tag = atom->tag; - int **bond_atom = atom->bond_atom; - int nbonds = atom->num_bond[i]; - - int m; - for (m = 0; m < nbonds; m++) { - if (n1 == tag[i] && n2 == bond_atom[i][m]) break; - if (n1 == bond_atom[i][m] && n2 == tag[i]) break; - } - if (m < nbonds) return m; - return -1; -} - -/* ---------------------------------------------------------------------- - find an angle with global end atoms n1 and n2 stored with local atom i - return -1 if don't find it - return angle index if do find it -------------------------------------------------------------------------- */ - -int FixShake::anglefind(int i, int n1, int n2) -{ - int **angle_atom1 = atom->angle_atom1; - int **angle_atom3 = atom->angle_atom3; - int nangles = atom->num_angle[i]; - - int m; - for (m = 0; m < nangles; m++) { - if (n1 == angle_atom1[i][m] && n2 == angle_atom3[i][m]) break; - if (n1 == angle_atom3[i][m] && n2 == angle_atom1[i][m]) break; - } - if (m < nangles) return m; - return -1; -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based arrays -------------------------------------------------------------------------- */ - -double FixShake::memory_usage() -{ - int nmax = atom->nmax; - double bytes = nmax * sizeof(int); - bytes += nmax*4 * sizeof(int); - bytes += nmax*3 * sizeof(int); - bytes += nmax*3 * sizeof(double); - bytes += maxvatom*6 * sizeof(double); - return bytes; -} - -/* ---------------------------------------------------------------------- - allocate local atom-based arrays -------------------------------------------------------------------------- */ - -void FixShake::grow_arrays(int nmax) -{ - memory->grow(shake_flag,nmax,"shake:shake_flag"); - memory->grow(shake_atom,nmax,4,"shake:shake_atom"); - memory->grow(shake_type,nmax,3,"shake:shake_type"); - memory->destroy(xshake); - memory->create(xshake,nmax,3,"shake:xshake"); -} - -/* ---------------------------------------------------------------------- - copy values within local atom-based arrays -------------------------------------------------------------------------- */ - -void FixShake::copy_arrays(int i, int j, int delflag) -{ - int flag = shake_flag[j] = shake_flag[i]; - if (flag == 1) { - shake_atom[j][0] = shake_atom[i][0]; - shake_atom[j][1] = shake_atom[i][1]; - shake_atom[j][2] = shake_atom[i][2]; - shake_type[j][0] = shake_type[i][0]; - shake_type[j][1] = shake_type[i][1]; - shake_type[j][2] = shake_type[i][2]; - } else if (flag == 2) { - shake_atom[j][0] = shake_atom[i][0]; - shake_atom[j][1] = shake_atom[i][1]; - shake_type[j][0] = shake_type[i][0]; - } else if (flag == 3) { - shake_atom[j][0] = shake_atom[i][0]; - shake_atom[j][1] = shake_atom[i][1]; - shake_atom[j][2] = shake_atom[i][2]; - shake_type[j][0] = shake_type[i][0]; - shake_type[j][1] = shake_type[i][1]; - } else if (flag == 4) { - shake_atom[j][0] = shake_atom[i][0]; - shake_atom[j][1] = shake_atom[i][1]; - shake_atom[j][2] = shake_atom[i][2]; - shake_atom[j][3] = shake_atom[i][3]; - shake_type[j][0] = shake_type[i][0]; - shake_type[j][1] = shake_type[i][1]; - shake_type[j][2] = shake_type[i][2]; - } -} - - -/* ---------------------------------------------------------------------- - initialize one atom's array values, called when atom is created -------------------------------------------------------------------------- */ - -void FixShake::set_arrays(int i) -{ - shake_flag[i] = 0; -} - -/* ---------------------------------------------------------------------- - update one atom's array values - called when molecule is created from fix gcmc -------------------------------------------------------------------------- */ - -void FixShake::update_arrays(int i, int atom_offset) -{ - int flag = shake_flag[i]; - - if (flag == 1) { - shake_atom[i][0] += atom_offset; - shake_atom[i][1] += atom_offset; - shake_atom[i][2] += atom_offset; - } else if (flag == 2) { - shake_atom[i][0] += atom_offset; - shake_atom[i][1] += atom_offset; - } else if (flag == 3) { - shake_atom[i][0] += atom_offset; - shake_atom[i][1] += atom_offset; - shake_atom[i][2] += atom_offset; - } else if (flag == 4) { - shake_atom[i][0] += atom_offset; - shake_atom[i][1] += atom_offset; - shake_atom[i][2] += atom_offset; - shake_atom[i][3] += atom_offset; - } -} - -/* ---------------------------------------------------------------------- - 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,vcm,quat ignored -------------------------------------------------------------------------- */ - -void FixShake::set_molecule(int nlocalprev, int tagprev, - double *xgeom, double *vcm, double *quat) -{ - int m,flag; - - int nlocal = atom->nlocal; - if (nlocalprev == nlocal) return; - - int *tag = atom->tag; - int **mol_shake_atom = onemol->shake_atom; - int **mol_shake_type = onemol->shake_type; - - for (int i = nlocalprev; i < nlocal; i++) { - m = tag[i] - tagprev-1; - - flag = shake_flag[i] = onemol->shake_flag[m]; - - if (flag == 1) { - shake_atom[i][0] = mol_shake_atom[m][0] + tagprev; - shake_atom[i][1] = mol_shake_atom[m][1] + tagprev; - shake_atom[i][2] = mol_shake_atom[m][2] + tagprev; - shake_type[i][0] = mol_shake_type[m][0]; - shake_type[i][1] = mol_shake_type[m][1]; - shake_type[i][2] = mol_shake_type[m][2]; - } else if (flag == 2) { - shake_atom[i][0] = mol_shake_atom[m][0] + tagprev; - shake_atom[i][1] = mol_shake_atom[m][1] + tagprev; - shake_type[i][0] = mol_shake_type[m][0]; - } else if (flag == 3) { - shake_atom[i][0] = mol_shake_atom[m][0] + tagprev; - shake_atom[i][1] = mol_shake_atom[m][1] + tagprev; - shake_atom[i][2] = mol_shake_atom[m][2] + tagprev; - shake_type[i][0] = mol_shake_type[m][0]; - shake_type[i][1] = mol_shake_type[m][1]; - } else if (flag == 4) { - shake_atom[i][0] = mol_shake_atom[m][0] + tagprev; - shake_atom[i][1] = mol_shake_atom[m][1] + tagprev; - shake_atom[i][2] = mol_shake_atom[m][2] + tagprev; - shake_atom[i][3] = mol_shake_atom[m][3] + tagprev; - shake_type[i][0] = mol_shake_type[m][0]; - shake_type[i][1] = mol_shake_type[m][1]; - shake_type[i][2] = mol_shake_type[m][2]; - } - } -} - -/* ---------------------------------------------------------------------- - pack values in local atom-based arrays for exchange with another proc -------------------------------------------------------------------------- */ - -int FixShake::pack_exchange(int i, double *buf) -{ - int m = 0; - buf[m++] = shake_flag[i]; - int flag = shake_flag[i]; - if (flag == 1) { - buf[m++] = shake_atom[i][0]; - buf[m++] = shake_atom[i][1]; - buf[m++] = shake_atom[i][2]; - buf[m++] = shake_type[i][0]; - buf[m++] = shake_type[i][1]; - buf[m++] = shake_type[i][2]; - } else if (flag == 2) { - buf[m++] = shake_atom[i][0]; - buf[m++] = shake_atom[i][1]; - buf[m++] = shake_type[i][0]; - } else if (flag == 3) { - buf[m++] = shake_atom[i][0]; - buf[m++] = shake_atom[i][1]; - buf[m++] = shake_atom[i][2]; - buf[m++] = shake_type[i][0]; - buf[m++] = shake_type[i][1]; - } else if (flag == 4) { - buf[m++] = shake_atom[i][0]; - buf[m++] = shake_atom[i][1]; - buf[m++] = shake_atom[i][2]; - buf[m++] = shake_atom[i][3]; - buf[m++] = shake_type[i][0]; - buf[m++] = shake_type[i][1]; - buf[m++] = shake_type[i][2]; - } - return m; -} - -/* ---------------------------------------------------------------------- - unpack values in local atom-based arrays from exchange with another proc -------------------------------------------------------------------------- */ - -int FixShake::unpack_exchange(int nlocal, double *buf) -{ - int m = 0; - int flag = shake_flag[nlocal] = static_cast (buf[m++]); - if (flag == 1) { - shake_atom[nlocal][0] = static_cast (buf[m++]); - shake_atom[nlocal][1] = static_cast (buf[m++]); - shake_atom[nlocal][2] = static_cast (buf[m++]); - shake_type[nlocal][0] = static_cast (buf[m++]); - shake_type[nlocal][1] = static_cast (buf[m++]); - shake_type[nlocal][2] = static_cast (buf[m++]); - } else if (flag == 2) { - shake_atom[nlocal][0] = static_cast (buf[m++]); - shake_atom[nlocal][1] = static_cast (buf[m++]); - shake_type[nlocal][0] = static_cast (buf[m++]); - } else if (flag == 3) { - shake_atom[nlocal][0] = static_cast (buf[m++]); - shake_atom[nlocal][1] = static_cast (buf[m++]); - shake_atom[nlocal][2] = static_cast (buf[m++]); - shake_type[nlocal][0] = static_cast (buf[m++]); - shake_type[nlocal][1] = static_cast (buf[m++]); - } else if (flag == 4) { - shake_atom[nlocal][0] = static_cast (buf[m++]); - shake_atom[nlocal][1] = static_cast (buf[m++]); - shake_atom[nlocal][2] = static_cast (buf[m++]); - shake_atom[nlocal][3] = static_cast (buf[m++]); - shake_type[nlocal][0] = static_cast (buf[m++]); - shake_type[nlocal][1] = static_cast (buf[m++]); - shake_type[nlocal][2] = static_cast (buf[m++]); - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int FixShake::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) -{ - int i,j,m; - double dx,dy,dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = xshake[j][0]; - buf[m++] = xshake[j][1]; - buf[m++] = xshake[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; - dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; - dz = pbc[2]*domain->zprd; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = xshake[j][0] + dx; - buf[m++] = xshake[j][1] + dy; - buf[m++] = xshake[j][2] + dz; - } - } - return 3; -} - -/* ---------------------------------------------------------------------- */ - -void FixShake::unpack_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - xshake[i][0] = buf[m++]; - xshake[i][1] = buf[m++]; - xshake[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixShake::reset_dt() -{ - if (strstr(update->integrate_style,"verlet")) { - dtv = update->dt; - dtfsq = update->dt * update->dt * force->ftm2v; - } else { - dtv = step_respa[0]; - dtf_innerhalf = 0.5 * step_respa[0] * force->ftm2v; - dtf_inner = step_respa[0] * force->ftm2v; - } -} - -/* ---------------------------------------------------------------------- - extract Molecule ptr -------------------------------------------------------------------------- */ - -void *FixShake::extract(const char *str, int &dim) -{ - dim = 0; - if (strcmp(str,"onemol") == 0) { - return onemol; - } - return NULL; -} diff --git a/src/fix_shake.h b/src/fix_shake.h deleted file mode 100644 index a1cfe0e2f4..0000000000 --- a/src/fix_shake.h +++ /dev/null @@ -1,242 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#ifdef FIX_CLASS - -FixStyle(shake,FixShake) - -#else - -#ifndef LMP_FIX_SHAKE_H -#define LMP_FIX_SHAKE_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixShake : public Fix { - public: - FixShake(class LAMMPS *, int, char **); - ~FixShake(); - int setmask(); - void init(); - void setup(int); - void pre_neighbor(); - void post_force(int); - void post_force_respa(int, int, int); - - double memory_usage(); - void grow_arrays(int); - void copy_arrays(int, int, int); - void set_arrays(int); - void update_arrays(int, int); - void set_molecule(int, int, double *, double *, double *); - - int pack_exchange(int, double *); - int unpack_exchange(int, double *); - int pack_comm(int, int *, double *, int, int *); - void unpack_comm(int, int, double *); - - int dof(int); - void reset_dt(); - void *extract(const char *, int &); - - private: - int me,nprocs; - double tolerance; // SHAKE tolerance - int max_iter; // max # of SHAKE iterations - int output_every; // SHAKE stat output every so often - bigint next_output; // timestep for next output - - // settings from input command - int *bond_flag,*angle_flag; // bond/angle types to constrain - int *type_flag; // constrain bonds to these types - double *mass_list; // constrain bonds to these masses - int nmass; // # of masses in mass_list - - double *bond_distance,*angle_distance; // constraint distances - - int ifix_respa; // rRESPA fix needed by SHAKE - int nlevels_respa; // copies of needed rRESPA variables - int *loop_respa; - double *step_respa; - - double **x,**v,**f; // local ptrs to atom class quantities - double *mass,*rmass; - int *type; - int nlocal; - // atom-based arrays - int *shake_flag; // 0 if atom not in SHAKE cluster - // 1 = size 3 angle cluster - // 2,3,4 = size of bond-only cluster - int **shake_atom; // global IDs of atoms in cluster - // central atom is 1st - // lowest global ID is 1st for size 2 - int **shake_type; // bondtype of each bond in cluster - // for angle cluster, 3rd value - // is angletype - double **xshake; // unconstrained atom coords - int *nshake; // count - - double dtv,dtfsq; // timesteps for trial move - double dtf_inner,dtf_innerhalf; // timesteps for rRESPA trial move - - int *list; // list of clusters to SHAKE - int nlist,maxlist; // size and max-size of list - - // stat quantities - int *b_count,*b_count_all; // counts for each bond type - double *b_ave,*b_max,*b_min; // ave/max/min dist for each bond type - double *b_ave_all,*b_max_all,*b_min_all; // MPI summing arrays - int *a_count,*a_count_all; // ditto for angle types - double *a_ave,*a_max,*a_min; - double *a_ave_all,*a_max_all,*a_min_all; - - // molecules added on-the-fly with SHAKE constraints - - class Molecule *onemol; - - void find_clusters(); - int masscheck(double); - void unconstrained_update(); - void unconstrained_update_respa(int); - void shake(int); - void shake3(int); - void shake4(int); - void shake3angle(int); - void stats(); - int bondfind(int, int, int); - int anglefind(int, int, int); - - // static variable for ring communication callback to access class data - // callback functions for ring communication - - static FixShake *fsptr; - static void ring_bonds(int, char *); - static void ring_nshake(int, char *); - static void ring_shake(int, char *); -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Cannot use fix shake with non-molecular system - -Your choice of atom style does not have bonds. - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: Invalid bond type index for fix shake - -Self-explanatory. Check the fix shake command in the input script. - -E: Invalid angle type index for fix shake - -Self-explanatory. - -E: Invalid atom type index for fix shake - -Atom types must range from 1 to Ntypes inclusive. - -E: Invalid atom mass for fix shake - -Mass specified in fix shake command must be > 0.0. - -E: Too many masses for fix shake - -The fix shake command cannot list more masses than there are atom -types. - -E: More than one fix shake - -Only one fix shake can be defined. - -E: Fix shake cannot be used with minimization - -Cannot use fix shake while doing an energy minimization since -it turns off bonds that should contribute to the energy. - -E: Shake fix must come before NPT/NPH fix - -NPT fix must be defined in input script after SHAKE fix, else the -SHAKE fix contribution to the pressure virial is incorrect. - -E: Bond potential must be defined for SHAKE - -Cannot use fix shake unless bond potential is defined. - -E: Angle potential must be defined for SHAKE - -When shaking angles, an angle_style potential must be used. - -E: Shake angles have different bond types - -All 3-atom angle-constrained SHAKE clusters specified by the fix shake -command that are the same angle type, must also have the same bond -types for the 2 bonds in the angle. - -E: Shake atoms %d %d missing on proc %d at step %ld - -The 2 atoms in a single shake cluster specified by the fix shake -command are not all accessible to a processor. This probably means -an atom has moved too far. - -E: Shake atoms %d %d %d missing on proc %d at step %ld - -The 3 atoms in a single shake cluster specified by the fix shake -command are not all accessible to a processor. This probably means -an atom has moved too far. - -E: Shake atoms %d %d %d %d missing on proc %d at step %ld - -The 4 atoms in a single shake cluster specified by the fix shake -command are not all accessible to a processor. This probably means -an atom has moved too far. - -E: Did not find fix shake partner info - -Could not find bond partners implied by fix shake command. This error -can be triggered if the delete_bonds command was used before fix -shake, and it removed bonds without resetting the 1-2, 1-3, 1-4 -weighting list via the special keyword. - -E: Shake cluster of more than 4 atoms - -A single cluster specified by the fix shake command can have no more -than 4 atoms. - -E: Shake clusters are connected - -A single cluster specified by the fix shake command must have a single -central atom with up to 3 other atoms bonded to it. - -W: Shake determinant < 0.0 - -The determinant of the quadratic equation being solved for a single -cluster specified by the fix shake command is numerically suspect. LAMMPS -will set it to 0.0 and continue. - -E: Shake determinant = 0.0 - -The determinant of the matrix being solved for a single cluster -specified by the fix shake command is numerically invalid. - -*/