[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[IMP-dev] removing modeldata



I wanted some easy coding so I made a quick pass at removing ModelData. The interesting changes are that Optimizer now has a couple of ModelData-esque helper methods to enable the ConjugateGradients and SteepestDescent optimizers to work with minimal changes (these allow easy access to the optimized attributes). Modifying the optimizers could make them more efficient (by grouping the attributes by particle), but that can wait.

In addition, due to some issues that emerged with the changes I modified and cleaned up MolecularDynamics a bit. The main change is it now skips particles which don't have optimizable x,y,z. This seems more consistent with the idea of an optimizer (it doesn' have to actually improve every optimizable attribute) and we need the optimizer to work when there are other sorts of particles floating around. It uses an IMP_LIST in order to simplify moving to reference counted particles if we do that.

Anyway, the tests pass. Check standards still bombs, so I don't know if it passes that.


Index: kernel/test/misc/test_restraint_sets.py
===================================================================
--- kernel/test/misc/test_restraint_sets.py	(revision 434)
+++ kernel/test/misc/test_restraint_sets.py	(working copy)
@@ -28,7 +28,6 @@
 
     def test_weights(self):
         """Test that sets can be weighted"""
-        model_data = self.model.get_model_data()
         e1 = self.model.evaluate(True)
         d1 = self.particles[0].get_derivative(IMP.FloatKey("x"))
         self.rset.set_weight(0.5)
Index: kernel/test/xml/test_xml.py
===================================================================
--- kernel/test/xml/test_xml.py	(revision 434)
+++ kernel/test/xml/test_xml.py	(working copy)
@@ -20,7 +20,6 @@
         """Check reading of XML files"""
         self.doc = IMP.xml_loader.load_imp_model(self.imp_model,
                                                  "xml/model.xml")
-        model_data = self.imp_model.get_model_data()
 
         # test particles
         num_particles = 0
Index: kernel/test/modeller/test_proximity.py
===================================================================
--- kernel/test/modeller/test_proximity.py	(revision 434)
+++ kernel/test/modeller/test_proximity.py	(working copy)
@@ -110,7 +110,6 @@
         os.unlink('out_proximity.pdb')
 
         # min distances
-        model_data = self.imp_model.get_model_data()
         for i in range(len(coords)):
             irad = self.particles[i].get_value(radius)
             for j in range(i+1,len(coords)):
Index: kernel/test/optimizers/test_cg_optimizer.py
===================================================================
--- kernel/test/optimizers/test_cg_optimizer.py	(revision 434)
+++ kernel/test/optimizers/test_cg_optimizer.py	(working copy)
@@ -53,7 +53,6 @@
         opt.set_model(model)
         opt.set_threshold(1e-5)
         e = opt.optimize(100)
-        model_data = model.get_model_data()
         for p in particles:
             val = p.get_value(IMP.FloatKey("x"))
             self.assertAlmostEqual(val, 1.0, places=1)
Index: kernel/test/optimizers/test_md_optimizer.py
===================================================================
--- kernel/test/optimizers/test_md_optimizer.py	(revision 434)
+++ kernel/test/optimizers/test_md_optimizer.py	(working copy)
@@ -99,17 +99,17 @@
                                lambda a: a + strength * delttm)
 
     def test_non_xyz(self):
-        """Should be unable to do MD on optimizable non-xyz attributes"""
+        """Should skip particles without xyz attributes"""
         p = IMP.Particle()
         self.model.add_particle(p)
         p.add_attribute(IMP.FloatKey("attr"), 0.0, True)
-        self.assertRaises(ValueError, self.md.optimize, 50)
+        self.md.optimize(100)
 
     def test_make_velocities(self):
         """Test that MD generates particle velocities"""
         self.md.optimize(0)
         keys = [IMP.FloatKey(x) for x in ("vx", "vy", "vz")]
-        for p in self.particles:
+        for p in self.model.get_particles():
             for key in keys:
                 self.assert_(p.has_attribute(key))
 
Index: kernel/include/IMP/ModelData.h
===================================================================
--- kernel/include/IMP/ModelData.h	(revision 434)
+++ kernel/include/IMP/ModelData.h	(working copy)
@@ -1,202 +0,0 @@
-/**
- *  \file ModelData.h    \brief Storage for all model particle data.
- *
- *  Copyright 2007-8 Sali Lab. All rights reserved.
- *
- */
-
-#ifndef __IMP_MODEL_DATA_H
-#define __IMP_MODEL_DATA_H
-
-#include "IMP_config.h"
-#include "base_types.h"
-
-#include <boost/iterator/filter_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <vector>
-
-namespace IMP
-{
-
-class Model;
-class Particle;
-class DerivativeAccumulator;
-
-
-//! The interface that optimizers use to access Particle attributes.
-/** All data for particles is stored in this structure. Float values
-    are all differentiable. Each Float value has an associated flag
-    saying whether it is optimized or not. Non-optimized values should
-    not be changed by the optimizer.
-
-    \note ModelData should only be used in writing Optimizers.
- */
-class IMPDLLEXPORT ModelData
-{
-  friend class Model;
-  friend class Particle;
-  friend class OptFloatIndexIterator;
-  friend class DerivativeAccumulator;
-  friend class std::auto_ptr<ModelData>;
-
-  struct FloatData
-  {
-    Float value_;
-    Float deriv_;
-    int stats_index_;
-    bool is_optimized_;
-  };
-
-  struct FloatIsOptimized
-  {
-    const ModelData *m_;
-  public:
-    FloatIsOptimized(const ModelData *m): m_(m){}
-    bool operator()(FloatIndex f) const {
-      return m_->get_is_optimized(f);
-    }
-  };
-
-  typedef boost::counting_iterator<FloatIndex, 
-    std::forward_iterator_tag, unsigned int> FloatIndexIterator;
-
-public:
-
-
-  //! Add particle float attribute (assumed differentiable) to the model.
-  /** The returned index can be used for obtaining and setting the attribute
-      value.
-      \param[in] value Initial value of the attribute.
-      \return index of a new float attribute.
-   */
-  FloatIndex add_float(const Float value);
-
-  //! Set particle attribute value.
-  /** \param[in] idx Index of the attribute.
-      \param[in] value Value the attribute should be given.
-   */
-  void set_value(const FloatIndex idx, const Float value);
-
-  //! Get particle float attribute (inline).
-  /** \param[in] idx Index of the particle float attribute.
-      \return value of particle float attribute.
-   */
-  Float get_value(const FloatIndex idx) const {
-    IMP_assert(idx.get_index() < float_data_.size(),
-               "Out of range index requested");
-    return float_data_[idx.get_index()].value_;
-  }
-
-  //! Get derivative of the given particle float attribute.
-  /** \param[in] idx Index of the particle float attribute.
-   */
-  Float get_deriv(const FloatIndex idx) const;
-
-  //! Indicate if the particle float attribute is to be optimized.
-  /** \param[in] idx Index of the particle float attribute.
-      \return True if particle float attribute is to be optimized.
-   */
-  bool get_is_optimized(const FloatIndex idx) const;
-
-  //! Set whether the particle float attribute is to be optimized.
-  /** \param[in] idx Index of the particle float attribute.
-      \param[in] is_optimized True if particle float attribute is to be
-                              optimized.
-   */
-  void set_is_optimized(const FloatIndex idx, bool is_optimized);
-
-  //! Set all derivatives to zero.
-  void zero_derivatives();
-
-  //! Add particle int attribute to the model.
-  /** The returned index can be used for obtaining and setting the
-      attribute value.
-      \param[in] value Initial value of the attribute.
-      \return index of a new int attribute.
-   */
-  IntIndex add_int(const Int value);
-
-  //! Set particle attribute value.
-  /** \param[in] idx Index of the attribute.
-      \param[in] value Value the attribute should be given.
-   */
-  void set_value(const IntIndex idx, const Int value);
-
-  //! Get particle int attribute (inline).
-  /** \param[in] idx Index of the particle int attribute.
-      \return value of particle float attribute.
-   */
-  Int get_value(const IntIndex idx) const {
-    IMP_assert(idx.get_index() < int_data_.size(),
-               "Out of range int requested");
-    return int_data_[idx.get_index()];
-  }
-
-  //! Add particle string attribute to the model.
-  /** The returned index can be used for obtaining and setting the attribute
-      value.
-      \param[in] value Initial value of the attribute.
-      \return index of a new string attribute.
-   */
-  StringIndex add_string(const String value);
-
-  //! Set particle attribute value.
-  /** \param[in] idx Index of the attribute.
-      \param[in] value Value the attribute should be given.
-   */
-  void set_value(const StringIndex idx, const String value);
-
-  //! Get particle string attribute (inline).
-  /** \param[in] idx Index of the particle string attribute.
-      \return value of particle string attribute.
-   */
-  String get_value(const StringIndex idx) const {
-   IMP_assert(idx.get_index() < string_data_.size(),
-               "Out of range string requested");
-    return string_data_[idx.get_index()];
-  }
-
-  typedef boost::filter_iterator<FloatIsOptimized,
-    FloatIndexIterator > OptimizedFloatIndexIterator;
-  OptimizedFloatIndexIterator optimized_float_indexes_begin() const {
-    return OptimizedFloatIndexIterator(FloatIsOptimized(this),
-                                       FloatIndexIterator(0),
-                                       FloatIndexIterator(float_data_.size()));
-  }
-  OptimizedFloatIndexIterator optimized_float_indexes_end() const {
-    return OptimizedFloatIndexIterator(FloatIsOptimized(this),
-                                       FloatIndexIterator(float_data_.size()),
-                                       FloatIndexIterator(float_data_.size()));
-  }
-
-
-  void show(std::ostream &out=std::cout) const;
-protected:
-  ModelData();
-  ~ModelData();
-
-  //! Add value to derivative (used by DerivativeAccumulator).
-  /** \param[in] idx Index of the particle float attribute.
-      \param[in] value Value to add to the float attribute derivative.
-   */
-  void add_to_deriv(const FloatIndex idx, const Float value);
-
-  //! particle variables and attributes
-  /** these are stored outside of particles to allow
-      restraints to get access them directly through
-      indexes rather than through particle dereferencing.
-   */
-  std::vector<FloatData> float_data_;
-
-  //! See float_data_.
-  std::vector<Int> int_data_;
-
-  //! See float_data_.
-  std::vector<String> string_data_;
-};
-
-IMP_OUTPUT_OPERATOR(ModelData);
-
-} // namespace IMP
-
-#endif  /* __IMP_MODEL_DATA_H */
Index: kernel/include/IMP/base_types.h
===================================================================
--- kernel/include/IMP/base_types.h	(revision 434)
+++ kernel/include/IMP/base_types.h	(working copy)
@@ -54,9 +54,6 @@
 struct ScoreStateTag {};
 struct OptimizerStateTag {};
 
