diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H
index a1db14173..cd9dd788c 100644
--- a/expui/BasisFactory.H
+++ b/expui/BasisFactory.H
@@ -141,6 +141,9 @@ namespace BasisClasses
     //! Get the current pseudo acceleration value
     Eigen::Vector3d currentAccel(double time);
 
+    //! Gravitational constant
+    double G = 1.0;
+
   public:
 
     //! The current pseudo acceleration
diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc
index e8b0ccde5..bf1d47355 100644
--- a/expui/BiorthBasis.cc
+++ b/expui/BiorthBasis.cc
@@ -409,6 +409,7 @@ namespace BasisClasses
     if (expcoef.rows()>0 && expcoef.cols()>0) expcoef.setZero();
     totalMass = 0.0;
     used = 0;
+    G = 1.0;
   }
   
   
@@ -422,6 +423,8 @@ namespace BasisClasses
     cf->time   = time;
     cf->normed = true;
 
+    G = cf->G;
+
     // Angular storage dimension
     int ldim = (lmax+1)*(lmax+2)/2;
 
@@ -477,6 +480,10 @@ namespace BasisClasses
     
     CoefClasses::SphStruct* cf = dynamic_cast<CoefClasses::SphStruct*>(coef.get());
 
+    // Set gravitational constant
+    //
+    G = cf->G;
+
     // Cache the current coefficient structure
     //
     coefret = coef;
@@ -675,7 +682,7 @@ namespace BasisClasses
     }
     
     double densfac = 1.0/(scale*scale*scale) * 0.25/M_PI;
-    double potlfac = 1.0/scale;
+    double potlfac = G/scale;	// Grav constant G is 1 by default
 
     return
       {den0 * densfac,		// 0
@@ -781,7 +788,7 @@ namespace BasisClasses
       }
     }
     
-    double potlfac = 1.0/scale;
+    double potlfac = G/scale;	// Grav constant G is 1 by default
 
     potr *= (-potlfac)/scale;
     pott *= (-potlfac);
@@ -1590,6 +1597,12 @@ namespace BasisClasses
     double tdens0, tdens, tpotl0, tpotl, tpotR, tpotz, tpotp;
     
     sl->accumulated_eval(R, z, phi, tpotl0, tpotl, tpotR, tpotz, tpotp);
+
+    tpotl0 *= G;		// Apply gravitational constant to
+    tpotl  *= G;		// potential and forces
+    tpotR  *= G;
+    tpotz  *= G;
+    tpotp  *= G;
     
     tdens = sl->accumulated_dens_eval(R, z, phi, tdens0);
 
@@ -1610,6 +1623,12 @@ namespace BasisClasses
     double tdens0, tdens, tpotl0, tpotl, tpotR, tpotz, tpotp;
 
     sl->accumulated_eval(R, z, phi, tpotl0, tpotl, tpotR, tpotz, tpotp);
+
+    tpotl0 *= G;		// Apply G to potential and forces
+    tpotl  *= G;
+    tpotR  *= G;
+    tpotz  *= G;
+    tpotp  *= G;
     
     tdens = sl->accumulated_dens_eval(R, z, phi, tdens0);
 
@@ -1637,7 +1656,8 @@ namespace BasisClasses
     double tpotx = tpotR*x/R - tpotp*y/R ;
     double tpoty = tpotR*y/R + tpotp*x/R ;
 
-    acc << tpotx, tpoty, tpotz;
+    // Apply G to forces on return
+    acc << tpotx*G, tpoty*G, tpotz*G;
   }
 
   // Evaluate in cylindrical coordinates
@@ -1648,6 +1668,12 @@ namespace BasisClasses
     sl->accumulated_eval(R, z, phi, tpotl0, tpotl, tpotR, tpotz, tpotp);
     tdens = sl->accumulated_dens_eval(R, z, phi, tdens0);
 
+    tpotl0 *= G;		// Apply G to potential and forces
+    tpotl  *= G;
+    tpotR  *= G;
+    tpotz  *= G;
+    tpotp  *= G;
+
     if (midplane) {
       height = sl->accumulated_midplane_eval(R, -colh*hcyl, colh*hcyl, phi);
       return
@@ -1682,6 +1708,8 @@ namespace BasisClasses
     cf->nmax   = nmax;
     cf->time   = time;
 
+    G = cf->G;
+
     Eigen::VectorXd cos1(nmax), sin1(nmax);
 
     // Initialize the values
@@ -1711,6 +1739,10 @@ namespace BasisClasses
 
     CoefClasses::CylStruct* cf = dynamic_cast<CoefClasses::CylStruct*>(coef.get());
 
+    // Set gravitational constant
+    //
+    G = cf->G;
+
     // Cache the current coefficient structure
     //
     coefret = coef;
@@ -1954,6 +1986,7 @@ namespace BasisClasses
     if (expcoef.rows()>0 && expcoef.cols()>0) expcoef.setZero();
     totalMass = 0.0;
     used = 0;
+    G = 1.0;
   }
   
   
@@ -1965,6 +1998,8 @@ namespace BasisClasses
     cf->nmax   = nmax;
     cf->time   = time;
 
+    G = cf->G;
+
     // Allocate the coefficient storage
     cf->store.resize((mmax+1)*nmax);
 
@@ -2012,6 +2047,10 @@ namespace BasisClasses
     CoefClasses::CylStruct* cf = dynamic_cast<CoefClasses::CylStruct*>(coef.get());
     auto & cc = *cf->coefs;
 
+    // Set gravitational constant
+    //
+    G = cf->G;
+
     // Cache the current coefficient structure
     //
     coefret = coef;
@@ -2120,9 +2159,10 @@ namespace BasisClasses
     if (R>ortho->getRtable() or fabs(z)>ortho->getRtable()) {
       double r2 = R*R + z*z;
       double r  = sqrt(r2);
-      pot0 = -totalMass/r;
-      rpot = -totalMass*R/(r*r2 + 10.0*std::numeric_limits<double>::min());
-      zpot = -totalMass*z/(r*r2 + 10.0*std::numeric_limits<double>::min());
+      // Apply G to potential and forces
+      pot0 = -G*totalMass/r;
+      rpot = -G*totalMass*R/(r*r2 + 10.0*std::numeric_limits<double>::min());
+      zpot = -G*totalMass*z/(r*r2 + 10.0*std::numeric_limits<double>::min());
       
       return {den0, den1, den0+den1, pot0, pot1, pot0+pot1, rpot, zpot, ppot};
     }
@@ -2195,11 +2235,11 @@ namespace BasisClasses
 
     den0 *= -1.0;
     den1 *= -1.0;
-    pot0 *= -1.0;
-    pot1 *= -1.0;
-    rpot *= -1.0;
-    zpot *= -1.0;
-    ppot *= -1.0;
+    pot0 *= -G;			// Apply G to potential and forces
+    pot1 *= -G;
+    rpot *= -G;
+    zpot *= -G;
+    ppot *= -G;
 
     return {den0, den1, den0+den1, pot0, pot1, pot0+pot1, rpot, zpot, ppot};
   }
@@ -2226,8 +2266,9 @@ namespace BasisClasses
       double r2 = R*R + z*z;
       double r  = sqrt(r2);
 
-      rpot = -totalMass*R/(r*r2 + 10.0*std::numeric_limits<double>::min());
-      zpot = -totalMass*z/(r*r2 + 10.0*std::numeric_limits<double>::min());
+      // Apply G to forces
+      rpot = -G*totalMass*R/(r*r2 + 10.0*std::numeric_limits<double>::min());
+      zpot = -G*totalMass*z/(r*r2 + 10.0*std::numeric_limits<double>::min());
       
       acc << rpot, zpot, ppot;
     }
@@ -2286,9 +2327,9 @@ namespace BasisClasses
       }
     }
 
-    rpot *= -1.0;
-    zpot *= -1.0;
-    ppot *= -1.0;
+    rpot *= -G;			// Apply G to forces
+    zpot *= -G;
+    ppot *= -G;
 
     double potx = rpot*x/R - ppot*y/R;
     double poty = rpot*y/R + ppot*x/R;
@@ -2756,6 +2797,7 @@ namespace BasisClasses
     if (expcoef.rows()>0 && expcoef.cols()>0) expcoef.setZero();
     totalMass = 0.0;
     used = 0;
+    G = 1.0;
   }
   
   
@@ -2767,6 +2809,8 @@ namespace BasisClasses
     cf->nmax   = nmax;
     cf->time   = time;
 
+    G = cf->G;
+
     // Allocate the coefficient storage
     cf->store.resize((mmax+1)*nmax);
 
@@ -2814,6 +2858,10 @@ namespace BasisClasses
     CoefClasses::CylStruct* cf = dynamic_cast<CoefClasses::CylStruct*>(coef.get());
     auto & cc = *cf->coefs;
 
+    // Set gravitational constant
+    //
+    G = cf->G;
+
     // Cache the current coefficient structure
     coefret = coef;
 
@@ -2971,10 +3019,10 @@ namespace BasisClasses
 
     den0 *= -1.0;
     den1 *= -1.0;
-    pot0 *= -1.0;
-    pot1 *= -1.0;
-    rpot *= -1.0;
-    ppot *= -1.0;
+    pot0 *= -G;			// Apply G to potential and forces
+    pot1 *= -G;
+    rpot *= -G;
+    ppot *= -G;
 
     return {den0, den1, den0+den1, pot0, pot1, pot0+pot1, rpot, zpot, ppot};
   }
@@ -3039,9 +3087,9 @@ namespace BasisClasses
       }
     }
 
-    rpot *= -1.0;
-    ppot *= -1.0;
-
+    // Apply G to forces
+    rpot *= -G;
+    ppot *= -G;
 
     double potx = rpot*x/R - ppot*y/R;
     double poty = rpot*y/R + ppot*x/R;
@@ -3242,6 +3290,7 @@ namespace BasisClasses
     expcoef.setZero();
     totalMass = 0.0;
     used = 0;
+    G = 1.0;
   }
   
   
@@ -3254,6 +3303,8 @@ namespace BasisClasses
     cf->nmaxz   = nmaxz;
     cf->time    = time;
 
+    G = cf->G;
+
     cf->allocate();
 
     *cf->coefs = expcoef;
@@ -3287,6 +3338,10 @@ namespace BasisClasses
     auto cf = dynamic_cast<CoefClasses::SlabStruct*>(coef.get());
     expcoef = *cf->coefs;
 
+    // Set gravitational constant
+    //
+    G = cf->G;
+
     // Cache the current coefficient structure
     //
     coefret = coef;
@@ -3527,7 +3582,8 @@ namespace BasisClasses
       }
     }
 
-    acc << accx.real(), accy.real(), accz.real();
+    // Apply G to forces on return
+    acc << G*accx.real(), G*accy.real(), G*accz.real();
   }
 
 
@@ -3538,7 +3594,8 @@ namespace BasisClasses
 
     auto [pot, den, frcx, frcy, frcz] = eval(x, y, z);
 
-    return {0, den, den, 0, pot, pot, frcx, frcy, frcz};
+    // Apply G to potential and forces on return
+    return {0, den, den, 0, pot*G, pot*G, frcx*G, frcy*G, frcz*G};
   }
 
   std::vector<double> Slab::cyl_eval(double R, double z, double phi)
@@ -3555,11 +3612,12 @@ namespace BasisClasses
     double potp = -frcx*sin(phi) + frcy*cos(phi);
     double potz =  frcz;
 
-    potR *= -1;
-    potp *= -1;
-    potz *= -1;
+    potR *= -G;			// Apply G to forces
+    potp *= -G;
+    potz *= -G;
 
-    return {0, den, den, 0, pot, pot, potR, potz, potp};
+    // Apply Go to potential on return
+    return {0, den, den, 0, pot*G, pot*G, potR, potz, potp};
   }
 
   std::vector<double> Slab::sph_eval(double r, double costh, double phi)
@@ -3577,11 +3635,12 @@ namespace BasisClasses
     double pott =  frcx*cos(phi)*costh + frcy*sin(phi)*costh - frcz*sinth;
     double potp = -frcx*sin(phi)       + frcy*cos(phi);
 
