Adding multiphase support, new stress

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<FixRHEOThermal *>(fixes[0]);
   }
+
+  if (pressure_flag) {
+    fixes = modify->get_fix_by_style("rheo/pressure");
+    fix_pressure = dynamic_cast<FixRHEOPressure *>(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;

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;

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;
 }

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 (surface_flag) {
-
+    compute_surface->compute_peratom();
-  if (shift_flag)
+    if (shift_flag)
-    compute_vshift->correct_surfaces();
+      compute_vshift->correct_surfaces();
+  }
 }
 
 /* ---------------------------------------------------------------------- */

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;
 }

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;

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"

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;

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 * dT_prefactor;
-          heatflow[i] += dT;
 
           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 * dT_prefactor;
-            heatflow[j] -= dT;
           }
         }
 
@@ -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
+          if (evflag) {
-            ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], dx[0], dx[1], dx[2]);
+            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");

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;

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

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();