-typedef Index<Int> IntIndex;
-typedef Index<Float> FloatIndex;
-typedef Index<String> StringIndex;
 typedef Index<ParticleTag> ParticleIndex;
 typedef Index<RestraintTag> RestraintIndex;
 typedef Index<ScoreStateTag> ScoreStateIndex;
@@ -73,13 +70,6 @@
 
 
 typedef std::vector<ParticleIndex> ParticleIndexes;
-typedef std::vector<FloatIndex> FloatIndexes;
-typedef std::vector<IntIndex> IntIndexes;
-typedef std::vector<StringIndex> StringIndexes;
-// typedefs for the particle variable and attribute indexes
-// typedef DataIndex<Float> FloatIndex;
-// typedef DataIndex<Int> IntIndex;
-// typedef DataIndex<String> StringIndex;
 
 template <class T> class Key;
 
Index: kernel/include/IMP/Model.h
===================================================================
--- kernel/include/IMP/Model.h	(revision 434)
+++ kernel/include/IMP/Model.h	(working copy)
@@ -20,7 +20,6 @@
 
 class Particle;
 class Restraint;
-class ModelData;
 class ScoreState;
 typedef std::vector<Restraint*> Restraints;
 typedef std::vector<ScoreState*> ScoreStates;
@@ -38,13 +37,6 @@
   Model();
   ~Model();
 
-  //! Get pointer to all model particle data.
-  /** \return pointer to all model particle data.
-   */
-  ModelData* get_model_data() const {
-    return model_data_.get();
-  }
-
   IMP_CONTAINER(Particle, particle, ParticleIndex);
   IMP_CONTAINER(ScoreState, score_state, ScoreStateIndex);
   IMP_CONTAINER(Restraint, restraint, RestraintIndex);
@@ -66,10 +58,6 @@
   VersionInfo get_version_info() const {
     return internal::kernel_version_info;
   }
-
-protected:
-  //! all of the data associated with the particles
-  std::auto_ptr<ModelData> model_data_;
 };
 
 
Index: kernel/include/IMP/optimizers/MolecularDynamics.h
===================================================================
--- kernel/include/IMP/optimizers/MolecularDynamics.h	(revision 434)
+++ kernel/include/IMP/optimizers/MolecularDynamics.h	(working copy)
@@ -21,6 +21,9 @@
     and a non-optimizable mass attribute; this optimizer assumes the score
     to be energy in kcal/mol, the xyz coordinates to be in angstroms, and
     the mass to be in AMU (g/mol).
+
+    Particles without optimized x,y,z and nonoptimized mass are skipped.
+
     \ingroup optimizer
  */
 class IMPDLLEXPORT MolecularDynamics : public Optimizer
@@ -34,30 +37,26 @@
   //! Set time step in fs
   void set_time_step(Float t) { time_step_ = t; }
 
+  IMP_LIST(private, Particle, particle, Particle*);
+
 protected:
   //! Perform a single dynamics step.
   virtual void step();
 
   //! Get the set of particles to use in this optimization.
-  /** Populates particles_, and gives each particle velocity attributes if it
-      does not already have them.
-      \param[in] model The model to optimize.
-      \exception InvalidStateException The model does not contain only
-                                       xyz particles.
+  /** Scans for particles which have the necessary attributes to be optimized. Particles
+      without optimized x,y,z and nonoptimized mass are skipped.
    */
-  void setup_particles(Model& model);
+  void setup_particles();
 
   //! Time step in fs
   Float time_step_;
 
   //! Keys of the xyz coordinates and mass
-  FloatKey xkey_, ykey_, zkey_, masskey_;
+  FloatKey cs_[3], masskey_;
 
   //! Keys of the xyz velocities
-  FloatKey vxkey_, vykey_, vzkey_;
-
-  //! Particles to optimize
-  std::vector<Particle *> particles_;
+  FloatKey vs_[3];
 };
 
 } // namespace IMP
Index: kernel/include/IMP/optimizers/ConjugateGradients.h
===================================================================
--- kernel/include/IMP/optimizers/ConjugateGradients.h	(revision 434)
+++ kernel/include/IMP/optimizers/ConjugateGradients.h	(working copy)
@@ -31,6 +31,14 @@
   //! Set the threshold for the minimum gradient
   void set_threshold(Float t){ threshold_=t;}
 private:
+
+  Float get_score(std::vector<FloatIndex> float_indices,
+                  std::vector<Float> &x, std::vector<Float> &dscore);
+  bool line_search(std::vector<Float> &x, std::vector<Float> &dx,
+                   float &alpha, const std::vector<FloatIndex> &float_indices,
+                   int &ifun, float &f, float &dg, float &dg1,
+                   int max_steps, const std::vector<Float> &search,
+                   const std::vector<Float> &estimate);
   Float threshold_;
 };
 