-    potr *= -1;
-    pott *= -1;
-    potp *= -1;
+    potr *= -G;			// Apply G to forces
+    pott *= -G;
+    potp *= -G;
     
-    return {0, den, den, 0, pot, pot, potr, pott, potp};
+    // Apply G to potential on return
+    return {0, den, den, 0, pot*G, pot*G, potr, pott, potp};
   }
 
 
@@ -3755,6 +3814,7 @@ namespace BasisClasses
     expcoef.setZero();
     totalMass = 0.0;
     used = 0;
+    G = 1.0;
   }
   
   
@@ -3767,6 +3827,8 @@ namespace BasisClasses
     cf->nmaxz   = nmaxz;
     cf->time    = time;
 
+    G = cf->G;
+
     cf->allocate();
 
     *cf->coefs = expcoef;
@@ -3800,6 +3862,10 @@ namespace BasisClasses
     auto cf = dynamic_cast<CoefClasses::CubeStruct*>(coef.get());
     expcoef = *cf->coefs;
 
+    // Set gravitational constant
+    //
+    G = cf->G;
+
     // Cache the cuurent coefficient structure
     //
     coefret = coef;
@@ -3889,7 +3955,8 @@ namespace BasisClasses
     double frcy = -frc(1).real();
     double frcz = -frc(2).real();
 
-    return {0, den1, den1, 0, pot1, pot1, frcx, frcy, frcz};
+    // Apply G to potential and forces on return
+    return {0, den1, den1, 0, pot1*G, pot1*G, frcx*G, frcy*G, frcz*G};
   }
 
   void Cube::computeAccel(double x, double y, double z,
@@ -3904,7 +3971,8 @@ namespace BasisClasses
     // Get the basis fields
     auto frc = ortho->get_force(expcoef, pos);
     
-    acc << -frc(0).real(), -frc(1).real(), -frc(2).real();
+    // Apply G to forces on return
+    acc << -G*frc(0).real(), -G*frc(1).real(), -G*frc(2).real();
   }
 
   std::vector<double> Cube::cyl_eval(double R, double z, double phi)
@@ -3920,10 +3988,12 @@ namespace BasisClasses
 
     // Get the basis fields
     double den1 = ortho->get_dens(expcoef, pos).real();
-    double pot1 = ortho->get_pot (expcoef, pos).real();
+    double pot1 = ortho->get_pot (expcoef, pos).real() * G;
 
-    auto frc = ortho->get_force(expcoef, pos);
+    auto frc = ortho->get_force(expcoef, pos) * G;
     
+    // Gravitational constant G applied to potenial and forces above
+
     double frcx = frc(0).real(), frcy = frc(1).real(), frcz = frc(2).real();
 
     double potR =  frcx*cos(phi) + frcy*sin(phi);
@@ -3951,9 +4021,11 @@ namespace BasisClasses
 
     // Get the basis fields
     double den1 = ortho->get_dens(expcoef, pos).real();
-    double pot1 = ortho->get_pot (expcoef, pos).real();
+    double pot1 = ortho->get_pot (expcoef, pos).real() * G;
 
-    auto frc = ortho->get_force(expcoef, pos);
+    auto frc = ortho->get_force(expcoef, pos) * G;
+
+    // Gravitational constant G applied to potential and forces above
     
     double frcx = frc(0).real();
     double frcy = frc(1).real();
@@ -3963,9 +4035,9 @@ namespace BasisClasses
     double pott =  frcx*cos(phi)*costh + frcy*sin(phi)*costh - frcz*sinth;
     double potp = -frcx*sin(phi)       + frcy*cos(phi);
 
-    potr *= -1;
-    pott *= -1;
-    potp *= -1;
+    potr *= -1.0;
+    pott *= -1.0;
+    potp *= -1.0;
     
     return {0, den1, den1, 0, pot1, pot1, potr, pott, potp};
   }
diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt
index 30102d192..a6bcab62f 100644
--- a/expui/CMakeLists.txt
+++ b/expui/CMakeLists.txt
@@ -1,4 +1,4 @@
-set(bin_PROGRAMS nativetoh5 h5compare viewcoefs h5power makecoefs testread)
+set(bin_PROGRAMS nativetoh5 h5compare viewcoefs h5power makecoefs testread testunits)
 
 set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp exputil
   ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES})
@@ -38,7 +38,7 @@ endif()
 set(expui_SOURCES BasisFactory.cc BiorthBasis.cc FieldBasis.cc
   CoefContainer.cc CoefStruct.cc FieldGenerator.cc expMSSA.cc
   Coefficients.cc KMeans.cc Centering.cc ParticleIterator.cc
-  Koopman.cc BiorthBess.cc SvdSignChoice.cc)
+  Koopman.cc BiorthBess.cc SvdSignChoice.cc UnitValidator.cc)
 add_library(expui ${expui_SOURCES})
 set_target_properties(expui PROPERTIES OUTPUT_NAME expui)
 target_include_directories(expui PUBLIC ${common_INCLUDE})
@@ -111,6 +111,7 @@ if(ENABLE_UTILS)
   add_executable(h5power    h5power.cc)
   add_executable(makecoefs  makecoefs.cc)
   add_executable(testread   testread.cc)
+  add_executable(testunits  testunits.cc)
 
   foreach(program ${bin_PROGRAMS})
     target_link_libraries(${program} expui exputil ${common_LINKLIB})
diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H
index b7ee8ac7b..106e65ae2 100644
--- a/expui/CoefStruct.H
+++ b/expui/CoefStruct.H
@@ -25,6 +25,9 @@ namespace CoefClasses
     //! Time stamp
     double time;
     
+    //! Gravitational constant
+    double G = 1.0;
+
     //! Coefficient data
     Eigen::VectorXcd store;
     
@@ -88,6 +91,9 @@ namespace CoefClasses
     //! Read-only access to center (no copy)
     double getTime() { return time; }
 
+    //! Set gravitational constant
+    void setGravConstant(double val) { G = val; }
+
   };
   
   //! Specialization of CoefStruct for spheres
diff --git a/expui/Coefficients.H b/expui/Coefficients.H
index 65f738eb7..98fbed4f6 100644
--- a/expui/Coefficients.H
+++ b/expui/Coefficients.H
@@ -6,8 +6,9 @@
 
 // Needed by member functions for writing parameters and stanzas
 #include <highfive/H5File.hpp>
+#include <highfive/H5Group.hpp>
 
-#include <Eigen/Eigen>	      // Store coefficient matrices and 2d grids
+#include <Eigen/Eigen>	    // Store coefficient matrices and 2d grids
 #include <unsupported/Eigen/CXX11/Tensor> // For 3d rectangular grids?
 
 // YAML support
@@ -16,6 +17,10 @@
 // The EXP native coefficient classes
 #include "CoefStruct.H"
 
