Index: kernel/include/IMP/internal/units.h
===================================================================
--- kernel/include/IMP/internal/units.h	(revision 0)
+++ kernel/include/IMP/internal/units.h	(revision 0)
@@ -0,0 +1,460 @@
+/**
+ *  \file units.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 "../utility.h"
+
+#include <sstream>
+#include <cmath>
+#include <iomanip>
+
+namespace IMP
+{
+
+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;
+
+  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;
+    /*for (int i=0; i< diff; ++i) {
+      v_*=10;
+    }
+    for (int i=0; i> diff; --i) {
+      v_/=10;
+      }*/
+    IMP_assert((get_normalized_value() -o.get_normalized_value())
+               <= .09 * 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:
+
+
+  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;
+  }
+
+  //! 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_);
+  }
+
+  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;
+}
+
+// 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<< " ";
+    }
+  }
+}
+
+
+//! 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.
+
+    The template parameters are
+    - EXP the shift in base 10
+    - M the power for meters
+    - KG the power for kilograms
+    - K the power for degrees 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;
+
+  // 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;}
+
+
+  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_){
+  }
+  //! 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();
+  }
+
+  //! 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 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_);
+  }
+
+  This operator+(This o) const {
+    return This(v_+o.v_);
+  }
+  This operator-(This o) const {
+    return This(v_-o.v_);
+  }
+
+  //! 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();
+  }
+};
+
+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()));
+}
+
+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()));
+}
+
+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;
+//typedef MKSUnit<-3, -1, 1, 0, -1> MilliPascalSeconds;
+
+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;
+}
+
+
+//! 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 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();}
+
+
+  //! Get the value ignoring the exponent
+  V get_exponential_value() const {return v_;}
+
+  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> JoulesPerKiloCalorie;
+extern IMPDLLEXPORT const ExponentialNumber<-4> KiloCaloriesPerJoule;
+
+
+/** \function convert_to_kcal Convert from MKS units to kcal units
+    \function convert_to_mks 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()
+                                      *KiloCaloriesPerJoule*NA);
+}
+
+template <int EXP>
+KCalPerMolUnit<EXP+9, 1, -1>
+convert_to_kcal(MKSUnit<EXP, 1, 1, 0, -2> j) {
+  std::cout << j.get_exponential_value() << "*"
+            << KiloCaloriesPerJoule << "*" 
+            << NA
+            << "*" << internal::ExponentialNumber<-10>(1) << std::endl;
+  internal::ExponentialNumber<EXP+9> en(j.get_exponential_value()
+                                        *KiloCaloriesPerJoule*NA
+                                        *internal::ExponentialNumber<-10>(1));
+  std::cout << "en is " << en << std::endl;
+  return KCalPerMolUnit<EXP+9, 1, -1>(en);
+}
+
+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()
+                                   *JoulesPerKiloCalorie/NA);
+}
+
+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()
+                                   *JoulesPerKiloCalorie/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;
+}
+
+} // internal
+
+}
+
+#endif
Index: kernel/include/IMP/internal/constants.h
===================================================================
--- kernel/include/IMP/internal/constants.h	(revision 0)
+++ kernel/include/IMP/internal/constants.h	(revision 0)
@@ -0,0 +1,37 @@
+/**
+ *  \file constants.h    \brief Various useful constants.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#ifndef __IMP_CONSTANTS_H
+#define __IMP_CONSTANTS_H
+
+#include "units.h"
+
+#include <cmath>
+
+namespace IMP
+{
+
+namespace internal
+{
+
+// make them doubles so we don't have to worry about digits
+
+//! Pi?
+static const double PI=M_PI;
+//! Avagadro's number from wikipedia
+extern IMPDLLEXPORT const ExponentialNumber<23> NA;
+//! Boltzmann constant in J/K
+extern IMPDLLEXPORT const MKSUnit<-23, 2, 1, -1, -2>  KB;
+
+//! the default temperature
+extern IMPDLLEXPORT const Kelvin T;
+
+}
+
+} // namespace IMP
+
+#endif  /* __IMP_KEY_H */
Index: kernel/src/internal/constants.cpp
===================================================================
--- kernel/src/internal/constants.cpp	(revision 0)
+++ kernel/src/internal/constants.cpp	(revision 0)
@@ -0,0 +1,22 @@
+#include "IMP/internal/constants.h"
+
+namespace IMP
+{
+
+namespace internal
+{
+
+const ExponentialNumber<23> NA(6.02214179);
+
+// definition of KB
+const MKSUnit<-23, 2, 1, -1, -2>  KB(1.3806503);
+
+const Kelvin T(297.15);
+
+
+const ExponentialNumber<3> JoulesPerKiloCalorie(4.1868);
+const ExponentialNumber<-4> KiloCaloriesPerJoule(2.388459);
+
+}
+
+}