Index: kernel/include/IMP/internal/AttributeTable.h
===================================================================
--- kernel/include/IMP/internal/AttributeTable.h	(revision 434)
+++ kernel/include/IMP/internal/AttributeTable.h	(working copy)
@@ -11,26 +11,33 @@
 #include "../base_types.h"
 #include "../utility.h"
 #include "../log.h"
-#include "../ModelData.h"
 
-#include <map>
+#include <boost/iterator/filter_iterator.hpp>
+#include <boost/iterator/counting_iterator.hpp>
+
 #include <vector>
 
 namespace IMP
 {
 
-class ModelData;
 
 namespace internal
 {
 
 /** \internal */
-template <class T>
+template <class T, class VT>
 class AttributeTable
 {
-  std::vector<Index<T> > map_;
+  struct Bin {
+    bool first;
+    VT second;
+    Bin(): first(false){}
+  };
+  typedef AttributeTable<T, VT> This;
+  typedef std::vector<Bin > Map;
+  Map map_;
 public:
-  typedef Index<T> Value;
+  typedef VT Value;
   typedef Key<T> Key;
   AttributeTable() {} 
   const Value get_value(Key k) const {
@@ -39,73 +46,99 @@
               << "\" not found in table.",
               IndexException((std::string("Invalid attribute \"")
                               + k.get_string() + "\" requested").c_str()));
-    return map_[k.get_index()];
+    return map_[k.get_index()].second;
   }
+
+  Value& get_value(Key k) {
+    IMP_check(contains(k),
+              "Attribute \"" << k.get_string()
+              << "\" not found in table.",
+              IndexException((std::string("Invalid attribute \"")
+                              + k.get_string() + "\" requested").c_str()));
+    return map_[k.get_index()].second;
+  }
   void insert(Key k, Value v);
   bool contains(Key k) const {
     IMP_check(k != Key(), "Can't search for default key",
               IndexException("Bad index"));
     return k.get_index() < map_.size()
-           && map_[k.get_index()] != Value();
+           && map_[k.get_index()].first;
   }
-  void show(std::ostream &out, const char *prefix="",
-            ModelData* md=NULL) const;
+  void show(std::ostream &out, const char *prefix="") const;
   std::vector<Key> get_keys() const;
+
+  class IsAttribute{
+    const This *map_;
+  public:
+    IsAttribute(): map_(NULL){}
+    IsAttribute(const This *map): map_(map) {}
+    bool operator()(Key k) const {
+      return map_->contains(k);
+    }
+  };
+  
+  typedef boost::counting_iterator<Key, boost::random_access_traversal_tag, std::size_t>
+  KeyIterator;
+  typedef boost::filter_iterator<IsAttribute, KeyIterator> AttributeKeyIterator;
+  
+  AttributeKeyIterator attribute_keys_begin() const {
+    return AttributeKeyIterator(IsAttribute(this),
+                                KeyIterator(Key(0U)),
+                                KeyIterator(Key(map_.size())));
+  }
+  AttributeKeyIterator attribute_keys_end() const {
+    return AttributeKeyIterator(IsAttribute(this),
+                                KeyIterator(Key(map_.size())),
+                                KeyIterator(Key(map_.size())));
+  }
+  
+
 };
 
-IMP_OUTPUT_OPERATOR_1(AttributeTable)
+IMP_OUTPUT_OPERATOR_2(AttributeTable)
 
 
 
 
-template <class T>
-inline void AttributeTable<T>::insert(Key k, Value v)
+template <class T, class VT>
+inline void AttributeTable<T, VT>::insert(Key k, Value v)
 {
   IMP_check(k != Key(),
             "Can't insert default key",
             IndexException("bad index"));
-  IMP_assert(v != Value(),
-             "Can't add attribute with no index");
   if (map_.size() <= k.get_index()) {
     map_.resize(k.get_index()+1);
   }
-  IMP_assert(map_[k.get_index()]== Value(),
+  IMP_assert(!map_[k.get_index()].first,
              "Trying to add attribute \"" << k.get_string()
              << "\" twice");
-  map_[k.get_index()]= v;
+  map_[k.get_index()].second= v;
+  map_[k.get_index()].first= true;
   IMP_assert(contains(k), "Something is broken");
 }
 
 
-template <class T>
-inline void AttributeTable<T>::show(std::ostream &out,
-                                    const char *prefix,
-                                    ModelData *md) const
+  template <class T, class VT>
+  inline void AttributeTable<T, VT>::show(std::ostream &out,
+                                          const char *prefix) const
 {
   for (unsigned int i=0; i< map_.size(); ++i) {
-    if (map_[i] != Value()) {
+    if (map_[i].first) {
       out << prefix;
       out << Key(i).get_string() << ": ";
-      if (md != NULL) {
-        out << md->get_value(map_[i]);
-      }
+      out << map_[i].second;
       out << std::endl;
     }
   }
 }
 
 
-template <class T>
-inline std::vector<typename AttributeTable<T>::Key> 
-AttributeTable<T>::get_keys() const
+template <class T, class VT>
+inline std::vector<typename AttributeTable<T, VT>::Key> 
+  AttributeTable<T, VT>::get_keys() const
 {
-  std::vector<Key> ret;
-   for (unsigned int i=0; i< map_.size(); ++i) {
-     if (map_[i] != Value()) {
-       ret.push_back(Key(i));
-     }
-   }
-   return ret;
+  std::vector<Key> ret(attribute_keys_begin(), attribute_keys_end());
+  return ret;
 }
 
 inline void show_attributes(std::ostream &out)
Index: kernel/include/IMP/Key.h
===================================================================
--- kernel/include/IMP/Key.h	(revision 434)
+++ kernel/include/IMP/Key.h	(working copy)
@@ -104,7 +104,7 @@
   };
 
   explicit Key(unsigned int i): str_(i) {
-    IMP_assert(data().rmap.size() > i, "There is no such attribute " << i);
+    //IMP_assert(data().rmap.size() > i, "There is no such attribute " << i);
   }
 
   //! Turn a key into a pretty string
@@ -143,6 +143,24 @@
   static unsigned int get_number_unique() {
     return data().map.size();
   }
+
+#ifndef SWIG
+  /** \todo These should be protected, I'll try to work how
+   */
+  This operator++() {
+    ++str_;
+    return *this;
+  }
+  This operator--() {
+    --str_;
+    return *this;
+  }
+  This operator+(int o) const {
+    This c=*this;
+    c.str_+= o;
+    return c;
+  }
+#endif
 };
 
 IMP_OUTPUT_OPERATOR_1(Key)
Index: kernel/include/IMP/SConscript
===================================================================
--- kernel/include/IMP/SConscript	(revision 434)
+++ kernel/include/IMP/SConscript	(working copy)
@@ -3,7 +3,7 @@
 
 files = ['base_types.h', 'random.h', 'Index.h', 'Model.h',
          'Particle.h', 'ScoreState.h', 'OptimizerState.h', 'IMP_config.h',
-         'log.h', 'DerivativeAccumulator.h', 'ModelData.h',
+         'log.h', 'DerivativeAccumulator.h', 
          'Key.h', 'utility.h', 'Restraint.h', 'Optimizer.h',
          'DecoratorBase.h', 'Vector3D.h', 'ScoreFuncParams.h',
          'UnaryFunction.h', 'PairScore.h', 'SingletonScore.h', 'macros.h',
Index: kernel/include/IMP/Particle.h
===================================================================
--- kernel/include/IMP/Particle.h	(revision 434)
+++ kernel/include/IMP/Particle.h	(working copy)
@@ -8,8 +8,6 @@
 #ifndef __IMP_PARTICLE_H
 #define __IMP_PARTICLE_H
 
-#include <map>
-
 #include "IMP_config.h"
 #include "base_types.h"
 #include "Model.h"
@@ -20,10 +18,30 @@
 #include "DerivativeAccumulator.h"
 #include "internal/ObjectPointer.h"
 
+#include <limits>
+
 namespace IMP
 {
-
-
+#ifndef SWIG
+namespace internal
+{
+  struct IMPDLLEXPORT FloatData {
+    bool is_optimized;
+    Float value;
+    Float derivative;
+    FloatData(Float v, bool opt): is_optimized(opt), value(v), derivative(0){}
+    FloatData():is_optimized(false), value(std::numeric_limits<Float>::max()),
+                derivative(std::numeric_limits<Float>::max()){}
+                                           
+  };
+  inline std::ostream &operator<<(std::ostream &out, const FloatData &d) {
+    out << d.value << ": " << d.derivative;
+    if (d.is_optimized) out << "(opt)";
+    return out;
+  }
+  
+}
+#endif
 class Model;
 
 //! Class to handle individual model particles.
@@ -184,10 +202,9 @@
 
   //! Return a vector containing all the FloatKeys for the Particle
   /**
-     This is for use in python mostly. C++ iterator will be added
-     at some point.
+     This is for use in python mostly. C++ users should use the iterators.
 
-     I would like to have a type-agnostic way of calling this 
+     \todo I would like to have a type-agnostic way of calling this 
      to be used to writing generic functions in python. The only
      ways I can think of doing this are to pass dummy arguments,
      which seems inelegant.
@@ -206,22 +223,40 @@
     return string_indexes_.get_keys();
   }
 
-protected:
+  //! An iterator through the keys of the float attributes of this particle
+  typedef internal::AttributeTable<Float, internal::FloatData>::AttributeKeyIterator
+    FloatKeyIterator;
+  //! Iterate through the keys of float attributes of the particle
+  FloatKeyIterator float_keys_begin() const {
+    return float_indexes_.attribute_keys_begin();
+  }
+  FloatKeyIterator float_keys_end() const {
+    return float_indexes_.attribute_keys_end();
+  }
 
-  template <class T>
-  T get_value_t(Index<T> k) const {
-    IMP_check(get_is_active(), "get_value called on inactive Particle",
-              InactiveParticleException());
-    return model_->get_model_data()->get_value(k);
+  //! An iterator through the keys of the int attributes of this particle
+  typedef internal::AttributeTable<Int, Int>::AttributeKeyIterator IntKeyIterator;
+  //! Iterate through the keys of int attributes of the particle
+  IntKeyIterator int_keys_begin() const {
+    return int_indexes_.attribute_keys_begin();
   }
+  IntKeyIterator int_keys_end() const {
+    return int_indexes_.attribute_keys_end();
+  }
 
-  template <class T>
-  void set_value_t(Index<T> k, const T&v) {
-    IMP_check(get_is_active(), "set_value called on inactive Particle",
-              InactiveParticleException());
-    model_->get_model_data()->set_value(k, v);
+  //! An iterator through the keys of the string attributes of this particle
+  typedef internal::AttributeTable<String, String>::AttributeKeyIterator StringKeyIterator;
+  //! Iterate through the keys of string attributes of the particle
+  StringKeyIterator string_keys_begin() const {
+    return string_indexes_.attribute_keys_begin();
   }
+  StringKeyIterator string_keys_end() const {
+    return string_indexes_.attribute_keys_end();
+  }
 
+protected:
+  void zero_derivatives();
+
   // Set pointer to model particle data.
   void set_model(Model *md, ParticleIndex pi);
 
@@ -232,11 +267,11 @@
   bool is_active_;
 
   // float attributes associated with the particle
-  internal::AttributeTable<Float> float_indexes_;
+  internal::AttributeTable<Float, internal::FloatData> float_indexes_;
   // int attributes associated with the particle
-  internal::AttributeTable<Int>  int_indexes_;
+  internal::AttributeTable<Int, Int>  int_indexes_;
   // string attributes associated with the particle
-  internal::AttributeTable<String>  string_indexes_;
+  internal::AttributeTable<String, String>  string_indexes_;
 
   ParticleIndex pi_;
 };
@@ -255,40 +290,52 @@
 
 inline Float Particle::get_value(FloatKey name) const
 {
-  return get_value_t(float_indexes_.get_value(name));
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
+  return float_indexes_.get_value(name).value;
 }
 
 inline Float Particle::get_derivative(FloatKey name) const
 {
-  return model_->get_model_data()->get_deriv(float_indexes_.get_value(name));
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
+  return float_indexes_.get_value(name).derivative;
 }
 
 inline void Particle::set_value(FloatKey name, Float value)
 {
-  set_value_t(float_indexes_.get_value(name), value);
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
+  IMP_assert(value==value, "Can't set value to NaN");
+  float_indexes_.get_value(name).value= value;
 }
 
 inline bool Particle::get_is_optimized(FloatKey name) const
 {
-  return model_->get_model_data()
-    ->get_is_optimized(float_indexes_.get_value(name));
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
+  return float_indexes_.get_value(name).is_optimized;
 }
 
 inline void Particle::set_is_optimized(FloatKey name, bool tf)
 {
-  model_->get_model_data()
-    ->set_is_optimized(float_indexes_.get_value(name), tf);
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
+  float_indexes_.get_value(name).is_optimized=tf;
 }
 
 inline void Particle::add_to_derivative(FloatKey name, Float value,
                                         const DerivativeAccumulator &da)
 {
-  return model_->get_model_data()->add_to_deriv(float_indexes_.get_value(name),
-                                                da(value));
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
+  float_indexes_.get_value(name).derivative+= da(value);
 }
 
 inline bool Particle::has_attribute(IntKey name) const
 {
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
   return int_indexes_.contains(name);
 }
 
@@ -296,17 +343,23 @@
 
 inline Int Particle::get_value(IntKey name) const
 {
-  return get_value_t(int_indexes_.get_value(name));
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
+  return int_indexes_.get_value(name);
 }
 
 
 inline void Particle::set_value(IntKey name, Int value)
 {
-  return set_value_t(int_indexes_.get_value(name), value);
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
+  int_indexes_.get_value(name)= value;
 }
 
 inline bool Particle::has_attribute(StringKey name) const
 {
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
   return string_indexes_.contains(name);
 }
 
@@ -314,12 +367,16 @@
 
 inline String Particle::get_value(StringKey name) const
 {
-  return get_value_t(string_indexes_.get_value(name));
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
+  return string_indexes_.get_value(name);
 }
 
 inline void Particle::set_value(StringKey name, String value)
 {
-  return set_value_t(string_indexes_.get_value(name), value);
+  IMP_check(get_is_active(), "Do not touch inactive particles",
+            ValueException("Don't touch inactive particles"));
+  string_indexes_.get_value(name)= value;
 }
 
 } // namespace IMP
Index: kernel/include/IMP/Optimizer.h
===================================================================
--- kernel/include/IMP/Optimizer.h	(revision 434)
+++ kernel/include/IMP/Optimizer.h	(working copy)
@@ -16,6 +16,7 @@
 #include "internal/Object.h"
 #include "utility.h"
 #include "Model.h"
+#include "Particle.h"
 #include "internal/ObjectPointer.h"
 
 namespace IMP
@@ -68,11 +69,133 @@
   //! Update optimizer state, should be called at each successful step
   void update_states();
 
-private:
+  //! An index to an optimized particle attribute
+  class FloatIndexInderator;
+  class FloatIndex{
+    friend class Optimizer;
+    friend class FloatIndexIterator;
+    Model::ParticleIterator p_;
+    Particle::FloatKeyIterator fk_;
+    FloatIndex(Model::ParticleIterator p): p_(p){}
+  public:
+    FloatIndex() {}
+  };
+
+
+  //! An interator through the optimized attributes
+  /**
+     The value type is an FloatIndex
+   */
+  class FloatIndexIterator {
+    typedef FloatIndexIterator This;
+    Model::ParticleIterator pe_;
+    mutable FloatIndex i_;
+
+    void search_valid() const {
+      while (i_.fk_ == (*i_.p_)->float_keys_end()
+             || !(*i_.p_)->get_is_optimized(*i_.fk_)) {
+        if (i_.fk_ == (*i_.p_)->float_keys_end()) {
+          ++i_.p_;
+          if (i_.p_== pe_) return;
+          else {
+            i_.fk_= (*i_.p_)->float_keys_begin();
+          }
+        } else {
+          ++i_.fk_;
+        }
+      }
+      IMP_assert(i_.p_ != pe_, "Should have just returned");
+      IMP_assert(i_.fk_ != (*i_.p_)->float_keys_end(),
+                 "Broken iterator end");
+      IMP_assert((*i_.p_)->get_is_optimized(*i_.fk_),
+                   "Why did the loop end?");
+    }
+    void find_next() const {
+      ++i_.fk_;
+      search_valid();
+    }
+  public:
+    FloatIndexIterator(Model::ParticleIterator pc,
+                       Model::ParticleIterator pe): pe_(pe), i_(pc){
+      if (pc != pe) {
+        i_.fk_= (*pc)->float_keys_begin();
+        search_valid();
+      }
+    }
+    typedef FloatIndex value_type;
+    typedef FloatIndex& reference;
+    typedef FloatIndex* pointer;
+    typedef std::forward_iterator_tag iterator_category;
+    typedef int difference_type;
+
+    const This &operator++() {
+      find_next();
+      return *this;
+    }
+    reference operator*() const {
+      IMP_assert((*i_.p_)->get_is_optimized(*i_.fk_),
+                 "The iterator is broken");
+      return i_;
+    }
+    pointer operator->() const {
+      return &i_;
+    }
+    bool operator==(const This &o) const {
+      if (i_.p_ != o.i_.p_) return false;
+      if (i_.p_== pe_) return o.i_.p_ ==pe_;
+      else return i_.fk_ == o.i_.fk_;
+    }
+    bool operator!=(const This &o) const {
+      return !operator==(o);
+    }
+  };
+  
+  //! Iterate through the optimized attributes
+  FloatIndexIterator float_indexes_begin() {
+    return FloatIndexIterator(model_->particles_begin(),
+                              model_->particles_end());
+  }
+
+  FloatIndexIterator float_indexes_end() {
+    return FloatIndexIterator(model_->particles_end(),
+                                  model_->particles_end());
+  }
+  
+  //! Set the value of an optimized attribute
+  /**
+     The attribute must be optimized or an ErrorException is thrown.
+   */
+  void set_value(FloatIndex fi, Float v) {
+    IMP_assert(fi.p_ != model_->particles_end(), 
+               "Out of range FloatIndex in Optimizer");
+    IMP_assert((*fi.p_)->get_is_optimized(*fi.fk_),
+               "Keep your mits off unoptimized attributes "
+               << (*fi.p_)->get_index() << " " << *fi.fk_ << std::endl);
+    (*fi.p_)->set_value(*fi.fk_, v);
+  }
+
+  //! Get the value of an optimized attribute
+  Float get_value(FloatIndex fi) const {
+    IMP_assert(fi.p_ != model_->particles_end(), 
+               "Out of range FloatIndex in Optimizer");
+    return (*fi.p_)->get_value(*fi.fk_);
+  }
+
+  //! Get the derivative of an optimized attribute
+  Float get_derivative(FloatIndex fi) const {
+    IMP_assert(fi.p_ != model_->particles_end(), 
+               "Out of range FloatIndex in Optimizer");
+    return (*fi.p_)->get_derivative(*fi.fk_);
+  }
+  
+  typedef std::vector<FloatIndex> FloatIndexes;
+
+ private:
   internal::ObjectPointer<Model, false> model_;
-
 };
 
+  
+
 } // namespace IMP
 
 #endif  /* __IMP_OPTIMIZER_H */
Index: kernel/include/IMP/macros.h
===================================================================
--- kernel/include/IMP/macros.h	(revision 434)
+++ kernel/include/IMP/macros.h	(working copy)
@@ -124,6 +124,17 @@
     return out;                                                         \
   }
 
