diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index 50bcb2a2d0..e450eaaf0b 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -14,7 +14,7 @@ /* ---------------------------------------------------------------------- Contributing authors: - Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) + Joel Clemmer (SNL) ----------------------------------------------------------------------- */ #include "compute_rheo_property_atom.h" @@ -29,6 +29,7 @@ #include "domain.h" #include "error.h" #include "fix_rheo.h" +#include "fix_rheo_pressure.h" #include "fix_rheo_thermal.h" #include "memory.h" #include "modify.h" @@ -44,7 +45,7 @@ using namespace RHEO_NS; /* ---------------------------------------------------------------------- */ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), fix_rheo(nullptr), fix_thermal(nullptr), compute_interface(nullptr), + Compute(lmp, narg, arg), fix_rheo(nullptr), fix_pressure(nullptr), fix_thermal(nullptr), compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr), compute_vshift(nullptr), compute_grad(nullptr), avec_index(nullptr), pack_choice(nullptr), col_index(nullptr) { @@ -55,7 +56,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a if (nvalues == 1) size_peratom_cols = 0; else size_peratom_cols = nvalues; - thermal_flag = interface_flag = surface_flag = shift_flag = 0; + pressure_flag = thermal_flag = interface_flag = surface_flag = shift_flag = 0; // parse input values // customize a new keyword by adding to if statement @@ -90,6 +91,9 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a col_index[i] = get_vector_index(arg[iarg]); } else if (strcmp(arg[iarg],"coordination") == 0) { pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination; + } else if (strcmp(arg[iarg],"pressure") == 0) { + pressure_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_pressure; } else if (strcmp(arg[iarg],"cv") == 0) { thermal_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv; @@ -155,6 +159,11 @@ void ComputeRHEOPropertyAtom::init() fixes = modify->get_fix_by_style("rheo/thermal"); fix_thermal = dynamic_cast(fixes[0]); } + + if (pressure_flag) { + fixes = modify->get_fix_by_style("rheo/pressure"); + fix_pressure = dynamic_cast(fixes[0]); + } } /* ---------------------------------------------------------------------- */ @@ -343,6 +352,22 @@ void ComputeRHEOPropertyAtom::pack_cv(int n) /* ---------------------------------------------------------------------- */ +void ComputeRHEOPropertyAtom::pack_pressure(int n) +{ + int *type = atom->type; + int *mask = atom->mask; + double *rho = atom->rho; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = fix_pressure->calc_pressure(rho[i], type[i]); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + void ComputeRHEOPropertyAtom::pack_shift_v(int n) { double **vshift = compute_vshift->vshift; diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h index bfae870ee5..f3596fbbf9 100644 --- a/src/RHEO/compute_rheo_property_atom.h +++ b/src/RHEO/compute_rheo_property_atom.h @@ -34,7 +34,7 @@ class ComputeRHEOPropertyAtom : public Compute { private: int nvalues, nmax; - int thermal_flag, interface_flag, surface_flag, shift_flag; + int pressure_flag, thermal_flag, interface_flag, surface_flag, shift_flag; int *avec_index; int *col_index; double *buf; @@ -53,12 +53,14 @@ class ComputeRHEOPropertyAtom : public Compute { void pack_cv(int); void pack_shift_v(int); void pack_gradv(int); + void pack_pressure(int); void pack_atom_style(int); int get_vector_index(char*); int get_tensor_index(char*); class FixRHEO *fix_rheo; + class FixRHEOPressure *fix_pressure; class FixRHEOThermal *fix_thermal; class ComputeRHEOInterface *compute_interface; class ComputeRHEOKernel *compute_kernel; diff --git a/src/RHEO/compute_rheo_rho_sum.cpp b/src/RHEO/compute_rheo_rho_sum.cpp index 0a2096a2b9..82d3aa4bc6 100644 --- a/src/RHEO/compute_rheo_rho_sum.cpp +++ b/src/RHEO/compute_rheo_rho_sum.cpp @@ -88,7 +88,7 @@ void ComputeRHEORhoSum::compute_peratom() // initialize arrays, local with quintic self-contribution, ghosts are zeroed for (i = 0; i < nlocal; i++) { w = compute_kernel->calc_w_quintic(i, i, 0.0, 0.0, 0.0, 0.0); - rho[i] += w * mass[type[i]]; + rho[i] = w * mass[type[i]]; } for (i = nlocal; i < nall; i++) rho[i] = 0.0; @@ -131,12 +131,11 @@ int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf, { int i, j, k, m; double *rho = atom->rho; - int *coordination = compute_kernel->coordination; m = 0; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = coordination[j]; + buf[m++] = rho[j]; } return m; } diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index d07f0d1a1f..f0f380f23a 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -279,6 +279,9 @@ void FixRHEO::setup(int /*vflag*/) error->one(FLERR, "Fix rheo/viscosity does not fully cover all atoms"); if (!t_coverage_flag) error->one(FLERR, "Fix rheo/thermal does not fully cover all atoms"); + + if (rhosum_flag) + compute_rhosum->compute_peratom(); } /* ---------------------------------------------------------------------- */ @@ -419,10 +422,11 @@ void FixRHEO::pre_force(int /*vflag*/) status[i] &= OPTIONSMASK; // Calculate surfaces, update status - if (surface_flag) compute_surface->compute_peratom(); - - if (shift_flag) - compute_vshift->correct_surfaces(); + if (surface_flag) { + compute_surface->compute_peratom(); + if (shift_flag) + compute_vshift->correct_surfaces(); + } } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index eac4b34046..8c523b2b35 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -31,14 +31,14 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum {NONE, LINEAR, CUBIC, TAITWATER}; +enum {NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL}; static constexpr double SEVENTH = 1.0 / 7.0; /* ---------------------------------------------------------------------- */ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), fix_rheo(nullptr), rho0(nullptr), csq(nullptr), rho0inv(nullptr), csqinv(nullptr), c_cubic(nullptr), pressure_style(nullptr) + Fix(lmp, narg, arg), fix_rheo(nullptr), rho0(nullptr), csq(nullptr), rho0inv(nullptr), csqinv(nullptr), c_cubic(nullptr), tpower(nullptr), pbackground(nullptr), pressure_style(nullptr) { if (narg < 4) error->all(FLERR,"Illegal fix command"); @@ -51,6 +51,9 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : int i, nlo, nhi; int n = atom->ntypes; memory->create(pressure_style, n + 1, "rheo:pressure_style"); + memory->create(c_cubic, n + 1, "rheo:c_cubic"); + memory->create(tpower, n + 1, "rheo:tpower"); + memory->create(pbackground, n + 1, "rheo:pbackground"); for (i = 1; i <= n; i++) pressure_style[i] = NONE; int iarg = 3; @@ -62,9 +65,21 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg + 1], "linear") == 0) { for (i = nlo; i <= nhi; i++) pressure_style[i] = LINEAR; - } else if (strcmp(arg[iarg + 1], "taitwater") == 0) { + } else if (strcmp(arg[iarg + 1], "tait/water") == 0) { for (i = nlo; i <= nhi; i++) pressure_style[i] = TAITWATER; + } else if (strcmp(arg[iarg + 1], "tait/general") == 0) { + if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure tait", error); + + double tpower_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + double pbackground_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 2; + + for (i = nlo; i <= nhi; i++) { + pressure_style[i] = TAITGENERAL; + tpower[i] = tpower_one; + pbackground[i] = pbackground_one; + } } else if (strcmp(arg[iarg + 1], "cubic") == 0) { if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure cubic", error); @@ -94,6 +109,8 @@ FixRHEOPressure::~FixRHEOPressure() memory->destroy(csqinv); memory->destroy(rho0inv); memory->destroy(c_cubic); + memory->destroy(tpower); + memory->destroy(pbackground); } /* ---------------------------------------------------------------------- */ @@ -204,6 +221,10 @@ double FixRHEOPressure::calc_pressure(double rho, int type) rho_ratio = rho * rho0inv[type]; rr3 = rho_ratio * rho_ratio * rho_ratio; p = csq[type] * rho0[type] * SEVENTH * (rr3 * rr3 * rho_ratio - 1.0); + } else if (pressure_style[type] == TAITGENERAL) { + rho_ratio = rho * rho0inv[type]; + p = csq[type] * rho0[type] * (pow(rho_ratio, tpower[type]) - 1.0) / tpower[type]; + p += pbackground[type]; } return p; } @@ -222,6 +243,11 @@ double FixRHEOPressure::calc_rho(double p, int type) rho = pow(7.0 * p + csq[type] * rho0[type], SEVENTH); rho *= pow(rho0[type], 6.0 * SEVENTH); rho *= pow(csq[type], -SEVENTH); + } else if (pressure_style[type] == TAITGENERAL) { + p -= pbackground[type]; + rho = pow(tpower[type] * p + csq[type] * rho0[type], 1.0 / tpower[type]); + rho *= pow(rho0[type], 1.0 - 1.0 / tpower[type]); + rho *= pow(csq[type], -1.0 / tpower[type]); } return rho; } diff --git a/src/RHEO/fix_rheo_pressure.h b/src/RHEO/fix_rheo_pressure.h index ee86c5e184..ca165b1ed5 100644 --- a/src/RHEO/fix_rheo_pressure.h +++ b/src/RHEO/fix_rheo_pressure.h @@ -38,7 +38,7 @@ class FixRHEOPressure : public Fix { double calc_rho(double, int); private: - double *c_cubic, *csq, *csqinv, *rho0, *rho0inv; + double *c_cubic, *csq, *csqinv, *rho0, *rho0inv, *tpower, *pbackground; int *pressure_style; class FixRHEO *fix_rheo; diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index dd2c6eddbd..e6c598418b 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -13,7 +13,7 @@ /* ---------------------------------------------------------------------- Contributing authors: - Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) + Joel Clemmer (SNL) ----------------------------------------------------------------------- */ #include "fix_rheo_thermal.h" diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp index 91799ccfd0..2fffa8b29c 100644 --- a/src/RHEO/fix_rheo_viscosity.cpp +++ b/src/RHEO/fix_rheo_viscosity.cpp @@ -52,6 +52,10 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) : int i, nlo, nhi; int n = atom->ntypes; memory->create(viscosity_style, n + 1, "rheo:viscosity_style"); + memory->create(eta, n + 1, "rheo:eta"); + memory->create(gd0, n + 1, "rheo:gd0"); + memory->create(K, n + 1, "rheo:K"); + memory->create(npow, n + 1, "rheo:npow"); for (i = 1; i <= n; i++) viscosity_style[i] = NONE; int iarg = 3; diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index bb3d5c3fda..339efed866 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -56,6 +56,7 @@ PairRHEO::PairRHEO(LAMMPS *lmp) : artificial_visc_flag = 0; rho_damp_flag = 0; thermal_flag = 0; + harmonic_means_flag = 0; comm_reverse = 3; } @@ -80,7 +81,7 @@ void PairRHEO::compute(int eflag, int vflag) int pair_force_flag, pair_rho_flag, pair_avisc_flag; int fluidi, fluidj; double xtmp, ytmp, ztmp, w, wp, Ti, Tj, dT, csq_ave, cs_ave; - double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj; + double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, eta_ave, kappa_ave,dT_prefactor; double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij; double *dWij, *dWji, *dW1ij, *dW1ji; double dx[3], du[3], dv[3], fv[3], dfp[3], fsolid[3], ft[3], vi[3], vj[3]; @@ -112,6 +113,7 @@ void PairRHEO::compute(int eflag, int vflag) int *type = atom->type; int *status = atom->status; tagint *tag = atom->tag; + double fnorm, ftang[3]; double **fp_store, *chi; if (compute_interface) { @@ -243,14 +245,19 @@ void PairRHEO::compute(int eflag, int vflag) // Thermal Evolution if (thermal_flag) { + if (harmonic_means_flag) { + kappa_ave = 2.0 * kappai * kappaj / (kappai + kappaj); + } else { + kappa_ave = 0.5 * (kappai * kappaj); + } + dT_prefactor = 2.0 * kappa_ave * (Ti - Tj) * rinv * rinv * voli * volj * 2.0 / (rhoi + rhoj); + dT = dot3(dx, dWij); - dT *= (kappai + kappaj) * (Ti - Tj) * rinv * rinv * voli * volj * 2.0 / (rhoi + rhoj); - heatflow[i] += dT; + heatflow[i] += dT * dT_prefactor; if (newton_pair || j < nlocal) { dT = dot3(dx, dWji); - dT *= (kappai + kappaj) * (Tj - Ti) * rinv * rinv * voli * volj * 2.0 / (rhoi + rhoj); - heatflow[j] -= dT; + heatflow[j] += dT * dT_prefactor; } } @@ -260,6 +267,12 @@ void PairRHEO::compute(int eflag, int vflag) fp_prefactor = voli * volj * (Pj + Pi); sub3(vi, vj, dv); + if (harmonic_means_flag) { + eta_ave = 2.0 * etai * etaj / (etai + etaj); + } else { + eta_ave = 0.5 * (etai * etaj); + } + //Add artificial viscous pressure if required if (artificial_visc_flag && pair_avisc_flag) { //Interpolate velocities to midpoint and use this difference for artificial viscosity @@ -283,7 +296,7 @@ void PairRHEO::compute(int eflag, int vflag) fv[a] = 0.0; for (b = 0; b < dim; b++) fv[a] += dv[a] * dx[b] * dWij[b]; - fv[a] *= (etai + etaj) * voli * volj * rinv * rinv; + fv[a] *= 2.0 * eta_ave * voli * volj * rinv * rinv; } add3(fv, dfp, ft); @@ -293,26 +306,38 @@ void PairRHEO::compute(int eflag, int vflag) f[i][1] += ft[1]; f[i][2] += ft[2]; - if (evflag) // Does not account for unbalanced forces - ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], dx[0], dx[1], dx[2]); + if (evflag) { + fnorm = dot3(ft, dx) * rinv * rinv * 0.5; + ftang[0] = ft[0] * 0.5 - dx[0] * fnorm; + ftang[1] = ft[1] * 0.5 - dx[1] * fnorm; + ftang[2] = ft[2] * 0.5 - dx[2] * fnorm; + ev_tally_nt(i, j, nlocal, newton_pair, 0.0, 0.0, fnorm, ftang[0], ftang[1], ftang[2], dx[0], dx[1], dx[2]); + } if (newton_pair || j < nlocal) { for (a = 0; a < dim; a ++) { fv[a] = 0.0; for (b = 0; b < dim; b++) fv[a] += (vi[a] - vj[a]) * dx[b] * dWji[b]; - fv[a] *= -(etai + etaj) * voli * volj * rinv * rinv; + fv[a] *= -2.0 * eta_ave * voli * volj * rinv * rinv; // flip sign here b/c -= at accummulator } scale3(fp_prefactor, dWji, dfp); - add3(fv, dfp, ft); add3(fsolid, ft, ft); f[j][0] -= ft[0]; f[j][1] -= ft[1]; f[j][2] -= ft[2]; + + if (evflag) { + fnorm = - dot3(ft, dx) * rinv * rinv * 0.5; + ftang[0] = ft[0] * 0.5 + dx[0] * fnorm; + ftang[1] = ft[1] * 0.5 + dx[1] * fnorm; + ftang[2] = ft[2] * 0.5 + dx[2] * fnorm; + ev_tally_nt(i, j, nlocal, newton_pair, 0.0, 0.0, fnorm, ftang[0], ftang[1], ftang[2], -dx[0], -dx[1], -dx[2]); + } } if (compute_interface) { @@ -360,7 +385,7 @@ void PairRHEO::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); if (compute_interface) { - comm->reverse_comm(this); + if (newton_pair) comm->reverse_comm(this); comm->forward_comm(this); } } @@ -404,6 +429,8 @@ void PairRHEO::settings(int narg, char **arg) artificial_visc_flag = 1; av = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg++; + } else if (strcmp(arg[iarg], "harmonic/means") == 0) { + harmonic_means_flag = 1; } else error->all(FLERR, "Illegal pair_style command, {}", arg[iarg]); iarg++; } @@ -469,7 +496,8 @@ void PairRHEO::setup() int n = atom->ntypes; memory->create(cs, n + 1, "rheo:cs"); - for (int i = 0; i <= n; i++) cs[i] = sqrt(csq[i]); + for (int i = 1; i <= n; i++) + cs[i] = sqrt(csq[i]); if (comm->ghost_velocity == 0) error->all(FLERR, "Pair RHEO requires ghost atoms store velocity"); diff --git a/src/RHEO/pair_rheo.h b/src/RHEO/pair_rheo.h index c43d450b8b..7a47927962 100644 --- a/src/RHEO/pair_rheo.h +++ b/src/RHEO/pair_rheo.h @@ -46,6 +46,8 @@ class PairRHEO : public Pair { int thermal_flag; int interface_flag; + int harmonic_means_flag; + void allocate(); class ComputeRHEOKernel *compute_kernel; diff --git a/src/pair.cpp b/src/pair.cpp index 5d789fbb9b..56a6283afa 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -1246,6 +1246,108 @@ void Pair::ev_tally_xyz(int i, int j, int nlocal, int newton_pair, } } + +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into global or per-atom accumulators + for virial, have delx,dely,delz and fnormal and ftangential +------------------------------------------------------------------------- */ + +void Pair::ev_tally_nt(int i, int j, int nlocal, int newton_pair, + double evdwl, double ecoul, double fn, + double ftx, double fty, double ftz, + double delx, double dely, double delz) +{ + double evdwlhalf,ecoulhalf,epairhalf,v[6]; + + if (eflag_either) { + if (eflag_global) { + if (newton_pair) { + eng_vdwl += evdwl; + eng_coul += ecoul; + } else { + evdwlhalf = 0.5*evdwl; + ecoulhalf = 0.5*ecoul; + if (i < nlocal) { + eng_vdwl += evdwlhalf; + eng_coul += ecoulhalf; + } + if (j < nlocal) { + eng_vdwl += evdwlhalf; + eng_coul += ecoulhalf; + } + } + } + if (eflag_atom) { + epairhalf = 0.5 * (evdwl + ecoul); + if (newton_pair || i < nlocal) eatom[i] += epairhalf; + if (newton_pair || j < nlocal) eatom[j] += epairhalf; + } + } + + if (vflag_either) { + v[0] = delx*delx*fn; + v[1] = dely*dely*fn; + v[2] = delz*delz*fn; + v[3] = delx*dely*fn; + v[4] = delx*delz*fn; + v[5] = dely*delz*fn; + + v[0] += delx*ftx; + v[1] += dely*fty; + v[2] += delz*ftz; + v[3] += delx*fty + dely*ftx; + v[4] += delx*ftz + delz*ftx; + v[5] += dely*ftz + delz*fty; + + if (vflag_global) { + if (newton_pair) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } else { + if (i < nlocal) { + virial[0] += 0.5*v[0]; + virial[1] += 0.5*v[1]; + virial[2] += 0.5*v[2]; + virial[3] += 0.5*v[3]; + virial[4] += 0.5*v[4]; + virial[5] += 0.5*v[5]; + } + if (j < nlocal) { + virial[0] += 0.5*v[0]; + virial[1] += 0.5*v[1]; + virial[2] += 0.5*v[2]; + virial[3] += 0.5*v[3]; + virial[4] += 0.5*v[4]; + virial[5] += 0.5*v[5]; + } + } + } + + if (vflag_atom) { + if (newton_pair || i < nlocal) { + vatom[i][0] += 0.5*v[0]; + vatom[i][1] += 0.5*v[1]; + vatom[i][2] += 0.5*v[2]; + vatom[i][3] += 0.5*v[3]; + vatom[i][4] += 0.5*v[4]; + vatom[i][5] += 0.5*v[5]; + } + if (newton_pair || j < nlocal) { + vatom[j][0] += 0.5*v[0]; + vatom[j][1] += 0.5*v[1]; + vatom[j][2] += 0.5*v[2]; + vatom[j][3] += 0.5*v[3]; + vatom[j][4] += 0.5*v[4]; + vatom[j][5] += 0.5*v[5]; + } + } + } +} + /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global or per-atom accumulators for virial, have delx,dely,delz and fx,fy,fz diff --git a/src/pair.h b/src/pair.h index 885a2c45ff..6533c7b124 100644 --- a/src/pair.h +++ b/src/pair.h @@ -295,6 +295,8 @@ class Pair : protected Pointers { void ev_tally_tip4p(int, int *, double *, double, double); void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double, double); + void ev_tally_nt(int, int, int, int, double, double, double, double, double, double, double, + double, double); void v_tally2(int, int, double, double *); void v_tally_tensor(int, int, int, int, double, double, double, double, double, double); void virial_fdotr_compute();