adding documentation and integration fix

src/.gitignore

@@ -421,14 +421,14 @@
 /atom_vec_full.h
 /atom_vec_full_hars.cpp
 /atom_vec_full_hars.h
-/atom_vec_granular.cpp
-/atom_vec_granular.h
 /atom_vec_molecular.cpp
 /atom_vec_molecular.h
 /atom_vec_oxdna.cpp
 /atom_vec_oxdna.h
 /atom_vec_peri.cpp
 /atom_vec_peri.h
+/atom_vec_sphere_temperature.cpp
+/atom_vec_sphere_temperature.h
 /atom_vec_template.cpp
 /atom_vec_template.h
 /body_nparticle.cpp
@@ -1510,6 +1510,8 @@
 /fix_smd_wall_surface.h
 /fix_srp.cpp
 /fix_srp.h
+/fix_temp_integrate.cpp
+/fix_temp_integrate.h
 /fix_tfmc.cpp
 /fix_tfmc.h
 /fix_ttm.cpp

src/GRANULAR/atom_vec_sphere_temperature.cpp

@@ -0,0 +1,160 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   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 "atom_vec_sphere_temperature.h"
+
+#include "atom.h"
+#include "error.h"
+#include "fix.h"
+#include "fix_adapt.h"
+#include "math_const.h"
+#include "modify.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+AtomVecSphereTemperature::AtomVecSphereTemperature(LAMMPS *lmp) : AtomVec(lmp)
+{
+  mass_type = PER_ATOM;
+  molecular = Atom::ATOMIC;
+
+  atom->sphere_flag = 1;
+  atom->radius_flag = atom->rmass_flag = atom->omega_flag = atom->torque_flag = 1;
+  atom->temperature_flag = atom->heatflux_flag = 1;
+
+  // strings with peratom variables to include in each AtomVec method
+  // strings cannot contain fields in corresponding AtomVec default strings
+  // order of fields in a string does not matter
+  // except: fields_data_atom & fields_data_vel must match data file
+
+  fields_grow = {"radius", "rmass", "omega", "torque", "temperature", "heatflux"};
+  fields_copy = {"radius", "rmass", "omega", "temperature"};
+  fields_comm_vel = {"omega", "temperature"};
+  fields_reverse = {"torque", "heatflux"};
+  fields_border = {"radius", "rmass", "temperature"};
+  fields_border_vel = {"radius", "rmass", "omega", "temperature"};
+  fields_exchange = {"radius", "rmass", "omega", "temperature"};
+  fields_restart = {"radius", "rmass", "omega", "temperature"};
+  fields_create = {"radius", "rmass", "omega", "temperature"};
+  fields_data_atom = {"id", "type", "radius", "rmass", "x", "temperature"};
+  fields_data_vel = {"id", "v", "omega"};
+}
+
+/* ----------------------------------------------------------------------
+   process sub-style args
+   optional arg = 0/1 for static/dynamic particle radii
+------------------------------------------------------------------------- */
+
+void AtomVecSphereTemperature::process_args(int narg, char **arg)
+{
+  if (narg != 0 && narg != 1) error->all(FLERR, "Illegal atom_style sphere command");
+
+  radvary = 0;
+  if (narg == 1) {
+    radvary = utils::numeric(FLERR, arg[0], true, lmp);
+    if (radvary < 0 || radvary > 1) error->all(FLERR, "Illegal atom_style sphere command");
+  }
+
+  // dynamic particle radius and mass must be communicated every step
+
+  if (radvary) {
+    fields_comm = {"radius", "rmass"};
+    fields_comm_vel = {"radius", "rmass", "omega"};
+  }
+
+  // delay setting up of fields until now
+
+  setup_fields();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecSphereTemperature::init()
+{
+  AtomVec::init();
+
+  // check if optional radvary setting should have been set to 1
+
+  for (int i = 0; i < modify->nfix; i++)
+    if (strcmp(modify->fix[i]->style, "adapt") == 0) {
+      auto fix = dynamic_cast<FixAdapt *>(modify->fix[i]);
+      if (fix->diamflag && radvary == 0)
+        error->all(FLERR, "Fix adapt changes particle radii but atom_style sphere is not dynamic");
+    }
+}
+
+/* ----------------------------------------------------------------------
+   set local copies of all grow ptrs used by this class, except defaults
+   needed in replicate when 2 atom classes exist and it calls pack_restart()
+------------------------------------------------------------------------- */
+
+void AtomVecSphereTemperature::grow_pointers()
+{
+  radius = atom->radius;
+  rmass = atom->rmass;
+  omega = atom->omega;
+}
+
+/* ----------------------------------------------------------------------
+   initialize non-zero atom quantities
+------------------------------------------------------------------------- */
+
+void AtomVecSphereTemperature::create_atom_post(int ilocal)
+{
+  radius[ilocal] = 0.5;
+  rmass[ilocal] = 4.0 * MY_PI / 3.0 * 0.5 * 0.5 * 0.5;
+}
+
+/* ----------------------------------------------------------------------
+   modify what AtomVec::data_atom() just unpacked
+   or initialize other atom quantities
+------------------------------------------------------------------------- */
+
+void AtomVecSphereTemperature::data_atom_post(int ilocal)
+{
+  radius_one = 0.5 * atom->radius[ilocal];
+  radius[ilocal] = radius_one;
+  if (radius_one > 0.0) rmass[ilocal] *= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
+
+  if (rmass[ilocal] <= 0.0) error->one(FLERR, "Invalid density in Atoms section of data file");
+
+  omega[ilocal][0] = 0.0;
+  omega[ilocal][1] = 0.0;
+  omega[ilocal][2] = 0.0;
+}
+
+/* ----------------------------------------------------------------------
+   modify values for AtomVec::pack_data() to pack
+------------------------------------------------------------------------- */
+
+void AtomVecSphereTemperature::pack_data_pre(int ilocal)
+{
+  radius_one = radius[ilocal];
+  rmass_one = rmass[ilocal];
+
+  radius[ilocal] *= 2.0;
+  if (radius_one != 0.0)
+    rmass[ilocal] = rmass_one / (4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one);
+}
+
+/* ----------------------------------------------------------------------
+   unmodify values packed by AtomVec::pack_data()
+------------------------------------------------------------------------- */
+
+void AtomVecSphereTemperature::pack_data_post(int ilocal)
+{
+  radius[ilocal] = radius_one;
+  rmass[ilocal] = rmass_one;
+}

src/GRANULAR/atom_vec_sphere_temperature.h

@@ -0,0 +1,51 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   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 ATOM_CLASS
+// clang-format off
+AtomStyle(sphere,AtomVecSphereTemperature);
+// clang-format on
+#else
+
+#ifndef LMP_ATOM_VEC_SPHERE_TEMPERATURE_H
+#define LMP_ATOM_VEC_SPHERE_TEMPERATURE_H
+
+#include "atom_vec.h"
+
+namespace LAMMPS_NS {
+
+class AtomVecSphereTemperature : public AtomVec {
+ public:
+  AtomVecSphereTemperature(class LAMMPS *);
+  void process_args(int, char **) override;
+  void init() override;
+
+  void grow_pointers() override;
+  void create_atom_post(int) override;
+  void data_atom_post(int) override;
+  void pack_data_pre(int) override;
+  void pack_data_post(int) override;
+
+ private:
+  double *radius, *rmass;
+  double **omega;
+  double *temperature, *heatflux;
+
+  int radvary;
+  double radius_one, rmass_one;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/GRANULAR/fix_temp_integrate.cpp

@@ -0,0 +1,142 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   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 "fix_temp_integrate.h"
+
+#include "atom.h"
+#include "error.h"
+#include "force.h"
+#include "memory.h"
+#include "respa.h"
+#include "update.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+enum {NONE, CONSTANT, TYPE};
+
+/* ---------------------------------------------------------------------- */
+
+FixTempIntegrate::FixTempIntegrate(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 4) error->all(FLERR,"Illegal fix command");
+
+  cp_style = NONE;
+
+  int ntypes = atom->ntypes;
+  int iarg = 3;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg+1],"constant") == 0) {
+      if (iarg+2 >= narg) error->all(FLERR,"Illegal fix command");
+      cp_style = CONSTANT;
+      cp = utils::numeric(FLERR,arg[iarg+2],false,lmp);
+      if (cp < 0.0) error->all(FLERR,"Illegal fix command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg+1],"type") == 0) {
+      if (iarg+1+ntypes >= narg) error->all(FLERR,"Illegal fix command");
+      cp_style = TYPE;
+      memory->create(cp_type,ntypes+1,"fix/temp/integrate:cp_type");
+      for (int i = 1; i <= ntypes; i++) {
+        cp_type[i] = utils::numeric(FLERR,arg[iarg+1+i],false,lmp);
+        if (cp_type[i] < 0.0) error->all(FLERR,"Illegal fix command");
+      }
+      iarg += 1+ntypes;
+    } else {
+      error->all(FLERR,"Illegal fix command");
+    }
+    iarg += 1;
+  }
+
+  if (cp_style == NONE)
+    error->all(FLERR, "Must specify specific heat in fix temp/integrate");
+  dynamic_group_allow = 1;
+  time_integrate = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixTempIntegrate::setmask()
+{
+  int mask = 0;
+  mask |= FINAL_INTEGRATE;
+  mask |= FINAL_INTEGRATE_RESPA;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTempIntegrate::init()
+{
+  dt = update->dt;
+
+  if (!atom->temperature_flag)
+    error->all(FLERR,"Fix temp/integrate requires atom style with temperature property");
+  if (!atom->heatflux_flag)
+    error->all(FLERR,"Fix temp/integrate requires atom style with heatflux property");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTempIntegrate::final_integrate()
+{
+  // update temperature of atoms in group
+
+  double *temperature = atom->temperature;
+  double *heatflux = atom->heatflux;
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  if (rmass) {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+        temperature[i] += dt * heatflux[i] / (calc_cp(i) * rmass[i]);
+      }
+  } else {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+        temperature[i] += dt * heatflux[i] / (calc_cp(i) * mass[type[i]]);
+      }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTempIntegrate::final_integrate_respa(int ilevel, int /*iloop*/)
+{
+  dt = update->dt;
+  final_integrate();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTempIntegrate::reset_dt()
+{
+  dt = update->dt;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixTempIntegrate::calc_cp(int i)
+{
+  if (cp_style == TYPE) {
+    return cp_type[atom->type[i]];
+  } else {
+    return cp;
+  }
+}

src/GRANULAR/fix_temp_integrate.h

@@ -0,0 +1,50 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   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 FIX_CLASS
+// clang-format off
+FixStyle(temp/integrate,FixTempIntegrate);
+// clang-format on
+#else
+
+#ifndef LMP_FIX_TEMP_INTEGRATE_H
+#define LMP_FIX_TEMP_INTEGRATE_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixTempIntegrate : public Fix {
+ public:
+  FixTempIntegrate(class LAMMPS *, int, char **);
+
+  int setmask() override;
+  void init() override;
+  void final_integrate() override;
+  void final_integrate_respa(int, int) override;
+  void reset_dt() override;
+
+ protected:
+  double dt;
+  double cp, *cp_type;
+  int cp_style;
+
+  int mass_require;
+
+  double calc_cp(int);
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/atom.cpp

@@ -123,6 +123,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   radius = rmass = nullptr;
   ellipsoid = line = tri = body = nullptr;
   quat = nullptr;
+  temperature = nullptr;
+  heatflux = nullptr;
 
   // molecular systems
 
@@ -410,6 +412,9 @@ void Atom::peratom_create()
   add_peratom("tri",&tri,INT,0);
   add_peratom("body",&body,INT,0);
 
+  add_peratom("temperature",&temperature,DOUBLE,0);
+  add_peratom("heatflux",&heatflux,DOUBLE,0);
+
   // BPM package
 
   add_peratom("quat",&quat,DOUBLE,4);
@@ -612,6 +617,7 @@ void Atom::set_atomflag_defaults()
   molecule_flag = molindex_flag = molatom_flag = 0;
   q_flag = mu_flag = 0;
   rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
+  temperature_flag = heatflux_flag = 0;
   vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
   cs_flag = csforce_flag = vforce_flag = ervelforce_flag = etag_flag = 0;
   rho_flag = esph_flag = cv_flag = vest_flag = 0;
@@ -2602,6 +2608,14 @@ length of the data area, and a short description.
      - double
      - 4
      - four quaternion components of the particles
+   * - temperature
+     - double
+     - 1
+     - temperature of the particles
+   * - heatflux
+     - double
+     - 1
+     - heatflux of the particles
    * - i_name
      - int
      - 1
@@ -2658,6 +2672,8 @@ void *Atom::extract(const char *name)
   if (strcmp(name,"tri") == 0) return (void *) tri;
   if (strcmp(name,"body") == 0) return (void *) body;
   if (strcmp(name,"quat") == 0) return (void *) quat;
+  if (strcmp(name,"temperature") == 0) return (void *) temperature;
+  if (strcmp(name,"heatflux") == 0) return (void *) heatflux;
 
   if (strcmp(name,"vfrac") == 0) return (void *) vfrac;
   if (strcmp(name,"s0") == 0) return (void *) s0;
@@ -2781,6 +2797,8 @@ int Atom::extract_datatype(const char *name)
   if (strcmp(name,"tri") == 0) return LAMMPS_INT;
   if (strcmp(name,"body") == 0) return LAMMPS_INT;
   if (strcmp(name,"quat") == 0) return LAMMPS_DOUBLE_2D;
+  if (strcmp(name,"temperature") == 0) return LAMMPS_DOUBLE;
+  if (strcmp(name,"heatflux") == 0) return LAMMPS_DOUBLE;
 
   if (strcmp(name,"vfrac") == 0) return LAMMPS_DOUBLE;
   if (strcmp(name,"s0") == 0) return LAMMPS_DOUBLE;

src/atom.h

@@ -80,6 +80,7 @@ class Atom : protected Pointers {
   double **omega, **angmom, **torque;
   int *ellipsoid, *line, *tri, *body;
   double **quat;
+  double *temperature, *heatflux;
 
   // molecular systems
 
@@ -182,6 +183,7 @@ class Atom : protected Pointers {
   int molecule_flag, molindex_flag, molatom_flag;
   int q_flag, mu_flag;
   int rmass_flag, radius_flag, omega_flag, torque_flag, angmom_flag, quat_flag;
+  int temperature_flag, heatflux_flag;
   int vfrac_flag, spin_flag, eradius_flag, ervel_flag, erforce_flag;
   int cs_flag, csforce_flag, vforce_flag, ervelforce_flag, etag_flag;
   int rho_flag, esph_flag, cv_flag, vest_flag;

src/set.cpp

@@ -47,7 +47,7 @@ enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
 enum{TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET,
      MOLECULE,X,Y,Z,VX,VY,VZ,CHARGE,MASS,SHAPE,LENGTH,TRI,
      DIPOLE,DIPOLE_RANDOM,SPIN,SPIN_RANDOM,QUAT,QUAT_RANDOM,
-     THETA,THETA_RANDOM,ANGMOM,OMEGA,
+     THETA,THETA_RANDOM,ANGMOM,OMEGA,TEMPERATURE,
      DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
      SPH_E,SPH_CV,SPH_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY,
      SMD_CONTACT_RADIUS,DPDTHETA,EPSILON,IVEC,DVEC,IARRAY,DARRAY};
@@ -389,6 +389,15 @@ void Set::command(int narg, char **arg)
       set(DENSITY);
       iarg += 2;
 
+    } else if (strcmp(arg[iarg],"temperature") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
+      if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1);
+      else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      if (!atom->temperature_flag)
+        error->all(FLERR,"Cannot set this attribute for this atom style");
+      set(TEMPERATURE);
+      iarg += 2;
+
     } else if (strcmp(arg[iarg],"volume") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
       if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1);
@@ -768,6 +777,7 @@ void Set::set(int keyword)
   case SHAPE:
   case DIAMETER:
   case DENSITY:
+  case TEMPERATURE:
   case QUAT:
   case IMAGE:
     if (modify->check_rigid_list_overlap(select))
@@ -1008,6 +1018,13 @@ void Set::set(int keyword)
       atom->omega[i][2] = zvalue;
     }
 
+    // set temperature of particle
+
+    else if (keyword == ANGMOM) {
+      if (dvalue < 0.0) error->one(FLERR,"Invalid temperature in set command");
+      atom->temperature[i] = dvalue;
+    }
+
     // reset any or all of 3 image flags
 
     else if (keyword == IMAGE) {