Completing pair and fix

This commit is contained in:
Joel Clemmer
2021-05-14 17:26:51 -06:00
parent fdf12c4abc
commit d939a16d3b
7 changed files with 149 additions and 236 deletions

View File

@ -83,6 +83,7 @@ Dump::Dump(LAMMPS *lmp, int /*narg*/, char **arg) : Pointers(lmp)
unit_flag = 0;
unit_count = 0;
delay_flag = 0;
write_header_flag = 1;
maxfiles = -1;
numfiles = 0;
@ -377,7 +378,7 @@ void Dump::write()
if (multiproc)
MPI_Allreduce(&bnme,&nheader,1,MPI_LMP_BIGINT,MPI_SUM,clustercomm);
if (filewriter) write_header(nheader);
if (filewriter && write_header_flag) write_header(nheader);
// insure buf is sized for packing and communicating
// use nmax to insure filewriter proc can receive info from others
@ -931,6 +932,13 @@ void Dump::modify_params(int narg, char **arg)
else delay_flag = 0;
iarg += 2;
} else if (strcmp(arg[iarg],"header") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) write_header_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) write_header_flag = 0;
else error->all(FLERR,"Illegal dump_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"every") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
int idump;

View File

@ -79,6 +79,7 @@ class Dump : protected Pointers {
int unit_flag; // 1 if dump should contain unit information
int unit_count; // # of times the unit information was written
int delay_flag; // 1 if delay output until delaystep
int write_header_flag; // 1 if write header, 0 if not
bigint delaystep;
int refreshflag; // 1 if dump_modify refresh specified

View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
#include <string.h>
#include "fix_pair_tracking.h"
#include "fix_pair_tracker.h"
#include "pair_tracker.h"
#include "atom.h"
#include "atom_vec.h"
@ -34,15 +34,15 @@ using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixPairTracking::FixPairTracking(LAMMPS *lmp, int narg, char **arg) :
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");
store_flag = 0;
local_flag = 1;
nvalues = narg - 4;
tmin = -1;
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery <= 0) error->all(FLERR,"Illegal fix pair/tracker command");
@ -52,45 +52,45 @@ FixPairTracking::FixPairTracking(LAMMPS *lmp, int narg, char **arg) :
pack_choice = new FnPtrPack[nvalues];
int i;
for (int iarg = 4; iarg < narg; iarg++) {
i = iarg - 4;
int iarg = 4;
int i = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"id1") == 0) {
pack_choice[i] = &FixPairTracking::pack_id1;
pack_choice[i] = &FixPairTracker::pack_id1;
} else if (strcmp(arg[iarg],"id2") == 0) {
pack_choice[i] = &FixPairTracking::pack_id2;
pack_choice[i] = &FixPairTracker::pack_id2;
} else if (strcmp(arg[iarg],"time/created") == 0) {
pack_choice[i] = &FixPairTracking::pack_time_created;
pack_choice[i] = &FixPairTracker::pack_time_created;
} else if (strcmp(arg[iarg],"time/broken") == 0) {
pack_choice[i] = &FixPairTracking::pack_time_broken;
pack_choice[i] = &FixPairTracker::pack_time_broken;
} else if (strcmp(arg[iarg],"time/total") == 0) {
pack_choice[i] = &FixPairTracking::pack_time_total;
pack_choice[i] = &FixPairTracker::pack_time_total;
} else if (strcmp(arg[iarg],"x") == 0) {
pack_choice[i] = &FixPairTracking::pack_x;
pack_choice[i] = &FixPairTracker::pack_x;
} else if (strcmp(arg[iarg],"y") == 0) {
pack_choice[i] = &FixPairTracking::pack_y;
pack_choice[i] = &FixPairTracker::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;
pack_choice[i] = &FixPairTracker::pack_z;
} else if (strcmp(arg[iarg],"rmin") == 0) {
pack_choice[i] = &FixPairTracking::pack_rmin;
pack_choice[i] = &FixPairTracker::pack_rmin;
} else if (strcmp(arg[iarg],"rave") == 0) {
pack_choice[i] = &FixPairTracking::pack_rave;
pack_choice[i] = &FixPairTracker::pack_rave;
} else if (strcmp(arg[iarg],"tmin") == 0) {
if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in fix pair/tracker command");
tmin = utils::numeric(FLERR,arg[iarg+1],false,lmp);
i -= 1;
nvalues -= 2;
iarg ++;
} else error->all(FLERR, "Invalid keyword in fix pair/tracker command");
iarg ++;
i ++;
}
nmax = 0;
@ -101,13 +101,8 @@ FixPairTracking::FixPairTracking(LAMMPS *lmp, int narg, char **arg) :
/* ---------------------------------------------------------------------- */
FixPairTracking::~FixPairTracking()
FixPairTracker::~FixPairTracker()
{
if (modify->nfix & store_flag == 1) {
modify->delete_fix(id_fix);
delete [] id_fix;
}
delete [] pack_choice;
memory->destroy(vector);
@ -116,7 +111,7 @@ FixPairTracking::~FixPairTracking()
/* ---------------------------------------------------------------------- */
int FixPairTracking::setmask()
int FixPairTracker::setmask()
{
int mask = 0;
mask |= POST_FORCE;
@ -125,93 +120,7 @@ int FixPairTracking::setmask()
/* ---------------------------------------------------------------------- */
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()
void FixPairTracker::init()
{
// Set size of array/vector
ncount = 0;
@ -223,30 +132,34 @@ void FixPairTracking::init()
/* ---------------------------------------------------------------------- */
void FixPairTracking::lost_contact(int i, int j, double n, double rs, double rm)
void FixPairTracker::lost_contact(int i, int j, double n, double rs, double rm)
{
if (ncount == nmax) reallocate(ncount);
if (update->ntimestep-n > tmin) {
if (ncount == nmax) reallocate(ncount);
index_i = i;
index_j = j;
rmin = rm;
rave = ra;
ntimestep = n;
index_i = i;
index_j = j;
// 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);
rmin = rm;
rsum = rs;
ntimestep = n;
// 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;
}
ncount += 1;
}
/* ---------------------------------------------------------------------- */
void FixPairTracking::post_force(int /*vflag*/)
void FixPairTracker::post_force(int /*vflag*/)
{
if (update->ntimestep % nevery == 0) {
size_local_rows = ncount;
@ -257,7 +170,7 @@ void FixPairTracking::post_force(int /*vflag*/)
/* ---------------------------------------------------------------------- */
void FixPairTracking::reallocate(int n)
void FixPairTracker::reallocate(int n)
{
// grow vector or array
while (nmax <= n) nmax += DELTA;
@ -275,7 +188,7 @@ void FixPairTracking::reallocate(int n)
memory usage of local data
------------------------------------------------------------------------- */
double FixPairTracking::memory_usage()
double FixPairTracker::memory_usage()
{
double bytes = nmax*nvalues * sizeof(double);
bytes += nmax*2 * sizeof(int);
@ -290,7 +203,7 @@ double FixPairTracking::memory_usage()
/* ---------------------------------------------------------------------- */
void FixPairTracking::pack_time_created(int n)
void FixPairTracker::pack_time_created(int n)
{
if (nvalues == 1)
vector[ncount] = ntimestep;
@ -300,7 +213,7 @@ void FixPairTracking::pack_time_created(int n)
/* ---------------------------------------------------------------------- */
void FixPairTracking::pack_time_broken(int n)
void FixPairTracker::pack_time_broken(int n)
{
if (nvalues == 1)
vector[ncount] = update->ntimestep;
@ -310,7 +223,7 @@ void FixPairTracking::pack_time_broken(int n)
/* ---------------------------------------------------------------------- */
void FixPairTracking::pack_time_total(int n)
void FixPairTracker::pack_time_total(int n)
{
if (nvalues == 1)
vector[ncount] = update->ntimestep-ntimestep;
@ -321,7 +234,7 @@ void FixPairTracking::pack_time_total(int n)
/* ---------------------------------------------------------------------- */
void FixPairTracking::pack_id1(int n)
void FixPairTracker::pack_id1(int n)
{
tagint *tag = atom->tag;
@ -333,7 +246,7 @@ void FixPairTracking::pack_id1(int n)
/* ---------------------------------------------------------------------- */
void FixPairTracking::pack_id2(int n)
void FixPairTracker::pack_id2(int n)
{
tagint *tag = atom->tag;
@ -345,7 +258,7 @@ void FixPairTracking::pack_id2(int n)
/* ---------------------------------------------------------------------- */
void FixPairTracking::pack_x(int n)
void FixPairTracker::pack_x(int n)
{
double lx_new = domain->xprd;
double **x = atom->x;
@ -359,7 +272,7 @@ void FixPairTracking::pack_x(int n)
/* ---------------------------------------------------------------------- */
void FixPairTracking::pack_y(int n)
void FixPairTracker::pack_y(int n)
{
double lx_new = domain->yprd;
double **x = atom->x;
@ -373,7 +286,7 @@ void FixPairTracking::pack_y(int n)
/* ---------------------------------------------------------------------- */
void FixPairTracking::pack_z(int n)
void FixPairTracker::pack_z(int n)
{
double lx_new = domain->zprd;
double **x = atom->x;
@ -384,48 +297,9 @@ void FixPairTracking::pack_z(int n)
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)
void FixPairTracker::pack_rmin(int n)
{
if (nvalues == 1)
vector[ncount] = rmin;
@ -435,7 +309,7 @@ void FixPairTracking::pack_rmin(int n)
/* ---------------------------------------------------------------------- */
void FixPairTracking::pack_rave(int n)
void FixPairTracker::pack_rave(int n)
{
if (nvalues == 1)
vector[ncount] = rsum/(update->ntimestep-ntimestep);

View File

@ -13,7 +13,7 @@
#ifdef FIX_CLASS
FixStyle(pair/tracking,FixPairTracking)
FixStyle(pair/tracker,FixPairTracker)
#else
@ -24,12 +24,11 @@ FixStyle(pair/tracking,FixPairTracking)
namespace LAMMPS_NS {
class FixPairTracking : public Fix {
class FixPairTracker : public Fix {
public:
FixPairTracking(class LAMMPS *, int, char **);
~FixPairTracking();
FixPairTracker(class LAMMPS *, int, char **);
~FixPairTracker();
int setmask();
void post_constructor();
void init();
void post_force(int);
double memory_usage();
@ -37,9 +36,10 @@ class FixPairTracking : public Fix {
private:
int nvalues, never;
int nmax;
int nmax, tmin;
int store_flag;
int index_i, index_j;
double rmin, rsum, ntimestep;
double *vector;
double **array;
@ -66,15 +66,9 @@ class FixPairTracking : public Fix {
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;
};
}

View File

@ -580,6 +580,10 @@ void PairHybrid::init_style()
}
}
// check beyond contact (set during pair coeff) before init style
for (istyle = 0; istyle < nstyles; istyle++)
if (styles[istyle]->beyond_contact) beyond_contact = 1;
// each sub-style makes its neighbor list request(s)
for (istyle = 0; istyle < nstyles; istyle++) styles[istyle]->init_style();

View File

@ -27,6 +27,7 @@
#include "fix.h"
#include "fix_dummy.h"
#include "fix_neigh_history.h"
#include "fix_pair_tracker.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
@ -36,6 +37,9 @@
using namespace LAMMPS_NS;
#define RLARGE 100
/* ---------------------------------------------------------------------- */
PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp)
@ -50,7 +54,7 @@ PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp)
finitecutflag = 0;
// create dummy fix as placeholder for FixNeigTRACKistory
// create dummy fix as placeholder for FixNeighHistory
// this is so final order of Modify:fix will conform to input script
fix_history = nullptr;
@ -81,7 +85,7 @@ PairTracker::~PairTracker()
void PairTracker::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum;
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz;
double radi,radj,radsum,rsq,r;
int *ilist,*jlist,*numneigh,**firstneigh;
@ -90,10 +94,11 @@ void PairTracker::compute(int eflag, int vflag)
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;
@ -106,13 +111,13 @@ void PairTracker::compute(int eflag, int vflag)
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];
@ -126,17 +131,18 @@ void PairTracker::compute(int eflag, int vflag)
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
data = &alldata[size_history*jj];
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
@ -147,7 +153,7 @@ void PairTracker::compute(int eflag, int vflag)
if (touch[jj] == 0) {
data[0] = update->ntimestep;
data[1] = 0.0;
data[2] = 0.0;
data[2] = RLARGE;
}
touch[jj] = 1;
@ -158,29 +164,40 @@ void PairTracker::compute(int eflag, int vflag)
}
}
} 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]);
}
touch[jj] = 0;
data[0] = 0.0; // initial time
data[1] = 0.0; // sum of r
data[2] = 0.0; // min of r
if (rsq >= radsum*radsum) {
} else {
data = &alldata[size_history*jj];
if (touch[jj] == 0) {
data[0] = update->ntimestep;
data[1] = 0.0;
data[2] = RLARGE;
}
touch[jj] = 1;
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;
if (updateflag) {
r = sqrt(rsq);
data[1] += r;
if(data[2] > r) data[2] = r;
}
}
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
@ -212,10 +229,10 @@ void PairTracker::allocate()
void PairTracker::settings(int narg, char **arg)
{
if (narg != 0 || narg != 1) error->all(FLERR,"Illegal pair_style command");
if (narg != 0 && narg != 1) error->all(FLERR,"Illegal pair_style command");
if (narg == 1) {
if (strcmp(arg[0], "radius")) finitecutflag = 1;
if (strcmp(arg[0], "radius") == 0) finitecutflag = 1;
else error->all(FLERR,"Illegal pair_style command");
}
}
@ -235,7 +252,7 @@ void PairTracker::coeff(int narg, char **arg)
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);
if (!finitecutflag) cut_one = utils::numeric(FLERR,arg[2],false,lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
@ -283,11 +300,26 @@ void PairTracker::init_style()
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 = (FixNeighHistory *) modify->fix[ifix];
fix_history->pair = this;
}
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
@ -325,7 +357,7 @@ void PairTracker::init_style()
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];
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");

View File

@ -37,7 +37,6 @@ class PairTracker : public Pair {
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);
@ -50,6 +49,7 @@ class PairTracker : public Pair {
double **cut;
double *onerad_dynamic,*onerad_frozen;
double *maxrad_dynamic,*maxrad_frozen;
int freeze_group_bit;
class FixDummy *fix_dummy;
class FixNeighHistory *fix_history;