netcf: define units for all variable where this is possible

src/NETCDF/dump_netcdf.cpp

@@ -101,6 +101,7 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) :
     int idim = 0;
     int ndims = 1;
     std::string mangled = earg[i];
+    Quantity quantity = Quantity::Unknown;
     bool constant = false;
 
     // name mangling
@@ -109,26 +110,32 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) :
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "coordinates";
+      quantity = Quantity::Distance;
     } else if ((mangled == "vx") || (mangled == "vy") || (mangled == "vz")) {
       idim = mangled[1] - 'x';
       ndims = 3;
       mangled = "velocities";
+      quantity = Quantity::Velocity;
     } else if ((mangled == "xs") || (mangled == "ys") || (mangled == "zs")) {
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "scaled_coordinates";
+      // no unit for scaled coordinates
     } else if ((mangled == "xu") || (mangled == "yu") || (mangled == "zu")) {
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "unwrapped_coordinates";
+      quantity = Quantity::Distance;
     } else if ((mangled == "fx") || (mangled == "fy") || (mangled == "fz")) {
       idim = mangled[1] - 'x';
       ndims = 3;
       mangled = "forces";
+      quantity = Quantity::Force;
     } else if ((mangled == "mux") || (mangled == "muy") || (mangled == "muz")) {
       idim = mangled[2] - 'x';
       ndims = 3;
       mangled = "mu";
+      quantity = Quantity::DipoleMoment;
     } else if (utils::strmatch(mangled, "^c_")) {
       std::size_t found = mangled.find('[');
       if (found != std::string::npos) {
@@ -175,6 +182,7 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) :
     perat[inc].constant = constant;
     perat[inc].ndumped = 0;
     perat[inc].field[idim] = i;
+    perat[inc].quantity = quantity;
   }
 
   n_buffer = 0;
@@ -341,7 +349,6 @@ void DumpNetCDF::openfile()
 
       int dims[NC_MAX_VAR_DIMS];
       size_t index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
-      float d[1];
 
       if (singlefile_opened) return;
       singlefile_opened = 1;
@@ -426,6 +433,11 @@ void DumpNetCDF::openfile()
             NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims, &perat[i].var), perat[i].name );
           }
         }
+
+        std::string unit = this->unit_of(perat[i].quantity);
+        if (!unit.empty()) {
+          NCERR( nc_put_att_text(ncid, perat[i].var, NC_UNITS_STR, unit.size(), unit.c_str()) );
+        }
       }
 
       // perframe variables
@@ -457,43 +469,18 @@ void DumpNetCDF::openfile()
       NCERR( nc_put_att_text(ncid, NC_GLOBAL, "program", 6, "LAMMPS") );
       NCERR( nc_put_att_text(ncid, NC_GLOBAL, "programVersion",strlen(lmp->version), lmp->version) );
 
-      // units
+      // units & scale
-      if (!strcmp(update->unit_style, "lj")) {
+      std::string unit = this->unit_of(Quantity::Time);
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 2, "lj") );
+      NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 2, "lj") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 2, "lj") );
-      } else if (!strcmp(update->unit_style, "real")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") );
-      } else if (!strcmp(update->unit_style, "metal")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 10, "picosecond") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") );
-      } else if (!strcmp(update->unit_style, "si")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 5, "meter") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 5, "meter") );
-      } else if (!strcmp(update->unit_style, "cgs")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 10, "centimeter") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 10, "centimeter") );
-      } else if (!strcmp(update->unit_style, "electron")) {
-        NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") );
-        NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 4, "Bohr") );
-        NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 4, "Bohr") );
-      } else {
-        error->all(FLERR,"Unsupported unit style: {}", update->unit_style);
-      }
 
-      NCERR( nc_put_att_text(ncid, cell_angles_var, NC_UNITS_STR,6, "degree") );
+      unit = this->unit_of(Quantity::Distance);
+      NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
+      NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
 
