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

Moving update/special/bonds into bond style to ensure correct fix ordering

Browse Source
This commit is contained in:
Joel Thomas Clemmer
2021-11-07 17:55:53 -07:00
parent f5626a2b9d
commit acb1c8e7f2
12 changed files with 161 additions and 152 deletions
Download Patch File Download Diff File Expand all files Collapse all files

116
src/BPM/bond_bpm.cpp
View File

@ -33,10 +33,11 @@ BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp)
{
id_fix_store_local = nullptr;
id_fix_prop_atom = nullptr;
id_fix_update = nullptr;
fix_store_local = nullptr;
fix_update_special_bonds = nullptr;
fix_bond_history = nullptr;
overlay_flag = 0;
prop_atom_flag = 0;
nvalues = 0;
output_data = nullptr;
@ -44,6 +45,12 @@ BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp)
r0_max_estimate = 0.0;
max_stretch = 1.0;
// create dummy fix as placeholder for FixUpdateSpecialBonds
// this is so final order of Modify:fix will conform to input script
id_fix_dummy = utils::strdup("BPM_DUMMY");
modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy));
}
/* ---------------------------------------------------------------------- */
@ -53,6 +60,13 @@ BondBPM::~BondBPM()
delete [] pack_choice;
delete [] id_fix_store_local;
delete [] id_fix_prop_atom;
if (id_fix_dummy) modify->delete_fix(id_fix_dummy);
if (id_fix_update) modify->delete_fix(id_fix_update);
delete [] id_fix_dummy;
delete [] id_fix_update;
memory->destroy(output_data);
}
@ -70,14 +84,38 @@ void BondBPM::init_style()
fix_store_local->nvalues = nvalues;
}
ifix = modify->find_fix_by_style("update/special/bonds");
if (ifix >= 0)
fix_update_special_bonds = (FixUpdateSpecialBonds *) modify->fix[ifix];
else {
fix_update_special_bonds = nullptr;
if (overlay_flag) {
if (force->special_lj[1] != 1.0)
error->all(FLERR, "Without fix update/special/bonds, BPM bond styles "
error->all(FLERR, "With overlay/pair, BPM bond styles "
"require special_bonds weight of 1.0 for first neighbors");
if (id_fix_update) {
modify->delete_fix(id_fix_update);
delete [] id_fix_update;
id_fix_update = nullptr;
}
} else {
// Require atoms know about all of their bonds and if they break
if (force->newton_bond)
error->all(FLERR,"Without overlay/pair, BPM bond sytles require Newton bond off");
// special lj must be 0 1 1 to censor pair forces between bonded particles
// special coulomb must be 1 1 1 to ensure all pairs are included in the
// neighbor list and 1-3 and 1-4 special bond lists are skipped
if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 ||
force->special_lj[3] != 1.0)
error->all(FLERR,"Without overlay/pair, BPM bond sytles requires special LJ weights = 0,1,1");
if (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 ||
force->special_coul[3] != 1.0)
error->all(FLERR,"Without overlay/pair, BPM bond sytles requires special Coulomb weights = 1,1,1");
if (id_fix_dummy) {
id_fix_update = utils::strdup("BPM_update_special_bonds");
fix_update_special_bonds = (FixUpdateSpecialBonds *) modify->replace_fix(id_fix_dummy,
fmt::format("{} all update/special/bonds", id_fix_update),1);
delete [] id_fix_dummy;
id_fix_dummy = nullptr;
}
}
if (force->angle || force->dihedral || force->improper)
@ -105,38 +143,44 @@ void BondBPM::settings(int narg, char **arg)
{
int iarg = 0;
while (iarg < narg) {
if (id_fix_store_local) {
if (strcmp(arg[iarg], "id1") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_id1;
} else if (strcmp(arg[iarg], "id2") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_id2;
} else if (strcmp(arg[iarg], "time") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_time;
} else if (strcmp(arg[iarg], "x") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_x;
} else if (strcmp(arg[iarg], "y") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_y;
} else if (strcmp(arg[iarg], "z") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_z;
} else if (strcmp(arg[iarg], "x/ref") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_x_ref;
prop_atom_flag = 1;
} else if (strcmp(arg[iarg], "y/ref") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_y_ref;
prop_atom_flag = 1;
} else if (strcmp(arg[iarg], "z/ref") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_z_ref;
prop_atom_flag = 1;
} else {
error->all(FLERR, "Illegal bond_style command");
}
} else if (strcmp(arg[iarg], "store/local") == 0) {
if (strcmp(arg[iarg], "store/local") == 0) {
id_fix_store_local = utils::strdup(arg[iarg+1]);
iarg ++;
nvalues = 0;
pack_choice = new FnPtrPack[narg - iarg - 1];
iarg += 2;
while (iarg < narg) {
if (strcmp(arg[iarg], "id1") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_id1;
} else if (strcmp(arg[iarg], "id2") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_id2;
} else if (strcmp(arg[iarg], "time") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_time;
} else if (strcmp(arg[iarg], "x") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_x;
} else if (strcmp(arg[iarg], "y") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_y;
} else if (strcmp(arg[iarg], "z") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_z;
} else if (strcmp(arg[iarg], "x/ref") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_x_ref;
prop_atom_flag = 1;
} else if (strcmp(arg[iarg], "y/ref") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_y_ref;
prop_atom_flag = 1;
} else if (strcmp(arg[iarg], "z/ref") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_z_ref;
prop_atom_flag = 1;
} else {
break;
}
iarg ++;
}
} else if (strcmp(arg[iarg], "overlay/pair") == 0) {
overlay_flag = 1;
iarg ++;
} else {
error->all(FLERR, "Illegal pair_style command");
}
iarg ++;
}
if (id_fix_store_local) {
@ -147,7 +191,7 @@ void BondBPM::settings(int narg, char **arg)
// Use store property to save reference positions as it can transfer to ghost atoms
if (prop_atom_flag == 1) {
id_fix_prop_atom = utils::strdup("BPM_PROPERTY_ATOM");
id_fix_prop_atom = utils::strdup("BPM_property_atom");
int ifix = modify->find_fix(id_fix_prop_atom);
if (ifix < 0) {
modify->add_fix(fmt::format("{} all property/atom "

13
src/BPM/bond_bpm.h
View File

@ -25,7 +25,7 @@ class BondBPM : public Bond {
virtual void compute(int, int) = 0;
virtual void coeff(int, char **) = 0;
virtual void init_style();
void settings(int, char **);
void settings(int, char **);
double equilibrium_distance(int);
void write_restart(FILE *){};
void read_restart(FILE *){};
@ -35,18 +35,19 @@ class BondBPM : public Bond {
protected:
double r0_max_estimate;
double max_stretch;
char *id_fix_dummy, *id_fix_update;
char *id_fix_store_local, *id_fix_prop_atom;
class FixStoreLocal *fix_store_local;
class FixUpdateSpecialBonds *fix_update_special_bonds;
class FixBondHistory *fix_bond_history;
class FixUpdateSpecialBonds *fix_update_special_bonds;
void process_broken(int, int);
void process_broken(int, int);
typedef void (BondBPM::*FnPtrPack)(int,int,int);
FnPtrPack *pack_choice; // ptrs to pack functions
double *output_data;
int prop_atom_flag, nvalues;
int prop_atom_flag, nvalues, overlay_flag;
int index_x_ref, index_y_ref, index_z_ref;
void pack_id1(int,int,int);
@ -57,7 +58,7 @@ class BondBPM : public Bond {
void pack_z(int,int,int);
void pack_x_ref(int,int,int);
void pack_y_ref(int,int,int);
void pack_z_ref(int,int,int);
void pack_z_ref(int,int,int);
};
} // namespace LAMMPS_NS

32
src/MISC/pair_tracker.cpp
View File

@ -385,33 +385,35 @@ void PairTracker::init_style()
"Pair tracker incompatible with granular pairstyles that extend beyond contact");
// check for FixFreeze and set freeze_group_bit
int ifix = modify->find_fix_by_style("^freeze");
if (ifix < 0) freeze_group_bit = 0;
else freeze_group_bit = modify->fix[ifix]->groupbit;
auto fixlist = modify->get_fix_by_style("^freeze");
if (fixlist.size() == 0)
freeze_group_bit = 0;
else if (fixlist.size() > 1)
error->all(FLERR, "Only one fix freeze command at a time allowed");
else
freeze_group_bit = fixlist.front()->groupbit;
// 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;
auto pours = modify->get_fix_by_style("^pour");
auto deps = modify->get_fix_by_style("^deposit");
// 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) {
for (auto &ipour : pours) {
itype = i;
onerad_dynamic[i] = *((double *) modify->fix[ipour]->extract("radius", itype));
double maxrad = *((double *) ipour->extract("radius", itype));
if (maxrad > 0.0) onerad_dynamic[i] = maxrad;
}
if (idep >= 0) {
for (auto &idep : deps) {
itype = i;
onerad_dynamic[i] = *((double *) modify->fix[idep]->extract("radius", itype));
double maxrad = *((double *) idep->extract("radius", itype));
if (maxrad > 0.0) onerad_dynamic[i] = maxrad;
}
}