diff --git a/src/fix_pair_tracker.cpp b/src/fix_pair_tracker.cpp new file mode 100755 index 0000000000..6030fc0881 --- /dev/null +++ b/src/fix_pair_tracker.cpp @@ -0,0 +1,444 @@ +/* ---------------------------------------------------------------------- + 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 +#include "fix_pair_tracking.h" +#include "pair_tracker.h" +#include "atom.h" +#include "atom_vec.h" +#include "update.h" +#include "force.h" +#include "pair.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "modify.h" +#include "error.h" +#include "group.h" +#include "domain.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define DELTA 100 + +/* ---------------------------------------------------------------------- */ + +FixPairTracking::FixPairTracking(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), + nvalues(0), + array(NULL), vector(NULL), pack_choice(NULL) +{ + if (narg < 3) error->all(FLERR,"Illegal fix pair/tracker command"); + store_flag = 0; + local_flag = 1; + nvalues = narg - 4; + + nevery = utils::inumeric(FLERR,arg[3],false,lmp); + if (nevery <= 0) error->all(FLERR,"Illegal fix pair/tracker command"); + + if (nvalues == 1) size_local_cols = 0; + else size_local_cols = nvalues; + + pack_choice = new FnPtrPack[nvalues]; + + int i; + for (int iarg = 4; iarg < narg; iarg++) { + i = iarg - 4; + + if (strcmp(arg[iarg],"id1") == 0) { + pack_choice[i] = &FixPairTracking::pack_id1; + } else if (strcmp(arg[iarg],"id2") == 0) { + pack_choice[i] = &FixPairTracking::pack_id2; + + } else if (strcmp(arg[iarg],"time/created") == 0) { + pack_choice[i] = &FixPairTracking::pack_time_created; + } else if (strcmp(arg[iarg],"time/broken") == 0) { + pack_choice[i] = &FixPairTracking::pack_time_broken; + } else if (strcmp(arg[iarg],"time/total") == 0) { + pack_choice[i] = &FixPairTracking::pack_time_total; + + } else if (strcmp(arg[iarg],"x") == 0) { + pack_choice[i] = &FixPairTracking::pack_x; + } else if (strcmp(arg[iarg],"y") == 0) { + pack_choice[i] = &FixPairTracking::pack_y; + } else if (strcmp(arg[iarg],"z") == 0) { + pack_choice[i] = &FixPairTracking::pack_z; + + } else if (strcmp(arg[iarg],"xstore") == 0) { + pack_choice[i] = &FixPairTracking::pack_xstore; + store_flag = 1; + } else if (strcmp(arg[iarg],"ystore") == 0) { + pack_choice[i] = &FixPairTracking::pack_ystore; + store_flag = 1; + } else if (strcmp(arg[iarg],"zstore") == 0) { + pack_choice[i] = &FixPairTracking::pack_zstore; + store_flag = 1; + + } else if (strcmp(arg[iarg],"rmin") == 0) { + pack_choice[i] = &FixPairTracking::pack_rmin; + } else if (strcmp(arg[iarg],"rave") == 0) { + pack_choice[i] = &FixPairTracking::pack_rave; + + } else error->all(FLERR, "Invalid keyword in fix pair/tracker command"); + } + + nmax = 0; + ncount = 0; + vector = NULL; + array = NULL; +} + +/* ---------------------------------------------------------------------- */ + +FixPairTracking::~FixPairTracking() +{ + if (modify->nfix & store_flag == 1) { + modify->delete_fix(id_fix); + delete [] id_fix; + } + + delete [] pack_choice; + + memory->destroy(vector); + memory->destroy(array); +} + +/* ---------------------------------------------------------------------- */ + +int FixPairTracking::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::post_constructor() +{ + //If use stored x,y,z values, use store property (can transfer to ghost atoms) to store positions + + if(store_flag == 1){ + + lx = domain->xprd; + ly = domain->yprd; + lz = domain->zprd; + + int nn = strlen(id) + strlen("_FIX_PROP_ATOM") + 1; + id_property_fix = new char[nn]; + strcpy(id_property_fix,id); + strcat(id_property_fix,"_FIX_PROP_ATOM"); + + int ifix = modify->find_fix(id_property_fix); + if (ifix < 0) { + + int n_x = strlen(id) + 4; + + char * lab1 = new char[n_x]; + strcpy(lab1, "d_x"); + strcat(lab1, id); + char * lab2 = new char[n_x]; + strcpy(lab2, "d_y"); + strcat(lab2, id); + char * lab3 = new char[n_x]; + strcpy(lab3, "d_z"); + strcat(lab3, id); + + char **newarg = new char*[8]; + newarg[0] = id_property_fix; + newarg[1] = group->names[igroup]; + newarg[2] = (char *) "property/atom"; + newarg[3] = (char *) lab1; + newarg[4] = (char *) lab2; + newarg[5] = (char *) lab3; + newarg[6] = (char *) "ghost"; + newarg[7] = (char *) "yes"; + + modify->add_fix(8,newarg); + //Needs ghost atoms to calculate CoM + + int type_flag; + int col_flag; + + strcpy(lab1, "x"); + strcat(lab1, id); + strcpy(lab2, "y"); + strcat(lab2, id); + strcpy(lab3, "z"); + strcat(lab3, id); + + index_x = atom->find_custom(lab1, type_flag, col_flag); + index_y = atom->find_custom(lab2, type_flag, col_flag); + index_z = atom->find_custom(lab3, type_flag, col_flag); + delete [] newarg; + delete [] lab1; + delete [] lab2; + delete [] lab3; + } + + ifix = modify->find_fix(id_fix); + if (ifix < 0) error->all(FLERR,"Could not find fix ID for fix broken/bond"); + if (modify->fix[ifix]->restart_reset) { + modify->fix[ifix]->restart_reset = 0; + } else { + + double *xi = atom->dvector[index_x]; + double *yi = atom->dvector[index_y]; + double *zi = atom->dvector[index_z]; + + double **xs = atom->x; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + xi[i] = xs[i][0]; + yi[i] = xs[i][1]; + zi[i] = xs[i][2]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::init() +{ + // Set size of array/vector + ncount = 0; + + if (ncount > nmax) { + reallocate(ncount);} + size_local_rows = ncount; +} + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::lost_contact(int i, int j, double n, double rs, double rm) +{ + if (ncount == nmax) reallocate(ncount); + + index_i = i; + index_j = j; + rmin = rm; + rave = ra; + ntimestep = n; + + // fill vector or array with local values + if (nvalues == 1) { + (this->*pack_choice[0])(0); + } else { + for (int n = 0; n < nvalues; n++) + (this->*pack_choice[n])(n); + } + + ncount += 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::post_force(int /*vflag*/) +{ + if (update->ntimestep % nevery == 0) { + size_local_rows = ncount; + ncount = 0; + } +} + + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::reallocate(int n) +{ + // grow vector or array + while (nmax <= n) nmax += DELTA; + + if (nvalues == 1) { + memory->grow(vector,nmax,"fix_broken_bonds:vector"); + vector_local = vector; + } else { + memory->grow(array,nmax,nvalues,"fix_broken_bonds:array"); + array_local = array; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double FixPairTracking::memory_usage() +{ + double bytes = nmax*nvalues * sizeof(double); + bytes += nmax*2 * sizeof(int); + return bytes; +} + +/* ---------------------------------------------------------------------- + one method for every keyword fix property/local can output + the atom property is packed into buf starting at n with stride nvalues + customize a new keyword by adding a method +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_time_created(int n) +{ + if (nvalues == 1) + vector[ncount] = ntimestep; + else + array[ncount][n] = ntimestep; +} + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_time_broken(int n) +{ + if (nvalues == 1) + vector[ncount] = update->ntimestep; + else + array[ncount][n] = update->ntimestep; +} + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_time_total(int n) +{ + if (nvalues == 1) + vector[ncount] = update->ntimestep-ntimestep; + else + array[ncount][n] = update->ntimestep-ntimestep; +} + + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_id1(int n) +{ + tagint *tag = atom->tag; + + if (nvalues == 1) + vector[ncount] = tag[index_i]; + else + array[ncount][n] = tag[index_i]; +} + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_id2(int n) +{ + tagint *tag = atom->tag; + + if (nvalues == 1) + vector[ncount] = tag[index_j]; + else + array[ncount][n] = tag[index_j]; +} + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_x(int n) +{ + double lx_new = domain->xprd; + double **x = atom->x; + + if (nvalues == 1) + vector[ncount] = (x[index_i][0] + x[index_j][0])/2; + else + array[ncount][n] = (x[index_i][0] + x[index_j][0])/2; +} + + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_y(int n) +{ + double lx_new = domain->yprd; + double **x = atom->x; + + if (nvalues == 1) + vector[ncount] = (x[index_i][1] + x[index_j][1])/2; + else + array[ncount][n] = (x[index_i][1] + x[index_j][1])/2; +} + + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_z(int n) +{ + double lx_new = domain->zprd; + double **x = atom->x; + + if (nvalues == 1) + vector[ncount] = (x[index_i][2] + x[index_j][2])/2; + else + array[ncount][n] = (x[index_i][2] + x[index_j][2])/2; +} + + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_xstore(int n) +{ + double *x = atom->dvector[index_x]; + + if (nvalues == 1) + vector[ncount] = (x[index_i] + x[index_j])/2; + else + array[ncount][n] = (x[index_i] + x[index_j])/2; +} + + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_ystore(int n) +{ + double *y = atom->dvector[index_y]; + + if (nvalues == 1) + vector[ncount] = (y[index_i] + y[index_j])/2; + else + array[ncount][n] = (y[index_i] + y[index_j])/2; +} + + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_zstore(int n) +{ + double *z = atom->dvector[index_z]; + + if (nvalues == 1) + vector[ncount] = (z[index_i] + z[index_j])/2; + else + array[ncount][n] = (z[index_i] + z[index_j])/2; +} + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_rmin(int n) +{ + if (nvalues == 1) + vector[ncount] = rmin; + else + array[ncount][n] = rmin; +} + +/* ---------------------------------------------------------------------- */ + +void FixPairTracking::pack_rave(int n) +{ + if (nvalues == 1) + vector[ncount] = rsum/(update->ntimestep-ntimestep); + else + array[ncount][n] = rsum/(update->ntimestep-ntimestep); +} diff --git a/src/fix_pair_tracker.h b/src/fix_pair_tracker.h new file mode 100755 index 0000000000..8f63bac6d9 --- /dev/null +++ b/src/fix_pair_tracker.h @@ -0,0 +1,119 @@ +/* ---------------------------------------------------------------------- + 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(pair/tracking,FixPairTracking) + +#else + +#ifndef LMP_FIX_PAIR_TRACKING_H +#define LMP_FIX_PAIR_TRACKING_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixPairTracking : public Fix { + public: + FixPairTracking(class LAMMPS *, int, char **); + ~FixPairTracking(); + int setmask(); + void post_constructor(); + void init(); + void post_force(int); + double memory_usage(); + void lost_contact(int, int, double, double, double); + + private: + int nvalues, never; + int nmax; + int store_flag; + int index_i, index_j; + + double *vector; + double **array; + + double lx; + double ly; + double lz; + + int ncount; + + void reallocate(int); + + typedef void (FixPairTracker::*FnPtrPack)(int); + FnPtrPack *pack_choice; // ptrs to pack functions + + void pack_id1(int); + void pack_id2(int); + + void pack_time_created(int); + void pack_time_broken(int); + void pack_time_total(int); + + void pack_x(int); + void pack_y(int); + void pack_z(int); + + void pack_xstore(int); + void pack_ystore(int); + void pack_zstore(int); + + void pack_rmin(int); + void pack_rave(int); + + char *id_fix; + int index_x, index_y, index_z; +}; + +} + +#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: Fix property/local cannot use these inputs together + +Only inputs that generate the same number of datums can be used +togther. E.g. bond and angle quantities cannot be mixed. + +E: Invalid keyword in compute property/local command + +Self-explanatory. + +E: Fix property/local does not (yet) work with atom_style template + +Self-explanatory. + +E: Fix property/local for property that isn't allocated + +Self-explanatory. + +E: No pair style is defined for compute property/local + +Self-explanatory. + +E: Pair style does not support compute property/local + +The pair style does not have a single() function, so it can +not be invoked by fix bond/swap. + +*/ diff --git a/src/pair_tracker.cpp b/src/pair_tracker.cpp new file mode 100644 index 0000000000..4f40c1bdfe --- /dev/null +++ b/src/pair_tracker.cpp @@ -0,0 +1,469 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://lammps.sandia.gov/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Leo Silbert (SNL), Gary Grest (SNL) +------------------------------------------------------------------------- */ + +#include "pair_tracker.h" + +#include +#include + +#include "atom.h" +#include "force.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "fix_dummy.h" +#include "fix_neigh_history.h" +#include "comm.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 1; + no_virial_fdotr_compute = 1; + history = 1; + size_history = 3; + + neighprev = 0; + nondefault_history_transfer = 1; + + finitecutflag = 0; + + // create dummy fix as placeholder for FixNeigTRACKistory + // this is so final order of Modify:fix will conform to input script + + fix_history = nullptr; + modify->add_fix("NEIGH_HISTORY_TRACK_DUMMY all DUMMY"); + fix_dummy = (FixDummy *) modify->fix[modify->nfix-1]; +} + +/* ---------------------------------------------------------------------- */ + +PairTracker::~PairTracker() +{ + if (!fix_history) modify->delete_fix("NEIGH_HISTORY_TRACK_DUMMY"); + else modify->delete_fix("NEIGH_HISTORY_TRACK"); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(cut); + + delete [] onerad_dynamic; + delete [] onerad_frozen; + delete [] maxrad_dynamic; + delete [] maxrad_frozen; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairTracker::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum; + double xtmp,ytmp,ztmp,delx,dely,delz; + double radi,radj,radsum,rsq,r; + int *ilist,*jlist,*numneigh,**firstneigh; + int *touch,**firsttouch; + double *data,*alldata,**firstdata; + + int updateflag = 1; + if (update->setupflag) updateflag = 0; + + + double **x = atom->x; + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + firsttouch = fix_history->firstflag; + firstdata = fix_history->firstvalue; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + if (finitecutflag) radi = radius[i]; + touch = firsttouch[i]; + alldata = firstdata[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + if (finitecutflag) { + radj = radius[j]; + radsum = radi + radj; + + if (rsq >= radsum*radsum) { + // unset non-touching neighbors + if(touch[jj] == 1) { + fix_pair_tracker->lost_contact(i, j, data[0], data[1], data[2]); + } + touch[jj] = 0; + data = &alldata[size_history*jj]; + data[0] = 0.0; // initial time + data[1] = 0.0; // sum of r + data[2] = 0.0; // min of r + + } else { + + data = &alldata[size_history*jj]; + if (touch[jj] == 0) { + data[0] = update->ntimestep; + data[1] = 0.0; + data[2] = 0.0; + } + touch[jj] = 1; + + if (updateflag) { + r = sqrt(rsq); + data[1] += r; + if(data[2] > r) data[2] = r; + } + } + } else { + + } + + if (rsq >= radsum*radsum) { + + + f[i][0] += fx; + f[i][1] += fy; + f[i][2] += fz; + + + if (newton_pair || j < nlocal) { + f[j][0] -= fx; + f[j][1] -= fy; + f[j][2] -= fz; + torque[j][0] -= radj*tor1; + torque[j][1] -= radj*tor2; + torque[j][2] -= radj*tor3; + } + + } + } + } +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairTracker::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cut,n+1,n+1,"pair:cut"); + + onerad_dynamic = new double[n+1]; + onerad_frozen = new double[n+1]; + maxrad_dynamic = new double[n+1]; + maxrad_frozen = new double[n+1]; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairTracker::settings(int narg, char **arg) +{ + if (narg != 0 || narg != 1) error->all(FLERR,"Illegal pair_style command"); + + if (narg == 1) { + if (strcmp(arg[0], "radius")) finitecutflag = 1; + else error->all(FLERR,"Illegal pair_style command"); + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairTracker::coeff(int narg, char **arg) +{ + if (narg > 2 && finitecutflag) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg != 3 && !finitecutflag) error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); + utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + + double cut_one = 0.0; + if (!finitecutflag) cut_one = utils::numeric(FLERR,arg[3],false,lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + setflag[i][j] = 1; + cut[i][j] = cut_one; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairTracker::init_style() +{ + int i; + + // error and warning checks + + if (!atom->radius_flag && finitecutflag) + error->all(FLERR,"Pair tracker requires atom attribute radius for finite cutoffs"); + + // need a history neigh list + + int irequest = neighbor->request(this,instance_me); + if (finitecutflag) neighbor->requests[irequest]->size = 1; + neighbor->requests[irequest]->history = 1; + + // if history is stored and first init, create Fix to store history + // it replaces FixDummy, created in the constructor + // this is so its order in the fix list is preserved + + if (fix_history == nullptr) { + char dnumstr[16]; + sprintf(dnumstr,"%d",size_history); + char **fixarg = new char*[4]; + fixarg[0] = (char *) "NEIGH_HISTORY_TRACK"; + fixarg[1] = (char *) "all"; + fixarg[2] = (char *) "NEIGH_HISTORY"; + fixarg[3] = dnumstr; + modify->replace_fix("NEIGH_HISTORY_TRACK_DUMMY",4,fixarg,1); + delete [] fixarg; + int ifix = modify->find_fix("NEIGH_HISTORY_TRACK"); + fix_history = (FixNeigTRACKistory *) modify->fix[ifix]; + fix_history->pair = this; + } + + if (finitecutflag) { + // set maxrad_dynamic and maxrad_frozen for each type + // include future FixPour and FixDeposit particles as dynamic + + int itype; + for (i = 1; i <= atom->ntypes; i++) { + onerad_dynamic[i] = onerad_frozen[i] = 0.0; + if (ipour >= 0) { + itype = i; + onerad_dynamic[i] = + *((double *) modify->fix[ipour]->extract("radius",itype)); + } + if (idep >= 0) { + itype = i; + onerad_dynamic[i] = + *((double *) modify->fix[idep]->extract("radius",itype)); + } + } + + double *radius = atom->radius; + int *mask = atom->mask; + int *type = atom->type; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & freeze_group_bit) + onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius[i]); + else + onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]); + + MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes, + MPI_DOUBLE,MPI_MAX,world); + MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes, + MPI_DOUBLE,MPI_MAX,world); + } + + int ifix = modify->find_fix("NEIGH_HISTORY_TRACK"); + if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID"); + fix_history = (FixNeigTRACKistory *) modify->fix[ifix]; + + ifix = modify->find_fix_by_style("pair/tracker"); + if(ifix < 0) error->all(FLERR,"Cannot use pair tracker without fix pair/tracker"); + fix_pair_tracker = (FixPairTracker *) modify->fix[ifix]; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairTracker::init_one(int i, int j) +{ + if (!allocated) allocate(); + + // always mix prefactors geometrically + + if (setflag[i][j] == 0) { + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + cut[j][i] = cut[i][j]; + + double cutoff; + + // if finite, cutoff = sum of max I,J radii for + // dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen + if (finitecutflag) { + cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; + cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); + cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); + } else { + cutoff = cut[i][j]; + } + return cutoff; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairTracker::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairTracker::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); + } + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairTracker::write_restart_settings(FILE *fp) +{ + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairTracker::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); + } + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +double PairTracker::single(int i, int j, int /*itype*/, int /*jtype*/, + double rsq, + double /*factor_coul*/, double /*factor_lj*/, + double &fforce) +{ + return 0.0; +} + +/* ---------------------------------------------------------------------- + transfer history during fix/neigh/history exchange + only needed if any history entries i-j are not just negative of j-i entries +------------------------------------------------------------------------- */ + +void PairTracker::transfer_history(double* source, double* target) +{ + for (int i = 0; i < size_history; i++) + target[i] = source[i]; +} + + +/* ---------------------------------------------------------------------- + self-interaction range of particle if finite particles +------------------------------------------------------------------------- */ + +double PairTracker::atom2cut(int i) +{ + double cut = atom->radius[i]*2; + return cut; +} + +/* ---------------------------------------------------------------------- + maximum interaction range for two finite particles +------------------------------------------------------------------------- */ + +double PairTracker::radii2cut(double r1, double r2) +{ + double cut = r1+r2; + return cut; +} + diff --git a/src/pair_tracker.h b/src/pair_tracker.h new file mode 100644 index 0000000000..43dcf54559 --- /dev/null +++ b/src/pair_tracker.h @@ -0,0 +1,101 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://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 PAIR_CLASS + +PairStyle(tracker,PairTracker) + +#else + +#ifndef LMP_PAIR_TRACKER_H +#define LMP_PAIR_TRACKER_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairTracker : public Pair { + public: + PairTracker(class LAMMPS *); + virtual ~PairTracker(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + void init_style(); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + void reset_dt(); + virtual double single(int, int, int, int, double, double, double, double &); + double atom2cut(int); + double radii2cut(double,double); + + protected: + int sizeflag; + int history; + int size_history; + int neighprev; + double **cut; + double *onerad_dynamic,*onerad_frozen; + double *maxrad_dynamic,*maxrad_frozen; + + class FixDummy *fix_dummy; + class FixNeighHistory *fix_history; + class FixPairTracker *fix_pair_tracker; + + void transfer_history(double*, double*); + void allocate(); +}; + +} + +#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: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair granular requires atom attributes radius, rmass + +The atom style defined does not have these attributes. + +E: Pair granular requires ghost atoms store velocity + +Use the comm_modify vel yes command to enable this. + +E: Could not find pair fix neigh history ID + +UNDOCUMENTED + +U: Pair granular with shear history requires newton pair off + +This is a current restriction of the implementation of pair +granular styles with history. + +U: Could not find pair fix ID + +A fix is created internally by the pair style to store shear +history information. You cannot delete it. + +*/