-      d[0] = update->dt;
+      NCERR( nc_put_att_text(ncid, cell_angles_var, NC_UNITS_STR, 6, "degree") );
-      NCERR( nc_put_att_float(ncid, time_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, d) );
+
-      d[0] = 1.0;
+      float scale[1] = {static_cast<float>(update->dt)};
-      NCERR( nc_put_att_float(ncid, cell_origin_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, d) );
+      NCERR( nc_put_att_float(ncid, time_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, scale) );
-      d[0] = 1.0;
-      NCERR( nc_put_att_float(ncid, cell_lengths_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, d) );
 
       /*
        * Finished with definition
@@ -905,4 +892,94 @@ void DumpNetCDF::ncerr(int err, const char *descr, int line)
   }
 }
 
+/* ---------------------------------------------------------------------- */
+
+std::string DumpNetCDF::unit_of(Quantity quantity) {
+  if (!strcmp(update->unit_style, "lj")) {
+    if (quantity == Quantity::Unknown) {
+      return "";
+    } else {
+      return "lj";
+    }
+  } else if (!strcmp(update->unit_style, "real")) {
+    switch (quantity) {
+    case Quantity::Unknown:
+      return "";
+    case Quantity::Time:
+      return "femtosecond";
+    case Quantity::Distance:
+      return "Angstrom";
+    case Quantity::Velocity:
+      return "Angstrom/femtosecond";
+    case Quantity::Force:
+      return "(Kcal/mol)/Angstrom)";
+    case Quantity::DipoleMoment:
+      return "e * Angstrom";
+    }
+  } else if (!strcmp(update->unit_style, "metal")) {
+    switch (quantity) {
+    case Quantity::Unknown:
+      return "";
+    case Quantity::Time:
+      return "picosecond";
+    case Quantity::Distance:
+      return "Angstrom";
+    case Quantity::Velocity:
+      return "Angstrom/picosecond";
+    case Quantity::Force:
+      return "eV/Angstrom";
+    case Quantity::DipoleMoment:
+      return "e * Angstrom";
+    }
+  } else if (!strcmp(update->unit_style, "si")) {
+    switch (quantity) {
+    case Quantity::Unknown:
+      return "";
+    case Quantity::Time:
+      return "second";
+    case Quantity::Distance:
+      return "meter";
+    case Quantity::Velocity:
+      return "meter/second";
+    case Quantity::Force:
+      return "Newton";
+    case Quantity::DipoleMoment:
+      return "Coulomb * meter";
+    }
+  } else if (!strcmp(update->unit_style, "cgs")) {
+    switch (quantity) {
+    case Quantity::Unknown:
+      return "";
+    case Quantity::Time:
+      return "second";
+    case Quantity::Distance:
+      return "centimeter";
+    case Quantity::Velocity:
+      return "centimeter/second";
+    case Quantity::Force:
+      return "dynes";
+    case Quantity::DipoleMoment:
+      return "statcoul * cm";
+    }
+  } else if (!strcmp(update->unit_style, "electron")) {
+    switch (quantity) {
+    case Quantity::Unknown:
+      return "";
+    case Quantity::Time:
+      return "second";
+    case Quantity::Distance:
+      return "centimeter";
+    case Quantity::Velocity:
+      return "centimeter/second";
+    case Quantity::Force:
+      return "Hartrees/Bohr";
+    case Quantity::DipoleMoment:
+      return "Debye";
+    }
+  }
+
+  error->all(FLERR, "Unsupported unit style: {}", update->unit_style);
+  return "";
+}
+
 #endif /* defined(LMP_HAS_NETCDF) */

src/NETCDF/dump_netcdf.h

