diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index 3afeb03e43..615d5ea740 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -103,7 +103,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a if (avec_index[i] < 0) error->all(FLERR, "Invalid keyword {} for atom style {} in compute rheo/property/atom command ", - atom->get_style(), arg[iarg]); + arg[iarg], atom->get_style()); pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style; if (strcmp(arg[iarg],"temperature") == 0) thermal_flag = 1; diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp index 88b33c09af..882977f3e9 100644 --- a/src/RHEO/compute_rheo_surface.cpp +++ b/src/RHEO/compute_rheo_surface.cpp @@ -152,16 +152,12 @@ void ComputeRHEOSurface::compute_peratom() grow_arrays(atom->nmax); size_t nbytes = nmax_store * sizeof(double); - memset(&divr, 0, nbytes); - memset(&rsurface, 0, nbytes); - memset(&nsurface, 0, 3 * nbytes); - memset(&gradC, 0, 3 * 3 * nbytes); - memset(&B, 0, 3 * 3 * nbytes); + memset(&divr[0], 0, nbytes); + memset(&rsurface[0], 0, nbytes); + memset(&nsurface[0][0], 0, dim * nbytes); + memset(&gradC[0][0], 0, dim * dim * nbytes); + memset(&B[0][0], 0, dim * dim * nbytes); - // Remove surface settings - int nall = nlocal + atom->nghost; - for (i = 0; i < nall; i++) - status[i] &= SURFACEMASK; // loop over neighbors to calculate the average orientation of neighbors for (ii = 0; ii < inum; ii++) { @@ -208,7 +204,7 @@ void ComputeRHEOSurface::compute_peratom() wp = compute_kernel->calc_dw_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq), dWij, dWji); - for (a = 0; a < dim; a++){ + for (a = 0; a < dim; a++) { divr[i] -= dWij[a] * dx[a] * Volj; gradC[i][a] += dWij[a] * Volj; } @@ -245,35 +241,32 @@ void ComputeRHEOSurface::compute_peratom() } } - // Find the free-surface - if (threshold_style == DIVR) { - for (i = 0; i < nall; i++) { - if (mask[i] & groupbit) { + // Remove surface settings and assign new values + int nall = nlocal + atom->nghost; + int test; + + for (i = 0; i < nall; i++) { + status[i] &= SURFACEMASK; + if (mask[i] & groupbit) { + if (threshold_style == DIVR) + test = divr[i] < threshold_divr; + else + test = coordination[i] < threshold_z; + + if (test) { + if (coordination[i] < threshold_splash) + status[i] |= STATUS_SPLASH; + else + status[i] |= STATUS_SURFACE; + rsurface[i] = 0.0; + } else { status[i] |= STATUS_BULK; rsurface[i] = cut; - if (divr[i] < threshold_divr) { - status[i] |= STATUS_SURFACE; - rsurface[i] = 0.0; - if (coordination[i] < threshold_splash) - status[i] |= STATUS_SPLASH; - } - } - } - } else { - for (i = 0; i < nall; i++) { - if (mask[i] & groupbit) { - status[i] |= STATUS_BULK; - rsurface[i] = cut; - if (coordination[i] < threshold_z) { - status[i] |= STATUS_SURFACE; - rsurface[i] = 0.0; - if (coordination[i] < threshold_splash) - status[i] |= STATUS_SPLASH; - } } } } + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; @@ -297,7 +290,8 @@ void ComputeRHEOSurface::compute_peratom() status[i] |= STATUS_LAYER; } - if (status[j] & STATUS_SURFACE) rsurface[i] = MIN(rsurface[i], sqrt(rsq)); + if (status[j] & STATUS_SURFACE) + rsurface[i] = MIN(rsurface[i], sqrt(rsq)); if (j < nlocal || newton) { @@ -306,7 +300,8 @@ void ComputeRHEOSurface::compute_peratom() status[j] |= STATUS_LAYER; } - if (status[i] & STATUS_SURFACE) rsurface[j] = MIN(rsurface[j], sqrt(rsq)); + if (status[i] & STATUS_SURFACE) + rsurface[j] = MIN(rsurface[j], sqrt(rsq)); } } } @@ -351,7 +346,8 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf) int i,a,b,k,j,m; int dim = domain->dimension; int *status = atom->status; - int temp; + int tmp1; + double tmp2; m = 0; for (i = 0; i < n; i++) { @@ -362,12 +358,13 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf) for (b = 0; b < dim; b ++) gradC[j][a * dim + b] += buf[m++]; } else if (comm_stage == 1) { - - temp = (int) buf[m++]; - if ((status[j] & STATUS_BULK) && (temp & STATUS_LAYER)) - status[j] = temp; - - rsurface[j] = MIN(rsurface[j], buf[m++]); + tmp1 = (int) buf[m++]; + if ((status[j] & STATUS_BULK) && (tmp1 & STATUS_LAYER)) { + status[j] &= SURFACEMASK; + status[j] |= STATUS_LAYER; + } + tmp2 = buf[m++]; + rsurface[j] = MIN(rsurface[j], tmp2); } } } diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 97b867179e..b1bc0fe0c0 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -133,7 +133,7 @@ FixRHEO::~FixRHEO() if (compute_kernel) modify->delete_compute("rheo_kernel"); if (compute_grad) modify->delete_compute("rheo_grad"); if (compute_interface) modify->delete_compute("rheo_interface"); - if (compute_surface) modify->delete_compute("compute_surface"); + if (compute_surface) modify->delete_compute("rheo_surface"); if (compute_rhosum) modify->delete_compute("rheo_rhosum"); if (compute_vshift) modify->delete_compute("rheo_vshift"); } @@ -211,7 +211,10 @@ void FixRHEO::setup_pre_force(int /*vflag*/) error->all(FLERR, "Fix RHEO must be defined before all other RHEO fixes"); // Calculate surfaces - if (surface_flag) compute_surface->compute_peratom(); + if (surface_flag) { + compute_kernel->compute_coordination(); + compute_surface->compute_peratom(); + } pre_force(0); } diff --git a/src/RHEO/pair_rheo_react.cpp b/src/RHEO/pair_rheo_react.cpp index 4846e077f6..d7412d5d0e 100644 --- a/src/RHEO/pair_rheo_react.cpp +++ b/src/RHEO/pair_rheo_react.cpp @@ -61,13 +61,13 @@ PairRHEOReact::PairRHEOReact(LAMMPS *lmp) : Pair(lmp), // between timesteps (fix property atom will handle callbacks) int tmp1, tmp2; - int index = atom->find_custom("react_nbond", tmp1, tmp2); - if (index == -1) { + index_nb = atom->find_custom("react_nbond", tmp1, tmp2); + if (index_nb == -1) { id_fix = utils::strdup("pair_rheo_react_fix_property_atom"); modify->add_fix(fmt::format("{} all property/atom i_react_nbond", id_fix)); - index = atom->find_custom("nbond", tmp1, tmp2); + index_nb = atom->find_custom("react_nbond", tmp1, tmp2); } - nbond = atom->ivector[index]; + nbond = atom->ivector[index_nb]; //Store non-persistent per atom quantities, intermediate @@ -79,9 +79,10 @@ PairRHEOReact::PairRHEOReact(LAMMPS *lmp) : Pair(lmp), PairRHEOReact::~PairRHEOReact() { - if (modify->nfix && fix_history) modify->delete_fix("NEIGH_HISTORY_RHEO_REACT"); - if (modify->nfix && fix_dummy) modify->delete_fix("NEIGH_HISTORY_RHEO_REACT_DUMMY"); - if (modify->nfix) modify->delete_fix("PROPERTY_ATOM_RHEO_REACT"); + if (modify->nfix && fix_history) modify->delete_fix("NEIGH_HISTORY_RHEO_REACT" + std::to_string(instance_me)); + if (modify->nfix && fix_dummy) modify->delete_fix("NEIGH_HISTORY_RHEO_REACT_DUMMY" + std::to_string(instance_me)); + if (modify->nfix) modify->delete_fix(id_fix); + delete[] id_fix; if (allocated) { memory->destroy(setflag); @@ -146,7 +147,7 @@ void PairRHEOReact::compute(int eflag, int vflag) } size_t nbytes = nmax_store * sizeof(int); - memset(&dbond, 0, nbytes); + memset(&dbond[0], 0, nbytes); // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { @@ -368,11 +369,11 @@ void PairRHEOReact::coeff(int narg, char **arg) void PairRHEOReact::init_style() { - int irequest = neighbor->request(this,instance_me); + int irequest = neighbor->request(this, instance_me); //neighbor->requests[irequest]->history = 1; if (fix_history == nullptr) { - auto cmd = fmt::format("NEIGH_HISTORY_RHEO_REACT {} all NEIGH_HISTORY {}", instance_me, size_history); + auto cmd = fmt::format("NEIGH_HISTORY_RHEO_REACT{} all NEIGH_HISTORY {}", instance_me, size_history); fix_history = dynamic_cast( modify->replace_fix("NEIGH_HISTORY_RHEO_REACT_DUMMY" + std::to_string(instance_me), cmd, 1)); fix_history->pair = this;