+//! Implement operator<< on class name, assuming it has two template arguments
+/** class name should also define the method std::ostream &show(std::ostream&)
+ */
+#define IMP_OUTPUT_OPERATOR_2(name) /** write to a stream*/             \
+  template <class L, class M>                                           \
+  inline std::ostream& operator<<(std::ostream &out, const name<L, M> &i) \
+  {                                                                     \
+    i.show(out);                                                        \
+    return out;                                                         \
+  }
+
 //! Implement operator<< on class name
 /** class name should also define the method std::ostream &show(std::ostream&)
  */
@@ -301,8 +312,8 @@
    \param[in] OnChanged Code to get executed when the container changes.
  */
 #define IMP_LIST_IMPL(Class, Ucname, lcname, Data, init, OnChanged)     \
-  IMP_CONTAINER_CORE_IMPL(Class, Ucname, lcname, Data, unsigned int,\
-                          init, OnChanged)\
+  IMP_CONTAINER_CORE_IMPL(Class, Ucname, lcname, Data, unsigned int,    \
+                          init, OnChanged)                              \
   /** \short Clear the contents of the container */                     \
   void Class::clear_##lcname##s(){                                      \
     lcname##_vector_.clear();                                           \
Index: kernel/src/Model.cpp
===================================================================
--- kernel/src/Model.cpp	(revision 434)
+++ kernel/src/Model.cpp	(working copy)
@@ -7,7 +7,6 @@
  */
 
 #include "IMP/Model.h"