@@ -40,12 +40,26 @@ class DumpNetCDF : public DumpCustom {
   virtual void write();
 
  private:
+   // type of quantity for per-atom values (used to get the unit)
+   enum Quantity {
+      Unknown = 0,
+
+      Time,
+      Distance,
+      Velocity,
+      Force,
+      DipoleMoment,
+   };
+   // get the name of the unit for the given quantity
+   std::string unit_of(Quantity quantity);
+
   // per-atoms quantities (positions, velocities, etc.)
   struct nc_perat_t {
     int dims;                        // number of dimensions
     int field[DUMP_NC_MAX_DIMS];     // field indices corresponding to the dim.
     char name[NC_FIELD_NAME_MAX];    // field name
     int var;                         // NetCDF variable
+    Quantity quantity;               // type of the quantity
 
     bool constant;    // is this property per file (not per frame)
     int ndumped;      // number of enties written for this prop.

src/NETCDF/dump_netcdf_mpiio.cpp

@@ -102,6 +102,7 @@ DumpNetCDFMPIIO::DumpNetCDFMPIIO(LAMMPS *lmp, int narg, char **arg) :
     int ndims = 1;
     std::string mangled = earg[i];
     bool constant = false;
+    Quantity quantity = Quantity::Unknown;
 
     // name mangling
     // in the AMBER specification
@@ -109,26 +110,32 @@ DumpNetCDFMPIIO::DumpNetCDFMPIIO(LAMMPS *lmp, int narg, char **arg) :
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "coordinates";
+      quantity = Quantity::Distance;
     } else if ((mangled == "vx") || (mangled == "vy") || (mangled == "vz")) {
       idim = mangled[1] - 'x';
       ndims = 3;
       mangled = "velocities";
+      quantity = Quantity::Velocity;
     } else if ((mangled == "xs") || (mangled == "ys") || (mangled == "zs")) {
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "scaled_coordinates";
+      // no quantity for scaled_coordinates
     } else if ((mangled == "xu") || (mangled == "yu") || (mangled == "zu")) {
       idim = mangled[0] - 'x';
       ndims = 3;
       mangled = "unwrapped_coordinates";
+      quantity = Quantity::Distance;
     } else if ((mangled == "fx") || (mangled == "fy") || (mangled == "fz")) {
       idim = mangled[1] - 'x';
       ndims = 3;
       mangled = "forces";
+      quantity = Quantity::Force;
     } else if ((mangled == "mux") || (mangled == "muy") || (mangled == "muz")) {
       idim = mangled[2] - 'x';
       ndims = 3;
       mangled = "mu";
+      quantity = Quantity::DipoleMoment;
     } else if (utils::strmatch(mangled, "^c_")) {
       std::size_t found = mangled.find('[');
       if (found != std::string::npos) {
@@ -173,6 +180,7 @@ DumpNetCDFMPIIO::DumpNetCDFMPIIO(LAMMPS *lmp, int narg, char **arg) :
     }
 
     perat[inc].field[idim] = i;
+    perat[inc].quantity = quantity;
   }
 
   n_buffer = 0;
@@ -414,6 +422,11 @@ void DumpNetCDFMPIIO::openfile()
         dims[2] = vector_dim[perat[i].dims];
         NCERRX( ncmpi_def_var(ncid, perat[i].name, xtype, 3, dims, &perat[i].var), perat[i].name );
       }
+
+      std::string unit = this->unit_of(perat[i].quantity);
+      if (!unit.empty()) {
+        NCERR( ncmpi_put_att_text(ncid, perat[i].var, NC_UNITS_STR, unit.size(), unit.c_str()) );
+      }
     }
 
     // perframe variables
@@ -441,43 +454,18 @@ void DumpNetCDFMPIIO::openfile()
     NCERR( ncmpi_put_att_text(ncid, NC_GLOBAL, "program", 6, "LAMMPS") );
     NCERR( ncmpi_put_att_text(ncid, NC_GLOBAL, "programVersion", strlen(lmp->version), lmp->version) );
 
-    // units
+    // units & scale
-    if (!strcmp(update->unit_style, "lj")) {
+    std::string unit = this->unit_of(Quantity::Time);
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 2, "lj") );
+    NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 2, "lj") );
+
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 2, "lj") );
+    unit = this->unit_of(Quantity::Distance);
-    } else if (!strcmp(update->unit_style, "real")) {
+    NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") );
+    NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, unit.size(), unit.c_str()) );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") );
-    } else if (!strcmp(update->unit_style, "metal")) {
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 10, "picosecond") );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") );
-    } else if (!strcmp(update->unit_style, "si")) {
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 5, "meter") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 5, "meter") );
-    } else if (!strcmp(update->unit_style, "cgs")) {
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 10, "centimeter") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 10, "centimeter") );
-    } else if (!strcmp(update->unit_style, "electron")) {
-      NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") );
-      NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 4, "Bohr") );
-      NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 4, "Bohr") );
-    } else {
-      error->all(FLERR,"Unsupported unit style: {}", update->unit_style);
-    }
 
     NCERR( ncmpi_put_att_text(ncid, cell_angles_var, NC_UNITS_STR, 6, "degree") );
 
