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

update to master

Browse Source
This commit is contained in:
jtclemm
2021-06-02 12:15:41 -06:00
parent 6f00e98945
commit 96ac2dc9f6
4 changed files with 1048 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

376
src/fix_pair_tracker.cpp Normal file
View File

@ -0,0 +1,376 @@
/* ----------------------------------------------------------------------
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 "fix_pair_tracker.h"
#include "atom.h"
#include "atom_vec.h"
#include "error.h"
#include "group.h"
#include "memory.h"
#include "modify.h"
#include "tokenizer.h"
#include "update.h"
#include <string.h>
using namespace LAMMPS_NS;
using namespace FixConst;
#define DELTA 1000
/* ---------------------------------------------------------------------- */
FixPairTracker::FixPairTracker(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");
local_flag = 1;
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery <= 0) error->all(FLERR,"Illegal fix pair/tracker command");
local_freq = nevery;
// If optional arguments included, this will be oversized
nvalues = narg - 4;
pack_choice = new FnPtrPack[nvalues];
tmin = -1;
type_filter = nullptr;
int iarg = 4;
nvalues = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"id1") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_id1;
} else if (strcmp(arg[iarg],"id2") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_id2;
} else if (strcmp(arg[iarg],"time/created") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_time_created;
} else if (strcmp(arg[iarg],"time/broken") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_time_broken;
} else if (strcmp(arg[iarg],"time/total") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_time_total;
} else if (strcmp(arg[iarg],"x") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_x;
} else if (strcmp(arg[iarg],"y") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_y;
} else if (strcmp(arg[iarg],"z") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_z;
} else if (strcmp(arg[iarg],"r/min") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_rmin;
} else if (strcmp(arg[iarg],"r/ave") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_rave;
} else if (strcmp(arg[iarg],"time/min") == 0) {
if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in fix pair/tracker command");
tmin = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg ++;
} else if (strcmp(arg[iarg],"type/include") == 0) {
if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in fix pair/tracker command");
int ntypes = atom->ntypes;
int i, j, itype, jtype, in, jn, infield, jnfield;
int inlo, inhi, jnlo, jnhi;
char *istr, *jstr;
if (! type_filter) {
memory->create(type_filter,ntypes+1,ntypes+1,"fix/pair/tracker:type_filter");
for (i = 0; i <= ntypes; i ++) {
for (j = 0; j <= ntypes; j ++) {
type_filter[i][j] = 0;
}
}
}
in = strlen(arg[iarg+1]) + 1;
istr = new char[in];
strcpy(istr,arg[iarg+1]);
std::vector<std::string> iwords = Tokenizer(istr, ",").as_vector();
infield = iwords.size();
jn = strlen(arg[iarg+2]) + 1;
jstr = new char[jn];
strcpy(jstr,arg[iarg+2]);
std::vector<std::string> jwords = Tokenizer(jstr, ",").as_vector();
jnfield = jwords.size();
for (i = 0; i < infield; i++) {
const char * ifield = iwords[i].c_str();
utils::bounds(FLERR,ifield,1,ntypes,inlo,inhi,error);
for (j = 0; j < jnfield; j++) {
const char * jfield = jwords[j].c_str();
utils::bounds(FLERR,jfield,1,ntypes,jnlo,jnhi,error);
for (itype = inlo; itype <= inhi; itype++) {
for (jtype = jnlo; jtype <= jnhi; jtype++) {
type_filter[itype][jtype] = 1;
type_filter[jtype][itype] = 1;
}
}
}
}
delete [] istr;
delete [] jstr;
iarg += 2;
} else error->all(FLERR, "Invalid keyword in fix pair/tracker command");
iarg ++;
}
if (nvalues == 1) size_local_cols = 0;
else size_local_cols = nvalues;
nmax = 0;
ncount = 0;
vector = NULL;
array = NULL;
}
/* ---------------------------------------------------------------------- */
FixPairTracker::~FixPairTracker()
{
delete [] pack_choice;
memory->destroy(vector);
memory->destroy(array);
memory->destroy(type_filter);
}
/* ---------------------------------------------------------------------- */
int FixPairTracker::setmask()
{
int mask = 0;
mask |= POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::init()
{
// Set size of array/vector
ncount = 0;
if (ncount > nmax)
reallocate(ncount);
size_local_rows = ncount;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::lost_contact(int i, int j, double time_tmp, double nstep_tmp, double rsum_tmp, double rmin_tmp)
{
double time = update->atime + (update->ntimestep-update->atimestep)*update->dt;
if ((time-time_tmp) < tmin) return;
if (type_filter) {
int *type = atom->type;
if (type_filter[type[i]][type[j]] == 0) return;
}
int *mask = atom->mask;
if (!(mask[i] & groupbit)) return;
if (!(mask[j] & groupbit)) return;
if (ncount == nmax) reallocate(ncount);
index_i = i;
index_j = j;
rmin = rmin_tmp;
rsum = rsum_tmp;
time_initial = time_tmp;
nstep_initial = nstep_tmp;
// fill vector or array with local values
if (nvalues == 1) {
(this->*pack_choice[0])(0);
} else {
for (int k = 0; k < nvalues; k++) {
(this->*pack_choice[k])(k);
}
}
ncount += 1;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::post_force(int /*vflag*/)
{
if (update->ntimestep % nevery == 0) {
size_local_rows = ncount;
ncount = 0;
}
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::reallocate(int n)
{
// grow vector or array
while (nmax <= n) nmax += DELTA;
if (nvalues == 1) {
memory->grow(vector,nmax,"fix_pair_tracker:vector");
vector_local = vector;
} else {
memory->grow(array,nmax,nvalues,"fix_pair_tracker:array");
array_local = array;
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double FixPairTracker::memory_usage()
{
double bytes = nmax*nvalues * sizeof(double);
bytes += nmax*2 * sizeof(int);
return bytes;
}
/* ----------------------------------------------------------------------
one method for every keyword fix pair/tracker can output
the atom property is packed into the local vector or array
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_time_created(int n)
{
if (nvalues == 1)
vector[ncount] = time_initial;
else
array[ncount][n] = time_initial;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_time_broken(int n)
{
double time = update->atime + (update->ntimestep-update->atimestep)*update->dt;
if (nvalues == 1)
vector[ncount] = time;
else
array[ncount][n] = time;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_time_total(int n)
{
double time = update->atime + (update->ntimestep-update->atimestep)*update->dt;
if (nvalues == 1)
vector[ncount] = time-time_initial;
else
array[ncount][n] = time-time_initial;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_id1(int n)
{
tagint *tag = atom->tag;
if (nvalues == 1)
vector[ncount] = tag[index_i];
else
array[ncount][n] = tag[index_i];
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_id2(int n)
{
tagint *tag = atom->tag;
if (nvalues == 1)
vector[ncount] = tag[index_j];
else
array[ncount][n] = tag[index_j];
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_x(int n)
{
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 FixPairTracker::pack_y(int n)
{
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 FixPairTracker::pack_z(int n)
{
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 FixPairTracker::pack_rmin(int n)
{
if (nvalues == 1)
vector[ncount] = rmin;
else
array[ncount][n] = rmin;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_rave(int n)
{
if (nvalues == 1)
vector[ncount] = rsum/(update->ntimestep-nstep_initial);
else
array[ncount][n] = rsum/(update->ntimestep-nstep_initial);
}

92
src/fix_pair_tracker.h Normal file
View File

@ -0,0 +1,92 @@
/* ----------------------------------------------------------------------
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/tracker,FixPairTracker)
#else
#ifndef LMP_FIX_PAIR_TRACKING_H
#define LMP_FIX_PAIR_TRACKING_H
#include "fix.h"
namespace LAMMPS_NS {
class FixPairTracker : public Fix {
public:
FixPairTracker(class LAMMPS *, int, char **);
~FixPairTracker();
int setmask();
void init();
void post_force(int);
double memory_usage();
void lost_contact(int, int, double, double, double, double);
private:
int nvalues, never;
int nmax;
int store_flag;
int index_i, index_j;
double tmin, rmin, rsum, time_initial, nstep_initial;
double *vector;
double **array;
int **type_filter;
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_rmin(int);
void pack_rave(int);
};
}
#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: Invalid keyword in fix pair/tracker command
Self-explanatory.
*/

489
src/pair_tracker.cpp Normal file
View File

@ -0,0 +1,489 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "pair_tracker.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "fix.h"
#include "fix_dummy.h"
#include "fix_neigh_history.h"
#include "fix_pair_tracker.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 1;
no_virial_fdotr_compute = 1;
neighprev = 0;
history = 1;
size_history = 4;
nondefault_history_transfer = 1;
finitecutflag = 0;
// create dummy fix as placeholder for FixNeighHistory
// 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,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,time;
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;
ev_init(eflag,vflag);
double **x = atom->x;
double *radius = atom->radius;
int *type = atom->type;
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];
itype = type[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;
r = sqrt(rsq);
if (finitecutflag) {
radj = radius[j];
radsum = radi + radj;
if (rsq >= radsum*radsum) {
data = &alldata[size_history*jj];
if(touch[jj] == 1) {
fix_pair_tracker->lost_contact(i, j, data[0], data[1], data[2], data[3]);
}
touch[jj] = 0;
data[0] = 0.0; // initial time
data[1] = 0.0; // initial timestep
data[2] = 0.0; // sum of r, may overflow
data[3] = 0.0; // min of r
} else {
data = &alldata[size_history*jj];
if (touch[jj] == 0) {
time = update->atime + (update->ntimestep-update->atimestep)*update->dt;
data[0] = time;
data[1] = (double) update->ntimestep;
data[2] = r;
data[3] = r;
} else if (updateflag) {
data[2] += r;
if(data[3] > r) data[3] = r;
}
touch[jj] = 1;
}
} else {
jtype = type[j];
if (rsq >= cutsq[itype][jtype]) {
data = &alldata[size_history*jj];
if(touch[jj] == 1) {
fix_pair_tracker->lost_contact(i, j, data[0], data[1], data[2], data[3]);
}
touch[jj] = 0;
data[0] = 0.0; // initial time
data[1] = 0.0; // initial timestep
data[2] = 0.0; // sum of r, may overflow
data[3] = 0.0; // min of r
} else {
data = &alldata[size_history*jj];
if (touch[jj] == 0) {
time = update->atime + (update->ntimestep-update->atimestep)*update->dt;
data[0] = time;
data[1] = (double) update->ntimestep;
data[2] = r;
data[3] = r;
} else if (updateflag) {
data[2] += r;
if(data[3] > r) data[3] = r;
}
touch[jj] = 1;
}
}
}
}
}
/* ----------------------------------------------------------------------
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], "finite") == 0) 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[2],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 = (FixNeighHistory *) modify->fix[ifix];
fix_history->pair = this;
fix_history->use_bit_flag = 0;
}
if (finitecutflag) {
if(force->pair->beyond_contact)
error->all(FLERR, "Pair tracker incompatible with granular pairstyles that extend beyond contact");
// check for FixPour and FixDeposit so can extract particle radii
int ipour;
for (ipour = 0; ipour < modify->nfix; ipour++)
if (strcmp(modify->fix[ipour]->style,"pour") == 0) break;
if (ipour == modify->nfix) ipour = -1;
int idep;
for (idep = 0; idep < modify->nfix; idep++)
if (strcmp(modify->fix[idep]->style,"deposit") == 0) break;
if (idep == modify->nfix) idep = -1;
// 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 = (FixNeighHistory *) 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];
// if finite, cutoff = sum of max I,J radii for
// dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen
double cutoff;
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;
}

91
src/pair_tracker.h Normal file
View File

@ -0,0 +1,91 @@
/* -*- 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 *);
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;
int freeze_group_bit;
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 tracker requires atom attribute radius for finite cutoffs
The atom style defined does not have these attributes.
E: Could not find pair fix neigh history ID
The associated fix neigh/history is missing
E: Cannot use pair tracker without fix pair/tracker
This pairstyle requires one to define a pair/tracker fix
*/