diff --git a/doc/src/pair_rheo_solid.rst b/doc/src/pair_rheo_solid.rst index b6ff6d809d..ee86992776 100644 --- a/doc/src/pair_rheo_solid.rst +++ b/doc/src/pair_rheo_solid.rst @@ -21,36 +21,88 @@ Examples Description """"""""""" -pair style... +Style *rheo/solid* is effectively a copy of pair style +:doc:`bpm/spring ` except it only applies forces +between solid RHEO particles, determined by checking the status of +each pair of neighboring particles before calculating forces. + +The style computes pairwise forces with the formula + +.. math:: + + F = k (r - r_c) + +where :math:`k` is a stiffness and :math:`r_c` is the cutoff length. +An additional damping force is also applied to interacting +particles. The force is proportional to the difference in the +normal velocity of particles + +.. math:: + + F_D = - \gamma w (\hat{r} \bullet \vec{v}) + +where :math:`\gamma` is the damping strength, :math:`\hat{r}` is the +radial normal vector, :math:`\vec{v}` is the velocity difference +between the two particles, and :math:`w` is a smoothing factor. +This smoothing factor is constructed such that damping forces go to zero +as particles come out of contact to avoid discontinuities. It is +given by + +.. math:: + + w = 1.0 - \left( \frac{r}{r_c} \right)^8 . + +The following coefficients must be defined for each pair of atom types +via the :doc:`pair_coeff ` command as in the examples +above, or in the data file or restart files read by the +:doc:`read_data ` or :doc:`read_restart ` +commands, or by mixing as described below: + +* :math:`k` (force/distance units) +* :math:`r_c` (distance units) +* :math:`\gamma` (force/velocity units) -* :math:`k` (force/distance units) -* :math:`\sigma` (distance units) -* :math:`\gamma` (force/velocity units) ---------- Mixing, shift, table, tail correction, restart, rRESPA info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" -This style does not support the :doc:`pair_modify ` -shift, table, and tail options. +For atom type pairs I,J and I != J, the A coefficient and cutoff +distance for this pair style can be mixed. A is always mixed via a +*geometric* rule. The cutoff is mixed according to the pair_modify +mix value. The default mix value is *geometric*\ . See the +"pair_modify" command for details. -This style does not write information to :doc:`binary restart files `. Thus, you need to re-specify the pair_style and -pair_coeff commands in an input script that reads a restart file. +This pair style does not support the :doc:`pair_modify ` +shift option, since the pair interaction goes to 0.0 at the cutoff. -This style can only be used via the *pair* keyword of the :doc:`run_style respa ` command. It does not support the *inner*, *middle*, *outer* keywords. +The :doc:`pair_modify ` table and tail options are not +relevant for this pair style. + +This pair style writes its information to :doc:`binary restart files +`, so pair_style and pair_coeff commands do not need to be +specified in an input script that reads a restart file. + +This pair style can only be used via the *pair* keyword of the +:doc:`run_style respa ` command. It does not support the +*inner*, *middle*, *outer* keywords. + +---------- Restrictions """""""""""" -This fix is part of the RHEO package. It is only enabled if -LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +This pair style is part of the BPM package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. Related commands """""""""""""""" :doc:`fix rheo `, -:doc:`pair bpm/spring `, +:doc:`fix rheo/thermal `, +:doc:`pair bpm/spring ` Default """"""" diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index 34554497ad..351cff1420 100644 --- a/src/BPM/bond_bpm.cpp +++ b/src/BPM/bond_bpm.cpp @@ -17,6 +17,7 @@ #include "comm.h" #include "domain.h" #include "error.h" +#include "fix.h" #include "fix_bond_history.h" #include "fix_store_local.h" #include "fix_update_special_bonds.h" @@ -39,10 +40,14 @@ BondBPM::BondBPM(LAMMPS *_lmp) : pack_choice(nullptr), output_data(nullptr) { overlay_flag = 0; + ignore_special_flag = 0; prop_atom_flag = 0; break_flag = 1; nvalues = 0; + nhistory = 0; + update_flag = 0; + r0_max_estimate = 0.0; max_stretch = 1.0; @@ -93,39 +98,46 @@ void BondBPM::init_style() fix_store_local->nvalues = nvalues; } - if (overlay_flag) { - if (force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || - force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) - error->all(FLERR, - "With overlay/pair yes, BPM bond styles require a value of 1.0 for all special_bonds weights"); - 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 && break_flag) - error->all(FLERR, "With overlay/pair no, or break yes, BPM bond styles require Newton bond off"); + if (!ignore_special_flag) { + if (overlay_flag) { + if (force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || + force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) + error->all(FLERR, + "With overlay/pair yes, BPM bond styles require a value of 1.0 for all special_bonds weights"); + 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 && break_flag) + error->all(FLERR, "With overlay/pair no, or break yes, BPM bond styles require Newton bond off"); - // special lj must be 0 1 1 to censor pair forces between bonded particles - if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) - error->all(FLERR, - "With overlay/pair no, BPM bond styles require special LJ weights = 0,1,1"); - // if bonds can break, 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 (break_flag && (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || - force->special_coul[3] != 1.0)) - error->all(FLERR, - "With overlay/pair no, and break yes, BPM bond styles requires special Coulomb weights = 1,1,1"); + // special lj must be 0 1 1 to censor pair forces between bonded particles + if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) + error->all(FLERR, + "With overlay/pair no, BPM bond styles require special LJ weights = 0,1,1"); + // if bonds can break, 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 (break_flag && (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || + force->special_coul[3] != 1.0)) + error->all(FLERR, + "With overlay/pair no, and break yes, BPM bond styles requires special Coulomb weights = 1,1,1"); - if (id_fix_dummy && break_flag) { - id_fix_update = utils::strdup("BPM_UPDATE_SPECIAL_BONDS"); - fix_update_special_bonds = dynamic_cast(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 (id_fix_dummy && break_flag) { + id_fix_update = utils::strdup("BPM_UPDATE_SPECIAL_BONDS"); + fix_update_special_bonds = dynamic_cast(modify->replace_fix( + id_fix_dummy, fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update), 1)); + delete[] id_fix_dummy; + id_fix_dummy = nullptr; + } } + + // special 1-3 and 1-4 weights must be 1 to prevent building 1-3 and 1-4 special bond lists + if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || force->special_coul[2] != 1.0 || + force->special_coul[3] != 1.0) + error->all(FLERR, "Bond style bpm requires 1-3 and 1-4 special weights of 1.0"); } if (force->angle || force->dihedral || force->improper) @@ -133,10 +145,16 @@ void BondBPM::init_style() if (atom->molecular == 2) error->all(FLERR, "Bond style bpm cannot be used with atom style template"); - // special 1-3 and 1-4 weights must be 1 to prevent building 1-3 and 1-4 special bond lists - if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || force->special_coul[2] != 1.0 || - force->special_coul[3] != 1.0) - error->all(FLERR, "Bond style bpm requires 1-3 and 1-4 special weights of 1.0"); + // find all instances of bond history to delete/shift data + // (bond hybrid may create multiple) + histories = modify->get_fix_by_style("BOND_HISTORY"); + n_histories = histories.size(); + + // If a bond type isn't set, must be using bond style hybrid + hybrid_flag = 0; + for (int i = 1; i <= atom->nbondtypes; i++) + if (!setflag[i]) hybrid_flag = 1; + fix_bond_history->setflag = setflag; } /* ---------------------------------------------------------------------- @@ -253,6 +271,14 @@ void BondBPM::settings(int narg, char **arg) } } } + + // Set up necessary history fix + if (!fix_bond_history) { + fix_bond_history = dynamic_cast(modify->replace_fix( + id_fix_dummy2, fmt::format("{} all BOND_HISTORY {} {}", id_fix_bond_history, update_flag, nhistory), 1)); + delete[] id_fix_dummy2; + id_fix_dummy2 = nullptr; + } } /* ---------------------------------------------------------------------- @@ -356,13 +382,16 @@ void BondBPM::process_broken(int i, int j) if (i < nlocal) { for (m = 0; m < num_bond[i]; m++) { - if (bond_atom[i][m] == tag[j]) { + if (bond_atom[i][m] == tag[j] && setflag[bond_type[i][m]]) { bond_type[i][m] = 0; n = num_bond[i]; bond_type[i][m] = bond_type[i][n - 1]; bond_atom[i][m] = bond_atom[i][n - 1]; - fix_bond_history->shift_history(i, m, n - 1); - fix_bond_history->delete_history(i, n - 1); + for (auto &ihistory: histories) { + auto fix_bond_history2 = dynamic_cast (ihistory); + fix_bond_history2->shift_history(i, m, n - 1); + fix_bond_history2->delete_history(i, n - 1); + } num_bond[i]--; break; } @@ -371,13 +400,16 @@ void BondBPM::process_broken(int i, int j) if (j < nlocal) { for (m = 0; m < num_bond[j]; m++) { - if (bond_atom[j][m] == tag[i]) { + if (bond_atom[j][m] == tag[i] && setflag[bond_type[j][m]]) { bond_type[j][m] = 0; n = num_bond[j]; bond_type[j][m] = bond_type[j][n - 1]; bond_atom[j][m] = bond_atom[j][n - 1]; - fix_bond_history->shift_history(j, m, n - 1); - fix_bond_history->delete_history(j, n - 1); + for (auto &ihistory: histories) { + auto fix_bond_history2 = dynamic_cast (ihistory); + fix_bond_history2->shift_history(j, m, n - 1); + fix_bond_history2->delete_history(j, n - 1); + } num_bond[j]--; break; } diff --git a/src/BPM/bond_bpm.h b/src/BPM/bond_bpm.h index 815b3b751f..28e4e7187e 100644 --- a/src/BPM/bond_bpm.h +++ b/src/BPM/bond_bpm.h @@ -16,8 +16,12 @@ #include "bond.h" +#include + namespace LAMMPS_NS { +class Fix; + class BondBPM : public Bond { public: BondBPM(class LAMMPS *); @@ -34,7 +38,7 @@ class BondBPM : public Bond { protected: double r0_max_estimate; double max_stretch; - int store_local_freq; + int store_local_freq, nhistory, update_flag, hybrid_flag; std::vector leftover_iarg; @@ -50,9 +54,12 @@ class BondBPM : public Bond { FnPtrPack *pack_choice; // ptrs to pack functions double *output_data; - int prop_atom_flag, nvalues, overlay_flag, break_flag; + int prop_atom_flag, nvalues, overlay_flag, break_flag, ignore_special_flag; int index_x_ref, index_y_ref, index_z_ref; + int n_histories; + std::vector histories; + void pack_id1(int, int, int); void pack_id2(int, int, int); void pack_time(int, int, int); diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp index ffb0d9521d..b904a2ac07 100644 --- a/src/BPM/bond_bpm_rotational.cpp +++ b/src/BPM/bond_bpm_rotational.cpp @@ -52,6 +52,9 @@ BondBPMRotational::BondBPMRotational(LAMMPS *_lmp) : smooth_flag = 1; normalize_flag = 0; + nhistory = 4; + id_fix_bond_history = utils::strdup("HISTORY_BPM_ROTATIONAL"); + single_extra = 7; svector = new double[7]; } @@ -458,6 +461,9 @@ void BondBPMRotational::compute(int eflag, int vflag) store_data(); } + if (hybrid_flag) + fix_bond_history->compress_history(); + int i1, i2, itmp, n, type; double r[3], r0[3], rhat[3]; double rsq, r0_mag, r_mag, r_mag_inv; @@ -563,6 +569,9 @@ void BondBPMRotational::compute(int eflag, int vflag) ev_tally_xyz(i1, i2, nlocal, newton_bond, 0.0, -force1on2[0] * smooth, -force1on2[1] * smooth, -force1on2[2] * smooth, r[0], r[1], r[2]); } + + if (hybrid_flag) + fix_bond_history->uncompress_history(); } /* ---------------------------------------------------------------------- */ @@ -652,14 +661,6 @@ void BondBPMRotational::init_style() if (domain->dimension == 2) error->warning(FLERR, "Bond style bpm/rotational not intended for 2d use"); - - if (!id_fix_bond_history) { - id_fix_bond_history = utils::strdup("HISTORY_BPM_ROTATIONAL"); - fix_bond_history = dynamic_cast(modify->replace_fix( - id_fix_dummy2, fmt::format("{} all BOND_HISTORY 0 4", id_fix_bond_history), 1)); - delete[] id_fix_dummy2; - id_fix_dummy2 = nullptr; - } } /* ---------------------------------------------------------------------- */ diff --git a/src/BPM/bond_bpm_spring.cpp b/src/BPM/bond_bpm_spring.cpp index 37b79f93fb..0c882ae6c6 100644 --- a/src/BPM/bond_bpm_spring.cpp +++ b/src/BPM/bond_bpm_spring.cpp @@ -23,6 +23,8 @@ #include "modify.h" #include "neighbor.h" +#include "update.h" + #include #include @@ -39,6 +41,9 @@ BondBPMSpring::BondBPMSpring(LAMMPS *_lmp) : smooth_flag = 1; normalize_flag = 0; + nhistory = 1; + id_fix_bond_history = utils::strdup("HISTORY_BPM_SPRING"); + single_extra = 1; svector = new double[1]; } @@ -137,6 +142,9 @@ void BondBPMSpring::compute(int eflag, int vflag) store_data(); } + if (hybrid_flag) + fix_bond_history->compress_history(); + int i1, i2, itmp, n, type; double delx, dely, delz, delvx, delvy, delvz; double e, rsq, r, r0, rinv, smooth, fbond, dot; @@ -226,6 +234,9 @@ void BondBPMSpring::compute(int eflag, int vflag) if (evflag) ev_tally(i1, i2, nlocal, newton_bond, 0.0, fbond, delx, dely, delz); } + + if (hybrid_flag) + fix_bond_history->uncompress_history(); } /* ---------------------------------------------------------------------- */ @@ -283,14 +294,6 @@ void BondBPMSpring::init_style() if (comm->ghost_velocity == 0) error->all(FLERR, "Bond bpm/spring requires ghost atoms store velocity"); - - if (!id_fix_bond_history) { - id_fix_bond_history = utils::strdup("HISTORY_BPM_SPRING"); - fix_bond_history = dynamic_cast(modify->replace_fix( - id_fix_dummy2, fmt::format("{} all BOND_HISTORY 0 1", id_fix_bond_history), 1)); - delete[] id_fix_dummy2; - id_fix_dummy2 = nullptr; - } } /* ---------------------------------------------------------------------- */ diff --git a/src/BPM/compute_nbond_atom.cpp b/src/BPM/compute_nbond_atom.cpp index af1cc2b9bc..b6e7b0139f 100644 --- a/src/BPM/compute_nbond_atom.cpp +++ b/src/BPM/compute_nbond_atom.cpp @@ -14,7 +14,9 @@ #include "compute_nbond_atom.h" #include "atom.h" +#include "atom_vec.h" #include "comm.h" +#include "error.h" #include "force.h" #include "memory.h" @@ -25,18 +27,20 @@ using namespace LAMMPS_NS; ComputeNBondAtom::ComputeNBondAtom(LAMMPS *_lmp, int narg, char **arg) : Compute(_lmp, narg, arg), nbond(nullptr) { - if (narg < 4) utils::missing_cmd_args(FLERR, "compute nbond/atom", error); + if (narg < 3) utils::missing_cmd_args(FLERR, "compute nbond/atom", error); + + if (atom->avec->bonds_allow == 0) + error->all(FLERR,"Compute nbond/atom used when bonds are not allowed"); btype = -1; - - iarg = 3; + int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg], "bond/type") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "compute nbond/atom bond/type", error); btype = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; } else { - error->all(FLERR, "Unknown compute nbond/type command {}", arg[iarg]); + error->all(FLERR, "Unknown compute nbond/atom command {}", arg[iarg]); } } diff --git a/src/RHEO/bond_rheo_shell.cpp b/src/RHEO/bond_rheo_shell.cpp index d366163459..0a3caa1b4f 100644 --- a/src/RHEO/bond_rheo_shell.cpp +++ b/src/RHEO/bond_rheo_shell.cpp @@ -25,6 +25,7 @@ #include "error.h" #include "fix_bond_history.h" #include "fix_rheo.h" +#include "fix_rheo_oxidation.h" #include "fix_store_local.h" #include "force.h" #include "memory.h" @@ -48,7 +49,12 @@ BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) : partial_flag = 1; comm_reverse = 1; - tform = rmax = -1; + nhistory = 2; + update_flag = 1; + id_fix_bond_history = utils::strdup("HISTORY_RHEO_SHELL"); + ignore_special_flag = 1; + + tform = -1; single_extra = 1; svector = new double[1]; @@ -155,12 +161,14 @@ void BondRHEOShell::store_data() void BondRHEOShell::compute(int eflag, int vflag) { - if (!fix_bond_history->stored_flag) { fix_bond_history->stored_flag = true; store_data(); } + if (hybrid_flag) + fix_bond_history->compress_history(); + int i1, i2, itmp, n, type; double delx, dely, delz, delvx, delvy, delvz; double e, rsq, r, r0, rinv, dr, fbond, dot, t; @@ -168,6 +176,7 @@ void BondRHEOShell::compute(int eflag, int vflag) ev_init(eflag, vflag); + double *rsurface = compute_surface->rsurface; double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -223,7 +232,7 @@ void BondRHEOShell::compute(int eflag, int vflag) if (t < tform) { // Check if eligible - if (r > rmax || !(status[i1] & STATUS_SURFACE) || !(status[i2] & STATUS_SURFACE)) { + if (r > rmax || rsurface[i1] > rsurf || rsurface[i2] > rsurf) { bondlist[n][2] = 0; process_ineligibility(i1, i2); continue; @@ -286,6 +295,9 @@ void BondRHEOShell::compute(int eflag, int vflag) // If it has bonds, no shifting if (nbond[i] != 0) status[i] |= STATUS_NO_SHIFT; } + + if (hybrid_flag) + fix_bond_history->uncompress_history(); } /* ---------------------------------------------------------------------- */ @@ -339,6 +351,8 @@ void BondRHEOShell::coeff(int narg, char **arg) void BondRHEOShell::init_style() { + BondBPM::init_style(); + if (comm->ghost_velocity == 0) error->all(FLERR, "Bond rheo/shell requires ghost atoms store velocity"); @@ -350,37 +364,12 @@ void BondRHEOShell::init_style() "Bond rheo/shell requires surface calculation in fix rheo"); compute_surface = fix_rheo->compute_surface; - if (fix_rheo->oxidation_fix_defined != 1) - error->all(FLERR, "Need to define fix rheo/oxdiation to use bond rheo/shell"); - // check consistency in values (copy?), swap conditions to rsurf + fixes = modify->get_fix_by_style("^rheo/oxidation$"); + if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo/oxidation to use bond rheo/shell"); + class FixRHEOOxidation *fix_rheo_oxidation = dynamic_cast(fixes[0]); - if (!id_fix_bond_history) { - id_fix_bond_history = utils::strdup("HISTORY_RHEO_SHELL"); - fix_bond_history = dynamic_cast(modify->replace_fix( - id_fix_dummy2, fmt::format("{} all BOND_HISTORY 1 2", id_fix_bond_history), 1)); - delete[] id_fix_dummy2; - id_fix_dummy2 = nullptr; - } - - // Reproduce standard functions of BondBPM, removing special restrictions - // Since this bond is intended to be created by fix rheo/oxidation, it - // ignores special statuses - - if (id_fix_store_local) { - auto ifix = modify->get_fix_by_id(id_fix_store_local); - if (!ifix) error->all(FLERR, "Cannot find fix STORE/LOCAL id {}", id_fix_store_local); - if (strcmp(ifix->style, "STORE/LOCAL") != 0) - error->all(FLERR, "Incorrect fix style matched, not STORE/LOCAL: {}", ifix->style); - fix_store_local = dynamic_cast(ifix); - fix_store_local->nvalues = nvalues; - } - - id_fix_update = nullptr; - - if (force->angle || force->dihedral || force->improper) - error->all(FLERR, "Bond style rheo/shell cannot be used with 3,4-body interactions"); - if (atom->molecular == 2) - error->all(FLERR, "Bond style rheo/shell cannot be used with atom style template"); + rsurf = fix_rheo_oxidation->rsurf; + rmax = fix_rheo_oxidation->cut; } /* ---------------------------------------------------------------------- */ @@ -397,20 +386,13 @@ void BondRHEOShell::settings(int narg, char **arg) tform = utils::numeric(FLERR, arg[iarg + 1], false, lmp); if (tform < 0.0) error->all(FLERR, "Illegal bond rheo/shell value for t/form, {}", tform); i += 1; - } else if (strcmp(arg[iarg], "r/max") == 0) { - if (iarg + 1 > narg) error->all(FLERR, "Illegal bond rheo/shell command, missing option for r/max"); - rmax = utils::numeric(FLERR, arg[iarg + 1], false, lmp); - if (rmax < 0.0) error->all(FLERR, "Illegal bond rheo/shell value for r/max, {}", rmax); - i += 1; } else { error->all(FLERR, "Illegal bond rheo/shell command, invalid argument {}", arg[iarg]); } } if (tform < 0.0) - error->all(FLERR, "Illegal bond rheo/shell command, must specify t/form"); - if (rmax < 0.0) - error->all(FLERR, "Illegal bond rheo/shell command, must specify r/max"); + error->all(FLERR, "Illegal bond rheo/shell command, must specify formation time"); } @@ -467,7 +449,6 @@ void BondRHEOShell::read_restart(FILE *fp) void BondRHEOShell::write_restart_settings(FILE *fp) { fwrite(&tform, sizeof(double), 1, fp); - fwrite(&rmax, sizeof(double), 1, fp); } /* ---------------------------------------------------------------------- @@ -478,10 +459,8 @@ void BondRHEOShell::read_restart_settings(FILE *fp) { if (comm->me == 0) { utils::sfread(FLERR, &tform, sizeof(double), 1, fp, nullptr, error); - utils::sfread(FLERR, &rmax, sizeof(double), 1, fp, nullptr, error); } MPI_Bcast(&tform, 1, MPI_DOUBLE, 0, world); - MPI_Bcast(&rmax, 1, MPI_DOUBLE, 0, world); } @@ -570,13 +549,16 @@ void BondRHEOShell::process_ineligibility(int i, int j) if (i < nlocal) { for (m = 0; m < num_bond[i]; m++) { - if (bond_atom[i][m] == tag[j]) { + if (bond_atom[i][m] == tag[j] && setflag[bond_type[i][m]]) { bond_type[i][m] = 0; n = num_bond[i]; bond_type[i][m] = bond_type[i][n - 1]; bond_atom[i][m] = bond_atom[i][n - 1]; - fix_bond_history->shift_history(i, m, n - 1); - fix_bond_history->delete_history(i, n - 1); + for (auto &ihistory: histories) { + auto fix_bond_history2 = dynamic_cast (ihistory); + fix_bond_history2->shift_history(i, m, n - 1); + fix_bond_history2->delete_history(i, n - 1); + } num_bond[i]--; break; } @@ -585,13 +567,16 @@ void BondRHEOShell::process_ineligibility(int i, int j) if (j < nlocal) { for (m = 0; m < num_bond[j]; m++) { - if (bond_atom[j][m] == tag[i]) { + if (bond_atom[j][m] == tag[i] && setflag[bond_type[j][m]]) { bond_type[j][m] = 0; n = num_bond[j]; bond_type[j][m] = bond_type[j][n - 1]; bond_atom[j][m] = bond_atom[j][n - 1]; - fix_bond_history->shift_history(j, m, n - 1); - fix_bond_history->delete_history(j, n - 1); + for (auto &ihistory: histories) { + auto fix_bond_history2 = dynamic_cast (ihistory); + fix_bond_history2->shift_history(j, m, n - 1); + fix_bond_history2->delete_history(j, n - 1); + } num_bond[j]--; break; } diff --git a/src/RHEO/bond_rheo_shell.h b/src/RHEO/bond_rheo_shell.h index a6a747a3f1..828f693ea3 100644 --- a/src/RHEO/bond_rheo_shell.h +++ b/src/RHEO/bond_rheo_shell.h @@ -43,7 +43,7 @@ class BondRHEOShell : public BondBPM { protected: double *k, *ecrit, *gamma; - double tform, rmax; + double tform, rmax, rsurf; int *dbond, *nbond; int index_nb, nmax_store; diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index bdc9951f90..61b83a4542 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -160,27 +160,30 @@ void ComputeRHEOPropertyAtom::init() error->all(FLERR, "Cannot request velocity shifting property without corresponding option in fix rheo"); if (thermal_flag && !(fix_rheo->thermal_flag)) error->all(FLERR, "Cannot request thermal property without fix rheo/thermal"); - if (shell_flag && !(fix_rheo->oxidation_flag)) - error->all(FLERR, "Cannot request number of shell bonds without fix rheo/oxidation"); compute_interface = fix_rheo->compute_interface; compute_kernel = fix_rheo->compute_kernel; compute_surface = fix_rheo->compute_surface; compute_vshift = fix_rheo->compute_vshift; compute_grad = fix_rheo->compute_grad; +} +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::setup() +{ if (thermal_flag) { - fixes = modify->get_fix_by_style("rheo/thermal"); + auto fixes = modify->get_fix_by_style("rheo/thermal"); fix_thermal = dynamic_cast(fixes[0]); } if (pressure_flag) { - fixes = modify->get_fix_by_style("rheo/pressure"); + auto fixes = modify->get_fix_by_style("rheo/pressure"); fix_pressure = dynamic_cast(fixes[0]); } if (shell_flag) { - fixes = modify->get_fix_by_style("rheo/oxidation"); + auto fixes = modify->get_fix_by_style("rheo/oxidation"); fix_oxidation = dynamic_cast(fixes[0]); } } diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h index a66ef4ece6..fc5c5d1bd0 100644 --- a/src/RHEO/compute_rheo_property_atom.h +++ b/src/RHEO/compute_rheo_property_atom.h @@ -29,6 +29,7 @@ class ComputeRHEOPropertyAtom : public Compute { ComputeRHEOPropertyAtom(class LAMMPS *, int, char **); ~ComputeRHEOPropertyAtom() override; void init() override; + void setup() override; void compute_peratom() override; double memory_usage() override; diff --git a/src/RHEO/fix_rheo_oxidation.h b/src/RHEO/fix_rheo_oxidation.h index 5242d24460..be95efbf2c 100644 --- a/src/RHEO/fix_rheo_oxidation.h +++ b/src/RHEO/fix_rheo_oxidation.h @@ -37,12 +37,13 @@ class FixRHEOOxidation : public Fix { void pre_force(int) override; void post_integrate() override; int *nbond; + double rsurf, cut; private: int btype, index_nb; - double rsurf, cut, cutsq; - class NeighList *list; + double cutsq; + class NeighList *list; class ComputeRHEOSurface *compute_surface; class FixRHEO *fix_rheo; }; diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index 7d9aff5424..a133428d4a 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -270,13 +270,8 @@ void FixRHEOThermal::init() auto req = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL); req->set_cutoff(cut_kernel); - // find instances of bond history to delete data - // skip history for shell, only exception + // find instances of bond history to delete/shift data histories = modify->get_fix_by_style("BOND_HISTORY"); - if (n_histories > 0) - for (int i = 0; i < histories.size(); i++) - if (strcmp(histories[i]->id, "HISTORY_RHEO_SHELL") == 0) - histories.erase(histories.begin() + i); n_histories = histories.size(); } } @@ -498,7 +493,7 @@ void FixRHEOThermal::break_bonds() int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - // Rapidly delete all bonds for local atoms that melt of a given type + // Delete all bonds for local atoms that melt of a given type for (int i = 0; i < nlocal; i++) { if (!(status[i] & STATUS_MELTING)) continue; for (m = (num_bond[i] - 1); m >= 0; m--) { @@ -543,7 +538,7 @@ void FixRHEOThermal::break_bonds() bondlist[n][2] = 0; // Delete bonds for non-melted local atoms (shifting) - if (i < nlocal) { + if (i < nlocal && !(status[i] & STATUS_MELTING)) { for (m = 0; m < num_bond[i]; m++) { if (bond_atom[i][m] == tag[j] && bond_type[i][m] == btype) { nmax = num_bond[i] - 1; @@ -562,7 +557,7 @@ void FixRHEOThermal::break_bonds() } } - if (j < nlocal) { + if (j < nlocal && !(status[j] & STATUS_MELTING)) { for (m = 0; m < num_bond[j]; m++) { if (bond_atom[j][m] == tag[i] && bond_type[j][m] == btype) { nmax = num_bond[j] - 1; diff --git a/src/RHEO/pair_rheo_solid.cpp b/src/RHEO/pair_rheo_solid.cpp index f6dcd95879..9b358420d2 100644 --- a/src/RHEO/pair_rheo_solid.cpp +++ b/src/RHEO/pair_rheo_solid.cpp @@ -327,6 +327,10 @@ double PairRHEOSolid::single(int i, int j, int itype, int jtype, double rsq, dou if (rsq > cutsq[itype][jtype]) return 0.0; + int *status = atom->status; + if (!(status[i] & STATUS_SOLID)) return 0.0; + if (!(status[j] & STATUS_SOLID)) return 0.0; + double **x = atom->x; double **v = atom->v; diff --git a/src/fix_bond_history.cpp b/src/fix_bond_history.cpp index 277df75085..461f91bd52 100644 --- a/src/fix_bond_history.cpp +++ b/src/fix_bond_history.cpp @@ -33,7 +33,7 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ FixBondHistory::FixBondHistory(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), bondstore(nullptr), id_fix(nullptr), id_array(nullptr) + Fix(lmp, narg, arg), bondstore(nullptr), bondtype_orig(nullptr), bondstore_comp(nullptr), bondstore_orig(nullptr), id_fix(nullptr), id_array(nullptr) { if (narg != 5) error->all(FLERR, "Illegal fix bond/history command"); @@ -53,7 +53,6 @@ FixBondHistory::FixBondHistory(LAMMPS *lmp, int narg, char **arg) : updated_bond_flag = 0; maxbond = 0; - allocate(); } /* ---------------------------------------------------------------------- */ @@ -65,6 +64,8 @@ FixBondHistory::~FixBondHistory() delete[] id_array; memory->destroy(bondstore); + memory->destroy(bondstore_comp); + memory->destroy(bondtype_orig); } /* ---------------------------------------------------------------------- */ @@ -135,6 +136,7 @@ void FixBondHistory::pre_exchange() int nlocal = atom->nlocal; tagint **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; int *num_bond = atom->num_bond; tagint *tag = atom->tag; @@ -142,12 +144,12 @@ void FixBondHistory::pre_exchange() i1 = bondlist[n][0]; i2 = bondlist[n][1]; - // skip bond if already broken - if (bondlist[n][2] <= 0) { continue; } + // skip bond if already broken or not allocated + if (bondlist[n][2] <= 0 || !setflag[bondlist[n][2]]) { continue; } if (i1 < nlocal) { for (m = 0; m < num_bond[i1]; m++) { - if (bond_atom[i1][m] == tag[i2]) { + if (bond_atom[i1][m] == tag[i2] && setflag[bond_type[i1][m]]) { for (idata = 0; idata < ndata; idata++) { stored[i1][m * ndata + idata] = bondstore[n][idata]; } @@ -157,7 +159,7 @@ void FixBondHistory::pre_exchange() if (i2 < nlocal) { for (m = 0; m < num_bond[i2]; m++) { - if (bond_atom[i2][m] == tag[i1]) { + if (bond_atom[i2][m] == tag[i1] && setflag[bond_type[i2][m]]) { for (idata = 0; idata < ndata; idata++) { stored[i2][m * ndata + idata] = bondstore[n][idata]; } @@ -179,17 +181,21 @@ void FixBondHistory::allocate() else maxbond = static_cast(LB_FACTOR * atom->nbonds / comm->nprocs); memory->create(bondstore, maxbond, ndata, "fix_bond_store:bondstore"); + if (hybrid_flag) { + memory->create(bondstore_comp, maxbond, ndata, "fix_bond_store:bondstore_comp"); + memory->create(bondtype_orig, maxbond, "fix_bond_store:bondtype_orig"); + } } /* ---------------------------------------------------------------------- */ void FixBondHistory::setup_post_neighbor() { - //Grow array if number of bonds has increased - while (neighbor->nbondlist >= maxbond) { - maxbond += DELTA; - memory->grow(bondstore, maxbond, ndata, "fix_bond_store:bondstore"); - } + hybrid_flag = 0; + for (int i = 1; i <= atom->nbondtypes; i++) + if (!setflag[i]) hybrid_flag = 1; + + if (maxbond == 0) allocate(); pre_exchange(); post_neighbor(); @@ -206,6 +212,10 @@ void FixBondHistory::post_neighbor() while (neighbor->nbondlist >= maxbond) { maxbond += DELTA; memory->grow(bondstore, maxbond, ndata, "fix_bond_store:bondstore"); + if (hybrid_flag) { + memory->grow(bondstore_comp, maxbond, ndata, "fix_bond_store:bondstore_comp"); + memory->grow(bondtype_orig, maxbond, "fix_bond_store:bondtype_orig"); + } } int i1, i2, n, m, idata; @@ -215,6 +225,7 @@ void FixBondHistory::post_neighbor() int nlocal = atom->nlocal; tagint **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; int *num_bond = atom->num_bond; tagint *tag = atom->tag; @@ -222,12 +233,12 @@ void FixBondHistory::post_neighbor() i1 = bondlist[n][0]; i2 = bondlist[n][1]; - // skip bond if already broken - if (bondlist[n][2] <= 0) { continue; } + // skip bond if already broken or not allocated + if (bondlist[n][2] <= 0 || !setflag[bondlist[n][2]]) { continue; } if (i1 < nlocal) { for (m = 0; m < num_bond[i1]; m++) { - if (bond_atom[i1][m] == tag[i2]) { + if (bond_atom[i1][m] == tag[i2] && setflag[bond_type[i1][m]]) { for (idata = 0; idata < ndata; idata++) { bondstore[n][idata] = stored[i1][m * ndata + idata]; } @@ -237,7 +248,7 @@ void FixBondHistory::post_neighbor() if (i2 < nlocal) { for (m = 0; m < num_bond[i2]; m++) { - if (bond_atom[i2][m] == tag[i1]) { + if (bond_atom[i2][m] == tag[i1] && setflag[bond_type[i2][m]]) { for (idata = 0; idata < ndata; idata++) { bondstore[n][idata] = stored[i2][m * ndata + idata]; } @@ -246,6 +257,12 @@ void FixBondHistory::post_neighbor() } } + if (hybrid_flag) { + nbondlist_orig = nbondlist; + for (n = 0; n < nbondlist; n++) + bondtype_orig[n] = bondlist[n][2]; + } + updated_bond_flag = 1; } @@ -294,6 +311,55 @@ void FixBondHistory::set_arrays(int i) for (int idata = 0; idata < ndata; idata++) stored[i][m * ndata + idata] = 0.0; } +/* ---------------------------------------------------------------------- + Compress history arrays, cutting out unused types, for bond hybrid +------------------------------------------------------------------------- */ + +void FixBondHistory::compress_history() +{ + // if this is a re-neighbor step or updating, compress bondstore + + int type; + int ncomp = 0; + if (update_flag || (neighbor->ago == 0)) { + for (int n = 0; n < nbondlist_orig; n++) { + type = bondtype_orig[n]; + if (type <= 0) continue; + if (setflag[type]) { + for (int m = 0; m < ndata; m++) + bondstore_comp[ncomp][m] = bondstore[n][m]; + ncomp += 1; + } + } + } + + // replace ptr to original array + bondstore_orig = bondstore; + bondstore = bondstore_comp; +} + +/* ---------------------------------------------------------------------- */ + +void FixBondHistory::uncompress_history() +{ + if (update_flag) { + int ncomp = 0; + for (int n = 0; n < nbondlist_orig; n++) { + type = bondtype_orig[n]; + + if (type <= 0) continue; + if (!setflag[type]) continue; + + for (int m = 0; m < ndata; m++) + bondstore[n][m] = bondstore_comp[ncomp][m]; + ncomp += 1; + } + } + + // restore ptr to original array + bondstore = bondstore_orig; +} + /* ---------------------------------------------------------------------- Delete bond by zeroing data ------------------------------------------------------------------------- */ diff --git a/src/fix_bond_history.h b/src/fix_bond_history.h index fafcf52bd9..8ee3132ab1 100644 --- a/src/fix_bond_history.h +++ b/src/fix_bond_history.h @@ -52,18 +52,29 @@ class FixBondHistory : public Fix { void check_cache(int, int); void clear_cache(); + // methods for bond style hybrid + void compress_history(); + void uncompress_history(); + // if data is temporarily stored while the bond_atom array // is being reordered, use map of vectors with pairs for keys // to enable quick look up std::map, std::vector> cached_histories; + int *setflag; // Set by BondBPM, which bond types are used double **bondstore; int stored_flag; protected: void allocate(); - int update_flag; //Flag whether history values can evolve + int hybrid_flag; + int nbondlist_orig; + int *bondtype_orig; + double **bondstore_comp; + double **bondstore_orig; + + int update_flag; // Flag whether history values can evolve int updated_bond_flag; int nbond, maxbond, ndata; int index;