-    d[0] = update->dt;
+    float scale[1] = {static_cast<float>(update->dt)};
-    NCERR( ncmpi_put_att_float(ncid, time_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, d) );
+    NCERR( ncmpi_put_att_float(ncid, time_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, scale) );
-    d[0] = 1.0;
-    NCERR( ncmpi_put_att_float(ncid, cell_origin_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, d) );
-    d[0] = 1.0;
-    NCERR( ncmpi_put_att_float(ncid, cell_lengths_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, d) );
 
     /*
      * Finished with definition
@@ -900,4 +888,84 @@ void DumpNetCDFMPIIO::ncerr(int err, const char *descr, int line)
   }
 }
 
+/* ---------------------------------------------------------------------- */
+
+std::string DumpNetCDF::unit_of(Quantity quantity) {
+  if (quantity == Quantity::Unknown) {
+    return "";
+  }
+
+  if (!strcmp(update->unit_style, "lj")) {
+    return "lj";
+  } else if (!strcmp(update->unit_style, "real")) {
+    switch (quantity) {
+    case Quantity::Time:
+      return "femtosecond";
+    case Quantity::Distance:
+      return "Angstrom";
+    case Quantity::Velocity:
+      return "Angstrom/femtosecond";
+    case Quantity::Force:
+      return "(Kcal/mol)/Angstrom)";
+    case Quantity::DipoleMoment:
+      return "e * Angstrom";
+    }
+  } else if (!strcmp(update->unit_style, "metal")) {
+    switch (quantity) {
+    case Quantity::Time:
+      return "picosecond";
+    case Quantity::Distance:
+      return "Angstrom";
+    case Quantity::Velocity:
+      return "Angstrom/picosecond";
+    case Quantity::Force:
+      return "eV/Angstrom";
+    case Quantity::DipoleMoment:
+      return "e * Angstrom";
+    }
+  } else if (!strcmp(update->unit_style, "si")) {
+    switch (quantity) {
+    case Quantity::Time:
+      return "second";
+    case Quantity::Distance:
+      return "meter";
+    case Quantity::Velocity:
+      return "meter/second";
+    case Quantity::Force:
+      return "Newton";
+    case Quantity::DipoleMoment:
+      return "Coulomb * meter";
+    }
+  } else if (!strcmp(update->unit_style, "cgs")) {
+    switch (quantity) {
+    case Quantity::Time:
+      return "second";
+    case Quantity::Distance:
+      return "centimeter";
+    case Quantity::Velocity:
+      return "centimeter/second";
+    case Quantity::Force:
+      return "dynes";
+    case Quantity::DipoleMoment:
+      return "statcoul * cm";
+    }
+  } else if (!strcmp(update->unit_style, "electron")) {
+    switch (quantity) {
+    case Quantity::Time:
+      return "second";
+    case Quantity::Distance:
+      return "centimeter";
+    case Quantity::Velocity:
+      return "centimeter/second";
+    case Quantity::Force:
+      return "Hartrees/Bohr";
+    case Quantity::DipoleMoment:
+      return "Debye";
+    }
+  }
+
+  error->all(FLERR, "Unsupported unit style: {}", update->unit_style);
+  return "";
+}
+
 #endif /* defined(LMP_HAS_PNETCDF) */

src/NETCDF/dump_netcdf_mpiio.h

@@ -40,12 +40,26 @@ class DumpNetCDFMPIIO : public DumpCustom {
   virtual void write();
 
  private:
+   // type of quantity for per-atom values (used to get the unit)
+   enum Quantity {
+      Unknown = 0,
+
+      Time,
+      Distance,
+      Velocity,
+      Force,
+      DipoleMoment,
+   };
+   // get the name of the unit for the given quantity
+   std::string unit_of(Quantity quantity);
+
   // per-atoms quantities (positions, velocities, etc.)
   struct nc_perat_t {
     int dims;                              // number of dimensions
     int field[DUMP_NC_MPIIO_MAX_DIMS];     // field indices corresponding to the dim.
     char name[NC_MPIIO_FIELD_NAME_MAX];    // field name
     int var;                               // NetCDF variable
+    Quantity quantity;                     // type of the quantity
   };
 
   typedef void (DumpNetCDFMPIIO::*funcptr_t)(void *);