Index: kernel/include/IMP/optimizers/ConjugateGradients.h
===================================================================
--- kernel/include/IMP/optimizers/ConjugateGradients.h	(revision 533)
+++ kernel/include/IMP/optimizers/ConjugateGradients.h	(working copy)
@@ -30,6 +30,9 @@
 
   //! Set the threshold for the minimum gradient
   void set_threshold(Float t){ threshold_=t;}
+
+  //! Limit how far anything can change each time step
+  void set_max_change(Float t){ max_change_=t;}
 private:
 
   Float get_score(std::vector<FloatIndex> float_indices,
@@ -40,6 +43,7 @@
                    int max_steps, const std::vector<Float> &search,
                    const std::vector<Float> &estimate);
   Float threshold_;
+  Float max_change_;
 };
 
 } // namespace IMP
Index: kernel/src/optimizers/ConjugateGradients.cpp
===================================================================
--- kernel/src/optimizers/ConjugateGradients.cpp	(revision 533)
+++ kernel/src/optimizers/ConjugateGradients.cpp	(working copy)
@@ -33,6 +33,15 @@
   int i, opt_var_cnt = float_indices.size();
   /* set model state */
   for (i = 0; i < opt_var_cnt; i++) {
+    IMP_check(x[i]==x[i], "Got NaN in CG",
+              ValueException("Got NaN in CG"));
+    if (std::abs(x[i] - get_value(float_indices[i])) > max_change_) {
+      if (x[i] < get_value(float_indices[i])) {
+        x[i]= get_value(float_indices[i]) - max_change_;
+      } else {
+        x[i]= get_value(float_indices[i]) + max_change_;
+      }
+    }
     set_value(float_indices[i], x[i]);
   }
 
@@ -42,6 +51,10 @@
   /* get derivatives */
   for (i = 0; i < opt_var_cnt; i++) {
     dscore[i] = get_derivative(float_indices[i]);
+    IMP_check(dscore[i]== dscore[i] && dscore[i]
+              != std::numeric_limits<Float>::infinity()
+              && dscore[i] != - std::numeric_limits<Float>::infinity(),
+              "Bad input to CG", ValueException("Bad input to CG"));
   }
   return score;
 }
@@ -202,6 +215,7 @@
 ConjugateGradients::ConjugateGradients()
 {
   threshold_=std::numeric_limits<Float>::epsilon();
+  max_change_= std::numeric_limits<Float>::max()/100.0;
 }
 
 
@@ -230,6 +244,9 @@
   // get initial state in x(n):
   for (i = 0; i < n; i++) {
     x[i] = get_value(float_indices[i]);
+    IMP_check(x[i]== x[i] && x[i] != std::numeric_limits<Float>::infinity()
+              && x[i] != - std::numeric_limits<Float>::infinity(),
+              "Bad input to CG", ValueException("Bad input to CG"));
   }
 
   // Initialize optimization variables