+// Unit validation
+#include <UnitValidator.H>
+
+// All coefficient classes are in this namespace
 namespace CoefClasses
 { 
   //! An index key
@@ -24,9 +29,31 @@ namespace CoefClasses
   using E2d = 
   Eigen::Matrix<std::complex<double>,
 		Eigen::Dynamic,	Eigen::Dynamic,	Eigen::ColMajor>;
-    
+  
   using E3d = Eigen::Tensor<std::complex<double>, 3>;
 
+  // A struct for a single unit metadatum
+  struct Unit
+  {
+    char name[16];		//! The name of the physical value in ascii
+    char unit[16];		//! The unit type in ascii
+    float value;			//! The value as a float
+
+    Unit() { std::memset(name, 0, 16); std::memset(unit, 0, 16); value=0.0f; }
+
+    Unit(std::string n, std::string u, float v) : value(v)
+    {
+      // Constant string size
+      const size_t sz=16;
+      // Zero fill
+      std::memset(name, 0, sz); std::memset(unit, 0, sz);
+      // Copy strings to char arrays
+      strncpy(name, n.c_str(), std::min(n.length()+1, sz));
+      strncpy(unit, u.c_str(), std::min(u.length()+1, sz));
+    }
+
+  };
+
   /** 
       Abstract class for any type of coefficient database
       
@@ -87,6 +114,12 @@ namespace CoefClasses
     //! Time offset for interpolation
     double deltaT;
 
+    //! Default units (EXP native)
+    std::vector<Unit> units {{"G", "none", 1.0f}};
+
+    //! Unit validator
+    static UnitValidator check;
+
   public:
     
     //! Constructor
@@ -198,6 +231,51 @@ namespace CoefClasses
     //! Set maximum grid interpolation offset
     void setDeltaT(double dT) { deltaT = dT; }
 
+    //! Set units
+    void setUnits
+    (const std::string name, const std::string unit, const float value);
+
+    //! Set units using a vector of tuples (name, unit, value)
+    void setUnits
+    (const std::vector<std::tuple<std::string, std::string, float>>& units);
+
+    //! Remove a unit by name
+    void removeUnits(const std::string name);
+
+    //! Write units to H5
+    void WriteH5Units(HighFive::File& file);
+
+    //! Read units from H5
+    void ReadH5Units(HighFive::File& file);
+
+    //! Get units
+    std::vector<std::tuple<std::string, std::string, float>> getUnits();
+
+    //! Get allowed unit types
+    const std::vector<std::string> getAllowedUnitTypes() const
+    {
+      return check.getAllowedTypes();
+    }
+    
+    //! Get allowed aliases for each type
+    const std::vector<std::string>
+    getAllowedTypeAliases(const std::string& type) const
+    {
+      return check.getAllowedTypeAliases(type);
+    }
+
+    //! Get allowed units for a given type
+    const std::vector<std::string>
+    getAllowedUnitNames(const std::string& type) const
+    {
+      auto ret = check.getAllowedUnits(type);
+      std::sort(ret.begin(), ret.end());
+      return ret;
+    }
+
+    //! Get gravitational constant
+    double getGravConstant();
+
     class CoefsError : public std::runtime_error
     {
     public:
diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 6f91a1129..b8b743473 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -5,6 +5,7 @@
 #include <sstream>
 #include <iomanip>
 #include <complex>
+#include <cstring>
 
 #include "config_exp.h"
 #include "localmpi.H"
@@ -17,9 +18,48 @@
 
 #include "Coefficients.H"
 
+// Create a helper function to describe the compound type
+HighFive::CompoundType create_compound_Unit() {
+  return {
+    {"name",  HighFive::AtomicType<char[16]>()},
+    {"unit",  HighFive::AtomicType<char[16]>()},
+    {"value", HighFive::AtomicType<float>()}
+  };
+}
+
+// Register the type with HighFive.
+HIGHFIVE_REGISTER_TYPE(CoefClasses::Unit, create_compound_Unit)
+
+
+// Helper ostream manipulator for debugging Unit info
+std::ostream& operator<< (std::ostream& out,
+			  const std::vector<CoefClasses::Unit>& t)
+{
+  out << std::string(48, '-') << std::endl
+      << std::setw(16) << "Name"
+      << std::setw(16) << "Unit"
+      << std::setw(16) << "Value"
+      << std::endl
+      << std::setw(16) << std::string(10, '-')
+      << std::setw(16) << std::string(10, '-')
+      << std::setw(16) << std::string(10, '-')
+      << std::endl;
+  for (auto p : t)
+    out << std::setw(16) << p.name
+	<< std::setw(16) << p.unit
+	<< std::setw(16) << p.value
+	<< std::endl;
+  out << std::string(48, '-') << std::endl;
+  
+  return out;
+}
+
 namespace CoefClasses
 {
-  
+
+  // Static instance of the unit validator
+  UnitValidator Coefs::check = UnitValidator();
+
   void Coefs::copyfields(std::shared_ptr<Coefs> p)
   {
     // These variables will copy data, not pointers
@@ -31,6 +71,115 @@ namespace CoefClasses
     p->times    = times;
   }
 
+  void Coefs::removeUnits(const std::string name)
+  {
+    // Lambda test for matching name
+    auto test = [&name](Unit& elem) -> bool{return elem.name==name; };
+
+    // Explanation: std::remove_if shifts elements to be removed to
+    // the end of the range and returns an iterator to the new logical
+    // end. units.erase() then removes the elements from that point to
+    // the end.
+    units.erase(std::remove_if(units.begin(), units.end(), test), units.end()); 
+  }
+  
+
+  void Coefs::setUnits
+  (const std::string Name, const std::string Unit, const float Value)
+  {
+    auto copy_s = [](std::string& s, char* c) -> void {
+      const size_t sz = 16;
+      strncpy(c, s.c_str(), std::min(s.length()+1, sz));
+    };
+      
+    // Run the type validation
+    //
+    auto [valid, name, unit] = check(Name, Unit);
+
+    if (not valid) {
+      throw std::runtime_error(std::string("Coefs::setUnits: Warning, type '")
+			       + Name + "' with unit '" + Unit +
+			       "' is incompatible or not recognized.");
+    }
+
+    // Check for existing unit and update
+    //
+    for (auto & p : units) {
+      std::string p_name(p.name);
+
+      if (name == p_name) {
+	copy_s(unit, &p.unit[0]);
+	p.value = Value;
+	return;
+      }
+    }
+
+    // No matches
+    //
+    units.push_back({name, unit, Value});
+  }
+
+  // Set units using a vector of tuples (name, unit, value)
+  // Converts to a list of tuples in Python
+  void Coefs::setUnits
+  (const std::vector<std::tuple<std::string, std::string, float>>& units)
+  {
+    for (auto p : units) {
+      setUnits(std::get<0>(p), std::get<1>(p), std::get<2>(p));
+    }
+  }
+
+  //! Get units
+  std::vector<std::tuple<std::string, std::string, float>> Coefs::getUnits()
+  {
+    std::vector<std::tuple<std::string, std::string, float>> ret;
+    for (auto p : units) ret.push_back({p.name, p.unit, p.value});
+    return ret;
+  }
+
+
+  //! Get gravitational constant
+  double Coefs::getGravConstant()
+  {
+    for (auto p : units) {
+      std::string G(p.name);
+      if (G == "G") return p.value;
+    }
+    // Default if "G" is removed or not set for some unforseen reason
+    return 1.0;
+  }
+
+  void Coefs::WriteH5Units(HighFive::File& file)
+  {
+    if (units.size() != 4) {
+      std::ostringstream sout;
+      sout << "---- Coefs::WriteH5Units: Warning, expected 4 units: "
+	   << "(length, mass, time, G) or (length, mass, velocity, G), etc. "
+	   << "I found " << units.size() << " units instead.  Please "
+	   << " provide a consistent unit set.";
+      throw std::runtime_error(sout.str());
+    }
+
+    HighFive::DataSet dataset = file.createDataSet("Units", units);
+
+    if (units.size() == 4) {
+      std::cout << "Coefs::WriteH5Units: wrote units to HDF5 file:" << std::endl
+		<< units << std::endl;
+    }
+  }
+
+  void Coefs::ReadH5Units(HighFive::File& file)
+  {
+    if (file.exist("Units")) {
+      HighFive::DataSet dataset = file.getDataSet("Units");
+      units = dataset.read<std::vector<Unit>>();
+      if (verbose and myid==0) {
+	std::cout << "Coefs::ReadH5Units: read units from HDF5 file:" << std::endl;
+	std::cout << units;
+      }
+    }
+  }
+
   std::tuple<Eigen::VectorXcd&, bool> Coefs::interpolate(double time)
   {
     bool onGrid = true;
@@ -82,6 +231,10 @@ namespace CoefClasses
     unsigned count;
     double scale;
     
+    // Check for units
+    ReadH5Units(file);
+    double G = getGravConstant();
+
     file.getAttribute("name"    ).read(name    );
     file.getAttribute("lmax"    ).read(Lmax    );
     file.getAttribute("nmax"    ).read(Nmax    );
@@ -157,6 +310,7 @@ namespace CoefClasses
       coef->scale = scale;
       coef->geom  = geometry;
       coef->id    = forceID;
+      coef->setGravConstant(G);
 
       coef->allocate();
       *coef->coefs = in;
@@ -181,8 +335,10 @@ namespace CoefClasses
       ret->coefs[v.first] =
 	std::dynamic_pointer_cast<SphStruct>(v.second->deepcopy());
 
-    ret->Lmax = Lmax;
-    ret->Nmax = Nmax;
+    ret->Lmax  = Lmax;
+    ret->Nmax  = Nmax;
+    ret->units = units;
+
 
     return ret;
   }
@@ -203,6 +359,7 @@ namespace CoefClasses
     ret->Mmax  = Mmax;
     ret->Nmax  = Nmax;
     ret->angle = angle;
+    ret->units = units;
 
     return ret;
   }
@@ -223,6 +380,7 @@ namespace CoefClasses
     ret->NmaxX  = NmaxX;
     ret->NmaxY  = NmaxY;
     ret->NmaxZ  = NmaxZ;
+    ret->units  = units;
 
     return ret;
   }
@@ -243,6 +401,7 @@ namespace CoefClasses
     ret->NmaxX  = NmaxX;
     ret->NmaxY  = NmaxY;
     ret->NmaxZ  = NmaxZ;
+    ret->units  = units;
 
     return ret;
   }
@@ -295,6 +454,10 @@ namespace CoefClasses
     unsigned count;
     double scale;
     
+    // Check for units
+    ReadH5Units(file);
+    double G = getGravConstant();
+
     file.getAttribute("name"    ).read(name    );
     file.getAttribute("nfld"    ).read(Nfld    );
     file.getAttribute("lmax"    ).read(Lmax    );
@@ -354,6 +517,7 @@ namespace CoefClasses
       coef->scale = scale;
       coef->geom  = geometry;
       coef->id    = fieldID;
+      coef->setGravConstant(G);
 
       coef->allocate();
       coef->store = in;
@@ -394,6 +558,10 @@ namespace CoefClasses
     unsigned count;
     double scale;
     
+    // Check for units
+    ReadH5Units(file);
+    double G = getGravConstant();
+
     file.getAttribute("name"    ).read(name    );
     file.getAttribute("nfld"    ).read(Nfld    );
     file.getAttribute("mmax"    ).read(Mmax    );
@@ -454,6 +622,7 @@ namespace CoefClasses
       coef->scale = scale;
       coef->geom  = geometry;
       coef->id    = fieldID;
+      coef->setGravConstant(G);
 
       coef->allocate();
       coef->store = in;
@@ -659,6 +828,8 @@ namespace CoefClasses
   
   void SphCoefs::WriteH5Params(HighFive::File& file)
   {
+    WriteH5Units(file);
+
     double scale = coefs.begin()->second->scale;
     
     std::string forceID(coefs.begin()->second->id);
@@ -841,6 +1012,10 @@ namespace CoefClasses
     unsigned count;
     std::string config;
     
+    ReadH5Units(file);
+    double G = getGravConstant();
+
+
     file.getAttribute("name"   ).read(name  );
     file.getAttribute("mmax"   ).read(Mmax  );
     file.getAttribute("nmax"   ).read(Nmax  );
@@ -913,6 +1088,7 @@ namespace CoefClasses
       coef->rot = rot;
       coef->assign(in, Mmax, Nmax);
       coef->time = Time;
+      coef->setGravConstant(G);
       
       coefs[roundTime(Time)] = coef;
     }
@@ -1068,6 +1244,8 @@ namespace CoefClasses
   
   void CylCoefs::WriteH5Params(HighFive::File& file)
   {
+    WriteH5Units(file);
+
     std::string forceID(coefs.begin()->second->id);
 
     file.createAttribute<int>("mmax", HighFive::DataSpace::From(Mmax)).write(Mmax);
@@ -1247,6 +1425,10 @@ namespace CoefClasses
     unsigned count;
     std::string config;
     
+    // Check for units
+    ReadH5Units(file);
+    double G = getGravConstant();
+
     file.getAttribute("name"   ).read(name  );
     file.getAttribute("nmaxx"  ).read(NmaxX );
     file.getAttribute("nmaxy"  ).read(NmaxY );
@@ -1288,6 +1470,7 @@ namespace CoefClasses
       
       coef->assign(dat);
       coef->time = Time;
+      coef->setGravConstant(G);
       
       coefs[roundTime(Time)] = coef;
     }
@@ -1395,6 +1578,8 @@ namespace CoefClasses
 
   void SlabCoefs::WriteH5Params(HighFive::File& file)
   {
+    WriteH5Units(file);
+
     std::string forceID(coefs.begin()->second->id);
 
     file.createAttribute<int>("nmaxx", HighFive::DataSpace::From(NmaxX)).write(NmaxX);
@@ -1600,6 +1785,10 @@ namespace CoefClasses
     unsigned count;
     std::string config;
     
+    // Check for units
+    ReadH5Units(file);
+    double G = getGravConstant();
+
     file.getAttribute("name"   ).read(name  );
     file.getAttribute("nmaxx"  ).read(NmaxX );
     file.getAttribute("nmaxy"  ).read(NmaxY );
@@ -1641,7 +1830,8 @@ namespace CoefClasses
       
       coef->assign(dat);
       coef->time = Time;
-      
+      coef->setGravConstant(G);
+    
       coefs[roundTime(Time)] = coef;
     }
 
@@ -1748,6 +1938,8 @@ namespace CoefClasses
 
   void CubeCoefs::WriteH5Params(HighFive::File& file)
   {
+    WriteH5Units(file);
+
     std::string forceID(coefs.begin()->second->id);
 
     file.createAttribute<int>("nmaxx", HighFive::DataSpace::From(NmaxX)).write(NmaxX);
@@ -2021,6 +2213,7 @@ namespace CoefClasses
       // Pack the data into the coefficient variable
       //
       auto coef = std::make_shared<TrajStruct>();
+
       coef->traj  = traj;
       coef->rank  = rank;
       coef->time  = times[n];
@@ -2098,6 +2291,8 @@ namespace CoefClasses
   
   void TrajectoryData::WriteH5Params(HighFive::File& file)
   {
+    WriteH5Units(file);
+
     int traj = coefs.begin()->second->traj;
     int rank = coefs.begin()->second->rank;
     
@@ -2455,6 +2650,8 @@ namespace CoefClasses
       HighFive::Attribute geom = h5file.getAttribute("geometry");
       geom.read(geometry);
       
+      // Now try to deduce the coefficient type
+      //
       try {
 	// Is the set a biorthogonal basis (has the forceID attribute)
 	// or general basis (fieldID attribute)?
@@ -2492,6 +2689,10 @@ namespace CoefClasses
 	throw std::runtime_error(msg + err.what());
       }
 	
+      // Attempt to red units
+      //
+      coefs->ReadH5Units(h5file);
+
       return coefs;
       
     } catch (HighFive::Exception& err) {
@@ -2686,6 +2887,10 @@ namespace CoefClasses
       //
       HighFive::File file(prefix, HighFive::File::ReadWrite);
       
+      // Attempt to read units
+      //
+      ReadH5Units(file);
+
       // Get the dataset
       HighFive::DataSet dataset = file.getDataSet("count");
       
@@ -2936,6 +3141,8 @@ namespace CoefClasses
     
     double scale = coefs.begin()->second->scale;
 
+    WriteH5Units(file);
+
     // Write the remaining parameters
     //
     file.createAttribute<int>   ("nfld",  HighFive::DataSpace::From(Nfld)  ).write(Nfld);
@@ -3057,6 +3264,8 @@ namespace CoefClasses
     std::string fieldID("polar velocity orthgonal function coefficients");
     file.createAttribute<std::string>("fieldID", HighFive::DataSpace::From(fieldID)).write(fieldID);
     
+    WriteH5Units(file);
+
     double scale = coefs.begin()->second->scale;
 
     // Write the remaining parameters
diff --git a/expui/UnitValidator.H b/expui/UnitValidator.H
new file mode 100644
index 000000000..18603edf9
--- /dev/null
+++ b/expui/UnitValidator.H
@@ -0,0 +1,73 @@
+#include <unordered_map>
+#include <vector>
+#include <string>
+#include <tuple>
+#include <map>
+
+/** Class to validate unit types and names The allowed types and units
+    can be extended as needed by adding new names and aliases to
+    UnitValidator::createAllowedUnitTypes() and
+    UnitValidator::createAllowedUnitNames() in UnitValidator.cc
+
+    Strategy:
+    -------
+    - Types and their aliases are stored in an std::unordered_map
+    - Units and their aliases are stored in an std::map of
+      std::unordered_maps
+    - The operator()(type, unit, value) first checks if the provided
+      type exists in the allowed types map.  If it does, it retrieves
+      the canonical type string which is used to look up the allowed
+      units for that type in the std::map.  Then, the unit type is
+      checked against the allowed units for that type.  If both type
+      and unit are valid, the canonical strings are returned.
+
+    Usage:
+    -----
+    - Create an instance of UnitValidator.  Currently, this is done in
+      the Coefficients class by default.
+    - Call the operator() on an instance with the type and unit
+      strings to be validated
+    - The operator() returns a tuple with:
+       + a boolean indicating if the type and unit are valid
+       + the canonical type string, or "unknown" if invalid
+       + the canonical unit string, or "unknown" if invalid
+*/
+class UnitValidator
+{
+private:
+
+  //! Store allowed unit types and their aliases
+  std::unordered_map<std::string, std::string> allowed_types;
+
+  //! Store allowed unit names and their aliases
+  std::map<std::string, std::unordered_map<std::string, std::string>> allowed_units;
+
+  //! Create the allowed types map (length, mass, time, velocity, G
+  //! [and aliases])
+  std::unordered_map<std::string, std::string> createAllowedUnitTypes();
+
+  //! Create the allowed units map
+  std::map<std::string, std::unordered_map<std::string, std::string>> createAllowedUnitNames();
+
+public:
+
+  //! Constructor
+  UnitValidator();
+
+  //! Destructor
+  ~UnitValidator() {}
+
+  //! Do the full check and return canonical strings
+  std::tuple<bool, std::string, std::string>
+  operator()(const std::string& type, const std::string& unit);
+
+  //! Get allowed unit types
+  const std::vector<std::string> getAllowedTypes() const;
+
+  //! Get allowed aliases for each type
+  const std::vector<std::string> getAllowedTypeAliases(const std::string& type) const;
+
+  //! Get allowed units for a given type and their aliases
+  const std::vector<std::string> getAllowedUnits(const std::string& type) const;
+
+};
diff --git a/expui/UnitValidator.cc b/expui/UnitValidator.cc
new file mode 100644
index 000000000..d3321494f
--- /dev/null
+++ b/expui/UnitValidator.cc
@@ -0,0 +1,235 @@
+#include <algorithm>
+#include "UnitValidator.H"
+
+// Constructor/initializer
+UnitValidator::UnitValidator()
+{
+  // Initialize the 'dictionaries'
+  allowed_types = createAllowedUnitTypes();
+  allowed_units = createAllowedUnitNames();
+}
+
+// Do the full check and return canonical strings in two steps.  First
+// check the type.  Then the unit within that type.  Returns a tuple of
+// (is_valid, canonical_type, canonical_unit).
+std::tuple<bool, std::string, std::string>
+UnitValidator::operator()(const std::string& type, const std::string& unit)
+{
+  // Check type first
+  if (allowed_types.count(type) > 0) {
+
+    // Get canonical type
+    std::string canonical_type = allowed_types.at(type);
+    
+    // Now check unit in the type category
+    if (allowed_units[canonical_type].count(unit) > 0) {
+
+      // Get canonical unit name
+      std::string canonical_unit = allowed_units[canonical_type].at(unit);
+
+      // Return successful final results
+      return {true, canonical_type, canonical_unit};
+    }
+  }
+
+  // If we get here, we have a type or unit that is not recognized
+  return {false, "unknown", "unknown"};
+}
+
+
+// Map aliases to their canonical (primary) unit types
+std::unordered_map<std::string, std::string>
+UnitValidator::createAllowedUnitTypes()
+{
+  std::unordered_map<std::string, std::string> allowed;
+
+  // These are the canonical names
+  //
+  allowed["length"]   = "length";
+  allowed["mass"]     = "mass";
+  allowed["time"]     = "time";
+  allowed["velocity"] = "velocity";
+  allowed["G"]        = "G";
+  
+  // These are recognized aliases aliases for the types
+  //
+  allowed["Length"]   = "length";
+  allowed["Len"]      = "length";
+  allowed["len"]      = "length";
+  allowed["l"]        = "length";
+  allowed["L"]        = "length";
+
+  allowed["Mass"]     = "mass";
+  allowed["m"]        = "mass";
+  allowed["M"]        = "mass";
+
+  allowed["Time"]     = "time";
+  allowed["t"]        = "time";
+  allowed["T"]        = "time";
+
+  allowed["vel"]      = "velocity";
+  allowed["Vel"]      = "velocity";
+  allowed["Velocity"] = "velocity";
+  allowed["v"]        = "velocity";
+  allowed["V"]        = "velocity";
+  
+  allowed["Grav"]     = "G";
+  allowed["grav"]     = "G";
+  allowed["grav_constant"] = "G";
+  allowed["Grav_constant"] = "G";
+  allowed["gravitational_constant"] = "G";
+  allowed["Gravitational_constant"] = "G";
+  
+  return allowed;
+}
+
+// Map aliases to their canonical (primary) unit names for each type
+// category
+std::map<std::string, std::unordered_map<std::string, std::string>>
+UnitValidator::createAllowedUnitNames()
+{
+  std::map<std::string, std::unordered_map<std::string, std::string>> allowed;
+
+  // Allow 'none' as a unit for all types
+  //
+  allowed["length"]["none"]     = "none";
+  allowed["mass"]["none"]       = "none";
+  allowed["time"]["none"]       = "none";
+  allowed["length"]["None"]     = "none";
+  allowed["mass"]["None"]       = "none";
+  allowed["velocity"]["none"]   = "none";
+
+  // Special non-units
+  //
+  allowed["G"][""]              = "none";
+  allowed["G"]["mixed"]         = "mixed";
+  allowed["G"]["none"]          = "none";
+  allowed["G"]["unitless"]      = "none";
+
+  // Canonical length units
+  //
+  allowed["length"]["m"]        = "m";
+  allowed["length"]["cm"]       = "cm";
+  allowed["length"]["km"]       = "km";
+  allowed["length"]["um"]       = "um";
+  allowed["length"]["nm"]       = "nm";
+  allowed["length"]["Angstrom"] = "Angstrom";
+  allowed["length"]["AU"]       = "AU";
+  allowed["length"]["ly"]       = "ly";
+  allowed["length"]["pc"]       = "pc";
+  allowed["length"]["kpc"]      = "kpc";
+  allowed["length"]["Mpc"]      = "Mpc";
+
+  // Standard astronomical mass units
+  //
+  allowed["mass"]["Msun"]       = "Msun";
+  allowed["mass"]["Mearth"]     = "Mearth";
+  allowed["mass"]["g"]          = "g";
+  allowed["mass"]["kg"]         = "kg";
+  
+  // Time units
+  //
+  allowed["time"]["s"]          = "s";
+  allowed["time"]["min"]        = "min";
+  allowed["time"]["hr"]         = "hr";
+  allowed["time"]["day"]        = "day";
+  allowed["time"]["yr"]         = "yr";
+  allowed["time"]["Myr"]        = "Myr";
+  allowed["time"]["Gyr"]        = "Gyr";
+	
+  // Velocity units
+  //
+  allowed["velocity"]["cm/s"]     = "cmm/s";
+  allowed["velocity"]["m/s"]      = "m/s";
+  allowed["velocity"]["km/s"]     = "km/s";
+  allowed["velocity"]["km/hr"]    = "km/hr";
+  allowed["velocity"]["km/min"]   = "km/min";
+  allowed["velocity"]["c"]        = "c";
+  allowed["velocity"]["none"]     = "none";
+  
+
+  // Length aliases
+  //
+  allowed["length"]["meter"]      = "m";
+  allowed["length"]["centimeter"] = "cm";
+  allowed["length"]["kilometer"]  = "km";
+  allowed["length"]["nanometer"]  = "nm";
+  allowed["length"]["micrometer"] = "um";
+  allowed["length"]["micron"]     = "um";
+  allowed["length"]["angstrom"]   = "Angstrom";
+  allowed["length"]["AA"]         = "Angstrom";
+  allowed["length"]["astronomical_unit"] = "AU";
+  allowed["length"]["au"]         = "AU";
+  allowed["length"]["light_year"] = "ly";
+  allowed["length"]["lyr"]        = "ly";
+  allowed["length"]["parsec"]     = "pc";
+  allowed["length"]["kiloparsec"] = "kpc";
+  allowed["length"]["megaparsec"] = "Mpc";
+
+
+  // Mass aliases
+  //
+  allowed["mass"]["solar_mass"] = "Msun";
+  allowed["mass"]["earth_mass"] = "Mearth";
+  allowed["mass"]["gram"]       = "g";
+  allowed["mass"]["kilograms"]  = "kg";
+
+  // Time aliases
+  //
+  allowed["time"]["second"]     = "s";
+  allowed["time"]["minute"]     = "min";
+  allowed["time"]["hour"]       = "hr";
+  allowed["time"]["year"]       = "yr";
+  allowed["time"]["None"]       = "none";
+
+  // Velocity aliases
+  //
+  allowed["velocity"]["meter_per_second"]      = "m/s";
+  allowed["velocity"]["centimeter_per_second"] = "cm/s";
+  allowed["velocity"]["cm_per_s"]              = "cm/s";
+  allowed["velocity"]["m_per_s"]               = "m/s";
+  allowed["velocity"]["km_per_s"]              = "km/s";
+  allowed["velocity"]["km_per_hr"]             = "km/hr";
+  allowed["velocity"]["km_per_min"]            = "km/min";
+  allowed["velocity"]["speed_of_light"]        = "c";
+
+
+  return allowed;
+}
+
+// Get allowed aliases for each type
+const std::vector<std::string>
+UnitValidator::getAllowedTypeAliases(const std::string& type) const
+{
+  // Build the alias vector
+  std::vector<std::string> aliases;
+
+  for (const auto& pair : allowed_types) {
+    if (pair.second == type) {
+      aliases.push_back(pair.first);
+    }
+  }
+
+  // Sort the vector
+  std::sort(aliases.begin(), aliases.end());
+
+  return aliases;
+}
+
+// Canonical names for allowed unit types
+const std::vector<std::string> UnitValidator::getAllowedTypes() const
+{
+  return {"mass", "length", "time", "velocity", "G"};
+}
+
+// Get allowed units for a given type (and their aliases)
+const std::vector<std::string> UnitValidator::getAllowedUnits(const std::string& type) const
+{
+  std::vector<std::string> units;
+  if (allowed_units.count(type) > 0) {
+    for (const auto& pair : allowed_units.at(type)) {
+      units.push_back(pair.first);
+    }
+  }
+  return units;
+}
diff --git a/expui/testunits.cc b/expui/testunits.cc
new file mode 100644
index 000000000..e1204c8e6
--- /dev/null
+++ b/expui/testunits.cc
@@ -0,0 +1,37 @@
+// This is a simple test program for the UnitValidator class
+
+#include <iostream>
+#include "UnitValidator.H"
+
+int main()
+{
+  UnitValidator check;
+
+  std::string type, unit;
+  std::cout << "Enter type and unit: ";
+  std::cin >> type >> unit;
+
+  bool valid;
+  std::string canonical_type, canonical_unit;
+  std::tie(valid, canonical_type, canonical_unit) = check(type, unit);
+
+  if (valid) {
+    std::cout << "The type '" << type << "' with unit '" << unit
+	      << "' is valid." << std::endl;
+    std::cout << "The canonical names are: Type='" << canonical_type
+      	      << "', Unit='" << canonical_unit << "'" << std::endl;
+  } else {
+    std::cout << "The type '" << type << "' with unit '" << unit
+	      << "' is not valid." << std::endl;
+  }
+
+  // G test
+  std::tie(valid, canonical_type, canonical_unit) = check("G", "");
+  if (valid) {
+    std::cout << "The type 'G' with units '' is valid." << std::endl;
+  } else {
+    std::cout << "The type 'G' with units '' is not valid." << std::endl;
+  }
+
+  return 0;
+}
diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc
index 81ccdc266..9b24879e4 100644
--- a/pyEXP/CoefWrappers.cc
+++ b/pyEXP/CoefWrappers.cc
@@ -202,6 +202,19 @@ void CoefficientClasses(py::module &m) {
     void zerodata() override {
       PYBIND11_OVERRIDE_PURE(void, Coefs, zerodata,);
     }
+
+    void setUnits(const std::string name, const std::string unit, const float value) {
+      PYBIND11_OVERRIDE(void, Coefs, setUnits, name, unit, value);
+    }
+
+    void setUnits(const std::vector<std::tuple<std::string, std::string, float>>& units) {
+      PYBIND11_OVERRIDE(void, Coefs, setUnits, units);
+    }
+
+    void removeUnits(const std::string name) {
+      PYBIND11_OVERRIDE(void, Coefs, removeUnits, name);
+    }
+
   };
 
   class PySphCoefs : public SphCoefs
@@ -237,7 +250,7 @@ void CoefficientClasses(py::module &m) {
 
     std::shared_ptr<CoefStruct> getCoefStruct(double time) override {
       PYBIND11_OVERRIDE(std::shared_ptr<CoefStruct>, SphCoefs, getCoefStruct,
-			time);
+                        time);
     }
 
     std::vector<double> Times() override {
@@ -261,11 +274,11 @@ void CoefficientClasses(py::module &m) {
     }
 
     void clear() override {
-      PYBIND11_OVERRIDE(void, SphCoefs,	clear,);
+      PYBIND11_OVERRIDE(void, SphCoefs, clear,);
     }
 
     void add(CoefStrPtr coef) override {
-      PYBIND11_OVERRIDE(void, SphCoefs,	add, coef);
+      PYBIND11_OVERRIDE(void, SphCoefs, add, coef);
     }
 
     std::vector<Key> makeKeys(Key k) override {
@@ -287,7 +300,7 @@ void CoefficientClasses(py::module &m) {
   {
   protected:
     void readNativeCoefs(const std::string& file, int stride, double tmin, double tmax) override {
-      PYBIND11_OVERRIDE(void, CylCoefs,	readNativeCoefs, file, stride, tmin, tmax);
+      PYBIND11_OVERRIDE(void, CylCoefs, readNativeCoefs, file, stride, tmin, tmax);
     }
 
     std::string getYAML() override {
@@ -316,7 +329,7 @@ void CoefficientClasses(py::module &m) {
 
     std::shared_ptr<CoefStruct> getCoefStruct(double time) override {
       PYBIND11_OVERRIDE(std::shared_ptr<CoefStruct>, CylCoefs, getCoefStruct,
-			time);
+                        time);
     }
 
     std::vector<double> Times() override {
@@ -336,7 +349,7 @@ void CoefficientClasses(py::module &m) {
     }
 
     bool CompareStanzas(std::shared_ptr<Coefs> check) override {
-      PYBIND11_OVERRIDE(bool, CylCoefs, CompareStanzas,	check);
+      PYBIND11_OVERRIDE(bool, CylCoefs, CompareStanzas, check);
     }
 
     void clear() override {
@@ -344,7 +357,7 @@ void CoefficientClasses(py::module &m) {
     }
 
     void add(CoefStrPtr coef) override {
-      PYBIND11_OVERRIDE(void, CylCoefs,	add, coef);
+      PYBIND11_OVERRIDE(void, CylCoefs, add, coef);
     }
 
     std::vector<Key> makeKeys(Key k) override {
@@ -394,7 +407,7 @@ void CoefficientClasses(py::module &m) {
 
     std::shared_ptr<CoefStruct> getCoefStruct(double time) override {
       PYBIND11_OVERRIDE(std::shared_ptr<CoefStruct>, SlabCoefs, getCoefStruct,
-			time);
+                        time);
     }
 
     std::vector<double> Times() override {
@@ -418,11 +431,11 @@ void CoefficientClasses(py::module &m) {
     }
 
     void clear() override {
-      PYBIND11_OVERRIDE(void, SlabCoefs,	clear,);
+      PYBIND11_OVERRIDE(void, SlabCoefs,        clear,);
     }
 
     void add(CoefStrPtr coef) override {
-      PYBIND11_OVERRIDE(void, SlabCoefs,	add, coef);
+      PYBIND11_OVERRIDE(void, SlabCoefs,        add, coef);
     }
 
     std::vector<Key> makeKeys(Key k) override {
@@ -473,7 +486,7 @@ void CoefficientClasses(py::module &m) {
 
     std::shared_ptr<CoefStruct> getCoefStruct(double time) override {
       PYBIND11_OVERRIDE(std::shared_ptr<CoefStruct>, CubeCoefs, getCoefStruct,
-			time);
+                        time);
     }
 
     std::vector<double> Times() override {
@@ -497,11 +510,11 @@ void CoefficientClasses(py::module &m) {
     }
 
     void clear() override {
-      PYBIND11_OVERRIDE(void, CubeCoefs,	clear,);
+      PYBIND11_OVERRIDE(void, CubeCoefs,        clear,);
     }
 
     void add(CoefStrPtr coef) override {
-      PYBIND11_OVERRIDE(void, CubeCoefs,	add, coef);
+      PYBIND11_OVERRIDE(void, CubeCoefs,        add, coef);
     }
 
     std::vector<Key> makeKeys(Key k) override {
@@ -552,7 +565,7 @@ void CoefficientClasses(py::module &m) {
 
     std::shared_ptr<CoefStruct> getCoefStruct(double time) override {
       PYBIND11_OVERRIDE(std::shared_ptr<CoefStruct>, TableData, getCoefStruct,
-			time);
+                        time);
     }
 
     std::vector<double> Times() override {
@@ -630,7 +643,7 @@ void CoefficientClasses(py::module &m) {
 
     std::shared_ptr<CoefStruct> getCoefStruct(double time) override {
       PYBIND11_OVERRIDE(std::shared_ptr<CoefStruct>, TrajectoryData, getCoefStruct,
-			time);
+                        time);
     }
 
     std::vector<double> Times() override {
@@ -678,7 +691,7 @@ void CoefficientClasses(py::module &m) {
   py::class_<CoefClasses::CoefStruct, std::shared_ptr<CoefClasses::CoefStruct>, PyCoefStruct>
     (m, "CoefStruct")
     .def(py::init<>(),
-	 R"(
+         R"(
          Base class coefficient data structure object
 
          Returns
@@ -686,7 +699,7 @@ void CoefficientClasses(py::module &m) {
          CoefStruct
          )")
     .def("create", &CoefStruct::create,
-	 R"(
+         R"(
         Initialize a coefficient zeroed structure from user supplied dimensions
 
         Returns
@@ -708,12 +721,12 @@ void CoefficientClasses(py::module &m) {
         while preserving your original coefficients.
         )")
     .def_readonly("geometry", &CoefStruct::geom,
-		  R"(
+                  R"(
                   str
                       geometry type
                   )")
     .def_readonly("time", &CoefStruct::time,
-		   R"(
+                   R"(
                    float
                        data's time stamp
                    )")
@@ -775,8 +788,28 @@ void CoefficientClasses(py::module &m) {
           --------
           setCenter : read-write access to the center data
           )")
+    .def("setGravConstant", &CoefStruct::setGravConstant,
+          py::arg("G"),
+          R"(
+          Set the gravitational constant
+  
+          Parameters
+          ----------
+          G  : float
+             gravitational constant, default is 1.0
+  
+          Returns
+          -------
+          None
+  
+          Notes
+          -----
+          The gravitational constant is used for field evaluation for
+          biorthogonal basis sets.  It will be set automatically when
+          reading EXP coefficient files.
+          )")
     .def("setCoefCenter",
-	 static_cast<void (CoefStruct::*)(Eigen::Vector3d&)>(&CoefStruct::setCenter),
+         static_cast<void (CoefStruct::*)(Eigen::Vector3d&)>(&CoefStruct::setCenter),
           py::arg("mat"),
           R"(
           Set the center vector
@@ -811,7 +844,7 @@ void CoefficientClasses(py::module &m) {
           setCoefRotation : read-write access to the rotation matrix
           )")
     .def("setCoefRotation",
-	 static_cast<void (CoefStruct::*)(RowMatrix3d&)>(&CoefStruct::setRotation),
+         static_cast<void (CoefStruct::*)(RowMatrix3d&)>(&CoefStruct::setRotation),
           py::arg("mat"),
           R"(
           Set the rotation matrix
@@ -845,7 +878,7 @@ void CoefficientClasses(py::module &m) {
         --------
         setCoefs : read-write access to Coefs
         )")
-  .def("setCoefs",		// Member function overload
+  .def("setCoefs",              // Member function overload
         static_cast<void (CoefStruct::*)(Eigen::VectorXcd&)>(&CoefStruct::setCoefs),
         py::arg("mat"),
         R"(
@@ -870,7 +903,7 @@ void CoefficientClasses(py::module &m) {
         --------
         getCoefs : read-only access to Coefs
         )")
-  .def("setCoefs",		// Member function overload
+  .def("setCoefs",              // Member function overload
         static_cast<Eigen::Ref<Eigen::VectorXcd>(CoefStruct::*)()>(&CoefStruct::setCoefs),
         R"(
         Read-write access to the underlying data store
@@ -897,7 +930,7 @@ void CoefficientClasses(py::module &m) {
     (m, "SphStruct")
     .def(py::init<>(), "Spherical coefficient data structure object")
     .def("assign", &SphStruct::assign,
-	      R"(
+              R"(
         Assign a coefficient matrix to CoefStruct.
 
         Parameters
@@ -918,7 +951,7 @@ void CoefficientClasses(py::module &m) {
     (m, "CylStruct")
     .def(py::init<>(), "Cylindrical coefficient data structure object")
     .def("assign", &CylStruct::assign,
-	      R"(
+              R"(
         Assign a coefficient matrix to CoefStruct.
 
         Parameters
@@ -939,7 +972,7 @@ void CoefficientClasses(py::module &m) {
     (m, "SlabStruct")
     .def(py::init<>(), "Slab coefficient data structure object")
     .def("assign", &SlabStruct::assign,
-	      R"(
+              R"(
         Assign a coefficient matrix to CoefStruct.
 
         Parameters
@@ -961,7 +994,7 @@ void CoefficientClasses(py::module &m) {
     (m, "CubeStruct")
     .def(py::init<>(), "Cube coefficient data structure object")
     .def("assign", &CubeStruct::assign,
-	      R"(
+              R"(
         Assign a coefficient matrix to CoefStruct.
 
         Parameters
@@ -983,7 +1016,7 @@ void CoefficientClasses(py::module &m) {
     (m, "TblStruct")
     .def(py::init<>(), "Multicolumn table data structure object")
     .def("assign", &TblStruct::assign,
-	      R"(
+              R"(
         Assign a coefficient matrix to CoefStruct.
 
         Parameters
@@ -1000,7 +1033,7 @@ void CoefficientClasses(py::module &m) {
     (m, "SphFldStruct")
     .def(py::init<>(), "Spherical field coefficient data structure object")
     .def("assign", &SphFldStruct::assign,
-	R"(
+        R"(
         Assign a coefficient matrix to CoefStruct.
 
         Parameters
@@ -1023,7 +1056,7 @@ void CoefficientClasses(py::module &m) {
     (m, "CylFldStruct")
     .def(py::init<>(), "Cylindrical field coefficient data structure object")
     .def("assign", &CylFldStruct::assign,
-	      R"(
+              R"(
         Assign a flattened coefficient array to CylFldStruct.
 
         Parameters
@@ -1068,7 +1101,7 @@ void CoefficientClasses(py::module &m) {
          py::arg("type"),
          py::arg("verbose"))
     .def("__call__",
-	 &CoefClasses::Coefs::getData,
+         &CoefClasses::Coefs::getData,
          R"(
          Return the flattened coefficient structure for the desired time.
 
@@ -1145,18 +1178,116 @@ void CoefficientClasses(py::module &m) {
          You will get a runtime error if the entry does not exist.
          )",py::arg("time"))
     .def("Times",
-            &CoefClasses::Coefs::Times,
-            R"(
-            Return a list of times for coefficient sets currently in the container
+         &CoefClasses::Coefs::Times,
+         R"(
+         Return a list of times for coefficient sets currently in the container
 
-            Returns
-            -------
-            list(float,...)
-                list of times
-            )")
-    .def("WriteH5Coefs",
-            &CoefClasses::Coefs::WriteH5Coefs,
-            R"(
+         Returns
+         -------
+         list(float,...)
+             list of times
+         )")
+    .def("removeUnits",
+         &CoefClasses::Coefs::removeUnits,
+         R"(
+         Remove a unit from the coefficient structure.
+
+         Parameters
+         ----------
+         name : str
+             the name of physical quantity (G, Length, Mass, Time, etc)
+
+         Returns
+         -------
+         None
+         )")
+    .def("setUnits",
+         py::overload_cast<const std::string, const std::string, const float>(&CoefClasses::Coefs::setUnits),
+         R"(
+         Set the units for the coefficient structure.
+
+         Parameters
+         ----------
+         name : str
+            the name of physical quantity (G, Length, Mass, Time, etc)
+         unit : str
+            the unit string (scalar, mixed, kpc, Msun, Myr, km/s etc.).
+            This field is optional and can be empty.
+         value : float
+            the default value of the multiples of the unit
+
+         Returns
+         -------
+         None
+         )", py::arg("name"), py::arg("unit")="", py::arg("value")=1.0)
+    .def("setUnits",
+         py::overload_cast<const std::vector<std::tuple<std::string, std::string, float>>&>(&CoefClasses::Coefs::setUnits),
+         R"(
+         Set the units for the coefficient structure.
+
+         Parameters
+         ----------
+         list((str,str,float))
+             list of (name, unit, value) tuples, where each tuple contains the coefficient name (str), its unit (str), and its value (float).
+
+         Returns
+         -------
+         None
+         )", py::arg("units"))
+    .def("getAllowedUnitNames",
+         &CoefClasses::Coefs::getAllowedUnitNames,
+         R"(
+         Get the allowed unit names for the a given unit type 
+
+         Parameters
+         ----------
+         type : str
+                the type of unit
+
+         Returns
+         -------
+         list(str,...)
+                list of allowed unit names
+
+        )", py::arg("type"))
+    .def("getAllowedUnitTypes",  &CoefClasses::Coefs::getAllowedUnitTypes,
+         R"(
+         Get the allowed unit types for the coefficient structure
+
+         Parameters
+         ----------
+         None
+
+         Returns
+         -------
+         list(str,...)
+                list of allowed unit types
+         )")
+    .def("getAllowedTypeAliases",
+	 &CoefClasses::Coefs::getAllowedTypeAliases,
+	 R"(
+         Get the allowed type aliases for the coefficient structure
+
+         Parameters
+         ----------
+	 type : str
+                the type of unit
+
+         Returns
+         -------
+         list(str)
+		list of allowed type aliases
+	 )", py::arg("type"))
+   .def("WriteH5Coefs",
+         [](CoefClasses::Coefs& self, const std::string& filename) {
+           if (self.getUnits().size()!=4) {
+             std::cout << "Coefs::WriteH5Coefs: please set units for your coefficient set using the `setUnit()` member," << std::endl
+                       << "                     one for each unit.  We suggest explicitly setting 'G', 'Length', 'Mass'," << std::endl
+                       << "                     'Time', or optionally 'Velocity' before writing HDF5 coefficients" << std::endl;
+           }
+           self.WriteH5Coefs(filename);
+         },
+         R"(
             Write the coefficients into an EXP HDF5 coefficient file with the given prefix name.
 
             Parameters
@@ -1174,7 +1305,8 @@ void CoefficientClasses(py::module &m) {
             coefficient file already exists.  This is a safety
             feature.  If you'd like a new version of this file, delete
             the old before this call.
-            )",py::arg("filename"))
+            )",
+         py::arg("filename"))
     .def("ExtendH5Coefs",
             &CoefClasses::Coefs::ExtendH5Coefs,
             R"(
@@ -1298,6 +1430,15 @@ void CoefficientClasses(py::module &m) {
          bool
              True if the data is identical, False otherwise.
          )")
+    .def("getUnits", &CoefClasses::Coefs::getUnits,
+         R"(
+         Get the units of the coefficient data
+
+         Returns
+         -------
+         list((str,str,float))
+             list of (name, unit, value) tuples, where each tuple contains the coefficient name (str), its unit (str), and its value (float).
+         )")
     .def_static("factory", &CoefClasses::Coefs::factory,
               R"(
               Deduce the type and read coefficients from a native or HDF5 file
@@ -1322,7 +1463,7 @@ void CoefficientClasses(py::module &m) {
             py::arg("tmin")=-std::numeric_limits<double>::max(),
             py::arg("tmax")= std::numeric_limits<double>::max())
     .def_static("makecoefs", &CoefClasses::Coefs::makecoefs,
-		R"(
+                R"(
                 make a new coefficient container instance compatible
 
                 Parameters
@@ -1374,7 +1515,7 @@ void CoefficientClasses(py::module &m) {
   py::class_<CoefClasses::SphCoefs, std::shared_ptr<CoefClasses::SphCoefs>, PySphCoefs, CoefClasses::Coefs>
     (m, "SphCoefs", "Container for spherical coefficients")
     .def(py::init<bool>(),
-	 R"(
+         R"(
          Construct a null SphCoefs object
 
          Parameters
@@ -1387,7 +1528,7 @@ void CoefficientClasses(py::module &m) {
          SphCoefs instance
          )")
     .def("__call__",
-	 &CoefClasses::SphCoefs::getMatrix,
+         &CoefClasses::SphCoefs::getMatrix,
          R"(
          Return the coefficient Matrix for the desired time.
 
@@ -1408,7 +1549,7 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"))
     .def("setMatrix",
-	 &CoefClasses::SphCoefs::setMatrix,
+         &CoefClasses::SphCoefs::setMatrix,
          R"(
          Enter and/or rewrite the coefficient matrix at the provided time
 
@@ -1425,13 +1566,13 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"), py::arg("mat"))
     .def("getAllCoefs",
-	 [](CoefClasses::SphCoefs& A)
-	 {
-	   Eigen::Tensor<std::complex<double>, 3> M = A.getAllCoefs(); // Need a copy here
-	   py::array_t<std::complex<double>> ret = make_ndarray3<std::complex<double>>(M);
-	   return ret;
-	 },
-	 R"(
+         [](CoefClasses::SphCoefs& A)
+         {
+           Eigen::Tensor<std::complex<double>, 3> M = A.getAllCoefs(); // Need a copy here
+           py::array_t<std::complex<double>> ret = make_ndarray3<std::complex<double>>(M);
+           return ret;
+         },
+         R"(
         Provide a 3-dimensional ndarray indexed by spherical index, radial index,
         and time index
 
@@ -1449,7 +1590,7 @@ void CoefficientClasses(py::module &m) {
   py::class_<CoefClasses::CylCoefs, std::shared_ptr<CoefClasses::CylCoefs>, PyCylCoefs, CoefClasses::Coefs>
     (m, "CylCoefs", "Container for cylindrical coefficients")
     .def(py::init<bool>(),
-	 R"(
+         R"(
          Construct a null CylCoefs object
 
          Parameters
@@ -1462,7 +1603,7 @@ void CoefficientClasses(py::module &m) {
          CylCoefs instance
          )")
     .def("__call__",
-	 &CoefClasses::CylCoefs::getMatrix,
+         &CoefClasses::CylCoefs::getMatrix,
          R"(
          Return the coefficient Matrix for the desired time.
 
@@ -1483,7 +1624,7 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"))
     .def("setMatrix",
-	 &CoefClasses::CylCoefs::setMatrix,
+         &CoefClasses::CylCoefs::setMatrix,
          R"(
          Enter and/or rewrite the coefficient matrix at the provided time
 
@@ -1500,13 +1641,13 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"), py::arg("mat"))
     .def("getAllCoefs",
-	 [](CoefClasses::CylCoefs& A)
-	 {
-	   Eigen::Tensor<std::complex<double>, 3> M = A.getAllCoefs(); // Need a copy here
-	   py::array_t<std::complex<double>> ret = make_ndarray3<std::complex<double>>(M);
-	   return ret;
-	 },
-	 R"(
+         [](CoefClasses::CylCoefs& A)
+         {
+           Eigen::Tensor<std::complex<double>, 3> M = A.getAllCoefs(); // Need a copy here
+           py::array_t<std::complex<double>> ret = make_ndarray3<std::complex<double>>(M);
+           return ret;
+         },
+         R"(
          Provide a 3-dimensional ndarray indexed by azimuthal index, radial index, and time index
 
          Returns
@@ -1515,11 +1656,11 @@ void CoefficientClasses(py::module &m) {
              3-dimensional numpy array containing the cylindrical coefficients
          )")
     .def("EvenOddPower",
-	 [](CoefClasses::CylCoefs& A, int nodd, int min, int max)
-	 {
-	   return A.EvenOddPower(nodd, min, max);
-	 },
-	 R"(
+         [](CoefClasses::CylCoefs& A, int nodd, int min, int max)
+         {
+           return A.EvenOddPower(nodd, min, max);
+         },
+         R"(
          Get cylindrical coefficient power separated into vertically even and odd contributions.
 
          Parameters
@@ -1544,14 +1685,14 @@ void CoefficientClasses(py::module &m) {
          argument if it is not explicitly set in your EXP::Cylinder
          configuration. If in doubt, use the default.
          )",
-	 py::arg("nodd")=-1, py::arg("min")=0,
-	 py::arg("max")=std::numeric_limits<int>::max());
+         py::arg("nodd")=-1, py::arg("min")=0,
+         py::arg("max")=std::numeric_limits<int>::max());
 
 
   py::class_<CoefClasses::SphFldCoefs, std::shared_ptr<CoefClasses::SphFldCoefs>, CoefClasses::Coefs>
     (m, "SphFldCoefs", "Container for spherical field coefficients")
     .def(py::init<bool>(),
-	       R"(
+               R"(
          Construct a null SphFldCoefs object
 
          Parameters
@@ -1564,12 +1705,12 @@ void CoefficientClasses(py::module &m) {
          SphFldCoefs instance
          )")
     .def("__call__",
-	 [](CoefClasses::SphFldCoefs& A, double time)
-	 {
-	   // Need a copy here
-	   auto M = A.getMatrix(time);
-	   return make_ndarray3<std::complex<double>>(M);
-	 },
+         [](CoefClasses::SphFldCoefs& A, double time)
+         {
+           // Need a copy here
+           auto M = A.getMatrix(time);
+           return make_ndarray3<std::complex<double>>(M);
+         },
          R"(
          Return the coefficient tensor for the desired time.
 
@@ -1590,12 +1731,12 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"))
     .def("setMatrix",
-	 [](CoefClasses::SphFldCoefs& A, double time,
-	    py::array_t<std::complex<double>> mat)
-	 {
-	   auto M = make_tensor3<std::complex<double>>(mat);
-	   A.setMatrix(time, M);
-	 },
+         [](CoefClasses::SphFldCoefs& A, double time,
+            py::array_t<std::complex<double>> mat)
+         {
+           auto M = make_tensor3<std::complex<double>>(mat);
+           A.setMatrix(time, M);
+         },
          R"(
          Enter and/or rewrite the coefficient tensor at the provided time
 
@@ -1612,13 +1753,13 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"), py::arg("mat"))
     .def("getAllCoefs",
-	 [](CoefClasses::SphFldCoefs& A)
-	 {
-	   Eigen::Tensor<std::complex<double>, 4> M = A.getAllCoefs(); // Need a copy here
-	   py::array_t<std::complex<double>> ret = make_ndarray4<std::complex<double>>(M);
-	   return ret;
-	 },
-	 R"(
+         [](CoefClasses::SphFldCoefs& A)
+         {
+           Eigen::Tensor<std::complex<double>, 4> M = A.getAllCoefs(); // Need a copy here
+           py::array_t<std::complex<double>> ret = make_ndarray4<std::complex<double>>(M);
+           return ret;
+         },
+         R"(
         Provide a 4-dimensional ndarray indexed by channel index, spherical index, radial index, and time index
 
         Returns
@@ -1636,7 +1777,7 @@ void CoefficientClasses(py::module &m) {
   py::class_<CoefClasses::CylFldCoefs, std::shared_ptr<CoefClasses::CylFldCoefs>, CoefClasses::Coefs>
     (m, "CylFldCoefs", "Container for cylindrical field coefficients")
     .def(py::init<bool>(),
-	 R"(
+         R"(
          Construct a null CylFldCoefs object
 
          Parameters
@@ -1649,11 +1790,11 @@ void CoefficientClasses(py::module &m) {
          CylFldCoefs instance
          )")
     .def("__call__",
-	 [](CoefClasses::CylFldCoefs& A, double time)
-	 {
-	   auto M = A.getMatrix(time); // Need a copy here
-	   return make_ndarray3<std::complex<double>>(M);
-	 },
+         [](CoefClasses::CylFldCoefs& A, double time)
+         {
+           auto M = A.getMatrix(time); // Need a copy here
+           return make_ndarray3<std::complex<double>>(M);
+         },
          R"(
          Return the coefficient tensor for the desired time.
 
@@ -1674,12 +1815,12 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"))
     .def("setMatrix",
-	 [](CoefClasses::CylFldCoefs& A, double time,
-	    py::array_t<std::complex<double>> mat)
-	 {
-	   auto M = make_tensor3<std::complex<double>>(mat);
-	   A.setMatrix(time, M);
-	 },
+         [](CoefClasses::CylFldCoefs& A, double time,
+            py::array_t<std::complex<double>> mat)
+         {
+           auto M = make_tensor3<std::complex<double>>(mat);
+           A.setMatrix(time, M);
+         },
          R"(
          Enter and/or rewrite the coefficient tensor at the provided time
 
@@ -1696,13 +1837,13 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"), py::arg("mat"))
     .def("getAllCoefs",
-	 [](CoefClasses::CylFldCoefs& A)
-	 {
-	   Eigen::Tensor<std::complex<double>, 4> M = A.getAllCoefs(); // Need a copy here
-	   py::array_t<std::complex<double>> ret = make_ndarray4<std::complex<double>>(M);
-	   return ret;
-	 },
-	 R"(
+         [](CoefClasses::CylFldCoefs& A)
+         {
+           Eigen::Tensor<std::complex<double>, 4> M = A.getAllCoefs(); // Need a copy here
+           py::array_t<std::complex<double>> ret = make_ndarray4<std::complex<double>>(M);
+           return ret;
+         },
+         R"(
         Provide a 4-dimensional ndarray indexed by channel index, spherical index, radial index, and time index
 
         Returns
@@ -1715,7 +1856,7 @@ void CoefficientClasses(py::module &m) {
   py::class_<CoefClasses::SlabCoefs, std::shared_ptr<CoefClasses::SlabCoefs>, PySlabCoefs, CoefClasses::Coefs>
     (m, "SlabCoefs", "Container for cube coefficients")
     .def(py::init<bool>(),
-	       R"(
+               R"(
          Construct a null SlabCoefs object
 
          Parameters
@@ -1728,7 +1869,7 @@ void CoefficientClasses(py::module &m) {
          SlabCoefs instance
          )")
     .def("__call__",
-	 &CoefClasses::SlabCoefs::getTensor,
+         &CoefClasses::SlabCoefs::getTensor,
          R"(
          Return the coefficient tensor for the desired time.
 
@@ -1749,7 +1890,7 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"))
     .def("setTensor",
-	 &CoefClasses::SlabCoefs::setTensor,
+         &CoefClasses::SlabCoefs::setTensor,
          R"(
          Enter and/or rewrite the coefficient tensor at the provided time
 
@@ -1766,13 +1907,13 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"), py::arg("tensor"))
     .def("getAllCoefs",
-	 [](CoefClasses::SlabCoefs& A)
-	 {
-	   Eigen::Tensor<std::complex<double>, 4> M = A.getAllCoefs(); // Need a copy here
-	   py::array_t<std::complex<double>> ret = make_ndarray4<std::complex<double>>(M);
-	   return ret;
-	 },
-	 R"(
+         [](CoefClasses::SlabCoefs& A)
+         {
+           Eigen::Tensor<std::complex<double>, 4> M = A.getAllCoefs(); // Need a copy here
+           py::array_t<std::complex<double>> ret = make_ndarray4<std::complex<double>>(M);
+           return ret;
+         },
+         R"(
          Provide a 4-dimensional ndarray indexed by nx, ny, nz, and time indices.
 
          Returns
@@ -1781,11 +1922,11 @@ void CoefficientClasses(py::module &m) {
              4-dimensional numpy array containing the slab coefficients
          )")
     .def("PowerDim",
-	 [](CoefClasses::SlabCoefs& A, std::string d, int min, int max)
-	 {
-	   return A.Power(d[0], min, max);
-	 },
-	 R"(
+         [](CoefClasses::SlabCoefs& A, std::string d, int min, int max)
+         {
+           return A.Power(d[0], min, max);
+         },
+         R"(
          Get power for the coefficient DB as a function of harmonic index for a
          given dimension.  This Power() member is equivalent to PowerDim('x').
 
@@ -1803,13 +1944,13 @@ void CoefficientClasses(py::module &m) {
          numpy.ndarray
              2-dimensional numpy array containing the power
          )", py::arg("d"), py::arg("min")=0,
-	 py::arg("max")=std::numeric_limits<int>::max());
+         py::arg("max")=std::numeric_limits<int>::max());
 
 
   py::class_<CoefClasses::CubeCoefs, std::shared_ptr<CoefClasses::CubeCoefs>, PyCubeCoefs, CoefClasses::Coefs>
     (m, "CubeCoefs", "Container for cube coefficients")
     .def(py::init<bool>(),
-	       R"(
+               R"(
          Construct a null CubeCoefs object
 
          Parameters
@@ -1822,7 +1963,7 @@ void CoefficientClasses(py::module &m) {
          CubeCoefs instance
          )")
     .def("__call__",
-	 &CoefClasses::CubeCoefs::getTensor,
+         &CoefClasses::CubeCoefs::getTensor,
          R"(
          Return the coefficient tensor for the desired time.
 
@@ -1843,7 +1984,7 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"))
     .def("setTensor",
-	 &CoefClasses::CubeCoefs::setTensor,
+         &CoefClasses::CubeCoefs::setTensor,
          R"(
          Enter and/or rewrite the coefficient tensor at the provided time
 
@@ -1860,13 +2001,13 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"), py::arg("tensor"))
     .def("getAllCoefs",
-	 [](CoefClasses::CubeCoefs& A)
-	 {
-	   Eigen::Tensor<std::complex<double>, 4> M = A.getAllCoefs(); // Need a copy here
-	   py::array_t<std::complex<double>> ret = make_ndarray4<std::complex<double>>(M);
-	   return ret;
-	 },
-	 R"(
+         [](CoefClasses::CubeCoefs& A)
+         {
+           Eigen::Tensor<std::complex<double>, 4> M = A.getAllCoefs(); // Need a copy here
+           py::array_t<std::complex<double>> ret = make_ndarray4<std::complex<double>>(M);
+           return ret;
+         },
+         R"(
          Provide a 4-dimensional ndarray indexed by nx, ny, nz, and time indices.
 
          Returns
@@ -1875,11 +2016,11 @@ void CoefficientClasses(py::module &m) {
              4-dimensional numpy array containing the cube coefficients
          )")
     .def("PowerDim",
-	 [](CoefClasses::CubeCoefs& A, std::string d, int min, int max)
-	 {
-	   return A.Power(d[0], min, max);
-	 },
-	 R"(
+         [](CoefClasses::CubeCoefs& A, std::string d, int min, int max)
+         {
+           return A.Power(d[0], min, max);
+         },
+         R"(
          Get power for the coefficient DB as a function of harmonic index for a
          given dimension.  This Power() member is equivalent to PowerDim('x').
 
@@ -1897,12 +2038,12 @@ void CoefficientClasses(py::module &m) {
          numpy.ndarray
              2-dimensional numpy array containing the power
          )", py::arg("d"), py::arg("min")=0,
-	 py::arg("max")=std::numeric_limits<int>::max());
+         py::arg("max")=std::numeric_limits<int>::max());
 
   py::class_<CoefClasses::TableData, std::shared_ptr<CoefClasses::TableData>, PyTableData, CoefClasses::Coefs>
     (m, "TableData", "Container for simple data tables with multiple columns")
     .def(py::init<bool>(),
-	 R"(
+         R"(
          Construct a null TableData object
 
          Parameters
@@ -1915,7 +2056,7 @@ void CoefficientClasses(py::module &m) {
          TableData instance
          )", py::arg("verbose")=true)
     .def(py::init<std::string&>(),
-	 R"(
+         R"(
          Construct a TableData object from a data file
 
          Parameters
@@ -1928,7 +2069,7 @@ void CoefficientClasses(py::module &m) {
          TableData instance
          )")
     .def(py::init<std::string&, bool>(),
-	 R"(
+         R"(
          Construct a TableData object from a data file
 
          Parameters
@@ -1943,7 +2084,7 @@ void CoefficientClasses(py::module &m) {
          TableData instance
          )", py::arg("filename"), py::arg("verbose")=true)
     .def(py::init<std::vector<double>&, std::vector<std::vector<double>>&, bool>(),
-	 R"(
+         R"(
          Construct a TableData object from data arrays
 
          Parameters
@@ -1960,7 +2101,7 @@ void CoefficientClasses(py::module &m) {
          TableData instance
          )", py::arg("time"), py::arg("array"), py::arg("verbose")=true)
     .def("getAllCoefs", &CoefClasses::TableData::getAllCoefs,
-	 R"(
+         R"(
          Return a 2-dimensional ndarray indexed by column and time
 
          Returns
@@ -1972,7 +2113,7 @@ void CoefficientClasses(py::module &m) {
   py::class_<CoefClasses::TrajectoryData, std::shared_ptr<CoefClasses::TrajectoryData>, PyTrajectoryData, CoefClasses::Coefs>
     (m, "TrajectoryData", "Container for trajectory/orbit data")
     .def(py::init<bool>(),
-	 R"(
+         R"(
          Construct a null TrajectoryData object
 
          Parameters
@@ -1985,7 +2126,7 @@ void CoefficientClasses(py::module &m) {
          TrajectoryData instance
          )", py::arg("verbose")=true)
     .def(py::init<std::string&>(),
-	 R"(
+         R"(
          Construct a TrajectoryData object from a data file
 
          Parameters
@@ -1998,7 +2139,7 @@ void CoefficientClasses(py::module &m) {
          TrajectoryData instance
          )")
     .def(py::init<std::string&, bool>(),
-	 R"(
+         R"(
          Construct a TrajectoryData object from a data file
 
          Parameters
@@ -2013,7 +2154,7 @@ void CoefficientClasses(py::module &m) {
          TrajectoryData instance
          )", py::arg("filename"), py::arg("verbose")=true)
     .def(py::init<std::vector<double>&, std::vector<Eigen::MatrixXd>&, bool>(),
-	 R"(
+         R"(
          Construct a TrajectoryData object from data arrays
 
          Parameters
@@ -2030,14 +2171,14 @@ void CoefficientClasses(py::module &m) {
          TrajectoryData instance
          )", py::arg("time"), py::arg("array"), py::arg("verbose")=true)
     .def("getAllCoefs",
-	 [](CoefClasses::TrajectoryData& A)
-	 {
-	   Eigen::Tensor<double, 3> M = A.getAllCoefs(); // Need a copy here
-	   py::array_t<double> ret = make_ndarray3<double>(M);
-	   return ret;
-	 },
+         [](CoefClasses::TrajectoryData& A)
+         {
+           Eigen::Tensor<double, 3> M = A.getAllCoefs(); // Need a copy here
+           py::array_t<double> ret = make_ndarray3<double>(M);
+           return ret;
+         },
 
-	 R"(
+         R"(
          Return a 3-dimensional ndarray indexed by column and time
 
          Returns
diff --git a/src/Cube.cc b/src/Cube.cc
index bbf59fb7f..d1a2e7fad 100644
--- a/src/Cube.cc
+++ b/src/Cube.cc
@@ -520,6 +520,11 @@ void Cube::dump_coefs_h5(const std::string& file)
     // Add the name attribute.  We only need this on the first call.
     cubeCoefs.setName(component->name);
 
+    // Add the default units
+    cubeCoefs.setUnits({{"length", "none", 1.0},
+			{"mass",   "none", 1.0},
+			{"time",   "none", 1.0}});
+
     // And the new coefficients and write the new HDF5
     cubeCoefs.clear();
     cubeCoefs.add(cur);
diff --git a/src/Cylinder.cc b/src/Cylinder.cc
index 67685e13a..b33e5132e 100644
--- a/src/Cylinder.cc
+++ b/src/Cylinder.cc
@@ -1587,6 +1587,11 @@ void Cylinder::dump_coefs_h5(const std::string& file)
     // Add the name attribute.  We only need this on the first call.
     cylCoefs.setName(component->name);
 
+    // Add the default units
+    cylCoefs.setUnits({{"length", "none", 1.0},
+		       {"mass",   "none", 1.0},
+		       {"time",   "none", 1.0}});
+
     // Add the new coefficients and write the new HDF5
     cylCoefs.clear();
     cylCoefs.add(cur);
diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc
index 976169325..a8749a258 100644
--- a/src/PolarBasis.cc
+++ b/src/PolarBasis.cc
@@ -1842,6 +1842,11 @@ void PolarBasis::dump_coefs_h5(const std::string& file)
     // Add the name attribute.  We only need this on the first call.
     cylCoefs.setName(component->name);
 
+    // Add the default units
+    cylCoefs.setUnits({{"length", "none", 1.0},
+		       {"mass",   "none", 1.0},
+		       {"time",   "none", 1.0}});
+
     // And the new coefficients and write the new HDF5
     cylCoefs.clear();
     cylCoefs.add(cur);
diff --git a/src/ShearSL.cc b/src/ShearSL.cc
index 12bfe70f0..038df984e 100644
--- a/src/ShearSL.cc
+++ b/src/ShearSL.cc
@@ -374,6 +374,11 @@ void ShearSL::dump_coefs_h5(const std::string& file)
     // Add the name attribute.  We only need this on the first call.
     slabCoefs.setName(component->name);
 
+    // Add the default units
+    slabCoefs.setUnits({{"length", "none", 1.0},
+			{"mass",   "none", 1.0},
+			{"time",   "none", 1.0}});
+
     // And the new coefficients and write the new HDF5
     slabCoefs.clear();
     slabCoefs.add(cur);
diff --git a/src/SlabSL.cc b/src/SlabSL.cc
index f3d4549d2..702aa481d 100644
--- a/src/SlabSL.cc
+++ b/src/SlabSL.cc
@@ -513,6 +513,11 @@ void SlabSL::dump_coefs_h5(const std::string& file)
     // Add the name attribute.  We only need this on the first call.
     slabCoefs.setName(component->name);
 
+    // Add the default units
+    slabCoefs.setUnits({{"length", "none", 1.0},
+			{"mass",   "none", 1.0},
+			{"time",   "none", 1.0}});
+
     // And the new coefficients and write the new HDF5
     slabCoefs.clear();
     slabCoefs.add(cur);
diff --git a/src/SphericalBasis.cc b/src/SphericalBasis.cc
index a520ff6d0..406e748ca 100644
--- a/src/SphericalBasis.cc
+++ b/src/SphericalBasis.cc
@@ -1926,6 +1926,11 @@ void SphericalBasis::dump_coefs_h5(const std::string& file)
     // Add the name attribute.  We only need this on the first call.
     sphCoefs.setName(component->name);
 
+    // Add the default units
+    sphCoefs.setUnits({{"length", "none", 1.0},
+		       {"mass",   "none", 1.0},
+		       {"time",   "none", 1.0}});
+
     // And the new coefficients and write the new HDF5
     sphCoefs.clear();
     sphCoefs.add(cur);
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 560f9e8b5..30dedb2d4 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -56,105 +56,107 @@ if(ENABLE_NBODY)
   add_test(NAME expExecuteTest
     COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/src/exp -v)
 
-  # Makes some spherical ICs using utils/ICs/gensph
-  add_test(NAME makeICTest
-    COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph -N 10000 -i SLGridSph.model
-    WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo)
-
-  # Runs those ICs using exp
-  add_test(NAME expNbodyTest
-    COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/src/exp config.yml
-    WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo)
-
-  set_tests_properties(expNbodyTest PROPERTIES DEPENDS makeICTest)
-
-  # Check OUTLOG file for a sane 2T/W mean and stdv
-  add_test(NAME expNbodyCheck2TW
-    COMMAND ${CMAKE_COMMAND} -E env
-    PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
-    ${PYTHON_EXECUTABLE} check.py
-    WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo)
-
-  set_tests_properties(expNbodyCheck2TW PROPERTIES DEPENDS expNbodyTest LABELS "long")
-
-  # This adds a coefficient read test using pyEXP only if
-  # expNbodyTest is run and pyEXP has been built
-  if(ENABLE_PYEXP)
-     # Read coefficient file with pyEXP
-     add_test(NAME pyEXPCoefReadTest
-              COMMAND ${CMAKE_COMMAND} -E env
-              PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
-              ${PYTHON_EXECUTABLE} readCoefs.py
-              WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefReadTest PROPERTIES DEPENDS expNbodyTest LABELS "long")
-
-     add_test(NAME pyEXPCoefMatrixTest
-              COMMAND ${CMAKE_COMMAND} -E env
-              PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
-              ${PYTHON_EXECUTABLE} changeCoefs.py
-              WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefMatrixTest PROPERTIES DEPENDS expNbodyTest LABELS "long")
-
-     add_test(NAME pyEXPCoefCreateTest
-              COMMAND ${CMAKE_COMMAND} -E env
-              PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
-              ${PYTHON_EXECUTABLE} createCoefs.py
-              WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quick")
+  if (ENABLE_UTILS)
+    # Makes some spherical ICs using utils/ICs/gensph
+    add_test(NAME makeICTest
+      COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph -N 10000 -i SLGridSph.model
+      WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo)
+
+    # Runs those ICs using exp
+    add_test(NAME expNbodyTest
+      COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/src/exp config.yml
+      WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo)
+
+    set_tests_properties(expNbodyTest PROPERTIES DEPENDS makeICTest)
+
+    # Check OUTLOG file for a sane 2T/W mean and stdv
+    add_test(NAME expNbodyCheck2TW
+      COMMAND ${CMAKE_COMMAND} -E env
+      PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
+      ${PYTHON_EXECUTABLE} check.py
+      WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo)
+
+    set_tests_properties(expNbodyCheck2TW PROPERTIES DEPENDS expNbodyTest LABELS "long")
+
+    # This adds a coefficient read test using pyEXP only if
+    # expNbodyTest is run and pyEXP has been built
+    if(ENABLE_PYEXP)
+       # Read coefficient file with pyEXP
+       add_test(NAME pyEXPCoefReadTest
+                COMMAND ${CMAKE_COMMAND} -E env
+                PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
+                ${PYTHON_EXECUTABLE} readCoefs.py
+                WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
+       set_tests_properties(pyEXPCoefReadTest PROPERTIES DEPENDS expNbodyTest LABELS "long")
+
+       add_test(NAME pyEXPCoefMatrixTest
+                COMMAND ${CMAKE_COMMAND} -E env
+                PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
+                ${PYTHON_EXECUTABLE} changeCoefs.py
+                WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
+       set_tests_properties(pyEXPCoefMatrixTest PROPERTIES DEPENDS expNbodyTest LABELS "long")
+
+       add_test(NAME pyEXPCoefCreateTest
+                COMMAND ${CMAKE_COMMAND} -E env
+                PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
+                ${PYTHON_EXECUTABLE} createCoefs.py
+                WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
+       set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quick")
+    endif()
+
+    # A separate test to remove the generated files if they all exist;
+    # perhaps there is a better way?
+    add_test(NAME removeTempFiles
+      COMMAND ${CMAKE_COMMAND} -E remove
+      config.run0.yml current.processor.rates.run0 new.bods
+      OUTLOG.run0 run0.levels SLGridSph.cache.run0 test.grid
+      outcoef.halo.run0 SLGridSph.cache.run0
+      WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo)
+
+    # Remove the temporary files
+    set_tests_properties(removeTempFiles PROPERTIES DEPENDS expNbodyCheck2TW
+      REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;run0.levels;SLGridSph.cache.run0;test.grid;"
+      )
+
+    # Makes some cube ICs using utils/ICs/cubeics
+    add_test(NAME makeCubeICTest
+      COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/cubeics -N 4000 -z -d 2,2,2
+      WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube)
+
+    # Runs those ICs using exp
+    add_test(NAME expCubeTest
+      COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/src/exp
+      WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube)
+
+    set_tests_properties(expCubeTest PROPERTIES DEPENDS makeCubeICTest)
+
+    # Check OUTLOG file for mean position
+    add_test(NAME expCubeCheckPos
+      COMMAND ${CMAKE_COMMAND} -E env
+      PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
+      ${PYTHON_EXECUTABLE} check.py
+      WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube)
+
+    set_tests_properties(expCubeCheckPos PROPERTIES DEPENDS expCubeTest)
+
+    # A separate test to remove the generated files if they all exist
+    add_test(NAME removeCubeFiles
+      COMMAND ${CMAKE_COMMAND} -E remove
+      config.runS.yml current.processor.rates.runS cube.bods
+      OUTLOG.runS runS.levels
+      WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube)
+
+    # Remove the temporary files
+    set_tests_properties(removeCubeFiles PROPERTIES DEPENDS expCubeCheckPos
+      REQUIRED_FILES "config.runS.yml;current.processor.rates.runS;cube.bods;OUTLOG.runS;runS.levels;")
+
+    # Set labels for pyEXP tests
+    set_tests_properties(expExecuteTest PROPERTIES LABELS "quick")
+    set_tests_properties(makeICTest expNbodyTest expNbodyCheck2TW
+    removeTempFiles makeCubeICTest expCubeTest removeCubeFiles
+    PROPERTIES LABELS "long")
   endif()
 
-  # A separate test to remove the generated files if they all exist;
-  # perhaps there is a better way?
-  add_test(NAME removeTempFiles
-    COMMAND ${CMAKE_COMMAND} -E remove
-    config.run0.yml current.processor.rates.run0 new.bods
-    OUTLOG.run0 run0.levels SLGridSph.cache.run0 test.grid
-    outcoef.halo.run0 SLGridSph.cache.run0
-    WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo)
-
-  # Remove the temporary files
-  set_tests_properties(removeTempFiles PROPERTIES DEPENDS expNbodyCheck2TW
-    REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;run0.levels;SLGridSph.cache.run0;test.grid;"
-    )
-
-  # Makes some cube ICs using utils/ICs/cubeics
-  add_test(NAME makeCubeICTest
-    COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/cubeics -N 4000 -z -d 2,2,2
-    WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube)
-
-  # Runs those ICs using exp
-  add_test(NAME expCubeTest
-    COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/src/exp
-    WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube)
-
-  set_tests_properties(expCubeTest PROPERTIES DEPENDS makeCubeICTest)
-
-  # Check OUTLOG file for mean position
-  add_test(NAME expCubeCheckPos
-    COMMAND ${CMAKE_COMMAND} -E env
-    PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
-    ${PYTHON_EXECUTABLE} check.py
-    WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube)
-
-  set_tests_properties(expCubeCheckPos PROPERTIES DEPENDS expCubeTest)
-
-  # A separate test to remove the generated files if they all exist
-  add_test(NAME removeCubeFiles
-    COMMAND ${CMAKE_COMMAND} -E remove
-    config.runS.yml current.processor.rates.runS cube.bods
-    OUTLOG.runS runS.levels
-    WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube)
-
-  # Remove the temporary files
-  set_tests_properties(removeCubeFiles PROPERTIES DEPENDS expCubeCheckPos
-    REQUIRED_FILES "config.runS.yml;current.processor.rates.runS;cube.bods;OUTLOG.runS;runS.levels;")
-
-  # Set labels for pyEXP tests
-  set_tests_properties(expExecuteTest PROPERTIES LABELS "quick")
-  set_tests_properties(makeICTest expNbodyTest expNbodyCheck2TW
-  removeTempFiles makeCubeICTest expCubeTest removeCubeFiles
-  PROPERTIES LABELS "long")
-
 endif()
 
 # Tests that should work for any configuration.  Nothing so far.  This