-#include "IMP/ModelData.h"
 #include "IMP/Particle.h"
 #include "IMP/log.h"
 #include "IMP/Restraint.h"
@@ -18,7 +17,7 @@
 {
 
 //! Constructor
-Model::Model(): model_data_(new ModelData())
+Model::Model()
 {
 }
 
@@ -44,9 +43,13 @@
           "Model evaluate (" << number_of_restraints() << " restraints):"
           << std::endl);
   // If calcualting derivatives, first set all derivatives to zero
-  if (calc_derivs)
-    model_data_->zero_derivatives();
+  if (calc_derivs) {
+    for (ParticleIterator pit = particles_begin(); pit != particles_end(); ++pit) {
+      (*pit)->zero_derivatives();
+    }
 
+  }
+
   IMP_LOG(VERBOSE,
           "Updating ScoreStates " << std::flush);
   for (ScoreStateIterator it = score_states_begin(); it != score_states_end();
Index: kernel/src/Particle.cpp
===================================================================
--- kernel/src/Particle.cpp	(revision 434)
+++ kernel/src/Particle.cpp	(working copy)
@@ -40,10 +40,7 @@
 {
   IMP_assert(model_ ,
              "Particle must be added to Model before an attributes are added");
-  float_indexes_.insert(name, model_->get_model_data()->add_float(value));
-
-  model_->get_model_data()->set_is_optimized(float_indexes_.get_value(name),
-                                             is_optimized);
+  float_indexes_.insert(name, internal::FloatData(value, is_optimized));
 }
 
 
@@ -51,7 +48,7 @@
 {
   IMP_assert(model_,
              "Particle must be added to Model before an attributes are added");
-  int_indexes_.insert(name, model_->get_model_data()->add_int(value));
+  int_indexes_.insert(name, value);
 }
 
 
@@ -59,11 +56,18 @@
 {
   IMP_assert(model_,
              "Particle must be added to Model before an attributes are added");
-  string_indexes_.insert(name, model_->get_model_data()->add_string(value));
+  string_indexes_.insert(name, value);
 }
 
 
+void Particle::zero_derivatives()
+{
+  for (FloatKeyIterator it= float_keys_begin(); it != float_keys_end(); ++it) {
+    float_indexes_.get_value(*it).derivative=0;
+  }
+}
 
+
 void Particle::show(std::ostream& out) const
 {
   char* inset = "  ";
@@ -78,13 +82,13 @@
 
   if (get_model() != NULL) {
     out << inset << inset << "float attributes:" << std::endl;
-    float_indexes_.show(out, "    ", get_model()->get_model_data());
+    float_indexes_.show(out, "    ");
 
     out << inset << inset << "int attributes:" << std::endl;
-    int_indexes_.show(out, "    ", get_model()->get_model_data());
+    int_indexes_.show(out, "    ");
 
     out << inset << inset << "string attributes:" << std::endl;
-    string_indexes_.show(out, "    ", get_model()->get_model_data());
+    string_indexes_.show(out, "    ");
   }
 }
 
Index: kernel/src/SConscript
===================================================================
--- kernel/src/SConscript	(revision 434)
+++ kernel/src/SConscript	(working copy)
@@ -17,7 +17,7 @@
 internal_files = SConscript('internal/SConscript')
 
 # Source files
-files = ['base_types.cpp', 'Model.cpp', 'ModelData.cpp',
+files = ['base_types.cpp', 'Model.cpp',
          'Particle.cpp', 'ScoreState.cpp',
          'OptimizerState.cpp', 'Log.cpp', 'Restraint.cpp', 'Optimizer.cpp',
          'BasicScoreFuncParams.cpp', 'random.cpp'
Index: kernel/src/Optimizer.cpp
===================================================================
--- kernel/src/Optimizer.cpp	(revision 434)
+++ kernel/src/Optimizer.cpp	(working copy)
@@ -12,20 +12,18 @@
 namespace IMP
 {
 
-//! Constructor
 Optimizer::Optimizer()
 {
   IMP_LOG(VERBOSE, "MEMORY: Optimizer created " << this << std::endl);
 }
 
 
-//! Destructor
 Optimizer::~Optimizer()
 {
   IMP_LOG(VERBOSE, "MEMORY: Optimizer destroyed " << this << std::endl);
 }
 
-//! Update optimizer state, at each successful step
+
 void Optimizer::update_states()
 {
   IMP_LOG(VERBOSE,
@@ -39,6 +37,7 @@
   IMP_LOG(VERBOSE, "done." << std::endl);
 }
 
+
 IMP_CONTAINER_IMPL(Optimizer, OptimizerState, optimizer_state,
                    OptimizerStateIndex, obj->set_optimizer(this));
 
Index: kernel/src/optimizers/ConjugateGradients.cpp
===================================================================
--- kernel/src/optimizers/ConjugateGradients.cpp	(revision 434)
+++ kernel/src/optimizers/ConjugateGradients.cpp	(working copy)
@@ -5,13 +5,12 @@
  *
  */
 
-#include <cmath>
-
 #include "IMP/log.h"
 #include "IMP/optimizers/ConjugateGradients.h"
 #include "IMP/Model.h"
-#include "IMP/ModelData.h"
+
 #include <limits>
+#include <cmath>
 
 namespace IMP
 {
@@ -27,22 +26,21 @@
     \param[out] dscore First derivatives for current state.
     \return The model score.
  */
-static Float get_score(Model &model, ModelData *model_data,
-                       std::vector<FloatIndex> float_indices,
-                       std::vector<Float> &x, std::vector<Float> &dscore)
+Float ConjugateGradients::get_score(std::vector<FloatIndex> float_indices,
+                                    std::vector<Float> &x, std::vector<Float> &dscore)
 {
   int i, opt_var_cnt = float_indices.size();
   /* set model state */
   for (i = 0; i < opt_var_cnt; i++) {
-    model_data->set_value(float_indices[i], x[i]);
+    set_value(float_indices[i], x[i]);
   }
 
   /* get score */
-  Float score = model.evaluate(true);
+  Float score = get_model()->evaluate(true);
 
   /* get derivatives */
   for (i = 0; i < opt_var_cnt; i++) {
-    dscore[i] = model_data->get_deriv(float_indices[i]);
+    dscore[i] = get_derivative(float_indices[i]);
   }
   return score;
 }
@@ -62,12 +60,11 @@
     \return true if the line search succeeded, false if max_steps was exceeded
             or a minimum could not be found.
  */
-static bool line_search(std::vector<Float> &x, std::vector<Float> &dx,
-                        float &alpha, Model &model, ModelData *model_data,
-                        const std::vector<FloatIndex> &float_indices,
-                        int &ifun, float &f, float &dg, float &dg1,
-                        int max_steps, const std::vector<Float> &search,
-                        const std::vector<Float> &estimate)
+bool ConjugateGradients::line_search(std::vector<Float> &x, std::vector<Float> &dx,
+                                     float &alpha, const std::vector<FloatIndex> &float_indices,
+                                     int &ifun, float &f, float &dg, float &dg1,
+                                     int max_steps, const std::vector<Float> &search,
+                                     const std::vector<Float> &estimate)
 {
   float ap, fp, dp, step, minf, u1, u2;
   int i, n, ncalls = ifun;
@@ -106,7 +103,7 @@
     }
 
     /* EVALUATE THE FUNCTION AT THE TRIAL POINT. */
-    f = get_score(model, model_data, float_indices, x, dx);
+    f = get_score(float_indices, x, dx);
 
     /* TEST IF THE MAXIMUM NUMBER OF FUNCTION CALLS HAVE BEEN USED. */
     if (++ifun > max_steps) {
@@ -213,17 +210,17 @@
 {
   std::vector<Float> x, dx;
   int i;
-  ModelData* model_data = get_model()->get_model_data();
+  //ModelData* model_data = get_model()->get_model_data();
 
-  FloatIndexes float_indices(model_data->optimized_float_indexes_begin(),
-                             model_data->optimized_float_indexes_end());
+  FloatIndexes float_indices(float_indexes_begin(),
+                             float_indexes_end());
   int n = float_indices.size();
 
   x.resize(n);
   dx.resize(n);
   // get initial state in x(n):
   for (i = 0; i < n; i++) {
-    x[i] = model_data->get_value(float_indices[i]);
+    x[i] = get_value(float_indices[i]);
   }
 
   // Initialize optimization variables
@@ -251,7 +248,7 @@
      whether a Beale restart is being done. nrst=n means that this
      iteration is a restart iteration. */
 g20:
-  f = get_score(*get_model(), model_data, float_indices, x, dx);
+  f = get_score(float_indices, x, dx);
   ifun++;
   nrst = n;
   // this is a gradient, not restart, direction:
@@ -295,7 +292,7 @@
   destimate = dx;
 
   /* Try to find a better score by linear search */
-  if (!line_search(x, dx, alpha, *get_model(), model_data, float_indices,
+  if (!line_search(x, dx, alpha, float_indices,
                    ifun, f, dg, dg1, max_steps, search, estimate)) {
     /* If the line search failed, it was either because the maximum number
        of iterations was exceeded, or the minimum could not be found */
@@ -417,7 +414,7 @@
 end:
   // If the 'best current estimate' is better than the current state, return
   // that:
-  bestf = get_score(*get_model(), model_data, float_indices, estimate,
+  bestf = get_score(float_indices, estimate,
                     destimate);
   if (bestf < f) {
     for (i = 0; i < n; i++) {
@@ -427,7 +424,7 @@
   }
   // Set final model state
   for (i = 0; i < n; i++) {
-    model_data->set_value(float_indices[i], x[i]);
+    set_value(float_indices[i], x[i]);
   }
   update_states();
   return f;
Index: kernel/src/optimizers/SteepestDescent.cpp
===================================================================
--- kernel/src/optimizers/SteepestDescent.cpp	(revision 434)
+++ kernel/src/optimizers/SteepestDescent.cpp	(working copy)
@@ -8,7 +8,6 @@
 #include "IMP/log.h"
 #include "IMP/optimizers/SteepestDescent.h"
 #include "IMP/Model.h"
-#include "IMP/ModelData.h"
 
 namespace IMP
 {
@@ -39,13 +38,12 @@
   std::vector<Float> temp_vals;
   std::vector<Float> temp_derivs;
   Float last_score, new_score = 0.0;
-  ModelData* model_data = get_model()->get_model_data();
 
   // set up the indexes
 
 
-  FloatIndexes float_indexes(model_data->optimized_float_indexes_begin(),
-                             model_data->optimized_float_indexes_end());
+  FloatIndexes float_indexes(float_indexes_begin(),
+                             float_indexes_end());
   int opt_var_cnt = float_indexes.size();
 
   Float current_step_size = step_size_;
@@ -67,8 +65,8 @@
     for (int i = 0; i < opt_var_cnt; i++) {
       FloatIndex fi = float_indexes[0];
 
-      temp_vals[i] = model_data->get_value(float_indexes[i]);
-      temp_derivs[i] = model_data->get_deriv(float_indexes[i]);
+      temp_vals[i] = get_value(float_indexes[i]);
+      temp_derivs[i] = get_derivative(float_indexes[i]);
     }
 
     bool done = false;
@@ -86,7 +84,7 @@
                 << temp_vals[i] - temp_derivs[i] * current_step_size << "  "
                 << temp_derivs[i]);
 
-        model_data->set_value(float_indexes[i],
+        set_value(float_indexes[i],
                               temp_vals[i] - temp_derivs[i]
                                              * current_step_size);
       }
Index: kernel/src/optimizers/MolecularDynamics.cpp
===================================================================
--- kernel/src/optimizers/MolecularDynamics.cpp	(revision 434)
+++ kernel/src/optimizers/MolecularDynamics.cpp	(working copy)
@@ -5,17 +5,25 @@
  *
  */
 
-#include <cmath>
-
 #include "IMP/log.h"
 #include "IMP/optimizers/MolecularDynamics.h"
+#include "IMP/decorators/XYZDecorator.h"
 
+#include <cmath>
+
 namespace IMP
 {
 
 //! Constructor
 MolecularDynamics::MolecularDynamics() : time_step_(4.0)
 {
+  cs_[0] = FloatKey("x");
+  cs_[1] = FloatKey("y");
+  cs_[2] = FloatKey("z");
+  masskey_ = FloatKey("mass");
+  vs_[0] = FloatKey("vx");
+  vs_[1] = FloatKey("vy");
+  vs_[2] = FloatKey("vz");
 }
 
 
@@ -24,51 +32,36 @@
 {
 }
 
+IMP_LIST_IMPL(MolecularDynamics, Particle, particle, Particle*,
+              {
+                if (0) std::cout << index;
+                for (unsigned int i=0; i< 3; ++i) {
+                  if (!obj->has_attribute(vs_[i])) {
+                    obj->add_attribute(vs_[i], 0.0, false);
+                  }
+                }
+              },);
 
+
 //! Get the set of particles to use in this optimization.
 /** Populates particles_, and gives each particle a velocity.
     \param[in] model The model to optimize.
     \exception InvalidStateException The model does not contain only
                                      xyz particles.
  */
-void MolecularDynamics::setup_particles(Model& model)
+void MolecularDynamics::setup_particles()
 {
-  xkey_ = FloatKey("x");
-  ykey_ = FloatKey("y");
-  zkey_ = FloatKey("z");
-  masskey_ = FloatKey("mass");
+  clear_particles();
 
-  int npart = model.number_of_particles();
-  int i;
-  for (i = 0; i < npart; ++i) {
-    Particle *p = model.get_particle(i);
-    if (p->has_attribute(xkey_) && p->get_is_optimized(xkey_)
-        && p->has_attribute(ykey_) && p->get_is_optimized(ykey_)
-        && p->has_attribute(zkey_) && p->get_is_optimized(zkey_)
+  for (unsigned int i = 0; i < get_model()->number_of_particles(); ++i) {
+    Particle *p = get_model()->get_particle(i);
+    if (p->has_attribute(cs_[0]) && p->get_is_optimized(cs_[2])
+        && p->has_attribute(cs_[1]) && p->get_is_optimized(cs_[2])
+        && p->has_attribute(cs_[2]) && p->get_is_optimized(cs_[2])
         && p->has_attribute(masskey_) && !p->get_is_optimized(masskey_)) {
-      particles_.push_back(p);
+      add_particle(p);
     }
   }
-  ModelData *md= model.get_model_data();
-  unsigned nopt = std::distance(md->optimized_float_indexes_begin(),
-                                md->optimized_float_indexes_end());
-
-  if (particles_.size() * 3 != nopt) {
-    throw InvalidStateException("Can only do MD on xyz particles with mass");
-  }
-
-  vxkey_ = FloatKey("vx");
-  vykey_ = FloatKey("vy");
-  vzkey_ = FloatKey("vz");
-  FloatKey *derivs[3] = { &vxkey_, &vykey_, &vzkey_ };
-  for (i = 0; i < npart; ++i) {
-    Particle *p = model.get_particle(i);
-    for (int dind = 0; dind < 3; ++dind) {
-      if (!p->has_attribute(*derivs[dind])) {
-        p->add_attribute(*derivs[dind], 0.0, false);
-      }
-    }
-  }
 }
 
 
@@ -80,27 +73,27 @@
   // in angstrom/fs/fs from raw derivatives
   static const Float deriv_to_acceleration = -4.1868e-4;
 
-  for (std::vector<Particle *>::iterator iter = particles_.begin();
-       iter != particles_.end(); ++iter) {
+  for (ParticleIterator iter = particles_begin();
+       iter != particles_end(); ++iter) {
     Particle *p = *iter;
-    Float x = p->get_value(xkey_);
-    Float y = p->get_value(ykey_);
-    Float z = p->get_value(zkey_);
+    Float x = p->get_value(cs_[0]);
+    Float y = p->get_value(cs_[1]);
+    Float z = p->get_value(cs_[2]);
     Float invmass = 1.0 / p->get_value(masskey_);
-    Float dvx = p->get_derivative(xkey_);
-    Float dvy = p->get_derivative(ykey_);
-    Float dvz = p->get_derivative(zkey_);
+    Float dvx = p->get_derivative(cs_[0]);
+    Float dvy = p->get_derivative(cs_[1]);
+    Float dvz = p->get_derivative(cs_[2]);
 
     // calculate velocity at t+(delta t/2) from that at t-(delta t/2)
-    Float vx = p->get_value(vxkey_);
-    Float vy = p->get_value(vykey_);
-    Float vz = p->get_value(vzkey_);
+    Float vx = p->get_value(vs_[0]);
+    Float vy = p->get_value(vs_[1]);
+    Float vz = p->get_value(vs_[2]);
     vx += dvx * deriv_to_acceleration * invmass * time_step_;
     vy += dvy * deriv_to_acceleration * invmass * time_step_;
     vz += dvz * deriv_to_acceleration * invmass * time_step_;
-    p->set_value(vxkey_, vx);
-    p->set_value(vykey_, vy);
-    p->set_value(vzkey_, vz);
+    p->set_value(vs_[0], vx);
+    p->set_value(vs_[1], vy);
+    p->set_value(vs_[2], vz);
 
     // get atomic shift
     Float dx = vx * time_step_;
@@ -111,9 +104,9 @@
     x += dx;
     y += dy;
     z += dz;
-    p->set_value(xkey_, x);
-    p->set_value(ykey_, y);
-    p->set_value(zkey_, z);
+    p->set_value(cs_[0], x);
+    p->set_value(cs_[1], y);
+    p->set_value(cs_[2], z);
   }
 }
 
@@ -125,7 +118,7 @@
 Float MolecularDynamics::optimize(unsigned int max_steps)
 {
   Model *model = get_model();
-  setup_particles(*model);
+  setup_particles();
 
   // get initial system score
   Float score = model->evaluate(true);
Index: kernel/src/ModelData.cpp
===================================================================
--- kernel/src/ModelData.cpp	(revision 434)
+++ kernel/src/ModelData.cpp	(working copy)
@@ -1,179 +0,0 @@
-/**
- *  \file ModelData.cpp  \brief Storage for all model particle data.
- *
- *  Copyright 2007-8 Sali Lab. All rights reserved.
- *
- */
-
-#include "IMP/ModelData.h"
-#include "IMP/log.h"
-
-namespace IMP
-{
-
-//! Constructor
-ModelData::ModelData()
-{
-}
-
-
-//! Destructor
-ModelData::~ModelData()
-{
-}
-
-
-//! Add particle float attribute (assumed differentiable) to the model.
-/** The returned index can be used for obtaining and setting the attribute
-    value.
-    \param[in] value Initial value of the attribute.
-    \return index of a new float attribute.
- */
-FloatIndex ModelData::add_float(const Float value)
-{
-  int size = float_data_.size();
-
-  float_data_.resize(size + 1);
-
-  float_data_[size].value_ = value;
-
-  return FloatIndex(size);
-}
-
-
-//! Set particle attribute value.
-/** \param[in] idx Index of the attribute.
-    \param[in] value Value the attribute should be given.
- */
-void ModelData::set_value(const FloatIndex idx, const Float value)
-{
-  IMP_assert(idx.get_index() < float_data_.size(),
-             "Out of range float requested");
-  float_data_[idx.get_index()].value_ = value;
-}
-
-
-
-//! Add value to derivative.
-/** \param[in] idx Index of the particle float attribute.
-    \param[in] value Value to add to the float attribute derivative.
- */
-void ModelData::add_to_deriv(const FloatIndex idx, Float value)
-{
-  IMP_assert(idx.get_index() < float_data_.size(),
-             "Out of range float requested");
-  float_data_[idx.get_index()].deriv_ += value;
-}
-
-
-//! Get derivative of the given particle float attribute.
-/** \param[in] idx Index of the particle float attribute.
- */
-Float ModelData::get_deriv(const FloatIndex idx) const
-{
-  IMP_assert(idx.get_index() < float_data_.size(),
-             "Out of range float requested");
-  return float_data_[idx.get_index()].deriv_;
-}
-
-
-//! Indicate if the particle float attribute is to be optimized.
-/** \return True if particle float attribute is to be optimized.
- */
-bool ModelData::get_is_optimized(const FloatIndex idx) const
-{
-  IMP_assert(idx.get_index() < float_data_.size(),
-             "Out of range float requested");
-  return float_data_[idx.get_index()].is_optimized_;
-}
-
-
-//! Set whether the particle float attribute is to be optimized.
-/** \param[in] is_optimized True if particle float attribute is to be optimized.
- */
-void ModelData::set_is_optimized(const FloatIndex idx, bool is_optimized)
-{
-  IMP_assert(idx.get_index() < float_data_.size(),
-             "Out of range float requested");
-  float_data_[idx.get_index()].is_optimized_ = is_optimized;
-}
-
-
-//! Set all derivatives to zero.
-void ModelData::zero_derivatives()
-{
-  for (size_t i = 0; i < float_data_.size(); i++) {
-    float_data_[i].deriv_ = (Float) 0.0;
-  }
-}
-
-
-//! Add particle int attribute to the model.
-/** The returned index can be used for obtaining and setting the
-    attribute value.
-
-    \param[in] value Initial value of the attribute.
-    \return index of a new int attribute.
- */
-IntIndex ModelData::add_int(const Int value)
-{
-  int size = int_data_.size();
-
-  int_data_.resize(size + 1);
-  int_data_[size] = value;
-
-  return IntIndex(size);
-}
-
-//! Set particle attribute value.
-/** \param[in] idx Index of the attribute.
-    \param[in] value Value the attribute should be given.
- */
-void ModelData::set_value(const IntIndex idx, const Int value)
-{
-  IMP_assert(idx.get_index() < int_data_.size(),
-             "Out of range int requested");
-  int_data_[idx.get_index()] = value;
-}
-
-
-//! Add particle string attribute to the model.
-/** The returned index can be used for obtaining and setting the attribute
-    value.
-    \param[in] value Initial value of the attribute.
-    \return index of a new string attribute.
- */
-StringIndex ModelData::add_string(const String value)
-{
-  int size = string_data_.size();
-
-  string_data_.resize(size + 1);
-
-  string_data_[size] = value;
-
-  return StringIndex(size);
-}
-
-//! Set particle attribute value.
-/** \param[in] idx Index of the attribute.
-    \param[in] value Value the attribute should be given.
- */
-void ModelData::set_value(const StringIndex idx, const String value)
-{
-  IMP_assert(idx.get_index() < string_data_.size(),
-             "Out of range string requested");
-  string_data_[idx.get_index()] = value;
-}
-
-
-void ModelData::show(std::ostream &out) const
-{
-  out << "Float variables: " << std::endl;
-  for (unsigned int i=0; i< float_data_.size(); ++i) {
-    out << "  " << float_data_[i].value_ << " " << float_data_[i].deriv_ 
-        << (float_data_[i].is_optimized_?" optimized" : " not optimized")
-        << std::endl;
-  }
-}
-
-}  // namespace IMP
Index: kernel/pyext/IMP/utils.py
===================================================================
--- kernel/pyext/IMP/utils.py	(revision 434)
+++ kernel/pyext/IMP/utils.py	(working copy)
@@ -16,7 +16,6 @@
             XYZParticle._xkey = IMP.FloatKey("x")
             XYZParticle._ykey = IMP.FloatKey("y")
             XYZParticle._zkey = IMP.FloatKey("z")
-        self.model_data = self.get_model().get_model_data()
         if x is not None:
             self.add_attribute(self._xkey, x, True)
         if y is not None:
Index: kernel/pyext/IMP.i
===================================================================
--- kernel/pyext/IMP.i	(revision 434)
+++ kernel/pyext/IMP.i	(working copy)
@@ -158,9 +158,6 @@
 %include "IMP/score_states/BondDecoratorListScoreState.h"
 
 namespace IMP {
-  %template(IntIndex) Index<Int>;
-  %template(FloatIndex) Index<Float>;
-  %template(StringIndex) Index<String>;
   %template(ParticleIndex) Index<ParticleTag>;
   %template(RestraintIndex) Index<RestraintTag>;
   %template(ScoreStateIndex) Index<ScoreStateTag>;
Index: impEM/EMFitRestraint.cpp
===================================================================
--- impEM/EMFitRestraint.cpp	(revision 434)
+++ impEM/EMFitRestraint.cpp	(working copy)
@@ -14,7 +14,6 @@
 {
 
   target_dens_map = &em_map_;
-  model_data = model_.get_model_data();
   scalefac = scale_;
   model_dens_map =   new SampledDensityMap(*em_map_.get_header());
 
Index: impEM/EMFitRestraint.h
===================================================================
--- impEM/EMFitRestraint.h	(revision 434)
+++ impEM/EMFitRestraint.h	(working copy)
@@ -35,7 +35,6 @@
   DensityMap *target_dens_map;
   SampledDensityMap *model_dens_map;
   // reference to the IMP environment
-  ModelData *model_data;
   float scalefac;
   int num_particles;
   IMPParticlesAccessPoint access_p;
Index: impEM/IMPParticlesAccessPoint.h
===================================================================
--- impEM/IMPParticlesAccessPoint.h	(revision 434)
+++ impEM/IMPParticlesAccessPoint.h	(working copy)
@@ -4,7 +4,6 @@
 #include <vector>
 #include <map>
 
-#include "IMP/ModelData.h"
 #include "IMP/Model.h"
 #include "IMP/Particle.h"