From 08d2dd26991dce65f9beee61039952bd3f81ca37 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 22 Dec 2023 10:46:44 -0700 Subject: [PATCH] Prototyping latent heat, other misc fixes --- src/.gitignore | 35 +-- src/RHEO/atom_vec_rheo_thermal.cpp | 46 ++-- src/RHEO/atom_vec_rheo_thermal.h | 2 +- src/RHEO/compute_rheo_grad.cpp | 69 +++-- src/RHEO/compute_rheo_grad.h | 4 +- src/RHEO/compute_rheo_kernel.cpp | 15 +- src/RHEO/compute_rheo_kernel.h | 1 + src/RHEO/compute_rheo_property_atom.cpp | 8 + src/RHEO/compute_rheo_vshift.cpp | 4 + src/RHEO/fix_rheo.cpp | 3 +- src/RHEO/fix_rheo_thermal.cpp | 139 ++++++++-- src/RHEO/fix_rheo_thermal.h | 6 +- src/RHEO/pair_rheo.cpp | 1 - src/RHEO/pair_rheo_solid.cpp | 351 ++++++++++++++++++++++++ src/RHEO/pair_rheo_solid.h | 51 ++++ src/set.cpp | 2 +- 16 files changed, 611 insertions(+), 126 deletions(-) create mode 100644 src/RHEO/pair_rheo_solid.cpp create mode 100644 src/RHEO/pair_rheo_solid.h diff --git a/src/.gitignore b/src/.gitignore index 1e634782dc..ccf8072922 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -197,40 +197,7 @@ /pair_tdpd.cpp /pair_tdpd.h -/atom_vec_rheo.cpp -/atom_vec_rheo.h -/atom_vec_rheo_thermal.cpp -/atom_vec_rheo_thermal.h -/compute_rheo_grad.cpp -/compute_rheo_grad.h -/compute_rheo_interface.cpp -/compute_rheo_interface.h -/compute_rheo_kernel.cpp -/compute_rheo_kernel.h -/compute_rheo_property_atom.cpp -/compute_rheo_property_atom.h -/compute_rheo_rho_sum.cpp -/compute_rheo_rho_sum.h -/compute_rheo_surface.cpp -/compute_rheo_surface.h -/compute_rheo_vshift.cpp -/compute_rheo_vshift.h -/fix_rheo.cpp -/fix_rheo.h -/fix_rheo_pressure.cpp -/fix_rheo_pressure.h -/fix_rheo_stress.cpp -/fix_rheo_stress.h -/fix_rheo_tension.cpp -/fix_rheo_tension.h -/fix_rheo_thermal.cpp -/fix_rheo_thermal.h -/fix_rheo_viscosity.cpp -/fix_rheo_viscosity.h -/pair_rheo.cpp -/pair_rheo.h -/pair_rheo_react.cpp -/pair_rheo_react.h +/*rheo* /compute_grid.cpp /compute_grid.h diff --git a/src/RHEO/atom_vec_rheo_thermal.cpp b/src/RHEO/atom_vec_rheo_thermal.cpp index de0c7fa5d7..4ecb7136a8 100644 --- a/src/RHEO/atom_vec_rheo_thermal.cpp +++ b/src/RHEO/atom_vec_rheo_thermal.cpp @@ -35,6 +35,7 @@ AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp) atom->status_flag = 1; atom->conductivity_flag = 1; atom->temperature_flag = 1; + atom->esph_flag = 1; atom->heatflow_flag = 1; atom->pressure_flag = 1; atom->rho_flag = 1; @@ -45,17 +46,17 @@ AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp) // order of fields in a string does not matter // except: fields_data_atom & fields_data_vel must match data file - fields_grow = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"}; - fields_copy = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"}; - fields_comm = {"status", "rho", "temperature"}; - fields_comm_vel = {"status", "rho", "temperature"}; + fields_grow = {"status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_copy = {"status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_comm = {"status", "rho", "esph"}; + fields_comm_vel = {"status", "rho", "esph"}; fields_reverse = {"drho", "heatflow"}; - fields_border = {"status", "rho", "temperature"}; - fields_border_vel = {"status", "rho", "temperature"}; - fields_exchange = {"status", "rho", "temperature"}; - fields_restart = {"status", "rho", "temperature"}; - fields_create = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"}; - fields_data_atom = {"id", "type", "status", "rho", "temperature", "x"}; + fields_border = {"status", "rho", "esph"}; + fields_border_vel = {"status", "rho", "esph"}; + fields_exchange = {"status", "rho", "esph"}; + fields_restart = {"status", "rho", "esph"}; + fields_create = {"status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_data_atom = {"id", "type", "status", "rho", "esph", "x"}; fields_data_vel = {"id", "v"}; setup_fields(); @@ -71,6 +72,7 @@ void AtomVecRHEOThermal::grow_pointers() status = atom->status; conductivity = atom->conductivity; temperature = atom->temperature; + esph = atom->esph; heatflow = atom->heatflow; pressure = atom->pressure; rho = atom->rho; @@ -98,6 +100,7 @@ void AtomVecRHEOThermal::data_atom_post(int ilocal) { drho[ilocal] = 0.0; heatflow[ilocal] = 0.0; + temperature[ilocal] = 0.0; pressure[ilocal] = 0.0; viscosity[ilocal] = 0.0; conductivity[ilocal] = 0.0; @@ -114,10 +117,11 @@ int AtomVecRHEOThermal::property_atom(const std::string &name) if (name == "rho") return 1; if (name == "drho") return 2; if (name == "temperature") return 3; - if (name == "heatflow") return 4; - if (name == "conductivity") return 5; - if (name == "pressure") return 6; - if (name == "viscosity") return 7; + if (name == "esph") return 4; + if (name == "heatflow") return 5; + if (name == "conductivity") return 6; + if (name == "pressure") return 7; + if (name == "viscosity") return 8; return -1; } @@ -167,7 +171,7 @@ void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues, } else if (index == 4) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) - buf[n] = heatflow[i]; + buf[n] = esph[i]; else buf[n] = 0.0; n += nvalues; @@ -175,7 +179,7 @@ void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues, } else if (index == 5) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) - buf[n] = conductivity[i]; + buf[n] = heatflow[i]; else buf[n] = 0.0; n += nvalues; @@ -183,12 +187,20 @@ void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues, } else if (index == 6) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) - buf[n] = pressure[i]; + buf[n] = conductivity[i]; else buf[n] = 0.0; n += nvalues; } } else if (index == 7) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = pressure[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 8) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) buf[n] = viscosity[i]; diff --git a/src/RHEO/atom_vec_rheo_thermal.h b/src/RHEO/atom_vec_rheo_thermal.h index 27c6c3c9b5..ad467f9de7 100644 --- a/src/RHEO/atom_vec_rheo_thermal.h +++ b/src/RHEO/atom_vec_rheo_thermal.h @@ -36,7 +36,7 @@ class AtomVecRHEOThermal : virtual public AtomVec { private: int *status; - double *conductivity, *temperature, *heatflow; + double *conductivity, *temperature, *heatflow, *esph; double *pressure, *rho, *drho, *viscosity; }; diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index b2ca0c9dc4..acfc01d793 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -42,16 +42,16 @@ enum{COMMGRAD, COMMFIELD}; ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), fix_rheo(nullptr), list(nullptr), rho0(nullptr), compute_interface(nullptr), compute_kernel(nullptr), - gradv(nullptr), gradr(nullptr), gradt(nullptr), gradn(nullptr) + gradv(nullptr), gradr(nullptr), grade(nullptr), gradn(nullptr) { if (narg < 4) error->all(FLERR,"Illegal compute rheo/grad command"); - velocity_flag = temperature_flag = rho_flag = eta_flag = 0; + velocity_flag = energy_flag = rho_flag = eta_flag = 0; for (int iarg = 3; iarg < narg; iarg++) { - if (strcmp(arg[iarg],"velocity") == 0) velocity_flag = 1; - else if (strcmp(arg[iarg],"rho") == 0) rho_flag = 1; - else if (strcmp(arg[iarg],"temperature") == 0) temperature_flag = 1; - else if (strcmp(arg[iarg],"viscosity") == 0) eta_flag = 1; + if (strcmp(arg[iarg], "velocity") == 0) velocity_flag = 1; + else if (strcmp(arg[iarg], "rho") == 0) rho_flag = 1; + else if (strcmp(arg[iarg], "energy") == 0) energy_flag = 1; + else if (strcmp(arg[iarg], "viscosity") == 0) eta_flag = 1; else error->all(FLERR, "Illegal compute rheo/grad command, {}", arg[iarg]); } @@ -72,7 +72,7 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) : comm_reverse += dim; } - if (temperature_flag) { + if (energy_flag) { ncomm_grad += dim; ncomm_field += 1; comm_reverse += dim; @@ -96,7 +96,7 @@ ComputeRHEOGrad::~ComputeRHEOGrad() { memory->destroy(gradv); memory->destroy(gradr); - memory->destroy(gradt); + memory->destroy(grade); memory->destroy(gradn); } @@ -130,7 +130,7 @@ void ComputeRHEOGrad::compute_peratom() int i, j, k, ii, jj, jnum, itype, jtype, a, b, fluidi, fluidj; double xtmp, ytmp, ztmp, delx, dely, delz; double rsq, imass, jmass; - double rhoi, rhoj, Voli, Volj, drho, dT, deta; + double rhoi, rhoj, Voli, Volj, drho, de, deta; double vi[3], vj[3], vij[3]; double wp, *dWij, *dWji; @@ -141,7 +141,7 @@ void ComputeRHEOGrad::compute_peratom() double **x = atom->x; double **v = atom->v; double *rho = atom->rho; - double *temperature = atom->temperature; + double *energy = atom->esph; double *viscosity = atom->viscosity; int *status = atom->status; int *type = atom->type; @@ -166,9 +166,9 @@ void ComputeRHEOGrad::compute_peratom() for (k = 0; k < dim; k++) gradr[i][k] = 0.0; } - if (temperature_flag) { + if (energy_flag) { for (k = 0; k < dim; k++) - gradt[i][k] = 0.0; + grade[i][k] = 0.0; } if (eta_flag) { for (k = 0; k < dim; k++) @@ -234,7 +234,7 @@ void ComputeRHEOGrad::compute_peratom() vij[2] = vi[2] - vj[2]; if (rho_flag) drho = rhoi - rhoj; - if (temperature_flag) dT = temperature[i] - temperature[j]; + if (energy_flag) de = energy[i] - energy[j]; if (eta_flag) deta = viscosity[i] - viscosity[j]; wp = compute_kernel->calc_dw(i, j, delx, dely, delz, sqrt(rsq)); @@ -250,8 +250,8 @@ void ComputeRHEOGrad::compute_peratom() if (rho_flag) // P,x P,y P,z gradr[i][a] -= drho * Volj * dWij[a]; - if (temperature_flag) // T,x T,y T,z - gradt[i][a] -= dT * Volj * dWij[a]; + if (energy_flag) // e,x e,y e,z + grade[i][a] -= de * Volj * dWij[a]; if (eta_flag) // n,x n,y n,z gradn[i][a] -= deta * Volj * dWij[a]; @@ -267,8 +267,8 @@ void ComputeRHEOGrad::compute_peratom() if (rho_flag) // P,x P,y P,z gradr[j][a] += drho * Voli * dWji[a]; - if (temperature_flag) // T,x T,y T,z - gradt[j][a] += dT * Voli * dWji[a]; + if (energy_flag) // e,x e,y e,z + grade[j][a] += de * Voli * dWji[a]; if (eta_flag) // n,x n,y n,z gradn[j][a] += deta * Voli * dWji[a]; @@ -308,7 +308,7 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf, int i,j,k,m; int *mask = atom->mask; double *rho = atom->rho; - double *temperature = atom->temperature; + double *energy = atom->esph; double **v = atom->v; int dim = domain->dimension; double *h_rate = domain->h_rate; @@ -335,9 +335,9 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf, for (k = 0; k < dim; k++) buf[m++] = gradr[j][k]; - if (temperature_flag) + if (energy_flag) for (k = 0; k < dim; k++) - buf[m++] = gradt[j][k]; + buf[m++] = grade[j][k]; if (eta_flag) for (k = 0; k < dim; k++) @@ -358,8 +358,8 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf, if (rho_flag) buf[m++] = rho[j]; - if (temperature_flag) - buf[m++] = temperature[j]; + if (energy_flag) + buf[m++] = energy[j]; } } return m; @@ -371,8 +371,7 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf) { int i, k, m, last; double *rho = atom->rho; - double *temperature = atom->temperature; - double **v = atom->v; + double *energy = atom->esph; double **v = atom->v; int dim = domain->dimension; m = 0; @@ -387,9 +386,9 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf) for (k = 0; k < dim; k++) gradr[i][k] = buf[m++]; - if (temperature_flag) + if (energy_flag) for (k = 0; k < dim; k++) - gradt[i][k] = buf[m++]; + grade[i][k] = buf[m++]; if (eta_flag) for (k = 0; k < dim; k++) @@ -403,8 +402,8 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf) if (rho_flag) rho[i] = buf[m++]; - if (temperature_flag) - temperature[i] = buf[m++]; + if (energy_flag) + energy[i] = buf[m++]; } } } @@ -427,9 +426,9 @@ int ComputeRHEOGrad::pack_reverse_comm(int n, int first, double *buf) for (k = 0; k < dim; k++) buf[m++] = gradr[i][k]; - if (temperature_flag) + if (energy_flag) for (k = 0; k < dim; k++) - buf[m++] = gradt[i][k]; + buf[m++] = grade[i][k]; if (eta_flag) for (k = 0; k < dim; k++) @@ -456,9 +455,9 @@ void ComputeRHEOGrad::unpack_reverse_comm(int n, int *list, double *buf) for (k = 0; k < dim; k++) gradr[j][k] += buf[m++]; - if (temperature_flag) + if (energy_flag) for (k = 0; k < dim; k++) - gradt[j][k] += buf[m++]; + grade[j][k] += buf[m++]; if (eta_flag) for (k = 0; k < dim; k++) @@ -477,8 +476,8 @@ void ComputeRHEOGrad::grow_arrays(int nmax) if (rho_flag) memory->grow(gradr, nmax, dim, "rheo:grad_rho"); - if (temperature_flag) - memory->grow(gradt, nmax, dim, "rheo:grad_temp"); + if (energy_flag) + memory->grow(grade, nmax, dim, "rheo:grad_energy"); if (eta_flag) memory->grow(gradn, nmax, dim, "rheo:grad_eta"); @@ -498,7 +497,7 @@ double ComputeRHEOGrad::memory_usage() if (rho_flag) bytes = (size_t) nmax_store * dim * sizeof(double); - if (temperature_flag) + if (energy_flag) bytes = (size_t) nmax_store * dim * sizeof(double); if (eta_flag) diff --git a/src/RHEO/compute_rheo_grad.h b/src/RHEO/compute_rheo_grad.h index 489f3c641d..2d663a5b07 100644 --- a/src/RHEO/compute_rheo_grad.h +++ b/src/RHEO/compute_rheo_grad.h @@ -40,7 +40,7 @@ class ComputeRHEOGrad : public Compute { void forward_fields(); double **gradv; double **gradr; - double **gradt; + double **grade; double **gradn; class FixRHEO *fix_rheo; @@ -48,7 +48,7 @@ class ComputeRHEOGrad : public Compute { int comm_stage, ncomm_grad, ncomm_field, nmax_store; double cut, cutsq, *rho0; - int velocity_flag, temperature_flag, rho_flag, eta_flag; + int velocity_flag, energy_flag, rho_flag, eta_flag; int interface_flag, remap_v_flag; class ComputeRHEOKernel *compute_kernel; diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index 52380a4337..6f58d79243 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -86,6 +86,8 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) : } comm_forward_save = comm_forward; + corrections_calculated = 0; + gsl_error_flag = 0; } /* ---------------------------------------------------------------------- */ @@ -152,6 +154,10 @@ int ComputeRHEOKernel::check_corrections(int i) if (coordination[i] < zmin) return 0; + // Skip if corrections not yet calculated + if (!corrections_calculated) + return 0; + return 1; } @@ -479,6 +485,7 @@ void ComputeRHEOKernel::compute_peratom() gsl_error_tags.clear(); if (kernel_style == QUINTIC) return; + corrections_calculated = 1; int i, j, ii, jj, inum, jnum, itype, g, a, b, gsl_error; double xtmp, ytmp, ztmp, r, rsq, w, vj, rhoj; @@ -530,11 +537,11 @@ void ComputeRHEOKernel::compute_peratom() if (rsq < hsq) { r = sqrt(rsq); - w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r); + w = calc_w_quintic(i, j, dx[0], dx[1], dx[2], r); rhoj = rho[j]; if (interface_flag) if (status[j] & PHASECHECK) - rhoj = compute_interface->correct_rho(j,i); + rhoj = compute_interface->correct_rho(j, i); vj = mass[type[j]] / rhoj; M += w * vj; @@ -578,12 +585,12 @@ void ComputeRHEOKernel::compute_peratom() if (rsq < hsq) { r = sqrt(rsq); - w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r); + w = calc_w_quintic(i, j, dx[0], dx[1], dx[2], r); rhoj = rho[j]; if (interface_flag) if (status[j] & PHASECHECK) - rhoj = compute_interface->correct_rho(j,i); + rhoj = compute_interface->correct_rho(j, i); vj = mass[type[j]] / rhoj; diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h index 2c9f4768e1..ed190c19ce 100644 --- a/src/RHEO/compute_rheo_kernel.h +++ b/src/RHEO/compute_rheo_kernel.h @@ -53,6 +53,7 @@ class ComputeRHEOKernel : public Compute { int gsl_error_flag; std::unordered_set gsl_error_tags; + int corrections_calculated; int kernel_style, zmin, dim, Mdim, ncor; int nmax_store; double h, hsq, hinv, hsqinv, pre_w, pre_wp; diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index e450eaaf0b..380ff398d8 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -104,6 +104,14 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a } else if (utils::strmatch(arg[iarg], "^grad/v")) { pack_choice[i] = &ComputeRHEOPropertyAtom::pack_gradv; col_index[i] = get_tensor_index(arg[iarg]); + } else if (strcmp(arg[iarg], "energy") == 0) { + avec_index[i] = atom->avec->property_atom("esph"); + if (avec_index[i] < 0) + error->all(FLERR, + "Invalid keyword {} for atom style {} in compute rheo/property/atom command ", + arg[iarg], atom->get_style()); + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style; + thermal_flag = 1; } else { avec_index[i] = atom->avec->property_atom(arg[iarg]); if (avec_index[i] < 0) diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 1f9314e99b..569c8569f7 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -32,6 +32,8 @@ #include "neigh_list.h" #include "neigh_request.h" +#include "update.h" + using namespace LAMMPS_NS; using namespace RHEO_NS; @@ -195,6 +197,7 @@ void ComputeRHEOVShift::compute_peratom() if (mask[i] & groupbit) { vmag = sqrt(vi[0] * vi[0] + vi[1] * vi[1] + vi[2] * vi[2]); prefactor = vmag * volj * dr; + vshift[i][0] += prefactor * dx[0]; vshift[i][1] += prefactor * dx[1]; vshift[i][2] += prefactor * dx[2]; @@ -204,6 +207,7 @@ void ComputeRHEOVShift::compute_peratom() if (mask[j] & groupbit) { vmag = sqrt(vj[0] * vj[0] + vj[1] * vj[1] + vj[2] * vj[2]); prefactor = vmag * voli * dr; + vshift[j][0] -= prefactor * dx[0]; vshift[j][1] -= prefactor * dx[1]; vshift[j][2] -= prefactor * dx[2]; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 0c74b2bad1..beba940174 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -164,7 +164,7 @@ void FixRHEO::post_constructor() compute_kernel->fix_rheo = this; std::string cmd = "rheo_grad all RHEO/GRAD velocity rho viscosity"; - if (thermal_flag) cmd += " temperature"; + if (thermal_flag) cmd += " energy"; compute_grad = dynamic_cast(modify->add_compute(cmd)); compute_grad->fix_rheo = this; @@ -401,6 +401,7 @@ void FixRHEO::pre_force(int /*vflag*/) compute_rhosum->compute_peratom(); compute_kernel->compute_peratom(); + if (interface_flag) { // Note on first setup, have no forces for pressure to reference compute_interface->compute_peratom(); diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index e6c598418b..7b61b9821e 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -47,7 +47,8 @@ enum {NONE, CONSTANT}; FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), fix_rheo(nullptr), compute_grad(nullptr), compute_vshift(nullptr), - Tc(nullptr), kappa(nullptr), cv(nullptr), Tc_style(nullptr), kappa_style(nullptr), cv_style(nullptr), + Tc(nullptr), kappa(nullptr), cv(nullptr), L(nullptr), + Tc_style(nullptr), kappa_style(nullptr), cv_style(nullptr), L_style(nullptr), fix_update_special_bonds(nullptr) { if (narg < 4) error->all(FLERR,"Illegal fix command"); @@ -62,13 +63,22 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : int i, nlo, nhi; int n = atom->ntypes; + memory->create(Tc_style, n + 1, "rheo:Tc_style"); memory->create(kappa_style, n + 1, "rheo:kappa_style"); memory->create(cv_style, n + 1, "rheo:cv_style"); + memory->create(L_style, n + 1, "rheo:L_style"); + + memory->create(Tc, n + 1, "rheo:Tc"); + memory->create(kappa, n + 1, "rheo:kappa"); + memory->create(cv, n + 1, "rheo:cv"); + memory->create(L, n + 1, "rheo:L"); + for (i = 1; i <= n; i++) { Tc_style[i] = NONE; kappa_style[i] = NONE; cv_style[i] = NONE; + L_style[i] = NONE; } int iarg = 3; @@ -81,9 +91,9 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg + 2], "constant") == 0) { if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal conductivity constant", error); - double kappa_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + double kappa_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); if (kappa_one < 0.0) error->all(FLERR, "The conductivity must be positive"); - iarg += 1; + iarg += 2; for (i = nlo; i <= nhi; i++) { kappa_style[i] = CONSTANT; @@ -102,9 +112,9 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg + 2], "constant") == 0) { if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal specific/heat constant", error); - double cv_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + double cv_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); if (cv_one < 0.0) error->all(FLERR, "The specific heat must be positive"); - iarg += 1; + iarg += 2; for (i = nlo; i <= nhi; i++) { cv_style[i] = CONSTANT; @@ -124,8 +134,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg + 2], "constant") == 0) { if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal Tfreeze constant", error); - double Tc_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 1; + double Tc_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 2; for (i = nlo; i <= nhi; i++) { Tc_style[i] = CONSTANT; @@ -136,6 +146,28 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Illegal fix command, {}", arg[iarg + 2]); } + iarg += 2; + } else if (strcmp(arg[iarg], "latent/heat") == 0) { + if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal latent/heat", error); + utils::bounds(FLERR, arg[iarg + 1], 1, n, nlo, nhi, error); + + // Cv arguments + if (strcmp(arg[iarg + 2], "constant") == 0) { + if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal latent/heat constant", error); + + double L_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + if (L_one < 0.0) error->all(FLERR, "The latent heat must be positive"); + iarg += 2; + + for (i = nlo; i <= nhi; i++) { + L_style[i] = CONSTANT; + L[i] = L_one; + } + + } else { + error->all(FLERR,"Illegal fix command, {}", arg[iarg + 2]); + } + iarg += 2; } else if (strcmp(arg[iarg], "react") == 0) { if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal react", error); @@ -155,9 +187,11 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : for (i = 1; i <= n; i++) { if (cv_style[i] == NONE) - error->all(FLERR,"Must specify specific/heat for atom type {} in fix/rheo/thermal", i); + error->all(FLERR, "Must specify specific/heat for atom type {} in fix/rheo/thermal", i); if (kappa_style[i] == NONE) - error->all(FLERR,"Must specify conductivity for atom type {} in fix/rheo/thermal", i); + error->all(FLERR, "Must specify conductivity for atom type {} in fix/rheo/thermal", i); + if (Tc_style[i] == NONE && L_style[i] != NONE) + error->all(FLERR, "Must specify critical temperature for atom type {} to use latent heat in fix rheo/thermal", i); } } @@ -173,9 +207,11 @@ FixRHEOThermal::~FixRHEOThermal() memory->destroy(cv_style); memory->destroy(Tc_style); memory->destroy(kappa_style); + memory->destroy(L_style); memory->destroy(cv); memory->destroy(Tc); memory->destroy(kappa); + memory->destroy(L); } /* ---------------------------------------------------------------------- */ @@ -210,6 +246,8 @@ void FixRHEOThermal::init() dtf = 0.5 * update->dt * force->ftm2v; + if (atom->esph_flag != 1) + error->all(FLERR,"fix rheo/thermal command requires atom property esph"); if (atom->temperature_flag != 1) error->all(FLERR,"fix rheo/thermal command requires atom property temperature"); if (atom->heatflow_flag != 1) @@ -267,9 +305,9 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/) int i, a; int *status = atom->status; - double *temperature = atom->temperature; - double **gradt = compute_grad->gradt; - double **vshift = compute_vshift->array_atom; + double *energy = atom->esph; + double **grade = compute_grad->grade; + double **vshift = compute_vshift->vshift; int nlocal = atom->nlocal; int dim = domain->dimension; @@ -279,9 +317,8 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/) for (i = 0; i < nlocal; i++) { if (status[i] & STATUS_NO_SHIFT) continue; - for (a = 0; a < dim; a++) - temperature[i] += dtv * vshift[i][a] * gradt[i][a]; + energy[i] += dtv * vshift[i][a] * grade[i][a]; } } @@ -290,29 +327,34 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/) void FixRHEOThermal::post_integrate() { int i, itype; - double cvi, Tci, Ti; + double cvi, Tci, Ti, Li; int *status = atom->status; + double *energy = atom->esph; double *temperature = atom->temperature; double *heatflow = atom->heatflow; - double *rho = atom->rho; int *type = atom->type; int n_melt = 0; int n_freeze = 0; - //Integrate temperature and check status + //Integrate energy and check status for (i = 0; i < atom->nlocal; i++) { if (status[i] & STATUS_NO_INTEGRATION) continue; itype = type[i]; - cvi = calc_cv(i, type[i]); - temperature[i] += dtf * heatflow[i] / cvi; + cvi = calc_cv(i, itype); + energy[i] += dtf * heatflow[i]; + temperature[i] = energy[i] / cvi; if (Tc_style[itype] != NONE) { Ti = temperature[i]; - if (Tc_style[itype] == CONSTANT) { - Tci = Tc[itype]; + Tci = calc_Tc(i, itype); + + if (L_style[itype] != NONE) { + Li = calc_L(i, itype); + if (Ti > Tci) Ti = MAX(Tci, (energy[i] - Li) / cvi); + temperature[i] = Ti; } if (Ti > Tci) { @@ -370,11 +412,39 @@ void FixRHEOThermal::post_neighbor() } /* ---------------------------------------------------------------------- + Calculate temperature In the future, update & forward evolving conductivity styles every timestep ------------------------------------------------------------------------- */ void FixRHEOThermal::pre_force(int /*vflag*/) { + int i, itype; + double cvi, Tci, Ti, Li; + + double *energy = atom->esph; + double *temperature = atom->temperature; + int *type = atom->type; + + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + // Calculate temperature + for (i = 0; i < nall; i++) { + itype = type[i]; + cvi = calc_cv(i, itype); + temperature[i] = energy[i] / cvi; + + if (Tc_style[itype] != NONE) { + Ti = temperature[i]; + Tci = calc_Tc(i, itype); + + if (L_style[itype] != NONE) { + Li = calc_L(i, itype); + if (Ti > Tci) Ti = MAX(Tci, (energy[i] - Li) / cvi); + temperature[i] = Ti; + } + } + } } /* ---------------------------------------------------------------------- */ @@ -382,17 +452,13 @@ void FixRHEOThermal::pre_force(int /*vflag*/) void FixRHEOThermal::final_integrate() { int *status = atom->status; - int *type = atom->type; - double *temperature = atom->temperature; + double *energy = atom->esph; double *heatflow = atom->heatflow; - double cvi; - //Integrate temperature and check status + //Integrate energy for (int i = 0; i < atom->nlocal; i++) { if (status[i] & STATUS_NO_INTEGRATION) continue; - - cvi = calc_cv(i, type[i]); - temperature[i] += dtf * heatflow[i] / cvi; + energy[i] += dtf * heatflow[i]; } } @@ -535,6 +601,23 @@ double FixRHEOThermal::calc_cv(int i, int itype) } } +/* ---------------------------------------------------------------------- */ + +double FixRHEOThermal::calc_Tc(int i, int itype) +{ + if (Tc_style[itype] == CONSTANT) { + return Tc[itype]; + } +} + +/* ---------------------------------------------------------------------- */ + +double FixRHEOThermal::calc_L(int i, int itype) +{ + if (L_style[itype] == CONSTANT) { + return L[itype]; + } +} /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/fix_rheo_thermal.h b/src/RHEO/fix_rheo_thermal.h index dc412d20b9..c4c26eef6a 100644 --- a/src/RHEO/fix_rheo_thermal.h +++ b/src/RHEO/fix_rheo_thermal.h @@ -43,12 +43,14 @@ class FixRHEOThermal : public Fix { void unpack_forward_comm(int, int, double *) override; void reset_dt() override; double calc_cv(int, int); + double calc_Tc(int, int); + double calc_L(int, int); private: - double *cv, *Tc, *kappa; + double *cv, *Tc, *kappa, *L; double dtf, dtv; double cut_kernel, cut_bond, cutsq_bond; - int *cv_style, *Tc_style, *kappa_style; + int *cv_style, *Tc_style, *kappa_style, *L_style; int btype; class NeighList *list; diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index 4ef0066e93..b07e914af1 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -96,7 +96,6 @@ void PairRHEO::compute(int eflag, int vflag) ev_init(eflag, vflag); double **gradv = compute_grad->gradv; - double **gradt = compute_grad->gradt; double **gradr = compute_grad->gradr; double **v = atom->v; double **x = atom->x; diff --git a/src/RHEO/pair_rheo_solid.cpp b/src/RHEO/pair_rheo_solid.cpp new file mode 100644 index 0000000000..1068e9b329 --- /dev/null +++ b/src/RHEO/pair_rheo_solid.cpp @@ -0,0 +1,351 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "pair_rheo_solid.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "fix_rheo.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" + +#include + +using namespace LAMMPS_NS; +using namespace RHEO_NS; + +/* ---------------------------------------------------------------------- */ + +PairRHEOSolid::PairRHEOSolid(LAMMPS *_lmp) : Pair(_lmp) +{ + writedata = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairRHEOSolid::~PairRHEOSolid() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(k); + memory->destroy(cut); + memory->destroy(gamma); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairRHEOSolid::compute(int eflag, int vflag) +{ + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair; + double r, rsq, rinv, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; + double vxtmp, vytmp, vztmp, delvx, delvy, delvz, dot, smooth; + + evdwl = 0.0; + if (eflag || vflag) + ev_setup(eflag, vflag); + else + evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + int *type = atom->type; + int *status = atom->status; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + double *special_lj = force->special_lj; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(status[i] & STATUS_SOLID)) continue; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + vxtmp = v[i][0]; + vytmp = v[i][1]; + vztmp = v[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + if (!(status[j] & STATUS_SOLID)) continue; + + factor_lj = special_lj[sbmask(j)]; + + if (factor_lj == 0) continue; + + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + rinv = 1.0 / r; + fpair = k[itype][jtype] * (cut[itype][jtype] - r); + + smooth = rsq / cutsq[itype][jtype]; + smooth *= smooth; + smooth *= smooth; + smooth = 1.0 - smooth; + delvx = vxtmp - v[j][0]; + delvy = vytmp - v[j][1]; + delvz = vztmp - v[j][2]; + dot = delx * delvx + dely * delvy + delz * delvz; + fpair -= gamma[itype][jtype] * dot * smooth * rinv; + + fpair *= factor_lj * rinv; + if (eflag) evdwl = 0.0; + + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; + + if (newton_pair || j < nlocal) { + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; + } + + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairRHEOSolid::allocate() +{ + allocated = 1; + const int np1 = atom->ntypes + 1; + + memory->create(setflag, np1, np1, "pair:setflag"); + for (int i = 1; i < np1; i++) + for (int j = i; j < np1; j++) setflag[i][j] = 0; + + memory->create(cutsq, np1, np1, "pair:cutsq"); + + memory->create(k, np1, np1, "pair:k"); + memory->create(cut, np1, np1, "pair:cut"); + memory->create(gamma, np1, np1, "pair:gamma"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairRHEOSolid::settings(int narg, char ** /*arg*/) +{ + if (narg != 0) error->all(FLERR, "Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairRHEOSolid::coeff(int narg, char **arg) +{ + if (narg != 5) error->all(FLERR, "Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + + double k_one = utils::numeric(FLERR, arg[2], false, lmp); + double cut_one = utils::numeric(FLERR, arg[3], false, lmp); + double gamma_one = utils::numeric(FLERR, arg[4], false, lmp); + + if (cut_one <= 0.0) error->all(FLERR, "Incorrect args for pair coefficients"); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { + k[i][j] = k_one; + cut[i][j] = cut_one; + gamma[i][j] = gamma_one; + + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairRHEOSolid::init_style() +{ + if (comm->ghost_velocity == 0) + error->all(FLERR,"Pair rheo/solid requires ghost atoms store velocity"); + + if (!atom->status_flag) + error->all(FLERR,"Pair rheo/solid requires atom_style rheo"); + + neighbor->add_request(this); +} + + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairRHEOSolid::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); + k[i][j] = mix_energy(k[i][i], k[j][j], cut[i][i], cut[j][j]); + gamma[i][j] = mix_energy(gamma[i][i], gamma[j][j], cut[i][i], cut[j][j]); + } + + cut[j][i] = cut[i][j]; + k[j][i] = k[i][j]; + gamma[j][i] = gamma[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairRHEOSolid::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i, j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j], sizeof(int), 1, fp); + if (setflag[i][j]) { + fwrite(&k[i][j], sizeof(double), 1, fp); + fwrite(&cut[i][j], sizeof(double), 1, fp); + fwrite(&gamma[i][j], sizeof(double), 1, fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairRHEOSolid::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i, j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); + if (setflag[i][j]) { + if (me == 0) { + utils::sfread(FLERR, &k[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &gamma[i][j], sizeof(double), 1, fp, nullptr, error); + } + MPI_Bcast(&k[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&gamma[i][j], 1, MPI_DOUBLE, 0, world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairRHEOSolid::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp, "%d %g %g %g\n", i, k[i][i], cut[i][i], gamma[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairRHEOSolid::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp, "%d %d %g %g %g\n", i, j, k[i][j], cut[i][j], gamma[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairRHEOSolid::single(int i, int j, int itype, int jtype, double rsq, double /*factor_coul*/, + double factor_lj, double &fforce) +{ + double fpair, r, rinv; + double delx, dely, delz, delvx, delvy, delvz, dot, smooth; + + if (rsq > cutsq[itype][jtype]) return 0.0; + + double **x = atom->x; + double **v = atom->v; + + r = sqrt(rsq); + rinv = 1.0 / r; + + fpair = k[itype][jtype] * (cut[itype][jtype] - r); + + smooth = rsq / cutsq[itype][jtype]; + smooth *= smooth; + smooth = 1.0 - smooth; + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + delvx = v[i][0] - v[j][0]; + delvy = v[i][1] - v[j][1]; + delvz = v[i][2] - v[j][2]; + dot = delx * delvx + dely * delvy + delz * delvz; + fpair -= gamma[itype][jtype] * dot * rinv * smooth; + + fpair *= factor_lj; + fforce = fpair; + + return 0.0; +} diff --git a/src/RHEO/pair_rheo_solid.h b/src/RHEO/pair_rheo_solid.h new file mode 100644 index 0000000000..66c2ac4bf1 --- /dev/null +++ b/src/RHEO/pair_rheo_solid.h @@ -0,0 +1,51 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(rheo/solid,PairRHEOSolid); +// clang-format on +#else + +#ifndef LMP_PAIR_RHEO_SOLID_H +#define LMP_PAIR_RHEO_SOLID_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairRHEOSolid : public Pair { + public: + PairRHEOSolid(class LAMMPS *); + ~PairRHEOSolid() override; + void compute(int, int) override; + void settings(int, char **) override; + void coeff(int, char **) override; + void init_style() override; + double init_one(int, int) override; + void write_restart(FILE *) override; + void read_restart(FILE *) override; + void write_data(FILE *) override; + void write_data_all(FILE *) override; + double single(int, int, int, int, double, double, double, double &) override; + + protected: + double **k, **cut, **gamma; + + void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/set.cpp b/src/set.cpp index 3e1058b048..e0d20f1dc7 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -901,7 +901,7 @@ void Set::set(int keyword) } else if (keyword == RHEO_STATUS) { - if (ivalue != 0 && ivalue !=2) + if (ivalue != 0 && ivalue != 2) error->one(FLERR,"Invalid value {} in set command for rheo/status", ivalue); atom->status[i] = ivalue; }