Index: kernel/include/IMP/internal/Unit.h
===================================================================
--- kernel/include/IMP/internal/Unit.h	(revision 0)
+++ kernel/include/IMP/internal/Unit.h	(revision 0)
@@ -0,0 +1,408 @@
+/**
+ *  \file Unit.h     \brief Classes to help with converting between units.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#ifndef __IMP_UNIT_H
+#define __IMP_UNIT_H
+
+
+#include "ExponentialNumber.h"
+#include "../utility.h"
+
+#include <boost/mpl/vector_c.hpp>
+#include <boost/mpl/clear.hpp>
+#include <boost/mpl/at.hpp>
+#include <boost/mpl/transform.hpp>
+#include <boost/mpl/equal.hpp>
+#include <boost/mpl/placeholders.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/mpl/divides.hpp>
+#include <boost/mpl/multiplies.hpp>
+#include <boost/mpl/plus.hpp>
+#include <boost/mpl/minus.hpp>
+#include <boost/mpl/front.hpp>
+#include <boost/mpl/pop_front.hpp>
+#include <boost/mpl/size.hpp>
+
+#include <iomanip>
+#include <cmath>
+
+
+namespace IMP
+{
+
+namespace internal
+{
+
+
+namespace unit
+{
+
+
+namespace internal
+{
+
+template <class Tag>
+std::string get_unit_name(int o) {
+  bool should_not_be_compiled;
+  return "";
+}
+
+// for stupid line length limits
+using namespace boost::mpl;
+using namespace boost::mpl::placeholders;
+
+
+/*
+  These classes allow the Units part of a unit to be manipulated
+ */
+template <class UA, class UB>
+struct Divide {
+  typedef typename transform<UA, UB, minus<_1, _2> >::type type;
+};
+
+template <class UA, class UB>
+struct Multiply {
+  typedef typename transform<UA, UB, plus<_1, _2> >::type type;
+};
+
+template <class UA>
+struct Sqrt { 
+  typedef typename transform<UA, divides<_1, int_<2> > >::type type;
+};
+
+
+template <class UA, class UB>
+struct Exponentiate {
+  typedef typename transform<UA, UB, multiplies<_1, _2> >::type type;
+};
+
+template <class InputVector, int SZ>
+struct DoNormalize {
+  typedef int type;
+};
+
+/*
+  The boost MPI multiply sorts of routines return a slightly different
+  type that we want so we have to change it into the exact type to get
+  certain things work work right.
+ */
+template <class InputVector>
+struct DoNormalize<InputVector, 0> {
+  typedef vector_c<int> type;
+};
+
+template <class InputVector>
+struct DoNormalize<InputVector, 1> {
+  typedef vector_c<int, at<InputVector, int_<0> >::type::value> type;
+};
+
+template <class InputVector>
+struct DoNormalize<InputVector, 2> {
+  typedef vector_c<int,
+                   at<InputVector, int_<0> >::type::value,
+                   at<InputVector, int_<1> >::type::value> type;
+};
+
+template <class InputVector>
+struct DoNormalize<InputVector, 3> {
+  typedef vector_c<int,
+                   at<InputVector, int_<0> >::type::value,
+                   at<InputVector, int_<1> >::type::value,
+                   at<InputVector, int_<2> >::type::value> type;
+};
+
+template <class InputVector>
+struct DoNormalize<InputVector, 4> {
+  typedef vector_c<int,
+                   at<InputVector, int_<0> >::type::value,
+                   at<InputVector, int_<1> >::type::value,
+                   at<InputVector, int_<2> >::type::value,
+                   at<InputVector, int_<3> >::type::value> type;
+};
+
+template <class InputVector>
+struct DoNormalize<InputVector, 5> {
+  typedef vector_c<int,
+                   at<InputVector, int_<0> >::type::value,
+                   at<InputVector, int_<1> >::type::value,
+                   at<InputVector, int_<2> >::type::value,
+                   at<InputVector, int_<3> >::type::value,
+                   at<InputVector, int_<4> >::type::value> type;
+};
+
+template <class InputVector>
+struct Normalize {
+  typedef typename DoNormalize<InputVector,
+                               size<InputVector>::type::value > ::type type;
+};
+
+// Recursive type to print out the unit names and powers
+template <class Tag, int O, int SZ, class Units>
+struct PrintUnits {
+  PrintUnits<Tag, O+1, SZ, Units> p_;
+
+  void operator()(std::ostream &out) const {
+    std::string str= get_unit_name<Tag>(O);
+    int e=  boost::mpl::at<Units, typename boost::mpl::int_<O> >::type::value;
+    if (e != 0)  {
+      out << " " << str;
+      if (e != 1) {
+        out << "^" << e;
+      }
+    }
+    p_(out);
+  }
+};
+
+//! Terminate the recursion
+template <class Tag, int O, class Units>
+struct PrintUnits<Tag, O, O, Units > {
+  void operator()(std::ostream &out) const{};
+};
+
+
+//! Specializaton for singleton units
+template <class Tag, int O>
+struct PrintUnits<Tag, O, O, boost::mpl::vector_c<int> > {
+  void operator()(std::ostream &out) const{
+    std::string str= get_unit_name<Tag>(O);
+    out << " " << str;
+  }
+};
+
+// recursive type to check if all the unit powers are zero
+template <int O, int SZ, class Units>
+struct IsNoUnits {
+  typedef IsNoUnits<O+1, SZ, Units> Next;
+  static const bool value= Next::value 
+    && !(boost::mpl::at<Units, typename boost::mpl::int_<O> >::type::value);
+};
+
+
+template <int O, class Units>
+struct IsNoUnits<O, O, Units > {
+  static const bool value =true;
+};
+
+
+
+} // unit::internal
+
+
+
+
+//! A base class for units
+/**
+   \internal
+ */
+template <class TagT, int EXPT, class UnitsT>
+class Unit {
+  template <class T, int E, class U>
+  friend class Unit;
+  typedef Unit<TagT, EXPT, UnitsT> This;
+  typedef ExponentialNumber<EXPT> V;
+  V v_;
+
+
+  bool is_default() const {
+    return v_== V();
+  }
+public:
+  typedef TagT Tag;
+  typedef UnitsT Units;
+  static const int EXP= EXPT;
+
+  explicit Unit(V v): v_(v){}
+  explicit Unit(double v): v_(v){}
+  //Unit(int v): v_(v){}
+  Unit(){}
+  template <int OEXP, class OUnits>
+  Unit(Unit<Tag, OEXP, OUnits> o): v_(o.v_) {
+    BOOST_STATIC_ASSERT((
+                         boost::mpl::equal<Units,OUnits>::type::value
+    ));
+  }
+
+  template <int OEXP>
+  Unit(Unit<Tag, OEXP, Units> o): v_(o.v_) {}
+
+  double to_scalar() const {
+    BOOST_STATIC_ASSERT((internal::IsNoUnits<0,
+                         boost::mpl::size<Units>::type::value, Units>::value));
+    return get_normalized_value();
+  }
+
+  IMP_COMPARISONS_1(v_);
+
+  //! Get the current value
+  /** \note The value returned is the value before it is multiplied
+      by the appropriate power of 10. That means, 1 Angstrom returns 1,
+      as does 1 Meter.
+   */
+  double get_value() const {
+    return v_.get_value();
+  }
+  //! Get the value with the proper exponent
+  /** 1 Angstrom returns 1e-10.
+   */
+  double get_normalized_value() const {
+    return v_.get_normalized_value();
+  }
+
+  V get_exponential_value() const {
+    return v_;
+  }
+
+
+  template <int OEXP, class OUnits>
+  Unit<Tag, OEXP+EXP, 
+       typename internal::Multiply<Units, OUnits>::type >
+  operator*(Unit<Tag, OEXP, OUnits> o) const {
+    return Unit<Tag, OEXP+EXP, 
+      typename internal::Multiply<Units, OUnits>::type >(v_*o.v_);
+  }
+
+  template <int OEXP, class OUnits>
+  Unit<Tag, EXP-OEXP, 
+       typename internal::Divide<Units, OUnits>::type >
+  operator/(Unit<Tag, OEXP, OUnits> o) const {
+    return Unit<Tag, EXP-OEXP,
+      typename internal::Divide<Units, OUnits>::type > (v_/o.v_);
+  }
+
+  This operator+(This o) const {
+    return This(v_+o.v_);
+  }
+  This operator-(This o) const {
+    return This(v_-o.v_);
+  }
+
+  This operator-() const {
+    return This(-v_);
+  }
+
+  void show(std::ostream &out) const {
+    out << v_;
+    internal::PrintUnits<Tag, 0,
+      boost::mpl::size<Units>::type::value, Units> p;
+    p(out);
+  }
+};
+
+template <class Tag, int EXP, class Units>
+std::ostream &operator<<(std::ostream &out,
+                         Unit<Tag, EXP, Units> o) {
+  o.show(out);
+  return out;
+}
+
+
+typedef boost::mpl::vector_c<int> SingletonUnit;
+
+
+
+/** need to be careful of integer division
+    \internal
+ */
+template <class Tag, int EXP, class Units>
+Unit<Tag, EXP/2, typename internal::Sqrt<Units>::type >
+sqrt(Unit<Tag, EXP, Units> o) {
+  return Unit<Tag, EXP/2, typename internal::Sqrt<Units>::type >
+    (ExponentialNumber<EXP/2>(std::sqrt(o.get_normalized_value())));
+}
+
+/** \internal
+ */
+template <class Tag, int EXP, class Units>
+Unit<Tag, EXP*2, typename internal::Multiply<Units, Units>::type >
+square(Unit<Tag, EXP, Units> o) {
+  return Unit<Tag, EXP*2, typename internal::Multiply<Units, Units>::type >
+    (::IMP::square(o.get_value()));
+}
+
+template <class Tag, int EXP, class Units>
+Unit<Tag, EXP, Units >
+operator*(Unit<Tag, EXP, Units> o, double d) {
+  return Unit<Tag, EXP, Units>(d*o.get_value());
+}
+
+template <class Tag, int EXP, class Units>
+Unit<Tag, EXP, Units >
+operator*(double d, Unit<Tag, EXP, Units> o) {
+  return Unit<Tag, EXP, Units>(d*o.get_value());
+}
+
+template <class Tag, int EXP, class Units>
+Unit<Tag, EXP, Units >
+operator/(Unit<Tag, EXP, Units> o, double d) {
+  return Unit<Tag, EXP, Units>(o.get_value()/d);
+}
+
+template <class Tag, int EXP, class Units>
+Unit<Tag, EXP, Units >
+operator/(double d, Unit<Tag, EXP, Units> o) {
+  return Unit<Tag, EXP, Units>(d/o.get_value());
+}
+
+// Multiply and divide Unit instantiations
+
+template <class U0, class U1>
+struct Divide {
+  BOOST_STATIC_ASSERT((boost::mpl::equal<typename U0::Tag,
+                       typename U1::Tag>::type::value));
+  typedef typename internal::Divide<typename U0::Units,
+                                    typename U1::Units>::type raw_units;
+  typedef typename internal::Normalize<raw_units>::type units;
+  typedef Unit<typename U0::Tag, U0::EXP - U1::EXP, units > type;
+
+};
+
+template <class U0, class U1>
+struct Multiply {
+  BOOST_STATIC_ASSERT((boost::mpl::equal<typename U0::Tag,
+                       typename U1::Tag>::type::value));
+  typedef typename internal::Multiply<typename U0::Units,
+                                      typename U1::Units>::type raw_units;
+  typedef typename internal::Normalize<raw_units>::type units;       
+  typedef Unit<typename U0::Tag, U0::EXP + U1::EXP,
+                         units > type;
+};
+
+template <class U, int E>
+struct Shift {
+  typedef Unit<typename U::Tag, U::EXP + E,
+               typename U::Units > type;
+
+};
+
+
+
+template <class U, class R, class A, int DEXP>
+struct Exchange {
+  BOOST_STATIC_ASSERT((boost::mpl::equal<typename U::Tag,
+                      typename R::Tag>::type::value));
+  BOOST_STATIC_ASSERT((boost::mpl::equal<typename U::Tag,
+                      typename A::Tag>::type::value));
+  typedef typename internal::Divide<typename U::Units,
+                                    typename R::Units>::type Div;
+  typedef typename internal::Multiply<Div, typename A::Units>::type Mul;
+  typedef typename internal::Normalize<Mul>::type units;       
+  typedef Unit<typename U::Tag,
+               U::EXP - R::EXP + A::EXP  +DEXP,
+               units> type;
+
+};
+
+} // unit
+
+} // internal
+
+
+} // IMP
+
+
+#endif
Index: kernel/include/IMP/internal/units.h
===================================================================
--- kernel/include/IMP/internal/units.h	(revision 589)
+++ kernel/include/IMP/internal/units.h	(working copy)
@@ -8,11 +8,12 @@
 #ifndef __IMP_UNITS_H
 #define __IMP_UNITS_H
 
