From 0d54f99fc0d7516011a503846a54ae56c9ebcbed Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Thu, 4 Jul 2024 17:34:18 +0200 Subject: [PATCH 1/9] Fixed a nasty bug when running the ipi interface with many MPI workers Looks like the migrate() call should happen before trajectory folding --- src/MISC/fix_ipi.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/MISC/fix_ipi.cpp b/src/MISC/fix_ipi.cpp index 87668b9192..e08727a4ed 100644 --- a/src/MISC/fix_ipi.cpp +++ b/src/MISC/fix_ipi.cpp @@ -373,6 +373,9 @@ void FixIPI::initial_integrate(int /*vflag*/) if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); + // move atoms to new processors via irregular() + // only needed if migrate_check() says an atom moves to far + if (irregular->migrate_check()) irregular->migrate_atoms(); if (domain->triclinic) domain->lamda2x(atom->nlocal); // ensures continuity of trajectories relative to the @@ -394,11 +397,6 @@ void FixIPI::initial_integrate(int /*vflag*/) } } } - // move atoms to new processors via irregular() - // only needed if migrate_check() says an atom moves to far - if (domain->triclinic) domain->x2lamda(atom->nlocal); - if (irregular->migrate_check()) irregular->migrate_atoms(); - if (domain->triclinic) domain->lamda2x(atom->nlocal); // check if kspace solver is used if (reset_flag && kspace_flag) { From 810698dc07ca91b7df46f352b7b90d07b72c00d3 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Fri, 5 Jul 2024 09:07:30 +0200 Subject: [PATCH 2/9] Skip neighbor folding at first call There may be issues due to the fact that the NL is initialized with the LAMMPS data file, that does not have to be the same as the starting config in i-PI --- src/MISC/fix_ipi.cpp | 9 ++++++++- src/MISC/fix_ipi.h | 1 + 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/MISC/fix_ipi.cpp b/src/MISC/fix_ipi.cpp index e08727a4ed..80666790e2 100644 --- a/src/MISC/fix_ipi.cpp +++ b/src/MISC/fix_ipi.cpp @@ -188,6 +188,7 @@ FixIPI::FixIPI(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), irregul master = (comm->me == 0) ? 1 : 0; inet = 1; reset_flag = 0; + firsttime = 1; int iarg = 5; while (iarg < narg) { @@ -253,6 +254,7 @@ void FixIPI::init() ipisock = 0; // TODO: should check for success in socket opening, // but the current open_socket routine dies brutally if unsuccessful + // tell lammps we have assigned a socket socketflag = 1; @@ -382,7 +384,11 @@ void FixIPI::initial_integrate(int /*vflag*/) // snapshot at neighbor list creation, minimizing the // number of neighbor list updates auto xhold = neighbor->get_xhold(); - if (xhold != NULL) { // don't wrap if xhold is not used in the NL + if (xhold != NULL && !firsttime) { + // don't wrap if xhold is not used in the NL, or the + // first call (because the NL is initialized from the + // data file that might have nothing to do with the + // current structure for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { auto delx = x[i][0] - xhold[i][0]; @@ -397,6 +403,7 @@ void FixIPI::initial_integrate(int /*vflag*/) } } } + firsttime = 0; // check if kspace solver is used if (reset_flag && kspace_flag) { diff --git a/src/MISC/fix_ipi.h b/src/MISC/fix_ipi.h index de954c8578..31bd59b2fa 100644 --- a/src/MISC/fix_ipi.h +++ b/src/MISC/fix_ipi.h @@ -42,6 +42,7 @@ class FixIPI : public Fix { long bsize; int kspace_flag; int reset_flag; + int firsttime; private: class Irregular *irregular; From e02ad1a3b208797cc9c5b3a08c24cabb5065b7a6 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 25 Jul 2024 17:02:22 -0600 Subject: [PATCH 3/9] Cleaning up SPH package, minor bug fixes --- src/SPH/compute_sph_e_atom.cpp | 30 ++++++------ src/SPH/compute_sph_rho_atom.cpp | 33 ++++++------- src/SPH/compute_sph_t_atom.cpp | 32 ++++++------- src/SPH/fix_sph.cpp | 32 ++++++++----- src/SPH/fix_sph_stationary.cpp | 27 ++++++----- src/SPH/fix_sph_stationary.h | 3 -- src/SPH/pair_sph_heatconduction.cpp | 33 ++++++++----- src/SPH/pair_sph_idealgas.cpp | 52 ++++++++++---------- src/SPH/pair_sph_idealgas.h | 1 - src/SPH/pair_sph_lj.cpp | 52 ++++++++++---------- src/SPH/pair_sph_lj.h | 2 - src/SPH/pair_sph_rhosum.cpp | 56 +++++++++------------- src/SPH/pair_sph_taitwater.cpp | 64 ++++++++++--------------- src/SPH/pair_sph_taitwater_morris.cpp | 68 +++++++++++---------------- 14 files changed, 227 insertions(+), 258 deletions(-) diff --git a/src/SPH/compute_sph_e_atom.cpp b/src/SPH/compute_sph_e_atom.cpp index 05e752c49b..8bb0622acd 100644 --- a/src/SPH/compute_sph_e_atom.cpp +++ b/src/SPH/compute_sph_e_atom.cpp @@ -13,13 +13,15 @@ ------------------------------------------------------------------------- */ #include "compute_sph_e_atom.h" -#include + #include "atom.h" -#include "update.h" -#include "modify.h" #include "comm.h" -#include "memory.h" #include "error.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include using namespace LAMMPS_NS; @@ -31,7 +33,7 @@ ComputeSPHEAtom::ComputeSPHEAtom(LAMMPS *lmp, int narg, char **arg) : if (narg != 3) error->all(FLERR,"Number of arguments for compute sph/e/atom command != 3"); if (atom->esph_flag != 1) - error->all(FLERR,"Compute sph/e/atom command requires atom_style sph)"); + error->all(FLERR,"Compute sph/e/atom requires atom attribut energy, e.g. in atom_style sph"); peratom_flag = 1; size_peratom_cols = 0; @@ -51,12 +53,11 @@ ComputeSPHEAtom::~ComputeSPHEAtom() void ComputeSPHEAtom::init() { - int count = 0; for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style,"evector/atom") == 0) count++; + if (strcmp(modify->compute[i]->style,"sph/e/atom") == 0) count++; if (count > 1 && comm->me == 0) - error->warning(FLERR,"More than one compute evector/atom"); + error->warning(FLERR,"More than one compute sph/e/atom"); } /* ---------------------------------------------------------------------- */ @@ -78,14 +79,13 @@ void ComputeSPHEAtom::compute_peratom() int *mask = atom->mask; int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - evector[i] = esph[i]; - } - else { - evector[i] = 0.0; - } + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + evector[i] = esph[i]; + } else { + evector[i] = 0.0; } + } } /* ---------------------------------------------------------------------- diff --git a/src/SPH/compute_sph_rho_atom.cpp b/src/SPH/compute_sph_rho_atom.cpp index 9c596bdfc3..54d1237368 100644 --- a/src/SPH/compute_sph_rho_atom.cpp +++ b/src/SPH/compute_sph_rho_atom.cpp @@ -13,13 +13,15 @@ ------------------------------------------------------------------------- */ #include "compute_sph_rho_atom.h" -#include + #include "atom.h" -#include "update.h" -#include "modify.h" #include "comm.h" -#include "memory.h" #include "error.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include using namespace LAMMPS_NS; @@ -30,8 +32,7 @@ ComputeSPHRhoAtom::ComputeSPHRhoAtom(LAMMPS *lmp, int narg, char **arg) : { if (narg != 3) error->all(FLERR,"Illegal compute sph/rho/atom command"); if (atom->rho_flag != 1) - error->all(FLERR,"Compute sph/rho/atom command requires atom_style sph"); - + error->all(FLERR, "Compute sph/rho/atom requires atom attribute density, e.g. in atom_style sph"); peratom_flag = 1; size_peratom_cols = 0; @@ -50,12 +51,11 @@ ComputeSPHRhoAtom::~ComputeSPHRhoAtom() void ComputeSPHRhoAtom::init() { - int count = 0; for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style,"rhoVector/atom") == 0) count++; + if (strcmp(modify->compute[i]->style,"sph/rho/atom") == 0) count++; if (count > 1 && comm->me == 0) - error->warning(FLERR,"More than one compute rhoVector/atom"); + error->warning(FLERR,"More than one compute sph/rho/atom"); } /* ---------------------------------------------------------------------- */ @@ -73,20 +73,17 @@ void ComputeSPHRhoAtom::compute_peratom() vector_atom = rhoVector; } - // compute kinetic energy for each atom in group - double *rho = atom->rho; int *mask = atom->mask; int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - rhoVector[i] = rho[i]; - } - else { - rhoVector[i] = 0.0; - } + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + rhoVector[i] = rho[i]; + } else { + rhoVector[i] = 0.0; } + } } /* ---------------------------------------------------------------------- diff --git a/src/SPH/compute_sph_t_atom.cpp b/src/SPH/compute_sph_t_atom.cpp index 4a83013434..e795c32a5e 100644 --- a/src/SPH/compute_sph_t_atom.cpp +++ b/src/SPH/compute_sph_t_atom.cpp @@ -13,13 +13,15 @@ ------------------------------------------------------------------------- */ #include "compute_sph_t_atom.h" -#include + #include "atom.h" -#include "update.h" -#include "modify.h" #include "comm.h" -#include "memory.h" #include "error.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include using namespace LAMMPS_NS; @@ -31,7 +33,7 @@ ComputeSPHTAtom::ComputeSPHTAtom(LAMMPS *lmp, int narg, char **arg) : if (narg != 3) error->all(FLERR,"Number of arguments for compute sph/t/atom command != 3"); if ((atom->esph_flag != 1) || (atom->cv_flag != 1)) - error->all(FLERR,"Compute sph/t/atom command requires atom_style sph"); + error->all(FLERR, "Compute sph/t/atom requires atom attributes energy and specific heat, e.g. in atom_style sph"); peratom_flag = 1; size_peratom_cols = 0; @@ -51,12 +53,11 @@ ComputeSPHTAtom::~ComputeSPHTAtom() void ComputeSPHTAtom::init() { - int count = 0; for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style,"meso/t/atom") == 0) count++; + if (strcmp(modify->compute[i]->style,"sph/t/atom") == 0) count++; if (count > 1 && comm->me == 0) - error->warning(FLERR,"More than one compute meso/t/atom"); + error->warning(FLERR,"More than one compute sph/t/atom"); } /* ---------------------------------------------------------------------- */ @@ -79,16 +80,15 @@ void ComputeSPHTAtom::compute_peratom() int *mask = atom->mask; int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - if (cv[i] > 0.0) { - tvector[i] = esph[i] / cv[i]; - } - } - else { - tvector[i] = 0.0; + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (cv[i] > 0.0) { + tvector[i] = esph[i] / cv[i]; } + } else { + tvector[i] = 0.0; } + } } /* ---------------------------------------------------------------------- diff --git a/src/SPH/fix_sph.cpp b/src/SPH/fix_sph.cpp index 0dafcee7c0..96479be928 100644 --- a/src/SPH/fix_sph.cpp +++ b/src/SPH/fix_sph.cpp @@ -13,10 +13,13 @@ ------------------------------------------------------------------------- */ #include "fix_sph.h" + #include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" #include "force.h" #include "update.h" -#include "error.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -24,11 +27,10 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ FixSPH::FixSPH(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) { - - if ((atom->esph_flag != 1) || (atom->rho_flag != 1)) - error->all(FLERR, - "Fix sph command requires atom_style with both energy and density"); + Fix(lmp, narg, arg) +{ + if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->vest_flag != 1)) + error->all(FLERR, "Fix sph requires atom attributes energy, density, and velocity estimates, e.g. in atom_style sph"); if (narg != 3) error->all(FLERR,"Illegal number of arguments for fix sph command"); @@ -38,7 +40,8 @@ FixSPH::FixSPH(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -int FixSPH::setmask() { +int FixSPH::setmask() +{ int mask = 0; mask |= INITIAL_INTEGRATE; mask |= FINAL_INTEGRATE; @@ -48,11 +51,14 @@ int FixSPH::setmask() { /* ---------------------------------------------------------------------- */ -void FixSPH::init() { +void FixSPH::init() +{ dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; } +/* ---------------------------------------------------------------------- */ + void FixSPH::setup_pre_force(int /*vflag*/) { // set vest equal to v @@ -76,7 +82,8 @@ void FixSPH::setup_pre_force(int /*vflag*/) allow for both per-type and per-atom mass ------------------------------------------------------------------------- */ -void FixSPH::initial_integrate(int /*vflag*/) { +void FixSPH::initial_integrate(int /*vflag*/) +{ // update v and x and rho and e of atoms in group double **x = atom->x; @@ -129,8 +136,8 @@ void FixSPH::initial_integrate(int /*vflag*/) { /* ---------------------------------------------------------------------- */ -void FixSPH::final_integrate() { - +void FixSPH::final_integrate() +{ // update v, rho, and e of atoms in group double **v = atom->v; @@ -169,7 +176,8 @@ void FixSPH::final_integrate() { /* ---------------------------------------------------------------------- */ -void FixSPH::reset_dt() { +void FixSPH::reset_dt() +{ dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; } diff --git a/src/SPH/fix_sph_stationary.cpp b/src/SPH/fix_sph_stationary.cpp index 673a7d50bf..877bb5b17b 100644 --- a/src/SPH/fix_sph_stationary.cpp +++ b/src/SPH/fix_sph_stationary.cpp @@ -13,10 +13,11 @@ ------------------------------------------------------------------------- */ #include "fix_sph_stationary.h" + #include "atom.h" +#include "error.h" #include "force.h" #include "update.h" -#include "error.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -24,11 +25,10 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ FixSPHStationary::FixSPHStationary(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) { - + Fix(lmp, narg, arg) +{ if ((atom->esph_flag != 1) || (atom->rho_flag != 1)) - error->all(FLERR, - "Fix sph/stationary command requires atom_style with both energy and density, e.g. meso"); + error->all(FLERR, "Fix sph/stationary requires atom attributes energy and density, e.g. in atom_style sph"); if (narg != 3) error->all(FLERR,"Illegal number of arguments for fix sph/stationary command"); @@ -38,7 +38,8 @@ FixSPHStationary::FixSPHStationary(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -int FixSPHStationary::setmask() { +int FixSPHStationary::setmask() +{ int mask = 0; mask |= INITIAL_INTEGRATE; mask |= FINAL_INTEGRATE; @@ -47,7 +48,8 @@ int FixSPHStationary::setmask() { /* ---------------------------------------------------------------------- */ -void FixSPHStationary::init() { +void FixSPHStationary::init() +{ dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; } @@ -56,8 +58,8 @@ void FixSPHStationary::init() { allow for both per-type and per-atom mass ------------------------------------------------------------------------- */ -void FixSPHStationary::initial_integrate(int /*vflag*/) { - +void FixSPHStationary::initial_integrate(int /*vflag*/) +{ double *rho = atom->rho; double *drho = atom->drho; double *esph = atom->esph; @@ -80,8 +82,8 @@ void FixSPHStationary::initial_integrate(int /*vflag*/) { /* ---------------------------------------------------------------------- */ -void FixSPHStationary::final_integrate() { - +void FixSPHStationary::final_integrate() +{ double *esph = atom->esph; double *desph = atom->desph; double *rho = atom->rho; @@ -101,7 +103,8 @@ void FixSPHStationary::final_integrate() { /* ---------------------------------------------------------------------- */ -void FixSPHStationary::reset_dt() { +void FixSPHStationary::reset_dt() +{ dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; } diff --git a/src/SPH/fix_sph_stationary.h b/src/SPH/fix_sph_stationary.h index 51e336a038..78ece5e5d4 100644 --- a/src/SPH/fix_sph_stationary.h +++ b/src/SPH/fix_sph_stationary.h @@ -33,9 +33,6 @@ class FixSPHStationary : public Fix { void final_integrate() override; void reset_dt() override; - private: - class NeighList *list; - protected: double dtv, dtf; double *step_respa; diff --git a/src/SPH/pair_sph_heatconduction.cpp b/src/SPH/pair_sph_heatconduction.cpp index cb27982fa3..145d9cadee 100644 --- a/src/SPH/pair_sph_heatconduction.cpp +++ b/src/SPH/pair_sph_heatconduction.cpp @@ -13,14 +13,15 @@ ------------------------------------------------------------------------- */ #include "pair_sph_heatconduction.h" -#include + #include "atom.h" +#include "domain.h" +#include "error.h" #include "force.h" #include "memory.h" -#include "error.h" #include "neigh_list.h" -#include "domain.h" +#include using namespace LAMMPS_NS; @@ -28,12 +29,16 @@ using namespace LAMMPS_NS; PairSPHHeatConduction::PairSPHHeatConduction(LAMMPS *lmp) : Pair(lmp) { + if ((atom->esph_flag != 1) || (atom->rho_flag != 1)) + error->all(FLERR, "Pair sph/heatconduction requires atom attributes energy and density, e.g. in atom_style sph"); + restartinfo = 0; } /* ---------------------------------------------------------------------- */ -PairSPHHeatConduction::~PairSPHHeatConduction() { +PairSPHHeatConduction::~PairSPHHeatConduction() +{ if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); @@ -44,7 +49,8 @@ PairSPHHeatConduction::~PairSPHHeatConduction() { /* ---------------------------------------------------------------------- */ -void PairSPHHeatConduction::compute(int eflag, int vflag) { +void PairSPHHeatConduction::compute(int eflag, int vflag) +{ int i, j, ii, jj, inum, jnum, itype, jtype; double xtmp, ytmp, ztmp, delx, dely, delz; @@ -124,7 +130,6 @@ void PairSPHHeatConduction::compute(int eflag, int vflag) { if (newton_pair || j < nlocal) { desph[j] -= deltaE; } - } } } @@ -134,7 +139,8 @@ void PairSPHHeatConduction::compute(int eflag, int vflag) { allocate all arrays ------------------------------------------------------------------------- */ -void PairSPHHeatConduction::allocate() { +void PairSPHHeatConduction::allocate() +{ allocated = 1; int n = atom->ntypes; @@ -152,7 +158,8 @@ void PairSPHHeatConduction::allocate() { global settings ------------------------------------------------------------------------- */ -void PairSPHHeatConduction::settings(int narg, char **/*arg*/) { +void PairSPHHeatConduction::settings(int narg, char **/*arg*/) +{ if (narg != 0) error->all(FLERR, "Illegal number of arguments for pair_style sph/heatconduction"); @@ -162,7 +169,8 @@ void PairSPHHeatConduction::settings(int narg, char **/*arg*/) { set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairSPHHeatConduction::coeff(int narg, char **arg) { +void PairSPHHeatConduction::coeff(int narg, char **arg) +{ if (narg != 4) error->all(FLERR,"Incorrect number of args for pair_style sph/heatconduction coefficients"); if (!allocated) @@ -178,7 +186,6 @@ void PairSPHHeatConduction::coeff(int narg, char **arg) { int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { - //printf("setting cut[%d][%d] = %f\n", i, j, cut_one); cut[i][j] = cut_one; alpha[i][j] = alpha_one; setflag[i][j] = 1; @@ -194,7 +201,8 @@ void PairSPHHeatConduction::coeff(int narg, char **arg) { init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ -double PairSPHHeatConduction::init_one(int i, int j) { +double PairSPHHeatConduction::init_one(int i, int j) +{ if (setflag[i][j] == 0) { error->all(FLERR,"All pair sph/heatconduction coeffs are not set"); @@ -209,7 +217,8 @@ double PairSPHHeatConduction::init_one(int i, int j) { /* ---------------------------------------------------------------------- */ double PairSPHHeatConduction::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, - double /*rsq*/, double /*factor_coul*/, double /*factor_lj*/, double &fforce) { + double /*rsq*/, double /*factor_coul*/, double /*factor_lj*/, double &fforce) +{ fforce = 0.0; return 0.0; diff --git a/src/SPH/pair_sph_idealgas.cpp b/src/SPH/pair_sph_idealgas.cpp index 114d323bbd..97cbda2ab3 100644 --- a/src/SPH/pair_sph_idealgas.cpp +++ b/src/SPH/pair_sph_idealgas.cpp @@ -13,14 +13,15 @@ ------------------------------------------------------------------------- */ #include "pair_sph_idealgas.h" -#include -#include "atom.h" -#include "force.h" -#include "neigh_list.h" -#include "memory.h" -#include "error.h" -#include "domain.h" +#include "atom.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" + +#include using namespace LAMMPS_NS; @@ -28,12 +29,17 @@ using namespace LAMMPS_NS; PairSPHIdealGas::PairSPHIdealGas(LAMMPS *lmp) : Pair(lmp) { + if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->vest_flag != 1)) + error->all(FLERR, "Pair sph/idealgas requires atom attributes energy, density, and velocity estimates, e.g. in atom_style sph"); + restartinfo = 0; + single_enable = 0; } /* ---------------------------------------------------------------------- */ -PairSPHIdealGas::~PairSPHIdealGas() { +PairSPHIdealGas::~PairSPHIdealGas() +{ if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); @@ -45,7 +51,8 @@ PairSPHIdealGas::~PairSPHIdealGas() { /* ---------------------------------------------------------------------- */ -void PairSPHIdealGas::compute(int eflag, int vflag) { +void PairSPHIdealGas::compute(int eflag, int vflag) +{ int i, j, ii, jj, inum, jnum, itype, jtype; double xtmp, ytmp, ztmp, delx, dely, delz, fpair; @@ -160,10 +167,6 @@ void PairSPHIdealGas::compute(int eflag, int vflag) { if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz); - - if (evflag) - ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, - delz); } } } @@ -175,7 +178,8 @@ void PairSPHIdealGas::compute(int eflag, int vflag) { allocate all arrays ------------------------------------------------------------------------- */ -void PairSPHIdealGas::allocate() { +void PairSPHIdealGas::allocate() +{ allocated = 1; int n = atom->ntypes; @@ -194,7 +198,8 @@ void PairSPHIdealGas::allocate() { global settings ------------------------------------------------------------------------- */ -void PairSPHIdealGas::settings(int narg, char **/*arg*/) { +void PairSPHIdealGas::settings(int narg, char **/*arg*/) +{ if (narg != 0) error->all(FLERR, "Illegal number of arguments for pair_style sph/idealgas"); @@ -204,7 +209,8 @@ void PairSPHIdealGas::settings(int narg, char **/*arg*/) { set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairSPHIdealGas::coeff(int narg, char **arg) { +void PairSPHIdealGas::coeff(int narg, char **arg) +{ if (narg != 4) error->all(FLERR,"Incorrect number of args for pair_style sph/idealgas coefficients"); if (!allocated) @@ -221,7 +227,6 @@ void PairSPHIdealGas::coeff(int narg, char **arg) { for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { viscosity[i][j] = viscosity_one; - //printf("setting cut[%d][%d] = %f\n", i, j, cut_one); cut[i][j] = cut_one; setflag[i][j] = 1; count++; @@ -236,8 +241,8 @@ void PairSPHIdealGas::coeff(int narg, char **arg) { init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ -double PairSPHIdealGas::init_one(int i, int j) { - +double PairSPHIdealGas::init_one(int i, int j) +{ if (setflag[i][j] == 0) { error->all(FLERR,"All pair sph/idealgas coeffs are not set"); } @@ -246,12 +251,3 @@ double PairSPHIdealGas::init_one(int i, int j) { return cut[i][j]; } - -/* ---------------------------------------------------------------------- */ - -double PairSPHIdealGas::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, - double /*rsq*/, double /*factor_coul*/, double /*factor_lj*/, double &fforce) { - fforce = 0.0; - - return 0.0; -} diff --git a/src/SPH/pair_sph_idealgas.h b/src/SPH/pair_sph_idealgas.h index 3d61d76616..c3c11ac5bd 100644 --- a/src/SPH/pair_sph_idealgas.h +++ b/src/SPH/pair_sph_idealgas.h @@ -32,7 +32,6 @@ class PairSPHIdealGas : public Pair { void settings(int, char **) override; void coeff(int, char **) override; double init_one(int, int) override; - double single(int, int, int, int, double, double, double, double &) override; protected: double **cut, **viscosity; diff --git a/src/SPH/pair_sph_lj.cpp b/src/SPH/pair_sph_lj.cpp index d59a3ad992..4bdefca1a6 100644 --- a/src/SPH/pair_sph_lj.cpp +++ b/src/SPH/pair_sph_lj.cpp @@ -13,14 +13,15 @@ ------------------------------------------------------------------------- */ #include "pair_sph_lj.h" -#include -#include "atom.h" -#include "force.h" -#include "neigh_list.h" -#include "memory.h" -#include "error.h" -#include "domain.h" +#include "atom.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" + +#include using namespace LAMMPS_NS; @@ -28,12 +29,17 @@ using namespace LAMMPS_NS; PairSPHLJ::PairSPHLJ(LAMMPS *lmp) : Pair(lmp) { + if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->cv_flag != 1) || (atom->vest_flag != 1)) + error->all(FLERR, "Pair sph/lj requires atom attributes energy, density, specific heat, and velocity estimates, e.g. in atom_style sph"); + restartinfo = 0; + single_enable = 0; } /* ---------------------------------------------------------------------- */ -PairSPHLJ::~PairSPHLJ() { +PairSPHLJ::~PairSPHLJ() +{ if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); @@ -45,7 +51,8 @@ PairSPHLJ::~PairSPHLJ() { /* ---------------------------------------------------------------------- */ -void PairSPHLJ::compute(int eflag, int vflag) { +void PairSPHLJ::compute(int eflag, int vflag) +{ int i, j, ii, jj, inum, jnum, itype, jtype; double xtmp, ytmp, ztmp, delx, dely, delz, fpair; @@ -182,7 +189,8 @@ void PairSPHLJ::compute(int eflag, int vflag) { allocate all arrays ------------------------------------------------------------------------- */ -void PairSPHLJ::allocate() { +void PairSPHLJ::allocate() +{ allocated = 1; int n = atom->ntypes; @@ -201,7 +209,8 @@ void PairSPHLJ::allocate() { global settings ------------------------------------------------------------------------- */ -void PairSPHLJ::settings(int narg, char **/*arg*/) { +void PairSPHLJ::settings(int narg, char **/*arg*/) +{ if (narg != 0) error->all(FLERR, "Illegal number of arguments for pair_style sph/lj"); @@ -211,7 +220,8 @@ void PairSPHLJ::settings(int narg, char **/*arg*/) { set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairSPHLJ::coeff(int narg, char **arg) { +void PairSPHLJ::coeff(int narg, char **arg) +{ if (narg != 4) error->all(FLERR, "Incorrect args for pair_style sph/lj coefficients"); @@ -229,7 +239,6 @@ void PairSPHLJ::coeff(int narg, char **arg) { for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { viscosity[i][j] = viscosity_one; - printf("setting cut[%d][%d] = %f\n", i, j, cut_one); cut[i][j] = cut_one; setflag[i][j] = 1; count++; @@ -244,8 +253,8 @@ void PairSPHLJ::coeff(int narg, char **arg) { init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ -double PairSPHLJ::init_one(int i, int j) { - +double PairSPHLJ::init_one(int i, int j) +{ if (setflag[i][j] == 0) { error->all(FLERR,"All pair sph/lj coeffs are not set"); } @@ -256,16 +265,6 @@ double PairSPHLJ::init_one(int i, int j) { return cut[i][j]; } -/* ---------------------------------------------------------------------- */ - -double PairSPHLJ::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, - double /*rsq*/, double /*factor_coul*/, double /*factor_lj*/, double &fforce) { - fforce = 0.0; - - return 0.0; -} - - /*double PairSPHLJ::LJEOS2(double rho, double e, double cv) { @@ -297,7 +296,8 @@ double PairSPHLJ::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, Journal of Chemical Physics 73 pp. 5401-5403 (1980) */ -void PairSPHLJ::LJEOS2(double rho, double e, double cv, double *p, double *c) { +void PairSPHLJ::LJEOS2(double rho, double e, double cv, double *p, double *c) +{ double T = e/cv; double beta = 1.0 / T; double beta_sqrt = sqrt(beta); diff --git a/src/SPH/pair_sph_lj.h b/src/SPH/pair_sph_lj.h index 47c34159e1..eb93a14de9 100644 --- a/src/SPH/pair_sph_lj.h +++ b/src/SPH/pair_sph_lj.h @@ -32,8 +32,6 @@ class PairSPHLJ : public Pair { void settings(int, char **) override; void coeff(int, char **) override; double init_one(int, int) override; - double single(int, int, int, int, double, double, double, double &) override; - //double LJEOS(int); void LJEOS2(double, double, double, double *, double *); protected: diff --git a/src/SPH/pair_sph_rhosum.cpp b/src/SPH/pair_sph_rhosum.cpp index 34adf3ce87..4fd96319ee 100644 --- a/src/SPH/pair_sph_rhosum.cpp +++ b/src/SPH/pair_sph_rhosum.cpp @@ -29,6 +29,9 @@ using namespace LAMMPS_NS; PairSPHRhoSum::PairSPHRhoSum(LAMMPS *lmp) : Pair(lmp) { + if (atom->rho_flag != 1) + error->all(FLERR, "Pair sph/rhosum requires atom attribute density, e.g. in atom_style sph"); + restartinfo = 0; // set comm size needed by this Pair @@ -39,11 +42,11 @@ PairSPHRhoSum::PairSPHRhoSum(LAMMPS *lmp) : Pair(lmp) /* ---------------------------------------------------------------------- */ -PairSPHRhoSum::~PairSPHRhoSum() { +PairSPHRhoSum::~PairSPHRhoSum() +{ if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); - memory->destroy(cut); } } @@ -52,14 +55,16 @@ PairSPHRhoSum::~PairSPHRhoSum() { init specific to this pair style ------------------------------------------------------------------------- */ -void PairSPHRhoSum::init_style() { +void PairSPHRhoSum::init_style() +{ // need a full neighbor list neighbor->add_request(this, NeighConst::REQ_FULL); } /* ---------------------------------------------------------------------- */ -void PairSPHRhoSum::compute(int eflag, int vflag) { +void PairSPHRhoSum::compute(int eflag, int vflag) +{ int i, j, ii, jj, jnum, itype, jtype; double xtmp, ytmp, ztmp, delx, dely, delz; double rsq, imass, h, ih, ihsq; @@ -75,25 +80,6 @@ void PairSPHRhoSum::compute(int eflag, int vflag) { int *type = atom->type; double *mass = atom->mass; - // check consistency of pair coefficients - - if (first) { - for (i = 1; i <= atom->ntypes; i++) { - for (j = 1; i <= atom->ntypes; i++) { - if (cutsq[i][j] > 0.0) { - if (!setflag[i][i] || !setflag[j][j]) { - if (comm->me == 0) { - printf( - "SPH particle types %d and %d interact, but not all of their single particle properties are set.\n", - i, j); - } - } - } - } - } - first = 0; - } - inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -186,7 +172,6 @@ void PairSPHRhoSum::compute(int eflag, int vflag) { rho[i] += mass[jtype] * wf; } - } } } @@ -200,7 +185,8 @@ void PairSPHRhoSum::compute(int eflag, int vflag) { allocate all arrays ------------------------------------------------------------------------- */ -void PairSPHRhoSum::allocate() { +void PairSPHRhoSum::allocate() +{ allocated = 1; int n = atom->ntypes; @@ -210,7 +196,6 @@ void PairSPHRhoSum::allocate() { setflag[i][j] = 0; memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); - memory->create(cut, n + 1, n + 1, "pair:cut"); } @@ -218,7 +203,8 @@ void PairSPHRhoSum::allocate() { global settings ------------------------------------------------------------------------- */ -void PairSPHRhoSum::settings(int narg, char **arg) { +void PairSPHRhoSum::settings(int narg, char **arg) +{ if (narg != 1) error->all(FLERR, "Illegal number of arguments for pair_style sph/rhosum"); @@ -229,7 +215,8 @@ void PairSPHRhoSum::settings(int narg, char **arg) { set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairSPHRhoSum::coeff(int narg, char **arg) { +void PairSPHRhoSum::coeff(int narg, char **arg) +{ if (narg != 3) error->all(FLERR,"Incorrect number of args for sph/rhosum coefficients"); if (!allocated) @@ -244,7 +231,6 @@ void PairSPHRhoSum::coeff(int narg, char **arg) { int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { - //printf("setting cut[%d][%d] = %f\n", i, j, cut_one); cut[i][j] = cut_one; setflag[i][j] = 1; count++; @@ -259,7 +245,8 @@ void PairSPHRhoSum::coeff(int narg, char **arg) { init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ -double PairSPHRhoSum::init_one(int i, int j) { +double PairSPHRhoSum::init_one(int i, int j) +{ if (setflag[i][j] == 0) { error->all(FLERR,"All pair sph/rhosum coeffs are not set"); } @@ -272,7 +259,8 @@ double PairSPHRhoSum::init_one(int i, int j) { /* ---------------------------------------------------------------------- */ double PairSPHRhoSum::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, double /*rsq*/, - double /*factor_coul*/, double /*factor_lj*/, double &fforce) { + double /*factor_coul*/, double /*factor_lj*/, double &fforce) +{ fforce = 0.0; return 0.0; @@ -281,7 +269,8 @@ double PairSPHRhoSum::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, /* ---------------------------------------------------------------------- */ int PairSPHRhoSum::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) { + int /*pbc_flag*/, int * /*pbc*/) +{ int i, j, m; double *rho = atom->rho; @@ -295,7 +284,8 @@ int PairSPHRhoSum::pack_forward_comm(int n, int *list, double *buf, /* ---------------------------------------------------------------------- */ -void PairSPHRhoSum::unpack_forward_comm(int n, int first, double *buf) { +void PairSPHRhoSum::unpack_forward_comm(int n, int first, double *buf) +{ int i, m, last; double *rho = atom->rho; diff --git a/src/SPH/pair_sph_taitwater.cpp b/src/SPH/pair_sph_taitwater.cpp index 9a5991718d..d97492de63 100644 --- a/src/SPH/pair_sph_taitwater.cpp +++ b/src/SPH/pair_sph_taitwater.cpp @@ -13,15 +13,16 @@ ------------------------------------------------------------------------- */ #include "pair_sph_taitwater.h" -#include -#include "atom.h" -#include "force.h" -#include "comm.h" -#include "neigh_list.h" -#include "memory.h" -#include "error.h" -#include "domain.h" +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" + +#include using namespace LAMMPS_NS; @@ -29,6 +30,9 @@ using namespace LAMMPS_NS; PairSPHTaitwater::PairSPHTaitwater(LAMMPS *lmp) : Pair(lmp) { + if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->vest_flag != 1)) + error->all(FLERR, "Pair sph/taitwater requires atom attributes energy, density, and velocity estimates, e.g. in atom_style sph"); + restartinfo = 0; single_enable = 0; first = 1; @@ -36,7 +40,8 @@ PairSPHTaitwater::PairSPHTaitwater(LAMMPS *lmp) : Pair(lmp) /* ---------------------------------------------------------------------- */ -PairSPHTaitwater::~PairSPHTaitwater() { +PairSPHTaitwater::~PairSPHTaitwater() +{ if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); @@ -51,7 +56,8 @@ PairSPHTaitwater::~PairSPHTaitwater() { /* ---------------------------------------------------------------------- */ -void PairSPHTaitwater::compute(int eflag, int vflag) { +void PairSPHTaitwater::compute(int eflag, int vflag) +{ int i, j, ii, jj, inum, jnum, itype, jtype; double xtmp, ytmp, ztmp, delx, dely, delz, fpair; @@ -72,25 +78,6 @@ void PairSPHTaitwater::compute(int eflag, int vflag) { int nlocal = atom->nlocal; int newton_pair = force->newton_pair; - // check consistency of pair coefficients - - if (first) { - for (i = 1; i <= atom->ntypes; i++) { - for (j = 1; i <= atom->ntypes; i++) { - if (cutsq[i][j] > 1.e-32) { - if (!setflag[i][i] || !setflag[j][j]) { - if (comm->me == 0) { - printf( - "SPH particle types %d and %d interact with cutoff=%g, but not all of their single particle properties are set.\n", - i, j, sqrt(cutsq[i][j])); - } - } - } - } - } - first = 0; - } - inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -201,7 +188,8 @@ void PairSPHTaitwater::compute(int eflag, int vflag) { allocate all arrays ------------------------------------------------------------------------- */ -void PairSPHTaitwater::allocate() { +void PairSPHTaitwater::allocate() +{ allocated = 1; int n = atom->ntypes; @@ -223,7 +211,8 @@ void PairSPHTaitwater::allocate() { global settings ------------------------------------------------------------------------- */ -void PairSPHTaitwater::settings(int narg, char **/*arg*/) { +void PairSPHTaitwater::settings(int narg, char **/*arg*/) +{ if (narg != 0) error->all(FLERR, "Illegal number of arguments for pair_style sph/taitwater"); @@ -233,7 +222,8 @@ void PairSPHTaitwater::settings(int narg, char **/*arg*/) { set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairSPHTaitwater::coeff(int narg, char **arg) { +void PairSPHTaitwater::coeff(int narg, char **arg) +{ if (narg != 6) error->all(FLERR, "Incorrect args for pair_style sph/taitwater coefficients"); @@ -257,14 +247,8 @@ void PairSPHTaitwater::coeff(int narg, char **arg) { B[i] = B_one; for (int j = MAX(jlo,i); j <= jhi; j++) { viscosity[i][j] = viscosity_one; - //printf("setting cut[%d][%d] = %f\n", i, j, cut_one); cut[i][j] = cut_one; - setflag[i][j] = 1; - - //cut[j][i] = cut[i][j]; - //viscosity[j][i] = viscosity[i][j]; - //setflag[j][i] = 1; count++; } } @@ -277,8 +261,8 @@ void PairSPHTaitwater::coeff(int narg, char **arg) { init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ -double PairSPHTaitwater::init_one(int i, int j) { - +double PairSPHTaitwater::init_one(int i, int j) +{ if (setflag[i][j] == 0) { error->all(FLERR,"All pair sph/taitwater coeffs are set"); } diff --git a/src/SPH/pair_sph_taitwater_morris.cpp b/src/SPH/pair_sph_taitwater_morris.cpp index 50fcd270f6..2209cc4ecf 100644 --- a/src/SPH/pair_sph_taitwater_morris.cpp +++ b/src/SPH/pair_sph_taitwater_morris.cpp @@ -13,15 +13,16 @@ ------------------------------------------------------------------------- */ #include "pair_sph_taitwater_morris.h" -#include -#include "atom.h" -#include "force.h" -#include "comm.h" -#include "neigh_list.h" -#include "memory.h" -#include "error.h" -#include "domain.h" +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" + +#include using namespace LAMMPS_NS; @@ -29,6 +30,9 @@ using namespace LAMMPS_NS; PairSPHTaitwaterMorris::PairSPHTaitwaterMorris(LAMMPS *lmp) : Pair(lmp) { + if ((atom->esph_flag != 1) || (atom->rho_flag != 1) || (atom->vest_flag != 1)) + error->all(FLERR, "Pair sph/taitwater/morris requires atom attributes energy, density, and velocity estimates, e.g. in atom_style sph"); + restartinfo = 0; first = 1; single_enable = 0; @@ -36,7 +40,8 @@ PairSPHTaitwaterMorris::PairSPHTaitwaterMorris(LAMMPS *lmp) : Pair(lmp) /* ---------------------------------------------------------------------- */ -PairSPHTaitwaterMorris::~PairSPHTaitwaterMorris() { +PairSPHTaitwaterMorris::~PairSPHTaitwaterMorris() +{ if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); @@ -51,7 +56,8 @@ PairSPHTaitwaterMorris::~PairSPHTaitwaterMorris() { /* ---------------------------------------------------------------------- */ -void PairSPHTaitwaterMorris::compute(int eflag, int vflag) { +void PairSPHTaitwaterMorris::compute(int eflag, int vflag) +{ int i, j, ii, jj, inum, jnum, itype, jtype; double xtmp, ytmp, ztmp, delx, dely, delz, fpair; @@ -72,25 +78,6 @@ void PairSPHTaitwaterMorris::compute(int eflag, int vflag) { int nlocal = atom->nlocal; int newton_pair = force->newton_pair; - // check consistency of pair coefficients - - if (first) { - for (i = 1; i <= atom->ntypes; i++) { - for (j = 1; i <= atom->ntypes; i++) { - if (cutsq[i][j] > 1.e-32) { - if (!setflag[i][i] || !setflag[j][j]) { - if (comm->me == 0) { - printf( - "SPH particle types %d and %d interact with cutoff=%g, but not all of their single particle properties are set.\n", - i, j, sqrt(cutsq[i][j])); - } - } - } - } - } - first = 0; - } - inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -152,9 +139,9 @@ void PairSPHTaitwaterMorris::compute(int eflag, int vflag) { fj = tmp * tmp * tmp; fj = B[jtype] * (fj * fj * tmp - 1.0) / (rho[j] * rho[j]); - velx=vxtmp - v[j][0]; - vely=vytmp - v[j][1]; - velz=vztmp - v[j][2]; + velx = vxtmp - v[j][0]; + vely = vytmp - v[j][1]; + velz = vztmp - v[j][2]; // dot product of velocity delta and distance vector delVdotDelR = delx * velx + dely * vely + delz * velz; @@ -169,8 +156,6 @@ void PairSPHTaitwaterMorris::compute(int eflag, int vflag) { fpair = -imass * jmass * (fi + fj) * wfd; deltaE = -0.5 *(fpair * delVdotDelR + fvisc * (velx*velx + vely*vely + velz*velz)); - // printf("testvar= %f, %f \n", delx, dely); - f[i][0] += delx * fpair + velx * fvisc; f[i][1] += dely * fpair + vely * fvisc; f[i][2] += delz * fpair + velz * fvisc; @@ -189,6 +174,7 @@ void PairSPHTaitwaterMorris::compute(int eflag, int vflag) { drho[j] += imass * delVdotDelR * wfd; } + // viscous forces do not contribute to virial if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz); } @@ -202,7 +188,8 @@ void PairSPHTaitwaterMorris::compute(int eflag, int vflag) { allocate all arrays ------------------------------------------------------------------------- */ -void PairSPHTaitwaterMorris::allocate() { +void PairSPHTaitwaterMorris::allocate() +{ allocated = 1; int n = atom->ntypes; @@ -224,7 +211,8 @@ void PairSPHTaitwaterMorris::allocate() { global settings ------------------------------------------------------------------------- */ -void PairSPHTaitwaterMorris::settings(int narg, char **/*arg*/) { +void PairSPHTaitwaterMorris::settings(int narg, char **/*arg*/) +{ if (narg != 0) error->all(FLERR, "Illegal number of arguments for pair_style sph/taitwater/morris"); @@ -234,7 +222,8 @@ void PairSPHTaitwaterMorris::settings(int narg, char **/*arg*/) { set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairSPHTaitwaterMorris::coeff(int narg, char **arg) { +void PairSPHTaitwaterMorris::coeff(int narg, char **arg) +{ if (narg != 6) error->all(FLERR, "Incorrect args for pair_style sph/taitwater/morris coefficients"); @@ -258,7 +247,6 @@ void PairSPHTaitwaterMorris::coeff(int narg, char **arg) { B[i] = B_one; for (int j = MAX(jlo,i); j <= jhi; j++) { viscosity[i][j] = viscosity_one; - //printf("setting cut[%d][%d] = %f\n", i, j, cut_one); cut[i][j] = cut_one; setflag[i][j] = 1; @@ -274,8 +262,8 @@ void PairSPHTaitwaterMorris::coeff(int narg, char **arg) { init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ -double PairSPHTaitwaterMorris::init_one(int i, int j) { - +double PairSPHTaitwaterMorris::init_one(int i, int j) +{ if (setflag[i][j] == 0) { error->all(FLERR,"All pair sph/taitwater/morris coeffs are not set"); } From 6abbbdba6c98f9ba77b56dc5820053dc6ec01310 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 25 Jul 2024 17:10:23 -0600 Subject: [PATCH 4/9] Patching vremap in SPH --- src/SPH/fix_sph.cpp | 42 +++++++++++++++++++++++++++++++++++++++--- src/SPH/fix_sph.h | 2 ++ 2 files changed, 41 insertions(+), 3 deletions(-) diff --git a/src/SPH/fix_sph.cpp b/src/SPH/fix_sph.cpp index 96479be928..876d7ac4fb 100644 --- a/src/SPH/fix_sph.cpp +++ b/src/SPH/fix_sph.cpp @@ -61,6 +61,10 @@ void FixSPH::init() void FixSPH::setup_pre_force(int /*vflag*/) { + remap_v_flag = domain->deform_vremap; + if (remap_v_flag && (!comm->ghost_velocity)) + error->all(FLERR, "Fix sph requires ghost atoms store velocity when deforming with remap v"); + // set vest equal to v double **v = atom->v; double **vest = atom->vest; @@ -119,9 +123,16 @@ void FixSPH::initial_integrate(int /*vflag*/) rho[i] += dtf * drho[i]; // ... and density // extrapolate velocity for use with velocity-dependent potentials, e.g. SPH - vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0]; - vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1]; - vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2]; + // if velocities are remapped, perform this extrapolation after communication + if (remap_v_flag) { + vest[i][0] = dtfm * f[i][0]; + vest[i][1] = dtfm * f[i][1]; + vest[i][2] = dtfm * f[i][2]; + } else { + vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0]; + vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1]; + vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2]; + } v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; @@ -136,6 +147,31 @@ void FixSPH::initial_integrate(int /*vflag*/) /* ---------------------------------------------------------------------- */ +void FixSPH::pre_force(int /*vflag*/) +{ + // if velocities are remapped, calculate estimates here + // note that vest currently stores dtfm * force + if (!remap_v_flag) return; + + double **v = atom->v; + double **vest = atom->vest; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) + nlocal = atom->nfirst; + + int nall = nlocal + atom->nghost; + for (int i = 0; i < nall; i++) { + if (mask[i] & groupbit) { + vest[i][0] += v[i][0]; + vest[i][1] += v[i][1]; + vest[i][2] += v[i][2]; + } + } +} + +/* ---------------------------------------------------------------------- */ + void FixSPH::final_integrate() { // update v, rho, and e of atoms in group diff --git a/src/SPH/fix_sph.h b/src/SPH/fix_sph.h index 415e38bd21..7676844dd9 100644 --- a/src/SPH/fix_sph.h +++ b/src/SPH/fix_sph.h @@ -30,6 +30,7 @@ class FixSPH : public Fix { int setmask() override; void init() override; void setup_pre_force(int) override; + void pre_force(int) override; void initial_integrate(int) override; void final_integrate() override; void reset_dt() override; @@ -41,6 +42,7 @@ class FixSPH : public Fix { double dtv, dtf; double *step_respa; int mass_require; + int remap_v_flag; class Pair *pair; }; From 49b377fc3db56bac056cff1009c55c896c23f644 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 25 Jul 2024 17:51:48 -0600 Subject: [PATCH 5/9] Adding errors to unpatched uses of vest with vremap --- src/DPD-MESO/fix_mvv_dpd.cpp | 5 +++++ src/DPD-MESO/fix_mvv_edpd.cpp | 5 +++++ src/DPD-MESO/fix_mvv_tdpd.cpp | 5 +++++ src/DPD-SMOOTH/fix_meso_move.cpp | 7 +++++++ src/DPD-SMOOTH/fix_rigid_meso.cpp | 5 +++++ src/MACHDYN/fix_smd_integrate_tlsph.cpp | 5 +++++ src/MACHDYN/fix_smd_integrate_ulsph.cpp | 5 +++++ 7 files changed, 37 insertions(+) diff --git a/src/DPD-MESO/fix_mvv_dpd.cpp b/src/DPD-MESO/fix_mvv_dpd.cpp index d51000b15b..6def3305d8 100644 --- a/src/DPD-MESO/fix_mvv_dpd.cpp +++ b/src/DPD-MESO/fix_mvv_dpd.cpp @@ -65,6 +65,11 @@ void FixMvvDPD::init() if (!atom->vest_flag) error->all(FLERR,"Fix mvv/dpd requires atom attribute vest e.g. from atom style mdpd"); + // Cannot use vremap since its effects aren't propagated to vest + // see RHEO or SPH packages for examples patches + if (domain->deform_vremap) + error->all(FLERR, "Fix mvv/dpd cannot be used with velocity remapping"); + if (!force->pair_match("^mdpd",0) && !force->pair_match("^dpd",0)) { if (force->pair_match("^hybrid",0)) { if (!(force->pair_match("^mdpd",0,1) || force->pair_match("^dpd",0),1)) { diff --git a/src/DPD-MESO/fix_mvv_edpd.cpp b/src/DPD-MESO/fix_mvv_edpd.cpp index 7ac0dc3de7..389a8c97bf 100644 --- a/src/DPD-MESO/fix_mvv_edpd.cpp +++ b/src/DPD-MESO/fix_mvv_edpd.cpp @@ -73,6 +73,11 @@ void FixMvvEDPD::init() { if (!atom->edpd_flag) error->all(FLERR,"Fix mvv/edpd requires atom style edpd"); + // Cannot use vremap since its effects aren't propagated to vest + // see RHEO or SPH packages for examples patches + if (domain->deform_vremap) + error->all(FLERR, "Fix mvv/edpd cannot be used with velocity remapping"); + if (!force->pair_match("^edpd",0)) { if (force->pair_match("^hybrid",0)) { if (!force->pair_match("^edpd",0,1)) { diff --git a/src/DPD-MESO/fix_mvv_tdpd.cpp b/src/DPD-MESO/fix_mvv_tdpd.cpp index 122fd56365..1babbb8c5f 100644 --- a/src/DPD-MESO/fix_mvv_tdpd.cpp +++ b/src/DPD-MESO/fix_mvv_tdpd.cpp @@ -71,6 +71,11 @@ void FixMvvTDPD::init() { if (!atom->tdpd_flag) error->all(FLERR,"Fix mvv/tdpd requires atom style tdpd"); + // Cannot use vremap since its effects aren't propagated to vest + // see RHEO or SPH packages for examples patches + if (domain->deform_vremap) + error->all(FLERR, "Fix mvv/tdpd cannot be used with velocity remapping"); + if (!force->pair_match("^tdpd",0)) { if (force->pair_match("^hybrid",0)) { if (!force->pair_match("^tdpd",0,1)) { diff --git a/src/DPD-SMOOTH/fix_meso_move.cpp b/src/DPD-SMOOTH/fix_meso_move.cpp index 2c3213c3cd..89a7814f6a 100644 --- a/src/DPD-SMOOTH/fix_meso_move.cpp +++ b/src/DPD-SMOOTH/fix_meso_move.cpp @@ -350,7 +350,14 @@ void FixMesoMove::init () { } void FixMesoMove::setup_pre_force (int /*vflag*/) { + + // Cannot use vremap since its effects aren't propagated to vest + // see RHEO or SPH packages for examples patches + if (domain->deform_vremap) + error->all(FLERR, "Fix meso/move cannot be used with velocity remapping"); + // set vest equal to v + double **v = atom->v; double **vest = atom->vest; int *mask = atom->mask; diff --git a/src/DPD-SMOOTH/fix_rigid_meso.cpp b/src/DPD-SMOOTH/fix_rigid_meso.cpp index e0ad501e02..bfceb3295e 100644 --- a/src/DPD-SMOOTH/fix_rigid_meso.cpp +++ b/src/DPD-SMOOTH/fix_rigid_meso.cpp @@ -92,6 +92,11 @@ void FixRigidMeso::setup (int vflag) { conjqm[ibody][2] *= 2.0; conjqm[ibody][3] *= 2.0; } + + // Cannot use vremap since its effects aren't propagated to vest + // see RHEO or SPH packages for examples patches + if (domain->deform_vremap) + error->all(FLERR, "Fix rigid/meso cannot be used with velocity remapping"); } /* ---------------------------------------------------------------------- diff --git a/src/MACHDYN/fix_smd_integrate_tlsph.cpp b/src/MACHDYN/fix_smd_integrate_tlsph.cpp index f138a3c387..97cce5524d 100644 --- a/src/MACHDYN/fix_smd_integrate_tlsph.cpp +++ b/src/MACHDYN/fix_smd_integrate_tlsph.cpp @@ -115,6 +115,11 @@ void FixSMDIntegrateTlsph::init() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; vlimitsq = vlimit * vlimit; + + // Cannot use vremap since its effects aren't propagated to vest + // see RHEO or SPH packages for examples patches + if (domain->deform_vremap) + error->all(FLERR, "Fix smd/integrate_tlsph cannot be used with velocity remapping"); } /* ---------------------------------------------------------------------- diff --git a/src/MACHDYN/fix_smd_integrate_ulsph.cpp b/src/MACHDYN/fix_smd_integrate_ulsph.cpp index 8ef5c65b67..395946ba20 100644 --- a/src/MACHDYN/fix_smd_integrate_ulsph.cpp +++ b/src/MACHDYN/fix_smd_integrate_ulsph.cpp @@ -144,6 +144,11 @@ void FixSMDIntegrateUlsph::init() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; vlimitsq = vlimit * vlimit; + + // Cannot use vremap since its effects aren't propagated to vest + // see RHEO or SPH packages for examples patches + if (domain->deform_vremap) + error->all(FLERR, "Fix smd/integrate_ulsph cannot be used with velocity remapping"); } /* ---------------------------------------------------------------------- From d43e87bce1acfb33b5b4bb913c5bfc89a0039b1a Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 25 Jul 2024 17:54:52 -0600 Subject: [PATCH 6/9] Missing word --- src/DPD-MESO/fix_mvv_dpd.cpp | 2 +- src/DPD-MESO/fix_mvv_edpd.cpp | 2 +- src/DPD-MESO/fix_mvv_tdpd.cpp | 2 +- src/DPD-SMOOTH/fix_meso_move.cpp | 2 +- src/DPD-SMOOTH/fix_rigid_meso.cpp | 2 +- src/MACHDYN/fix_smd_integrate_tlsph.cpp | 4 ++-- src/MACHDYN/fix_smd_integrate_ulsph.cpp | 4 ++-- 7 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/DPD-MESO/fix_mvv_dpd.cpp b/src/DPD-MESO/fix_mvv_dpd.cpp index 6def3305d8..5e91b4bad9 100644 --- a/src/DPD-MESO/fix_mvv_dpd.cpp +++ b/src/DPD-MESO/fix_mvv_dpd.cpp @@ -66,7 +66,7 @@ void FixMvvDPD::init() error->all(FLERR,"Fix mvv/dpd requires atom attribute vest e.g. from atom style mdpd"); // Cannot use vremap since its effects aren't propagated to vest - // see RHEO or SPH packages for examples patches + // see RHEO or SPH packages for examples of patches if (domain->deform_vremap) error->all(FLERR, "Fix mvv/dpd cannot be used with velocity remapping"); diff --git a/src/DPD-MESO/fix_mvv_edpd.cpp b/src/DPD-MESO/fix_mvv_edpd.cpp index 389a8c97bf..d34e9b5e5a 100644 --- a/src/DPD-MESO/fix_mvv_edpd.cpp +++ b/src/DPD-MESO/fix_mvv_edpd.cpp @@ -74,7 +74,7 @@ void FixMvvEDPD::init() if (!atom->edpd_flag) error->all(FLERR,"Fix mvv/edpd requires atom style edpd"); // Cannot use vremap since its effects aren't propagated to vest - // see RHEO or SPH packages for examples patches + // see RHEO or SPH packages for examples of patches if (domain->deform_vremap) error->all(FLERR, "Fix mvv/edpd cannot be used with velocity remapping"); diff --git a/src/DPD-MESO/fix_mvv_tdpd.cpp b/src/DPD-MESO/fix_mvv_tdpd.cpp index 1babbb8c5f..a34c86b6dd 100644 --- a/src/DPD-MESO/fix_mvv_tdpd.cpp +++ b/src/DPD-MESO/fix_mvv_tdpd.cpp @@ -72,7 +72,7 @@ void FixMvvTDPD::init() if (!atom->tdpd_flag) error->all(FLERR,"Fix mvv/tdpd requires atom style tdpd"); // Cannot use vremap since its effects aren't propagated to vest - // see RHEO or SPH packages for examples patches + // see RHEO or SPH packages for examples of patches if (domain->deform_vremap) error->all(FLERR, "Fix mvv/tdpd cannot be used with velocity remapping"); diff --git a/src/DPD-SMOOTH/fix_meso_move.cpp b/src/DPD-SMOOTH/fix_meso_move.cpp index 89a7814f6a..c9d4a41400 100644 --- a/src/DPD-SMOOTH/fix_meso_move.cpp +++ b/src/DPD-SMOOTH/fix_meso_move.cpp @@ -352,7 +352,7 @@ void FixMesoMove::init () { void FixMesoMove::setup_pre_force (int /*vflag*/) { // Cannot use vremap since its effects aren't propagated to vest - // see RHEO or SPH packages for examples patches + // see RHEO or SPH packages for examples of patches if (domain->deform_vremap) error->all(FLERR, "Fix meso/move cannot be used with velocity remapping"); diff --git a/src/DPD-SMOOTH/fix_rigid_meso.cpp b/src/DPD-SMOOTH/fix_rigid_meso.cpp index bfceb3295e..54fd72f924 100644 --- a/src/DPD-SMOOTH/fix_rigid_meso.cpp +++ b/src/DPD-SMOOTH/fix_rigid_meso.cpp @@ -94,7 +94,7 @@ void FixRigidMeso::setup (int vflag) { } // Cannot use vremap since its effects aren't propagated to vest - // see RHEO or SPH packages for examples patches + // see RHEO or SPH packages for examples of patches if (domain->deform_vremap) error->all(FLERR, "Fix rigid/meso cannot be used with velocity remapping"); } diff --git a/src/MACHDYN/fix_smd_integrate_tlsph.cpp b/src/MACHDYN/fix_smd_integrate_tlsph.cpp index 97cce5524d..e00c88bc22 100644 --- a/src/MACHDYN/fix_smd_integrate_tlsph.cpp +++ b/src/MACHDYN/fix_smd_integrate_tlsph.cpp @@ -116,8 +116,8 @@ void FixSMDIntegrateTlsph::init() { dtf = 0.5 * update->dt * force->ftm2v; vlimitsq = vlimit * vlimit; - // Cannot use vremap since its effects aren't propagated to vest - // see RHEO or SPH packages for examples patches + // Cannot use vremap since its effects aren't propagated to vest + // see RHEO or SPH packages for examples of patches if (domain->deform_vremap) error->all(FLERR, "Fix smd/integrate_tlsph cannot be used with velocity remapping"); } diff --git a/src/MACHDYN/fix_smd_integrate_ulsph.cpp b/src/MACHDYN/fix_smd_integrate_ulsph.cpp index 395946ba20..c6530c989b 100644 --- a/src/MACHDYN/fix_smd_integrate_ulsph.cpp +++ b/src/MACHDYN/fix_smd_integrate_ulsph.cpp @@ -145,8 +145,8 @@ void FixSMDIntegrateUlsph::init() { dtf = 0.5 * update->dt * force->ftm2v; vlimitsq = vlimit * vlimit; - // Cannot use vremap since its effects aren't propagated to vest - // see RHEO or SPH packages for examples patches + // Cannot use vremap since its effects aren't propagated to vest + // see RHEO or SPH packages for examples of patches if (domain->deform_vremap) error->all(FLERR, "Fix smd/integrate_ulsph cannot be used with velocity remapping"); } From 3deffb0dfd96fb7b57d735762da56bc6d71694ad Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 25 Jul 2024 18:10:27 -0600 Subject: [PATCH 7/9] typo --- src/SPH/compute_sph_e_atom.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SPH/compute_sph_e_atom.cpp b/src/SPH/compute_sph_e_atom.cpp index 8bb0622acd..b98fd97479 100644 --- a/src/SPH/compute_sph_e_atom.cpp +++ b/src/SPH/compute_sph_e_atom.cpp @@ -33,7 +33,7 @@ ComputeSPHEAtom::ComputeSPHEAtom(LAMMPS *lmp, int narg, char **arg) : if (narg != 3) error->all(FLERR,"Number of arguments for compute sph/e/atom command != 3"); if (atom->esph_flag != 1) - error->all(FLERR,"Compute sph/e/atom requires atom attribut energy, e.g. in atom_style sph"); + error->all(FLERR,"Compute sph/e/atom requires atom attribute energy, e.g. in atom_style sph"); peratom_flag = 1; size_peratom_cols = 0; From 68a6bc069330a2e6583af78e67f93be70e1c4bcf Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 25 Jul 2024 18:15:34 -0600 Subject: [PATCH 8/9] Adding missing header file --- src/DPD-MESO/fix_mvv_dpd.cpp | 1 + src/DPD-MESO/fix_mvv_edpd.cpp | 1 + src/DPD-MESO/fix_mvv_tdpd.cpp | 1 + src/DPD-SMOOTH/fix_rigid_meso.cpp | 3 ++- src/MACHDYN/fix_smd_integrate_tlsph.cpp | 15 +++++++++------ src/MACHDYN/fix_smd_integrate_ulsph.cpp | 15 +++++++++------ 6 files changed, 23 insertions(+), 13 deletions(-) diff --git a/src/DPD-MESO/fix_mvv_dpd.cpp b/src/DPD-MESO/fix_mvv_dpd.cpp index 5e91b4bad9..6b5440902b 100644 --- a/src/DPD-MESO/fix_mvv_dpd.cpp +++ b/src/DPD-MESO/fix_mvv_dpd.cpp @@ -24,6 +24,7 @@ #include "fix_mvv_dpd.h" #include "atom.h" +#include "domain.h" #include "error.h" #include "force.h" #include "update.h" diff --git a/src/DPD-MESO/fix_mvv_edpd.cpp b/src/DPD-MESO/fix_mvv_edpd.cpp index d34e9b5e5a..32a1608d6e 100644 --- a/src/DPD-MESO/fix_mvv_edpd.cpp +++ b/src/DPD-MESO/fix_mvv_edpd.cpp @@ -33,6 +33,7 @@ #include "fix_mvv_edpd.h" #include "atom.h" +#include "domain.h" #include "error.h" #include "force.h" #include "update.h" diff --git a/src/DPD-MESO/fix_mvv_tdpd.cpp b/src/DPD-MESO/fix_mvv_tdpd.cpp index a34c86b6dd..d4e2d8a5c9 100644 --- a/src/DPD-MESO/fix_mvv_tdpd.cpp +++ b/src/DPD-MESO/fix_mvv_tdpd.cpp @@ -29,6 +29,7 @@ #include "fix_mvv_tdpd.h" #include "atom.h" +#include "domain.h" #include "error.h" #include "force.h" #include "update.h" diff --git a/src/DPD-SMOOTH/fix_rigid_meso.cpp b/src/DPD-SMOOTH/fix_rigid_meso.cpp index 54fd72f924..3e8b07e6dc 100644 --- a/src/DPD-SMOOTH/fix_rigid_meso.cpp +++ b/src/DPD-SMOOTH/fix_rigid_meso.cpp @@ -29,11 +29,12 @@ ------------------------------------------------------------------------- */ #include "fix_rigid_meso.h" + #include "math_extra.h" #include "atom.h" #include "domain.h" -#include "memory.h" #include "error.h" +#include "memory.h" using namespace LAMMPS_NS; using namespace FixConst; diff --git a/src/MACHDYN/fix_smd_integrate_tlsph.cpp b/src/MACHDYN/fix_smd_integrate_tlsph.cpp index e00c88bc22..ee617c1577 100644 --- a/src/MACHDYN/fix_smd_integrate_tlsph.cpp +++ b/src/MACHDYN/fix_smd_integrate_tlsph.cpp @@ -24,15 +24,18 @@ ------------------------------------------------------------------------- */ #include "fix_smd_integrate_tlsph.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "pair.h" +#include "update.h" + #include #include #include -#include "atom.h" -#include "force.h" -#include "update.h" -#include "error.h" -#include "pair.h" -#include "comm.h" using namespace Eigen; using namespace LAMMPS_NS; diff --git a/src/MACHDYN/fix_smd_integrate_ulsph.cpp b/src/MACHDYN/fix_smd_integrate_ulsph.cpp index c6530c989b..6e2934c2ec 100644 --- a/src/MACHDYN/fix_smd_integrate_ulsph.cpp +++ b/src/MACHDYN/fix_smd_integrate_ulsph.cpp @@ -24,15 +24,18 @@ ------------------------------------------------------------------------- */ #include "fix_smd_integrate_ulsph.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "pair.h" +#include "update.h" + #include #include #include -#include "atom.h" -#include "comm.h" -#include "force.h" -#include "update.h" -#include "error.h" -#include "pair.h" using namespace Eigen; using namespace LAMMPS_NS; From a972a282a792e4ad6cfd23e746cf7289b4c59b80 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 26 Jul 2024 17:43:04 -0600 Subject: [PATCH 9/9] Adding remap v restriction to doc files --- doc/src/fix_meso_move.rst | 5 +++++ doc/src/fix_mvv_dpd.rst | 5 +++++ doc/src/fix_rigid_meso.rst | 5 +++++ doc/src/fix_smd_integrate_tlsph.rst | 5 +++++ doc/src/fix_smd_integrate_ulsph.rst | 5 +++++ src/DPD-SMOOTH/fix_rigid_meso.cpp | 2 +- 6 files changed, 26 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_meso_move.rst b/doc/src/fix_meso_move.rst index 55d54b2107..281a405ab0 100644 --- a/doc/src/fix_meso_move.rst +++ b/doc/src/fix_meso_move.rst @@ -247,6 +247,11 @@ defined by the :doc:`atom_style sph ` command. All particles in the group must be mesoscopic SPH/SDPD particles. +versionchanged:: TBD + +This fix is incompatible with deformation controls that remap velocity, +for instance the *remap v* option of :doc:`fix deform `. + Related commands """""""""""""""" diff --git a/doc/src/fix_mvv_dpd.rst b/doc/src/fix_mvv_dpd.rst index ff5b169f97..740785d562 100644 --- a/doc/src/fix_mvv_dpd.rst +++ b/doc/src/fix_mvv_dpd.rst @@ -97,6 +97,11 @@ These fixes are part of the DPD-MESO package. They are only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +versionchanged:: TBD + +This fix is incompatible with deformation controls that remap velocity, +for instance the *remap v* option of :doc:`fix deform `. + Related commands """""""""""""""" diff --git a/doc/src/fix_rigid_meso.rst b/doc/src/fix_rigid_meso.rst index c8c994fd26..933369788e 100644 --- a/doc/src/fix_rigid_meso.rst +++ b/doc/src/fix_rigid_meso.rst @@ -353,6 +353,11 @@ defined by the :doc:`atom_style sph ` command. All particles in the group must be mesoscopic SPH/SDPD particles. +versionchanged:: TBD + +This fix is incompatible with deformation controls that remap velocity, +for instance the *remap v* option of :doc:`fix deform `. + Related commands """""""""""""""" diff --git a/doc/src/fix_smd_integrate_tlsph.rst b/doc/src/fix_smd_integrate_tlsph.rst index 515b30d4ff..e714c91a17 100644 --- a/doc/src/fix_smd_integrate_tlsph.rst +++ b/doc/src/fix_smd_integrate_tlsph.rst @@ -53,6 +53,11 @@ Restrictions This fix is part of the MACHDYN package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +versionchanged:: TBD + +This fix is incompatible with deformation controls that remap velocity, +for instance the *remap v* option of :doc:`fix deform `. + Related commands """""""""""""""" diff --git a/doc/src/fix_smd_integrate_ulsph.rst b/doc/src/fix_smd_integrate_ulsph.rst index 17dfdb7b17..60c185db5f 100644 --- a/doc/src/fix_smd_integrate_ulsph.rst +++ b/doc/src/fix_smd_integrate_ulsph.rst @@ -61,6 +61,11 @@ Restrictions This fix is part of the MACHDYN package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +versionchanged:: TBD + +This fix is incompatible with deformation controls that remap velocity, +for instance the *remap v* option of :doc:`fix deform `. + Related commands """""""""""""""" diff --git a/src/DPD-SMOOTH/fix_rigid_meso.cpp b/src/DPD-SMOOTH/fix_rigid_meso.cpp index 3e8b07e6dc..38b9d0a09c 100644 --- a/src/DPD-SMOOTH/fix_rigid_meso.cpp +++ b/src/DPD-SMOOTH/fix_rigid_meso.cpp @@ -30,10 +30,10 @@ #include "fix_rigid_meso.h" -#include "math_extra.h" #include "atom.h" #include "domain.h" #include "error.h" +#include "math_extra.h" #include "memory.h" using namespace LAMMPS_NS;