-#include "../utility.h"
+#include "Unit.h"
+#include "../base_types.h"
 
+
 #include <sstream>
 #include <cmath>
-#include <iomanip>
 #include <limits>
 
 namespace IMP
@@ -21,478 +22,211 @@
 namespace internal
 {
 
-/** This helper class implements a floating point number
-    with an extra exponent. The actual value is the number 
-    times 10^EXP. I needed this since I was running out of bits
-    in the exponent when converting between vastly different units.
-    It also means that you can have a type for both nanometer
-    and meter without much pain. Since the exponent is stored
-    as a template argument, there shouldn't be any runtime overhead.
- */
-template <int EXP>
-class ExponentialNumber
-{
-  typedef ExponentialNumber<EXP> This;
+// want it for conversions
+extern IMPDLLEXPORT const unit::ExponentialNumber<3> JOULES_PER_KILOCALORIE;
+extern IMPDLLEXPORT const unit::ExponentialNumber<23> NA;
 
-  template <int OEXP>
-  void copy_from(ExponentialNumber<OEXP> o) {
-    const int diff = OEXP-EXP;
-    v_=o.v_;
-    float factor= std::pow(10.0, static_cast<double>(diff));
-    v_*=factor;
-    IMP_assert(std::abs(get_normalized_value() -o.get_normalized_value())
-               <= .09 * std::abs(std::max(get_normalized_value(),
-                                          o.get_normalized_value())),
-               "Mismatch after scaling");
-  }
-  bool is_default() const {
-    return v_== std::numeric_limits<double>::infinity();
-  }
-  // should be private
-  double v_;
-  template <int E>
-  friend class ExponentialNumber;
-public:
 
+namespace unit
+{
 
-  ExponentialNumber(): v_(std::numeric_limits<double>::infinity()) {}
-  template <int OEXP>
-  ExponentialNumber(ExponentialNumber<OEXP> o) {
-    copy_from(o);
-  }
-  ExponentialNumber(double d): v_(d){
-  }
-  template <int OEXP>
-  This &operator=(ExponentialNumber<OEXP> o) {
-    copy_from(o);
-    return *this;
-  }
+namespace internal
+{
 
-  //! Compute the value with the exponent
-  double get_normalized_value() const {
-    double ret=v_;
-    return ret*std::pow(10.0, EXP);
-  }
-  //! Get the stored value ignoring the exponent
-  double get_value() const {
-    return v_;
-  }
-  template <int OEXP>
-  ExponentialNumber<EXP> operator+(ExponentialNumber<OEXP> o) const {
-    return This(v_+ This(o).v_);
-  }
-  template <int OEXP>
-  ExponentialNumber<EXP> operator-(ExponentialNumber<OEXP> o) const {
-    return This(v_- This(o).v_);
-  }
+struct AtomsPerMol {};
 
-  template <int OEXP>
-  ExponentialNumber<EXP-OEXP> operator/(ExponentialNumber<OEXP> o) const {
-    return ExponentialNumber<EXP-OEXP>(v_/o.v_);
-  }
 
-  template <int OEXP>
-  ExponentialNumber<EXP+OEXP> operator*(ExponentialNumber<OEXP> o) const {
-    return ExponentialNumber<OEXP+EXP>(v_*o.v_);
-  }
+struct MDEnergyTag;
 
-  This operator-() const {
-    return This(-v_);
-  }
+template <>
+inline std::string get_unit_name<MDEnergyTag>(int o) {
+  std::string os[]= {"Cal/Mol"};
+  return os[0];
+}
 
-  void show(std::ostream &out) const {
-    std::ios::fmtflags of = out.flags();
-    out << setiosflags(std::ios::fixed) << v_;
-    if (EXP != 0) out << "e" << EXP;
-    out.flags(of);
-    /** \todo I should restore the io flags.
-     */
-  }
 
-  IMP_COMPARISONS_1(v_);
-};
+struct MDDerivativeTag;
 
-template <int E>
-std::ostream &operator<<(std::ostream &out, ExponentialNumber<E> o)
-{
-  o.show(out);
-  return out;
+template <>
+inline std::string get_unit_name<MDDerivativeTag>(int o) {
+  std::string os[]= {"Cal/(A Mol)"};
+  return os[0];
 }
 
-// helper function generate unit string names
-inline void handle_unit(std::ostringstream &out, std::string name, int m)
-{
-  if (m>0) {
-    out << name;
-    if (m>1) {
-      out << "^" << m<< " ";
-    }
-  }
-}
 
+struct MKSTag{};
 
-//! A class to represent a quantity in the MKS system
-/** The correct units are preserved with arithmetic operations
-    (and only operations with compatible units are allowed). The EXP
-    parameter is there to allow representation of quantities in cm or 
-    other prefixed units.
+template <>
+inline std::string get_unit_name<MKSTag>(int o) {
+  std::string os[]= {"kg", "m", "s", "k", "Cal"};
+  return os[o];
+}
 
-    The template parameters are
-    - EXP the shift in base 10
-    - M the power for meters
-    - KG the power for kilograms
-    - K the power for Kelvin
-    - S the power for seconds
- */
-template <int EXP, int M, int KG, int K, int S>
-class MKSUnit
-{
-  typedef MKSUnit<EXP, M, KG, K, S> This;
-  typedef ExponentialNumber<EXP> V;
 
-  template <int OE, int OM, int OKG, int OK, int OS>
-  friend class MKSUnit;
+struct DaltonTag;
 
-  // for debugging since the debugger doesn't always show
-  // template arguments
-  int get_exponent() const {return EXP;}
-  int get_meters() const {return M;}
-  int get_kilograms() const {return KG;}
-  int get_kelvin() const {return K;}
-  int get_seconds() const {return S;}
+template <>
+inline std::string get_unit_name<DaltonTag>(int o) {
+  std::string os[]= {"Da"};
+  return os[0];
+}
 
+} // ns unit::internal
 
-  V v_;
 
-public:
-  explicit MKSUnit(V v): v_(v) {}
-  MKSUnit(){}
-  MKSUnit(double v): v_(v) {
-    if (v < -std::numeric_limits<double>::max()) {
-      std::cout << get_exponent() << get_meters()
-                << get_kilograms() << get_kelvin() << get_seconds();
-    }
-  }
-  MKSUnit(float v): v_(static_cast<double>(v)) {}
-  MKSUnit(int v): v_(static_cast<double>(v)) {}
 
-  template <int OEXP>
-  MKSUnit(MKSUnit<OEXP, M, KG, K, S> o): v_(o.v_) {}
+typedef boost::mpl::vector_c<int, 1,0,0,0,0> Mass;
+typedef boost::mpl::vector_c<int, 0,1,0,0,0> Length;
+typedef boost::mpl::vector_c<int, 0,0,1,0,0> Time;
+typedef boost::mpl::vector_c<int, 0,0,0,1,0> Temperature;
+typedef boost::mpl::vector_c<int, 0,0,0,0,1> HeatEnergy;
 
-  //! MKS units can be compared if their types match
-  template <int OEXP>
-  bool operator<(MKSUnit<OEXP, M, KG, K, S> o) const {
-    return v_ < o.v_;
-  }
+typedef boost::mpl::vector_c<int, 1,2,-2,0,0> Energy;
+typedef boost::mpl::vector_c<int, 1,1,-2,0,0> Force;
+typedef boost::mpl::vector_c<int, 1,-1,-2,0,0> Pressure;
+typedef boost::mpl::vector_c<int, 0,-1,0,0,1> HeatEnergyDerivative;
 
-  //! MKS units can be compared if their types match
-  template <int OEXP>
-  bool operator<=(MKSUnit<OEXP, M, KG, K, S> o) const {
-    return v_ <= o.v_;
-  }
-  //! MKS units can be compared if their types match
-  template <int OEXP>
-  bool operator>=(MKSUnit<OEXP, M, KG, K, S> o) const {
-    return v_ >= o.v_;
-  }
-  //! only works with matching types
-  template <int OEXP>
-  bool operator>(MKSUnit<OEXP, M, KG, K, S> o) const {
-    return v_ > o.v_;
-  }
 
-  void show(std::ostream &out) const {
-    out << v_ << get_name();
-  }
+typedef Unit<internal::MKSTag, 0, Length> Meter;
+typedef Unit<internal::MKSTag, 2, Length> Centimeter;
+typedef Unit<internal::MKSTag, 0, Mass> Kilogram;
+typedef Unit<internal::MKSTag, -3, Mass> Gram;
+typedef Unit<internal::MKSTag, 0, Time> Second;
+typedef Unit<internal::MKSTag, 0, Energy> Joule;
+typedef Unit<internal::MKSTag, 0, Temperature> Kelvin;
+typedef Shift<Meter, -10>::type Angstrom;
+typedef Multiply<Angstrom, Angstrom>::type SquareAngstrom;
+typedef Shift<Joule, -15>::type Femtojoule;
+typedef Shift<Joule, -12>::type Picojoule;
+typedef Divide<Joule, Kelvin>::type JoulePerKelvin;
+typedef Multiply<Meter, Meter>::type SquareMeter;
+typedef Divide<SquareMeter, Second>::type SquareMeterPerSecond;
+typedef Shift<SquareMeterPerSecond, -4>::type SquareCentimeterPerSecond;
+typedef Unit<internal::MKSTag, -15, Time> Femtosecond;
+typedef Unit<internal::MKSTag, -15, Force> Femtonewton;
+typedef Unit<internal::MKSTag, -12,Force> Piconewton;
+typedef Unit<internal::MKSTag, -18, Mass> Femtogram;
+typedef Divide<Piconewton, Angstrom>::type
+PiconewtonPerAngstrom;
+typedef Unit<internal::MKSTag, 0, Pressure> Pascal;
+typedef Unit<internal::MKSTag, 0, HeatEnergy> Kilocalorie;
+typedef Unit<internal::MKSTag, -21, HeatEnergy> YoctoKilocalorie;
+typedef Divide<Kilocalorie, Meter>::type KilocaloriePerMeter;
+typedef Divide<Kilocalorie, Angstrom>::type KilocaloriePerAngstrom;
+// Calorie is ambiguous given our naming convention
+typedef Shift<KilocaloriePerMeter, -21>::type YoctoKilocaloriePerMeter;
+typedef Shift<KilocaloriePerAngstrom, -21>::type YoctoKilocaloriePerAngstrom;
+typedef Multiply<Multiply<Centimeter, Centimeter>::type,
+                      Centimeter>::type CubicCentimeter;
+typedef Divide<Gram, CubicCentimeter>::type GramPerCubicCentimeter;
 
-  //! Get the current value
-  /** \note The value returned is the value before it is multiplied
-      by the appropriate power of 10. That means, 1 Angstrom returns 1,
-      as does 1 Meter.
-   */
-  double get_value() const {
-    return v_.get_value();
-  }
-  //! Get the value with the proper exponent
-  /** 1 Angstrom returns 1e-10.
-   */
-  double get_normalized_value() const {
-    return v_.get_normalized_value();
-  }
 
-  V get_exponential_value() const {
-    return v_;
-  }
 
-  template <int OEXP, int OM, int OKG, int OK, int OS>
-  MKSUnit<OEXP+EXP, OM+M, OKG+KG, OK+K, OS+S>
-  operator*(MKSUnit<OEXP, OM, OKG, OK, OS> o) const {
-    return MKSUnit<OEXP+EXP, OM+M, OKG+KG, OK+K, OS+S>(v_*o.v_);
-  }
+template <int EXP, class Units>
+typename Exchange<Unit<internal::MKSTag, EXP, Units>,
+        Kilocalorie, Joule, 4>::type
+convert_Cal_to_J(Unit<internal::MKSTag, EXP, Units> i) {
+  typedef typename Exchange<Unit<internal::MKSTag, EXP, Units>,
+    Kilocalorie, Joule, 4>::type Return;
+  return Return(i.get_exponential_value()
+                * JOULES_PER_KILOCALORIE);
+}
 
-  template <int OEXP, int OM, int OKG, int OK, int OS>
-  MKSUnit<EXP-OEXP, M-OM, KG-OKG, K-OK, S-OS>
-  operator/(MKSUnit<OEXP, OM, OKG, OK, OS> o) const {
-    return MKSUnit<EXP-OEXP, M-OM, KG-OKG, K-OK, S-OS>(v_/o.v_);
-  }
+template <int EXP, class Units>
+typename Exchange<Unit<internal::MKSTag, EXP, Units>,
+        Joule, Kilocalorie, -3>::type
+convert_J_to_Cal(Unit<internal::MKSTag, EXP, Units> i) {
+  typedef typename Exchange<Unit<internal::MKSTag, EXP, Units>,
+    Joule, Kilocalorie, -3>::type Return;
+  std::cout << i.get_exponential_value()
+    / JOULES_PER_KILOCALORIE << std::endl;
+  return Return(i.get_exponential_value()
+                / JOULES_PER_KILOCALORIE);
+}
 
-  This operator+(This o) const {
-    return This(v_+o.v_);
-  }
-  This operator-(This o) const {
-    return This(v_-o.v_);
-  }
 
-  This operator-() const {
-    return This(-v_);
-  }
+//! Marker constant to handle coversions
+extern const IMPDLLEXPORT internal::AtomsPerMol ATOMS_PER_MOL;
 
-  //! Get a string version of the name of the units
-  std::string get_name() const {
-    int m=M;
-    int kg=KG;
-    int k=K;
-    int s=S;
-    int j=0;
 
-    while (kg >=1 && m > 1 && s<=-2) {
-      m-=2;
-      kg-=1;
-      s+=2;
-      j+=1;
-    }
-    std::ostringstream name;
-    handle_unit(name, "m", m);
-    handle_unit(name, "kg", kg);
-    handle_unit(name, "K", k);
-    handle_unit(name, "s", s);
-    handle_unit(name, "J", j);
 
-    if (m < 0 || kg < 0 || k < 0 || s < 0 || j < 0) {
-      name << "/";
-      handle_unit(name, "m", -m);
-      handle_unit(name, "kg", -kg);
-      handle_unit(name, "K", -k);
-      handle_unit(name, "s", -s);
-      handle_unit(name, "J", -j);
-    }
-    return name.str();
-  }
-};
+typedef Unit<internal::MDEnergyTag, 0, SingletonUnit> KilocaloriePerMol;
+typedef Unit<internal::MDDerivativeTag, 0, SingletonUnit>
+KilocaloriePerAngstromPerMol;
 
-template <int EXP, int M, int KG, int K, int S>
-MKSUnit<EXP/2, M/2, KG/2, K/2, S/2>
-sqrt(MKSUnit<EXP, M, KG, K, S> o)
-{
-  return MKSUnit<EXP/2, M/2, KG/2, K/2, S/2>(std::sqrt(o.get_value()));
+inline Unit<internal::MKSTag, -24, HeatEnergy >
+operator/(KilocaloriePerMol k,
+                                                   internal::AtomsPerMol) {
+  return Unit<internal::MKSTag, -24, HeatEnergy >(k.get_value()*1.661);
 }
 
-template <int EXP, int M, int KG, int K, int S>
-MKSUnit<EXP*2, M*2, KG*2, K*2, S*2>
-square(MKSUnit<EXP, M, KG, K, S> o)
-{
-  return MKSUnit<EXP*2, M*2, KG*2, K*2, S*2>(IMP::square(o.get_value()));
+inline Unit<internal::MKSTag, -14, HeatEnergyDerivative >
+operator/(KilocaloriePerAngstromPerMol k,
+          internal::AtomsPerMol) {
+  return Unit<internal::MKSTag, -14,
+    HeatEnergyDerivative >(k.get_value()*1.661);
 }
 
-//! They have to be the same EXP to make the return type make sense
-template <int EXP, int M, int KG, int K, int S>
-MKSUnit<EXP, M, KG, K, S>
-min(MKSUnit<EXP, M, KG, K, S> a, MKSUnit<EXP, M, KG, K, S> b) {
-  if (a<b) return a;
-  else return b;
+template <int EXP>
+KilocaloriePerMol operator*(Unit<internal::MKSTag, EXP, HeatEnergy> k,
+                            internal::AtomsPerMol) {
+  return KilocaloriePerMol(k.get_exponential_value()
+                           *NA);
 }
 
-//! They have to be the same EXP to make the return type make sense
-template <int EXP, int M, int KG, int K, int S>
-MKSUnit<EXP, M, KG, K, S>
-max(MKSUnit<EXP, M, KG, K, S> a, MKSUnit<EXP, M, KG, K, S> b) {
-  if (a>b) return a;
-  else return b;
+template <int EXP>
+KilocaloriePerAngstromPerMol
+operator*(Unit<internal::MKSTag, EXP,
+          HeatEnergyDerivative> k,
+          internal::AtomsPerMol) {
+  return KilocaloriePerAngstromPerMol(k.get_exponential_value()
+                                      *NA*ExponentialNumber<-10>(1));
 }
 
-typedef MKSUnit<-10, 1, 0, 0, 0> Angstrom;
-typedef MKSUnit<-20, 2, 0, 0, 0> SquaredAngstrom;
-typedef MKSUnit<0, 2, 1, 0, -2> Joule;
-typedef MKSUnit<-15, 2, 1, 0, -2> FemtoJoule;
-typedef MKSUnit<-12, 2, 1, 0, -2> PicoJoule;
-typedef MKSUnit<0, 0, 0, 1, 0> Kelvin;
-typedef MKSUnit<0,2,1,-1,-2> JoulePerKelvin;
-typedef MKSUnit<-15, 0,0,0,1> FemtoSecond;
-typedef MKSUnit<-4,2,0,0,-1> Cm2PerSecond;
-typedef MKSUnit<-15,1,1,0,-2> FemtoNewton;
-typedef MKSUnit<-12,1,1,0,-2> PicoNewton;
-typedef MKSUnit<-2, 0, 1, 0, -2> PicoNewtonPerAngstrom;
-typedef MKSUnit<0,1,0,0,0> Meter;
-typedef MKSUnit<0,0,0,0,0> Scalar;
-
-template <int EXP, int M, int KG, int K, int S>
-std::ostream &operator<<(std::ostream &out, MKSUnit<EXP, M, KG, K, S> u)
-{
-  u.show(out);
-  return out;
+template <int EXP>
+KilocaloriePerMol operator*(internal::AtomsPerMol,
+                            Unit<internal::MKSTag, EXP,HeatEnergy> k) {
+  return operator*(k, ATOMS_PER_MOL);
 }
 
 
-//! Represent a value with units in kcal and angstroms per mol
-/** This class is just for representing energies and derivatives
-    for MD-type computations and so is a bit more limited than
-    the MKS unit type.
- */
-template <int EXP, int KCal, int A>
-class KCalPerMolUnit
-{
-  typedef KCalPerMolUnit<EXP, KCal, A> This;
-  template <int E, int K, int AA>
-  friend class KCalPerMolUnit;
-  typedef ExponentialNumber<EXP> V;
-  V v_;
 
-  int get_exponent() const {return EXP;}
-  int get_kilocalories() const {return KCal;}
-  int get_angstroms() const {return A;}
-public:
-  explicit KCalPerMolUnit(ExponentialNumber<EXP> e): v_(e) {}
-  KCalPerMolUnit(float f): v_(f) {
-    if (0 < -std::numeric_limits<double>::max()) {
-      std::cout << get_exponent() << get_angstroms()
-                << get_kilocalories();
-    }
-  }
+template <int EXP>
+KilocaloriePerAngstromPerMol
+operator*(internal::AtomsPerMol,
+          Unit<internal::MKSTag, EXP,
+          HeatEnergyDerivative >
+          /*typename YoctoKilocaloriePerAngstrom::Units>*/ k) {
+  return operator*(k, ATOMS_PER_MOL);
+}
 
-  template <int OEXP>
-  KCalPerMolUnit(KCalPerMolUnit<OEXP, KCal, A> o):v_(o.v_) {
-  }
 
-  template <int OEXP>
-  KCalPerMolUnit(MKSUnit<OEXP, 0,0,0,0> o): v_(o.v_) {}
 
-  //! Get the value ignoring the exponent
-  double get_value() const {return v_.get_value();}
 
-  This operator-() const {
-    return This(-v_);
-  }
+// define Daltons
 
-  //! Get the value ignoring the exponent
-  V get_exponential_value() const {return v_;}
+typedef Unit<internal::DaltonTag, 0, SingletonUnit> Dalton;
+typedef Unit<internal::DaltonTag, 3, SingletonUnit> Kilodalton;
 
-  void show(std::ostream &out) const {
-    out << v_ << get_name();
-  }
-
-  //! Get a string representing the units
-  std::string get_name() const {
-    int exp=EXP;
-    int kcal=KCal;
-    int a=A;
-
-    std::ostringstream name;
-    handle_unit(name, "kCal", kcal);
-    handle_unit(name, "A", a);
-
-    if (kcal < 0 || a < 0) {
-      name << "/";
-      handle_unit(name, "kCal", -kcal);
-      handle_unit(name, "A", -a);
-    }
-    return name.str();
-  }
-
-  template <int OEXP, int OKC, int OA>
-  KCalPerMolUnit<OEXP+EXP, OKC+KCal, OA+A>
-  operator*(KCalPerMolUnit<OEXP, OKC, OA> o) const {
-    return KCalPerMolUnit<OEXP+EXP, OKC+KCal, OA+A>(v_*o.v_);
-  }
-
-  template <int OEXP, int OKC, int OA>
-  KCalPerMolUnit<EXP-OEXP, KCal-OKC, A-OA>
-  operator/(KCalPerMolUnit<OEXP, OKC, OA> o) const {
-    return KCalPerMolUnit<EXP-OEXP, KCal-OKC, A-OA>(v_/o.v_);
-  }
-
-  template <int OEXP>
-  This operator+(KCalPerMolUnit<OEXP, KCal, A> o) const {
-    return This(v_+o.v_);
-  }
-  template <int OEXP>
-  This operator-(KCalPerMolUnit<OEXP, KCal, A> o) const {
-    return This(v_-o.v_);
-  }
-
-  template <int OEXP>
-  bool operator<(KCalPerMolUnit<OEXP, KCal, A> o) const {
-    return v_ < o.v_;
-  }
-  template <int OEXP>
-  bool operator>(KCalPerMolUnit<OEXP, KCal, A> o) const {
-    return v_ > o.v_;
-  }
-};
-
-
-
-// want it for conversions
-extern IMPDLLEXPORT const ExponentialNumber<23> NA;
-
-extern IMPDLLEXPORT const ExponentialNumber<3> JOULES_PER_KCAL;
-extern IMPDLLEXPORT const ExponentialNumber<-4> KCALS_PER_JOULE;
-
-
-/** \fn convert_to_kcal
-    \brief Convert from MKS units to kcal units
-
-    \fn convert_to_mks
-    \brief Convert from kcal units to MKS units.
-
-    These don't yet work for arbitrary units since it is kind of
-    annoying to get right and there are not equivalents for everything.
- */
 template <int EXP>
-KCalPerMolUnit<EXP+19, 1, 0>
-convert_to_kcal(MKSUnit<EXP, 2, 1, 0, -2> j) {
-  return KCalPerMolUnit<EXP+19, 1, 0>(j.get_exponential_value()
-                                      *KCALS_PER_JOULE*NA);
+Unit<internal::MKSTag, EXP-26, Mass>
+convert_to_mks(Unit<internal::DaltonTag, EXP, SingletonUnit> d) {
+  return Unit<internal::MKSTag, EXP-26, Mass>(d.get_value()/NA.get_value());
 }
 
 template <int EXP>
-KCalPerMolUnit<EXP+9, 1, -1>
-convert_to_kcal(MKSUnit<EXP, 1, 1, 0, -2> j) {
-  internal::ExponentialNumber<EXP+9> en(j.get_exponential_value()
-                                        *KCALS_PER_JOULE*NA
-                                        *internal::ExponentialNumber<-10>(1));
-  return KCalPerMolUnit<EXP+9, 1, -1>(en);
+Unit<internal::DaltonTag, EXP+26, SingletonUnit>
+convert_to_Dalton(Unit<internal::MKSTag, EXP, Mass> d) {
+  return Unit<internal::DaltonTag,
+    EXP+26, SingletonUnit>(d.get_value()*NA.get_value());
 }
 
-template <int EXP>
-MKSUnit<EXP-20, 2,1,0,-2>
-convert_to_mks(KCalPerMolUnit<EXP, 1, 0> o) {
-  return MKSUnit<EXP-20, 2,1,0,-2>(o.get_exponential_value()
-                                   *JOULES_PER_KCAL/NA);
-}
+} // namespace unit
 
-template <int EXP>
-MKSUnit<EXP-10, 1,1,0,-2>
-convert_to_mks(KCalPerMolUnit<EXP, 1, -1> o) {
-  return MKSUnit<EXP-10, 1,1,0,-2>(o.get_exponential_value()
-                                   *JOULES_PER_KCAL/NA
-                                   *internal::ExponentialNumber<10>(1));
-}
-
-//! The unit for MD energies
-typedef KCalPerMolUnit<0, 1, 0> KCalPerMol;
-//! The unit for MD-compatible derivatives
-typedef KCalPerMolUnit<0, 1, -1> KCalPerAMol;
-
-
-template <int EXP, int KCal, int A>
-std::ostream &operator<<(std::ostream &out, KCalPerMolUnit<EXP, KCal, A> u)
-{
-  u.show(out);
-  return out;
-}
-
 } // namespace internal
 
+namespace unit= internal::unit;
+
 } // namespace IMP
 
 #endif  /* __IMP_UNITS_H */
Index: kernel/include/IMP/internal/ExponentialNumber.h
===================================================================
--- kernel/include/IMP/internal/ExponentialNumber.h	(revision 0)
+++ kernel/include/IMP/internal/ExponentialNumber.h	(revision 0)
@@ -0,0 +1,129 @@
+/**
+ *  \file ExponentialNumber.h
+ *  \brief Classes to add compile time exponents to numbers.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#ifndef __IMP_EXPONENTIAL_NUMBER_H
+#define __IMP_EXPONENTIAL_NUMBER_H
+
+#include "../macros.h"
+
+#include <iostream>
+#include <cmath>
+#include <iomanip>
+#include <limits>
+
+namespace IMP
+{
+
+namespace internal
+{
+
+namespace unit
+{
+
+/** This helper class implements a floating point number
+    with an extra exponent. The actual value is the number 
+    times 10^EXP. I needed this since I was running out of bits
+    in the exponent when converting between vastly different units.
+    It also means that you can have a type for both nanometer
+    and meter without much pain. Since the exponent is stored
+    as a template argument, there shouldn't be any runtime overhead.
+ */
+template <int EXP>
+class ExponentialNumber
+{
+  double v_;
+
+  typedef ExponentialNumber<EXP> This;
+
+  template <int OEXP>
+  void copy_from(ExponentialNumber<OEXP> o) {
+    const int diff = OEXP-EXP;
+    v_=o.v_;
+    float factor= std::pow(10.0, static_cast<double>(diff));
+    v_*=factor;
+  }
+  bool is_default() const {
+    return v_== std::numeric_limits<double>::infinity();
+  }
+
+  template <int E>
+  friend class ExponentialNumber;
+public:
+
+
+  ExponentialNumber(): v_(std::numeric_limits<double>::infinity()) {}
+  template <int OEXP>
+  ExponentialNumber(ExponentialNumber<OEXP> o) {
+    copy_from(o);
+  }
+  explicit ExponentialNumber(double d): v_(d){
+  }
+  template <int OEXP>
+  This &operator=(ExponentialNumber<OEXP> o) {
+    copy_from(o);
+    return *this;
+  }
+
+  //! Compute the value with the exponent
+  double get_normalized_value() const {
+    double ret=v_;
+    return ret*std::pow(10.0, EXP);
+  }
+  //! Get the stored value ignoring the exponent
+  double get_value() const {
+    return v_;
+  }
+  template <int OEXP>
+  ExponentialNumber<EXP> operator+(ExponentialNumber<OEXP> o) const {
+    return This(v_+ This(o).v_);
+  }
+  template <int OEXP>
+  ExponentialNumber<EXP> operator-(ExponentialNumber<OEXP> o) const {
+    return This(v_- This(o).v_);
+  }
+
+  template <int OEXP>
+  ExponentialNumber<EXP-OEXP> operator/(ExponentialNumber<OEXP> o) const {
+    return ExponentialNumber<EXP-OEXP>(v_/o.v_);
+  }
+
+  template <int OEXP>
+  ExponentialNumber<EXP+OEXP> operator*(ExponentialNumber<OEXP> o) const {
+    return ExponentialNumber<OEXP+EXP>(v_*o.v_);
+  }
+
+  This operator-() const {
+    return This(-v_);
+  }
+
+  void show(std::ostream &out) const {
+    std::ios::fmtflags of = out.flags();
+    out << setiosflags(std::ios::fixed) << v_;
+    if (EXP != 0) out << "e" << EXP;
+    out.flags(of);
+    /** \todo I should restore the io flags.
+     */
+  }
+
+  IMP_COMPARISONS_1(v_);
+};
+
+template <int E>
+std::ostream &operator<<(std::ostream &out, ExponentialNumber<E> o)
+{
+  o.show(out);
+  return out;
+}
+
+} // namespace unit
+
+} // namespace internal
+
+} // namespace IMP
+
+#endif  /* __IMP_EXPONENTIAL_NUMBER_H */
Index: kernel/include/IMP/internal/constants.h
===================================================================
--- kernel/include/IMP/internal/constants.h	(revision 589)
+++ kernel/include/IMP/internal/constants.h	(working copy)
@@ -21,16 +21,20 @@
 // make them doubles so we don't have to worry about digits
 
 //! Avagadro's number
-extern IMPDLLEXPORT const ExponentialNumber<23> NA;
+extern IMPDLLEXPORT const unit::ExponentialNumber<23> NA;
+
 //! Boltzmann constant in J/K
-extern IMPDLLEXPORT const MKSUnit<-23, 2, 1, -1, -2>  KB;
+extern IMPDLLEXPORT const
+unit::Shift<unit::Divide<unit::Joule, unit::Kelvin>::type, -23>::type KB;
 
 //! the default temperature
-extern IMPDLLEXPORT const Kelvin T;
+extern IMPDLLEXPORT const unit::Kelvin T;
 
 //! Pi
 static const double PI = 3.1415926535897931;
 
+extern IMPDLLEXPORT const unit::ExponentialNumber<3> JOULES_PER_KILOCALORIE;
+
 } // namespace internal
 
 } // namespace IMP
Index: kernel/src/internal/constants.cpp
===================================================================
--- kernel/src/internal/constants.cpp	(revision 589)
+++ kernel/src/internal/constants.cpp	(working copy)
@@ -13,16 +13,15 @@
 namespace internal
 {
 
-const ExponentialNumber<23> NA(6.02214179);
-
+const unit::ExponentialNumber<23> NA(6.02214179);
 // definition of KB
-const MKSUnit<-23, 2, 1, -1, -2>  KB(1.3806503);
+const unit::Shift<unit::Divide<unit::Joule, unit::Kelvin>::type, -23>::type
+KB(1.3806503);
 
-const Kelvin T(297.15);
+const unit::Kelvin T(297.15);
 
 
-const ExponentialNumber<3> JOULES_PER_KCAL(4.1868);
-const ExponentialNumber<-4> KCALS_PER_JOULE(2.388459);
+const unit::ExponentialNumber<3> JOULES_PER_KILOCALORIE(4.1868);
 
 } // namespace internal
 
Index: kernel/include/IMP/optimizers/BrownianDynamics.h
===================================================================
--- kernel/include/IMP/optimizers/BrownianDynamics.h	(revision 589)
+++ kernel/include/IMP/optimizers/BrownianDynamics.h	(working copy)
@@ -40,17 +40,20 @@
   BrownianDynamics(FloatKey dkey= FloatKey("D"));
   virtual ~BrownianDynamics();
 
-  IMP_OPTIMIZER(internal::kernel_version_info)
+  IMP_OPTIMIZER(internal::kernel_version_info);
 
+  //! Simulate until the given time in fs
+  void simulate(float time_in_fs);
+
   //! Set the max step that should be allowed in angstroms
   /** The timestep is shrunk if any particles are moved more than this */
   void set_max_step(Float a) {
-    set_max_step(internal::Angstrom(a));
+    set_max_step(unit::Angstrom(a));
   }
 
   //! Set time step in fs
   void set_time_step(Float t) {
-    set_time_step(internal::FemtoSecond(t));
+    set_time_step(unit::Femtosecond(t));
   }
 
   //! Time step in fs
@@ -64,33 +67,37 @@
   }
 
 
-  void set_temperature(Float t) { set_temperature(internal::Kelvin(t)); }
+  void set_temperature(Float t) { set_temperature(unit::Kelvin(t)); }
   Float get_temperature() const { return get_temperature_units().get_value(); }
 
   FloatKey get_d_key() const { return dkey_; }
 
-  //! Estimate the diffusion coeffient from the mass in kD
-  /** This method does not do anything fancy, but is OK for large
-      globular proteins.
+  //! Estimate the radius of a protein from the mass
+  /** Proteins are assumed to be spherical. The density is estimated 
+      using the formula from 
+      <http://www.proteinscience.org/cgi/content/full/13/10/2825>
 
-      Note that it depends on temperature and so can't be static.
+      The formula is:
+      density= 1.410+ 0.145 exp(-M/13) g/cm^3
    */
-  Float estimate_radius_from_mass(Float mass_in_kd) const {
-    return estimate_radius_from_mass_units(mass_in_kd).get_value();
+  static Float estimate_radius_from_mass(Float mass_in_kd) {
+    return
+ estimate_radius_from_mass_units(unit::Kilodalton(mass_in_kd)).get_value();
   }
 
   //! Return the expected distance moved for a particle with a given D
   /** The units on D are cm^2/sec and the return has units of Angstroms.
    */
   Float compute_sigma_from_D(Float D) const {
-    return compute_sigma_from_D(internal::Cm2PerSecond(D)).get_value();
+    unit::SquareCentimeterPerSecond du(D);
+    return compute_sigma_from_D(du).get_value();
   }
 
   //! Estimate the diffusion coefficient from the radius in Angstroms
   /** This depends on the temperature.
    */
   Float estimate_D_from_radius(Float radius_in_angstroms) const {
-    internal::Angstrom r(radius_in_angstroms);
+    unit::Angstrom r(radius_in_angstroms);
     return estimate_D_from_radius(r).get_value();
   }
 
@@ -98,7 +105,8 @@
   /** This value is in KCal/A/mol
    */
   Float get_force_scale_from_D(Float D) const {
-    return get_force_scale_from_D(internal::Cm2PerSecond(D)).get_value();
+    unit::SquareCentimeterPerSecond du(D);
+    return get_force_scale_from_D(du).get_value();
   }
 
   //! Get the current time in femtoseconds
@@ -108,7 +116,7 @@
 
   //! Set the current time in femtoseconds
   void set_current_time(Float fs) {
-    set_current_time( internal::FemtoSecond(fs));
+    set_current_time( unit::Femtosecond(fs));
   }
 
   //! Return a histogram of timesteps
@@ -122,55 +130,62 @@
   IMP_LIST(private, Particle, particle, Particle*);
 
 protected:
+
   //! Perform a single dynamics step.
   //  \return true if the initial step size was OK
   bool step();
 
+  void take_step();
+
   /** Propose a single step, this will be accepted if all moves are
       small enough.
       \return true if it should be accepted.
    */
   bool propose_step(std::vector<Vector3D> &proposal);
 
-  internal::FemtoJoule kt() const;
+  unit::Femtojoule kt() const;
 
   void setup_particles();
 
-  void set_max_step(internal::Angstrom a) {
+  void set_max_step(unit::Angstrom a) {
     max_change_= a;
   }
 
-  void set_time_step(internal::FemtoSecond t);
+  void set_time_step(unit::Femtosecond t);
 
-  internal::FemtoSecond get_time_step_units() const {return max_dt_;}
+  unit::Femtosecond get_time_step_units() const {return max_dt_;}
 
-  internal::FemtoSecond get_last_time_step_units() const {return cur_dt_;}
+  unit::Femtosecond get_last_time_step_units() const {return cur_dt_;}
 
-  void set_temperature(internal::Kelvin t) { T_=t;}
+  void set_temperature(unit::Kelvin t) { T_=t;}
 
-  internal::Kelvin get_temperature_units() const { return T_;}
+  unit::Kelvin get_temperature_units() const { return T_;}
 
-  internal::Angstrom estimate_radius_from_mass_units(Float mass_in_kd) const;
+  static unit::Angstrom
+    estimate_radius_from_mass_units(unit::Kilodalton mass_in_kd);
 
-  internal::Angstrom compute_sigma_from_D(internal::Cm2PerSecond D) const;
+  unit::Angstrom
+    compute_sigma_from_D(unit::SquareCentimeterPerSecond D) const;
 
-  internal::Cm2PerSecond
-  estimate_D_from_radius(internal::Angstrom radius) const;
+  unit::SquareCentimeterPerSecond
+    estimate_D_from_radius(unit::Angstrom radius) const;
 
-  internal::KCalPerAMol get_force_scale_from_D(internal::Cm2PerSecond D) const;
+  unit::KilocaloriePerAngstromPerMol
+    get_force_scale_from_D(unit::SquareCentimeterPerSecond D) const;
 
-  internal::FemtoSecond get_current_time_units() const {
+  unit::Femtosecond get_current_time_units() const {
     return cur_time_;
   }
 
   //! Set the current time in femtoseconds
-  void set_current_time(internal::FemtoSecond fs) {
+  void set_current_time(unit::Femtosecond fs) {
     cur_time_= fs;
   }
 
-  internal::Angstrom max_change_;
-  internal::FemtoSecond max_dt_, cur_dt_, cur_time_;
-  internal::Kelvin T_;
+  unsigned int num_const_dt_;
+  unit::Angstrom max_change_;
+  unit::Femtosecond max_dt_, cur_dt_, cur_time_;
+  unit::Kelvin T_;
 
   FloatKey dkey_;
 
Index: kernel/include/IMP/unary_functions/WormLikeChain.h
===================================================================
--- kernel/include/IMP/unary_functions/WormLikeChain.h	(revision 589)
+++ kernel/include/IMP/unary_functions/WormLikeChain.h	(working copy)
@@ -31,7 +31,11 @@
   /** \param[in] l_max maximum length of the chain in angstroms
       \param[in] lp persistence length in angstroms
    */
-  WormLikeChain(Float l_max, Float lp) : lmax_(l_max), lp_(lp) {}
+  WormLikeChain(Float l_max, Float lp) : lmax_(l_max), lp_(lp) {
+    IMP_check(l_max > lp, "The persistance length should be less "
+              << "than the total length for this model",
+              ValueException(""));
+  }
 
   virtual ~WormLikeChain() {}
 
@@ -40,19 +44,20 @@
       \return Energy in kcal/mol
    */
   virtual Float evaluate(Float lf) {
-    static const internal::PicoJoule zero=eval(internal::Angstrom(0));
-    internal::Angstrom l(lf);
-    if (l < internal::Angstrom(0)) l=internal::Angstrom(0);
-    internal::PicoJoule ret;
+    static const unit::Picojoule zero=eval(unit::Angstrom(0));
+    unit::Angstrom l(lf);
+    if (l < unit::Angstrom(0)) l=unit::Angstrom(0);
+    unit::Picojoule ret;
     if (l < cutoff()) {
       ret= (eval(l) - zero);
     } else {
-      internal::PicoJoule springterm=(l-cutoff())*cderiv(cutoff());
+      unit::Picojoule springterm=(l-cutoff())*cderiv(cutoff());
       ret= (eval(cutoff())+ springterm -zero);
     }
-    /*std::cout << "Return is " << ret <<" " << l << " " << lp_ << " " 
-      << lmax_ << std::endl;*/
-    return internal::KCalPerMol(convert_to_kcal(ret)).get_value();
+    std::cout << "Return is " << ret <<" " << l << " " << lp_ << " " 
+              << lmax_ << std::endl;
+    unit::YoctoKilocalorie zc= convert_J_to_Cal(ret);
+    return (zc*unit::ATOMS_PER_MOL).get_value();
   }
 
   //! Calculate the WormLikeChain energy given the length
@@ -61,9 +66,9 @@
       \return Score
    */
   virtual Float evaluate_deriv(Float fl, Float& deriv) {
-    internal::Angstrom l(fl);
-    if (l < internal::Angstrom(0)) l=internal::Angstrom(0);
-    internal::PicoNewton doubled;
+    unit::Angstrom l(fl);
+    if (l < unit::Angstrom(0)) l=unit::Angstrom(0);
+    unit::Piconewton doubled;
     if (l < cutoff()) {
       doubled= cderiv(l);
     } else {
@@ -74,8 +79,9 @@
     }
     //std::cout << "Force is " << doubled << " for length " << l << std::endl;
     // convert from picoNewton
-    deriv = internal::KCalPerAMol(internal::convert_to_kcal(doubled))
-      .get_value();
+    unit::YoctoKilocaloriePerAngstrom du= unit::convert_J_to_Cal(doubled);
+
+    deriv = (du*unit::ATOMS_PER_MOL).get_value();
     //std::cout << "Which converts to " << d << std::endl;
     return evaluate(fl);
   }
@@ -86,28 +92,28 @@
 
 protected:
   //! \note named to avoid clash with 'deriv' argument
-  internal::PicoNewton cderiv(internal::Angstrom l) const {
-    internal::PicoNewton pn= internal::KB*internal::T
-      /lp_*(internal::Scalar(.25)/ square(internal::Scalar(1)-l/lmax_)
-            -internal::Scalar(.25)+l/lmax_);
+  unit::Piconewton cderiv(unit::Angstrom l) const {
+    unit::Piconewton pn= internal::KB*internal::T
+      /lp_*(.25/ square(1.0-(l/lmax_).get_normalized_value())
+            -.25+(l/lmax_).to_scalar());
      return pn;
   }
 
-  internal::PicoJoule eval(internal::Angstrom m) const {
-    internal::PicoJoule J
-      =  internal::KB*internal::T/lp_*(internal::Scalar(.25)*square(lmax_)
+  unit::Picojoule eval(unit::Angstrom m) const {
+    unit::Picojoule J
+      =  internal::KB*internal::T/lp_*(.25*square(lmax_)
                                        /(lmax_-m)
-                                       -m*internal::Scalar(.25) 
-                                       +internal::Scalar(.5)*square(m)
+                                       -m*.25
+                                       +.5*square(m)
                                        /lmax_);
     return J;
   }
 
-  internal::Angstrom cutoff() const {
-    return internal::Scalar(.99)*lmax_;
+  unit::Angstrom cutoff() const {
+    return .99*lmax_;
   }
 
-  internal::Angstrom lmax_, lp_;
+  unit::Angstrom lmax_, lp_;
 };
 
 } // namespace IMP
Index: kernel/src/optimizers/BrownianDynamics.cpp
===================================================================
--- kernel/src/optimizers/BrownianDynamics.cpp	(revision 589)
+++ kernel/src/optimizers/BrownianDynamics.cpp	(working copy)
@@ -17,31 +17,34 @@
 
 namespace IMP
 {
-typedef internal::MKSUnit<-3, -1, 1, 0, -1> MilliPascalSeconds;
+typedef
+unit::Shift<unit::Multiply<unit::Pascal, 
+                           unit::Second>::type, 
+            -3>::type MillipascalSecond;
 
-static MilliPascalSeconds eta(internal::Kelvin T)
+static MillipascalSecond eta(unit::Kelvin T)
 {
-  const std::pair<internal::Kelvin, MilliPascalSeconds> points[]
-    ={ std::make_pair(internal::Kelvin(273+10.0),
-                      MilliPascalSeconds(1.308)),
-       std::make_pair(internal::Kelvin(273+20.0),
-                      MilliPascalSeconds(1.003)),
-       std::make_pair(internal::Kelvin(273+30.0),
-                      MilliPascalSeconds(0.7978)),
-       std::make_pair(internal::Kelvin(273+40.0),
-                      MilliPascalSeconds(0.6531)),
-       std::make_pair(internal::Kelvin(273+50.0),
-                      MilliPascalSeconds(0.5471)),
-       std::make_pair(internal::Kelvin(273+60.0),
-                      MilliPascalSeconds(0.4668)),
-       std::make_pair(internal::Kelvin(273+70.0),
-                      MilliPascalSeconds(0.4044)),
-       std::make_pair(internal::Kelvin(273+80.0),
-                      MilliPascalSeconds(0.3550)),
-       std::make_pair(internal::Kelvin(273+90.0),
-                      MilliPascalSeconds(0.3150)),
-       std::make_pair(internal::Kelvin(273+100.0),
-                      MilliPascalSeconds(0.2822))};
+  const std::pair<unit::Kelvin, MillipascalSecond> points[]
+    ={ std::make_pair(unit::Kelvin(273+10.0),
+                      MillipascalSecond(1.308)),
+       std::make_pair(unit::Kelvin(273+20.0),
+                      MillipascalSecond(1.003)),
+       std::make_pair(unit::Kelvin(273+30.0),
+                      MillipascalSecond(0.7978)),
+       std::make_pair(unit::Kelvin(273+40.0),
+                      MillipascalSecond(0.6531)),
+       std::make_pair(unit::Kelvin(273+50.0),
+                      MillipascalSecond(0.5471)),
+       std::make_pair(unit::Kelvin(273+60.0),
+                      MillipascalSecond(0.4668)),
+       std::make_pair(unit::Kelvin(273+70.0),
+                      MillipascalSecond(0.4044)),
+       std::make_pair(unit::Kelvin(273+80.0),
+                      MillipascalSecond(0.3550)),
+       std::make_pair(unit::Kelvin(273+90.0),
+                      MillipascalSecond(0.3150)),
+       std::make_pair(unit::Kelvin(273+100.0),
+                      MillipascalSecond(0.2822))};
 
   const unsigned int npoints= sizeof(points)/sizeof(std::pair<float,float>);
   if (T < points[0].first) {
@@ -49,10 +52,11 @@
   } else {
     for (unsigned int i=1; i< npoints; ++i) {
       if (points[i].first > T) {
-        internal::Scalar f= (T - points[i-1].first)/(points[i].first
-                                                     - points[i-1].first);
-        MilliPascalSeconds ret= 
-          (internal::Scalar(1.0)-f) *points[i-1].second + f*points[i].second;
+        float f= ((T - points[i-1].first)
+                  /(points[i].first - points[i-1].first))
+          .get_normalized_value();
+        MillipascalSecond ret= 
+          (1.0-f) *points[i-1].second + f*points[i].second;
         return ret;
       }
     }
@@ -74,11 +78,10 @@
 {
 }
 
-IMP_LIST_IMPL(BrownianDynamics, Particle, particle, Particle*,
-              {if (0) std::cout << obj << index;},);
+IMP_LIST_IMPL(BrownianDynamics, Particle, particle, Particle*,,);
 
 
-void BrownianDynamics::set_time_step(internal::FemtoSecond t)
+void BrownianDynamics::set_time_step(unit::Femtosecond t)
 {
   time_steps_.clear();
   max_dt_ = t;
@@ -118,7 +121,7 @@
   std::vector<Vector3D> new_pos(number_of_particles());
   bool noshrink=true;
   while (!propose_step(new_pos)) {
-    cur_dt_= cur_dt_/internal::Scalar(2.0);
+    cur_dt_= cur_dt_/2.0;
     noshrink=false;
   }
 
@@ -135,13 +138,6 @@
 //! Perform a single dynamics step.
 bool BrownianDynamics::propose_step(std::vector<Vector3D>& new_pos)
 {
-  // Assuming score is in kcal/mol, its derivatives in kcal/mol/angstrom,
-  // and mass is in g/mol, conversion factor necessary to get accelerations
-  // in angstrom/fs/fs from raw derivatives
-  // we want 
-  //const double eta= 10;// kg/(m s) from ~1 cg/(cm s)
-  //const double pi=M_PI;
-  //const double kTnu= 1.3806504*T_.get_value();
   FloatKeys xyzk=XYZDecorator::get_xyz_keys();
 
   IMP_LOG(VERBOSE, "dt is " << cur_dt_ << std::endl);
@@ -150,18 +146,41 @@
   for (unsigned int i=0; i< number_of_particles(); ++i) {
     Particle *p= get_particle(i);
     XYZDecorator d(p);
+    IMP_IF_CHECK(CHEAP) {
+      for (unsigned int j=0; j< 3; ++j) {
+        // GDB 6.6 prints infinity as 0 on 64 bit machines. Grumble.
+        /*int szf= sizeof(Float);
+        int szi= sizeof(int);
+        Float one=1;*/
+        Float mx= std::numeric_limits<Float>::max();
+        Float c= d.get_coordinate(j);
+        bool ba= (c != c);
+        bool bb = (c >= mx);
+        bool bc= -d.get_coordinate(j) >= std::numeric_limits<Float>::max();
+        if (ba || bb || bc ) {
+          IMP_WARN("Bad value for coordinate in Brownian Dynamics on "
+                   << "particle " << p->get_index() << std::endl);
+          throw ValueException("Bad coordinate value");
+        }
+      }
+    }
+
     //double xi= 6*pi*eta*radius; // kg/s
-    internal::Cm2PerSecond D(p->get_value(dkey_));
-    internal::Angstrom sigma(compute_sigma_from_D(D));
-#ifndef NDEBUG
-    internal::Angstrom osigma(sqrt(internal::Scalar(2.0)*D*cur_dt_));
-#endif
-    IMP_assert(sigma
-               - osigma
-               <= internal::Scalar(.01)* sigma,
-               "Sigma computations don't match " << sigma 
-               << " " 
-               << sqrt(internal::Scalar(2.0)*D*cur_dt_));
+    unit::SquareCentimeterPerSecond D(p->get_value(dkey_));
+    IMP_check(D.get_value() > 0
+              && D.get_value() < std::numeric_limits<Float>::max(),
+              "Bad diffusion coefficient on particle " << p->get_index(),
+              ValueException(""));
+    unit::Angstrom sigma(compute_sigma_from_D(D));
+    IMP_IF_CHECK(EXPENSIVE) {
+      unit::Angstrom osigma(sqrt(2.0*D*cur_dt_));
+      IMP_check(sigma - osigma
+                 <= .01* sigma,
+                 "Sigma computations don't match " << sigma 
+                 << " " 
+                << sqrt(2.0*D*cur_dt_),
+                ErrorException());
+    }
     IMP_LOG(VERBOSE, p->get_index() << ": sigma is " 
             << sigma << std::endl);
     boost::normal_distribution<double> mrng(0, sigma.get_value());
@@ -171,25 +190,22 @@
 
     //std::cout << p->get_index() << std::endl;
 
-    internal::Angstrom delta[3];
+    unit::Angstrom delta[3];
 
     for (unsigned j = 0; j < 3; ++j) {
-      // force originally in kcal/mol/A
-      // derive* 2.390e-4 gives it in J/(mol A)
-      // that * e10 gives it in J/(mol m)
-      // that / NA gives it in J/m
-      internal::KCalPerAMol force(-d.get_coordinate_derivative(j));
-      //*4.1868e+13/NA; old conversion
-      internal::FemtoNewton nforce= convert_to_mks(force);
-      internal::Angstrom R(sampler());
-      internal::Angstrom force_term=nforce*cur_dt_*D/kt();
+      unit::KilocaloriePerAngstromPerMol
+        force( -d.get_coordinate_derivative(j));
+      unit::Femtonewton nforce
+        = unit::convert_Cal_to_J(force/unit::ATOMS_PER_MOL);
+      unit::Angstrom R(sampler());
+      unit::Angstrom force_term(nforce*cur_dt_*D/kt());
       //std::cout << "Force term is " << force_term << " and R is "
       //<< R << std::endl;
-      internal::Angstrom dr= force_term +  R; //std::sqrt(2*kb*T_*ldt/xi) *
+      unit::Angstrom dr= force_term +  R; //std::sqrt(2*kb*T_*ldt/xi) *
       // get back from meters
       delta[j]=dr;
     }
-    //internal::Angstrom max_motion= internal::Scalar(4)*sigma;
+    //unit::Angstrom max_motion= unit::Scalar(4)*sigma;
     /*std::cout << "delta is " << delta << " mag is " 
       << delta.get_magnitude() << " sigma " << sigma << std::endl;*/
 
@@ -200,12 +216,12 @@
             << ", " << d.get_coordinate_derivative(1)
             << ", " << d.get_coordinate_derivative(2) << "]" << std::endl);
 
-    internal::SquaredAngstrom t= square(delta[0]);
+    unit::SquareAngstrom t= square(delta[0]);
 
-    internal::SquaredAngstrom magnitude2=square(delta[0])+square(delta[1])
+    unit::SquareAngstrom magnitude2=square(delta[0])+square(delta[1])
       +square(delta[2]);
 
-    //internal::SquaredAngstrom max_motion2= square(max_motion);
+    //unit::SquaredAngstrom max_motion2= square(max_motion);
     if (magnitude2 > square(max_change_)) {
       return false;
     }
@@ -223,85 +239,112 @@
  */
 Float BrownianDynamics::optimize(unsigned int max_steps)
 {
+ IMP_check(get_model() != NULL, "Must set model before calling optimize",
+           ValueException(""));
   setup_particles();
-  IMP_LOG(SILENT, "Running brownian dynamics on " << get_particles().size() 
+  IMP_LOG(TERSE, "Running brownian dynamics on " << get_particles().size() 
           << " particles with a step of " << cur_dt_ << std::endl);
-  unsigned int num_const_dt=0;
+  setup_particles();
   for (unsigned int i = 0; i < max_steps; ++i) {
-    update_states();
-    get_model()->evaluate(true);
-    if (step()) {
-      ++num_const_dt;
-    }
+    take_step();
+  }
+  return get_model()->evaluate(false);
+}
 
-    cur_time_= cur_time_+cur_dt_;
-    if (num_const_dt == 10) {
-      num_const_dt=0;
-      cur_dt_= std::min(cur_dt_*internal::Scalar(2), max_dt_);
+
+void BrownianDynamics::take_step() {
+  update_states();
+  get_model()->evaluate(true);
+  if (step()) ++ num_const_dt_;
+
+  cur_time_= cur_time_+cur_dt_;
+  if (num_const_dt_ == 10) {
+    num_const_dt_=0;
+    cur_dt_= std::min(cur_dt_*2.0, max_dt_);
+  }
+  {
+    unit::Femtosecond lb(max_dt_/2.0);
+    std::vector<int>::size_type d=0;
+    while (lb > cur_dt_) {
+      lb = lb/2.0;
+      ++d;
     }
-    {
-      internal::FemtoSecond lb(max_dt_/internal::Scalar(2.0));
-      std::vector<int>::size_type d=0;
-      while (lb > cur_dt_) {
-        lb = lb/internal::Scalar(2.0);
-        ++d;
-      }
-      time_steps_.resize(std::max(time_steps_.size(), d+1), 0);
-      ++time_steps_[d];
-    }
+    time_steps_.resize(std::max(time_steps_.size(), d+1), 0);
+    ++time_steps_[d];
   }
-  return get_model()->evaluate(false);
 }
 
+void BrownianDynamics::simulate(float max_time) {
+  IMP_check(get_model() != NULL, "Must set model before calling simulate",
+            ValueException(""));
+  setup_particles();
+  unit::Femtosecond mt(max_time);
+  num_const_dt_=0;
+  while (cur_time_ < mt) {
+    take_step();
+  }
+}
 
 
-internal::FemtoJoule BrownianDynamics::kt() const
+unit::Femtojoule BrownianDynamics::kt() const
 {
-  return internal::KB*T_;
+  return unit::Femtojoule(internal::KB*T_);
 }
 
 
-internal::Angstrom
-BrownianDynamics::estimate_radius_from_mass_units(Float mass_in_kd) const
+unit::Angstrom
+BrownianDynamics
+::estimate_radius_from_mass_units(unit::Kilodalton mass_in_kd)
 {
-  /* Data points used:
-     thyroglobulin 669kD 85A
-     catalase 232kD 52A
-     aldolase 158kD 48A
-   */
+  //unit::KG kg= convert_to_mks(mass_in_kd)
+  unit::Kilodalton kd(mass_in_kd);
+  unit::GramPerCubicCentimeter p= unit::GramPerCubicCentimeter(1.410) 
+    + unit::GramPerCubicCentimeter(0.145)* exp(-mass_in_kd.get_value()/13);
 
-  internal::Angstrom r(std::pow(mass_in_kd, .3333f)*10);
-  return r;
+  unit::Kilogram m= convert_to_mks(kd);
+  // m/(4/3 pi r^3) =p
+  // r= (m/p)/(4/(3 pi))^(1/3)
+  typedef unit::Multiply<unit::Angstrom,
+    unit::Angstrom>::type SquareAngstrom;
+  typedef unit::Multiply<unit::SquareAngstrom,
+    unit::Angstrom>::type CubicAngstrom;
+
+  CubicAngstrom v((m/p)*(4.0/(3.0*internal::PI)));
+
+  return unit::Angstrom(std::pow(v.get_value(), .3333));
 }
 
-internal::Cm2PerSecond
-BrownianDynamics::estimate_D_from_radius(internal::Angstrom r) const
+unit::SquareCentimeterPerSecond
+BrownianDynamics::estimate_D_from_radius(unit::Angstrom r) const
 {
-  MilliPascalSeconds e=eta(T_);
-  internal::MKSUnit<-13, 0, 1, 0, -1> etar( e*r);
+  MillipascalSecond e=eta(T_);
+  //unit::MKSUnit<-13, 0, 1, 0, -1> etar( e*r);
   /*std::cout << e << " " << etar << " " << kt << std::endl;
-  std::cout << "scalar etar " << (internal::Scalar(6*internal::PI)*etar) 
+  std::cout << "scalar etar " << (unit::Scalar(6*unit::PI)*etar) 
             << std::endl;
-  std::cout << "ret pre conv " << (kt/(internal::Scalar(6* internal::PI)*etar))
+  std::cout << "ret pre conv " << (kt/(unit::Scalar(6* unit::PI)*etar))
   << std::endl;*/
-  internal::Cm2PerSecond ret( kt()/(internal::Scalar(6* internal::PI)*etar));
+  unit::SquareCentimeterPerSecond ret(kt()/(6.0* internal::PI*e*r));
   //std::cout << "ret " << ret << std::endl;
   return ret;
 }
 
-internal::Angstrom
-BrownianDynamics::compute_sigma_from_D(internal::Cm2PerSecond D) const
+unit::Angstrom
+BrownianDynamics
+::compute_sigma_from_D(unit::SquareCentimeterPerSecond D) const
 {
-  return sqrt(internal::Scalar(2.0)*D*cur_dt_); //6*xi*kb*T_;
+  return sqrt(2.0*D*cur_dt_); //6*xi*kb*T_;
 }
 
-internal::KCalPerAMol BrownianDynamics
-::get_force_scale_from_D(internal::Cm2PerSecond D) const
+unit::KilocaloriePerAngstromPerMol BrownianDynamics
+::get_force_scale_from_D(unit::SquareCentimeterPerSecond D) const
 {
   // force motion is f*cur_dt_*D/kT
   // sigma*kT/(cur_dt_*D)
-  return convert_to_kcal(internal::Angstrom(compute_sigma_from_D(D))*kt()
-                         /(max_dt_*D));
+  unit::Angstrom s=compute_sigma_from_D(D);
+  unit::Piconewton pn= s*kt()/(max_dt_*D);
+  unit::YoctoKilocaloriePerAngstrom yc=unit::convert_J_to_Cal(pn);
+  return unit::operator*(unit::ATOMS_PER_MOL,yc);
 }