diff --git a/src/Action.cpp b/src/Action.cpp
index df502c02bf0f62684a71f742dbbfc8675a51544a..bbd8e6791d9fde1d3c412d51bded249ce4f64ffe 100644
--- a/src/Action.cpp
+++ b/src/Action.cpp
@@ -185,7 +185,17 @@ void Action::warning( const std::string & msg ){
   log.printf("WARNING for action %s with label %s : %s \n", name.c_str(), label.c_str(), msg.c_str() ); 
 }
 
-
+void Action::calculateFromPDB( const PDB& pdb ){
+  activate();
+  for(Dependencies::iterator p=after.begin();p!=after.end();++p){
+     ActionWithValue*av=dynamic_cast<ActionWithValue*>(*p);
+     if(av){ av->clearInputForces(); av->clearDerivatives(); }
+     (*p)->readAtomsFromPDB( pdb ); 
+     (*p)->calculate();
+  }
+  readAtomsFromPDB( pdb );
+  calculate();
+}
 
 
 
diff --git a/src/Action.h b/src/Action.h
index 7919a2fd93fe2d250dfbcc6f21386805cab29f47..b5432059fc35846e60cc8f7caa3cdf8e47b616ee 100644
--- a/src/Action.h
+++ b/src/Action.h
@@ -7,6 +7,7 @@
 #include "Value.h"
 #include "Tools.h"
 #include "Log.h"
+#include "PDB.h"
 
 namespace PLMD{
 
@@ -203,6 +204,14 @@ public:
 
   FILE *fopen(const char *path, const char *mode);
   int   fclose(FILE*fp);
+
+/// Calculate the action given a pdb file as input.  This is used to initialize
+/// things like distance from a point in CV map space given a pdb as an input file
+  void calculateFromPDB( const PDB& pdb );
+/// This is overwritten in ActionAtomistic so that we can read 
+/// the atoms from the pdb input file rather than taking them from the 
+/// MD code
+  virtual void readAtomsFromPDB( const PDB& pdb ){}
 };
 
 /////////////////////
diff --git a/src/ActionAtomistic.cpp b/src/ActionAtomistic.cpp
index ce8a1226482f130412cf53e7178b5258eeb8b732..65d59ba47eb21f56eccad5a8f2718ba1c278412f 100644
--- a/src/ActionAtomistic.cpp
+++ b/src/ActionAtomistic.cpp
@@ -174,19 +174,12 @@ void ActionAtomistic::applyForces(){
   atoms.forceOnEnergy+=forceOnEnergy;
 }
 
-void ActionAtomistic::readAndCalculate( const PDB& pdb ){
+void ActionAtomistic::readAtomsFromPDB( const PDB& pdb ){
   for(unsigned j=0;j<indexes.size();j++){
       if( indexes[j].index()>pdb.size() ) error("there are not enough atoms in the input pdb file");
-      if( pdb.getAtomNumbers()[j].index()!=j ) error("there are atoms missing in the pdb file");  
+      if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) error("there are atoms missing in the pdb file");  
       positions[j]=pdb.getPositions()[indexes[j].index()];
   }
   for(unsigned j=0;j<indexes.size();j++) charges[j]=pdb.getBeta()[indexes[j].index()];
   for(unsigned j=0;j<indexes.size();j++) masses[j]=pdb.getOccupancy()[indexes[j].index()];
-  prepare(); calculate();
 }
-
-
-
-
-
-
diff --git a/src/ActionAtomistic.h b/src/ActionAtomistic.h
index bf2ef85c212e4b304764706b11672022852967e3..cc64f87cffdff43bfa4235832d88e63e4abdee54 100644
--- a/src/ActionAtomistic.h
+++ b/src/ActionAtomistic.h
@@ -100,7 +100,7 @@ public:
   const std::set<AtomNumber> & getUnique()const;
 /// Read in an input file containing atom positions and calculate the action for the atomic 
 /// configuration therin
-  void readAndCalculate( const PDB& pdb );
+  void readAtomsFromPDB( const PDB& pdb );
 };
 
 inline
diff --git a/src/ActionWithDistribution.cpp b/src/ActionWithDistribution.cpp
index 3f1fd141b690d3ad8e8e24fea80225e7ea8eaf3b..5c21fd43d13fd50d56b862e7ca7662373d46c021 100644
--- a/src/ActionWithDistribution.cpp
+++ b/src/ActionWithDistribution.cpp
@@ -17,7 +17,6 @@ ActionWithDistribution::ActionWithDistribution(const ActionOptions&ao):
   Action(ao),
   read(false),
   all_values(true),
-  use_field(false),
   serial(false),
   updateFreq(0),
   lastUpdate(0),
@@ -45,11 +44,10 @@ ActionWithDistribution::~ActionWithDistribution(){
 }
 
 void ActionWithDistribution::addField( std::string key, Field* ff ){
-  all_values=false; use_field=true; 
-  std::string freport, fieldin; 
-  if( key.length()==0 ) fieldin="yes";
-  else parse(key,fieldin);
+  plumed_assert( key.length()!=0 );
+  std::string fieldin; parse(key,fieldin);
   if( fieldin.length()==0 ) return ;
+  all_values=false;  std::string freport;
   myfield=ff; myfield->read( fieldin, getNumberOfFunctionsInDistribution(), freport );
   if( !myfield->check() ){
      log.printf("ERROR for keyword FIELD in action %s with label %s : %s \n \n", getName().c_str(), getLabel().c_str(), ( myfield->errorMessage() ).c_str() );
@@ -100,9 +98,8 @@ void ActionWithDistribution::requestDistribution(){
 
   if(all_values){
       error("No function has been specified");
-  } else if( !use_field ){
-      plumed_massert( functions.size()==final_values.size(), "number of functions does not match number of values" );
-  }
+  } 
+  plumed_massert( functions.size()==final_values.size(), "number of functions does not match number of values" );
   // This sets up the dynamic list that holds what we are calculating
   for(unsigned i=0;i<getNumberOfFunctionsInDistribution();++i){ members.addIndexToList(i); }
   members.activateAll(); members.updateActiveMembers();
@@ -117,12 +114,13 @@ void ActionWithDistribution::prepare(){
     ActionWithValue* av=dynamic_cast<ActionWithValue*>(this);
     for(unsigned i=0;i<functions.size();++i) functions[i]->setNumberOfDerivatives( (av->getPntrToComponent(i))->getNumberOfDerivatives() );
     // Setup the buffers for mpi gather
-    if( use_field ){
+    if( myfield ){
       std::vector<unsigned> cv_sizes( getNumberOfFunctionsInDistribution() ); 
       for(unsigned i=0;i<cv_sizes.size();++i){ cv_sizes[i]=getThisFunctionsNumberOfDerivatives(i); }
       myfield->resizeBaseQuantityBuffers( cv_sizes ); 
       myfield->resizeDerivatives( getNumberOfFieldDerivatives() );
-    } else {
+    } 
+    if( functions.size()!=0 ){
       unsigned bufsize=0;
       for(unsigned i=0;i<functions.size();++i) bufsize+=functions[i]->requiredBufferSpace();
       buffer.resize( bufsize ); 
@@ -137,12 +135,13 @@ void ActionWithDistribution::prepare(){
     ActionWithValue* av=dynamic_cast<ActionWithValue*>(this);
     for(unsigned i=0;i<functions.size();++i) functions[i]->setNumberOfDerivatives( (av->getPntrToComponent(i))->getNumberOfDerivatives() );
     // Setup the buffers for mpi gather
-    if( use_field ){
+    if( myfield ){
       std::vector<unsigned> cv_sizes( getNumberOfFunctionsInDistribution() ); unsigned kk;
       for(unsigned i=0;i<cv_sizes.size();++i){ cv_sizes[i]=getThisFunctionsNumberOfDerivatives(i); }
       myfield->resizeBaseQuantityBuffers( cv_sizes ); 
       myfield->resizeDerivatives( getNumberOfFieldDerivatives() );
-    } else {
+    } 
+    if( functions.size()!=0 ){
       unsigned bufsize=0;
       for(unsigned i=0;i<functions.size();++i) bufsize+=functions[i]->requiredBufferSpace();
       buffer.resize( bufsize ); 
diff --git a/src/ActionWithDistribution.h b/src/ActionWithDistribution.h
index 62edfeb4f7a5c7d7c60b6d75fe9bd1c972f0f40c..3a14c49770686854a1a4db39cb71aa323c8a03d8 100644
--- a/src/ActionWithDistribution.h
+++ b/src/ActionWithDistribution.h
@@ -18,8 +18,6 @@ private:
   bool read;
 /// This tells us we are calculating all values (not doing anything to the distribution)
   bool all_values;
-/// This tells us we are using the field cvs
-  bool use_field;
 /// Do all calculations in serial
   bool serial;
 /// Everything for controlling the updating of neighbor lists
@@ -68,8 +66,6 @@ public:
   virtual void prepare();
 /// Calculate the values of the object
   void calculate();
-/// Are we using distributions 
-  bool usingDistributionFunctions() const;
 /// Overwrite this in your inherited actions if neighbor list update is more complicated
 /// than just calculating everything and seeing whats big.
   virtual void prepareForNeighborListUpdate(){};
@@ -84,8 +80,6 @@ public:
   virtual bool isPeriodic(const unsigned nn)=0;
 /// What are the domains of the base quantities
   virtual void retrieveDomain( const unsigned nn, double& min, double& max);
-/// Get the number of derivatives for this action
-//  virtual unsigned getNumberOfDerivatives()=0;
 /// Get the number of functions from which we are calculating the distribtuion
   virtual unsigned getNumberOfFunctionsInDistribution()=0;
 /// Calculate one of the functions in the distribution
@@ -115,11 +109,6 @@ bool ActionWithDistribution::getSerial() const {
   return serial;
 }
 
-inline
-bool ActionWithDistribution::usingDistributionFunctions() const {
-  return ( !all_values && !use_field );
-}
-
 inline
 bool ActionWithDistribution::isTimeForNeighborListUpdate() const {
   return reduceAtNextStep;
diff --git a/src/Atoms.h b/src/Atoms.h
index 1d1aa5a9b1426e4a6e069f1a4408c4ec3a0165b4..22d1b7456e0f9c25d8ef76761e8770504bc86062 100644
--- a/src/Atoms.h
+++ b/src/Atoms.h
@@ -152,6 +152,7 @@ public:
   void readBinary(std::istream&);
   double getKBoltzmann()const;
   double getMDKBoltzmann()const;
+  bool usingNaturalUnits()const;
   void setNaturalUnits(bool n){naturalUnits=n;};
   void setMDNaturalUnits(bool n){MDnaturalUnits=n;};
 };
@@ -171,6 +172,11 @@ ActionWithVirtualAtom* Atoms::getVirtualAtomsAction(AtomNumber i)const{
   return virtualAtomsActions[i.index()-getNatoms()];
 }
 
+inline
+bool Atoms::usingNaturalUnits() const {
+  return naturalUnits;
+}
+
 
 
 }
diff --git a/src/ColvarRMSD.cpp b/src/ColvarRMSD.cpp
index 4dc139bd49c615d66d2affd94d05b0839af3e0b4..288cf2e882895c2ea4011c2626e592c7dc9f1380 100644
--- a/src/ColvarRMSD.cpp
+++ b/src/ColvarRMSD.cpp
@@ -53,7 +53,7 @@ PLUMED_REGISTER_ACTION(ColvarRMSD,"RMSD")
 
 void ColvarRMSD::registerKeywords(Keywords& keys){
   Colvar::registerKeywords(keys);
-  keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure.  The atoms involved in the CV are specified in the pdb file");
+  keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV. " + PDB::documentation() );
   keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed.  Should be OPTIMAL or SIMPLE.");
 }
 
@@ -72,9 +72,8 @@ PLUMED_COLVAR_INIT(ao),rmsd(log)
   addValueWithDerivatives(); setNotPeriodic();
   PDB pdb;
 
-  // read everything in ang and transform to nm
-
-  pdb.read(reference,0.1/atoms.getUnits().length);
+  // read everything in ang and transform to nm if we are not in natural units
+  pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().length);
 
   rmsd.set(pdb,type);
 
diff --git a/src/ColvarTarget.cpp b/src/ColvarTarget.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dd6dc01e492d48bb07c36c75ff04ded0459b2b27
--- /dev/null
+++ b/src/ColvarTarget.cpp
@@ -0,0 +1,68 @@
+#include "Function.h"
+#include "PlumedMain.h"
+#include "ActionRegister.h"
+#include "PDB.h"
+#include "TargetDist.h"
+#include "Atoms.h"
+
+using namespace std;
+
+namespace PLMD {
+
+//+PLUMEDOC FUNCTION TARGET
+/**
+
+The pythagorean distance from a particular structure measured in the space defined by some 
+set of collective variables
+
+*/
+//+ENDPLUMEDOC
+
+class TargetFrame : public Function {
+private:
+  TargetDist target;
+  std::vector<double> derivs;
+public:
+  TargetFrame(const ActionOptions&);
+  virtual void calculate();
+  static void registerKeywords(Keywords& keys );
+};
+
+PLUMED_REGISTER_ACTION(TargetFrame,"TARGET")
+
+void TargetFrame::registerKeywords(Keywords& keys){
+  Function::registerKeywords(keys);
+  keys.use("ARG");
+  keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure. " + PDB::documentation() );
+  keys.add("optional","REFERENCE_VEC","the vector of values for the CVs at the reference point (if you use this you don't need REFERENCE)");
+}
+
+TargetFrame::TargetFrame(const ActionOptions&ao):
+Action(ao),
+Function(ao),
+target(log)
+{
+  std::vector<double> targ;
+  parseVector("REFERENCE_VEC",targ);
+  if( targ.size()!=0 ){
+    target.read( targ, getArguments() );
+  } else {
+    string reference;
+    parse("REFERENCE",reference);
+    PDB pdb; 
+    pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/plumed.getAtoms().getUnits().length);
+    printf("Read pdb file with %d atoms inside\n",pdb.size());
+    target.read( pdb, getArguments() );
+  }
+  checkRead();
+  derivs.resize( getNumberOfArguments() );
+  addValueWithDerivatives(); setNotPeriodic();
+}
+
+void TargetFrame::calculate(){
+  double r=target.calculate( derivs );
+  setValue(r);
+  for(unsigned i=0;i<derivs.size();i++) setDerivative(i,derivs[i]);
+}
+
+}
diff --git a/src/DistributionFunctions.h b/src/DistributionFunctions.h
index 5658d37f2918b764596f5a0a7604c66c9c1c2105..91dc1084c53e57af9b024523371ec7ffeb8210ed 100644
--- a/src/DistributionFunctions.h
+++ b/src/DistributionFunctions.h
@@ -159,6 +159,7 @@ public:
 
 class within : public DistributionFunction {
 private:
+  bool usenorm;
   double a,b,sigma;
   HistogramBead hist;
 public:
diff --git a/src/DistributionWithin.cpp b/src/DistributionWithin.cpp
index 13c2d35949dd8397427aea2814c193b3600ab5d8..b3f89758035e3b50d7efe2f573cc1fac11a5a9d6 100644
--- a/src/DistributionWithin.cpp
+++ b/src/DistributionWithin.cpp
@@ -4,22 +4,29 @@ namespace PLMD {
 
 std::string within::documentation(){
   std::ostringstream ostr;
-  ostr<<"To make this quantity continuous it is calculated using "<<HistogramBead::documentation(false);
+  ostr<<"To make this quantity continuous it is calculated using "<<HistogramBead::documentation(false)<<". ";
+  ostr<<"Adding the NORM flag allows you to calculate the fraction of colvars in the particular range rather than the total number ";
+  ostr<<"(N.B. this option should probably not used if you are using neighbor lists and the derivatives of the WITHIN value).";
   return ostr.str();
 }
 
 within::within( const std::string& parameters ) :
-DistributionFunction(parameters)
+DistributionFunction(parameters),
+usenorm(false)
 { 
   std::string errormsg;
+  std::vector<std::string> data=Tools::getWords(parameters);
+  Tools::parseFlag(data,"NORM",usenorm);
   hist.set( parameters, "", errormsg );
   if( errormsg.size()!=0 ) error( errormsg );
   addAccumulator( true );
+  addAccumulator( false );
 }
 
 std::string within::message(){
   std::ostringstream ostr;
-  ostr<<"number of values "<<hist.description(); 
+  if(usenorm) ostr<<"fraction of values "<<hist.description();
+  else ostr<<"number of values "<<hist.description(); 
   return ostr.str();
 }
 
@@ -38,11 +45,17 @@ void within::calculate( Value* value_in, std::vector<Value>& aux ){
   copyValue( 0, value_in ); 
   double df, f; f=hist.calculate( value_in->get() , df );
   chainRule(0, df); setValue(0, f);
+  setValue( 1, 1.0 );
 }
 
 void within::finish( Value* value_out ){
   extractDerivatives( 0, value_out );
-  value_out->set( getPntrToAccumulator(0)->get() );
+  if( usenorm ) {
+    double nval=getPntrToAccumulator(1)->get();
+    value_out->chainRule(1.0/nval); value_out->set(getPntrToAccumulator(0)->get()/nval);
+  } else {
+    value_out->set( getPntrToAccumulator(0)->get() );
+  }
 }
 
 }
diff --git a/src/FieldBias.cpp b/src/FieldBias.cpp
index 729d7fc8aad361e480117b810809aa224ae994df..56ced3add23e06a45c648e8bebabf616f04dc558 100644
--- a/src/FieldBias.cpp
+++ b/src/FieldBias.cpp
@@ -36,8 +36,8 @@ FieldBias::FieldBias(const ActionOptions&ao):
   std::string ll; parse("FIELD",ll);
   ActionWithDistribution* field=plumed.getActionSet().selectWithLabel<ActionWithDistribution*>(ll);
   if(field){
-     if( field->usingDistributionFunctions() ) error("action " + ll + " calculates a colvar and not a field"); 
      myfield=field->getField();
+     if(!myfield) error("input action " + ll + " does not calcualte a field");
      apply_action=dynamic_cast<ActionWithValue*>( field );
   } else {
      apply_action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(ll);
diff --git a/src/HistogramBead.cpp b/src/HistogramBead.cpp
index 94d559edef833a0dfa784f342a495ad585004997..1c27d797bb59f7969ec28a9b4f5a654d42c3e9b4 100644
--- a/src/HistogramBead.cpp
+++ b/src/HistogramBead.cpp
@@ -32,7 +32,8 @@ std::string HistogramBead::histodocs() {
   ostr<<"The range is divided into a discete number of bins and the number of values that fall within each bin is calculated using ";
   ostr<<"\\f$ w(r)=\\int_a^b \\frac{1}{\\sqrt{2\\pi}\\sigma} \\exp\\left( -\\frac{(r'-r)^2}{2\\sigma^2} \\right) \\textrm{d}r' \\f$";
   ostr<<"where \\f$ \\sigma=(b-a)k \\f$.  The particular range of interest and number of bins are specified using ";
-  ostr<<"(NBINS=\\f$n\\f$ LOWER=\\f$a\\f$ UPPER=\\f$b\\f$ SMEAR=\\f$x\\f$). If the SMEAR keyword is not present then by default \\f$k=0.5\\f$";
+  ostr<<"(NBINS=\\f$n\\f$ LOWER=\\f$a\\f$ UPPER=\\f$b\\f$ SMEAR=\\f$x\\f$). If the SMEAR keyword is not present then by default \\f$k=0.5\\f$. ";
+  ostr<<"You can calculate a normalized histogram using the NORM flag (N.B. Don't use this if you are using derivatives of the histogram and neighbor lists)";
   return ostr.str();
 }
 
@@ -51,12 +52,15 @@ void HistogramBead::generateBins( const std::string& params, const std::string&
   plumed_massert(range[0]<range[1],"Range specification is dubious"); 
   bool found_b=Tools::parse(data,dd+"SMEAR",smear);
   if(!found_b){ Tools::convert(0.5,smear); }  
+  bool usenorm=false; std::string normstr;
+  if(dd=="") bool found_n=Tools::parseFlag(data,dd+"NORM",usenorm);
+  if(usenorm && dd==""){ normstr="NORM"; } else { normstr=""; } 
 
   std::string lb,ub; double delr = ( range[1]-range[0] ) / static_cast<double>( nbins );
   for(unsigned i=0;i<nbins;++i){
      Tools::convert( range[0]+i*delr, lb );
      Tools::convert( range[0]+(i+1)*delr, ub );
-     bins.push_back( dd + "LOWER=" + lb + " " + dd + "UPPER=" + ub + " " + dd + "SMEAR=" + smear );
+     bins.push_back( dd + "LOWER=" + lb + " " + dd + "UPPER=" + ub + " " + dd + "SMEAR=" + smear + " " + normstr );
   }
   plumed_assert(bins.size()==nbins);
   if( dd.empty() ) plumed_massert(data.empty(),"Error reading histogram"); 
@@ -76,6 +80,7 @@ void HistogramBead::set( const std::string& params, const std::string& dd, std::
   
   smear=0.5; bool found_b=Tools::parse(data,dd+"SMEAR",smear);
   width=smear*(highb-lowb); init=true;
+  bool usenorm; bool found_n=Tools::parseFlag(data,dd+"NORM",usenorm);
   if( dd.empty() ){ if( !data.empty()) errormsg="Error reading within"; }
 }
 
@@ -88,5 +93,6 @@ void HistogramBead::printKeywords( Log& log ) const {
   hkeys.add("compulsory","LOWER","the lower boundary for this particular bin");
   hkeys.add("compulsory","UPPER","the upper boundary for this particular bin");
   hkeys.add("compulsory","SMEAR","0.5","the ammount to smear the Gaussian for each value in the distribution"); 
+  hkeys.addFlag("NORM",false,"normalize the histogram according to the number of values we are histograming");
   hkeys.print( log );
 }
diff --git a/src/MultiColvar.cpp b/src/MultiColvar.cpp
index a1c92dffd735983b445d1fbc95e183010fb99931..d020de23ccd6cba5e2f6c88e85afdbe3f515e414 100644
--- a/src/MultiColvar.cpp
+++ b/src/MultiColvar.cpp
@@ -59,7 +59,7 @@ ActionWithValue(ao),
 ActionWithDistribution(ao),
 usepbc(true),
 readatoms(false),
-setperiods(false),
+//setperiods(false),
 needsCentralAtomPosition(false)
 {
   if( keywords.style("NOPBC", "flag") ){ 
@@ -194,7 +194,7 @@ void MultiColvar::readAtoms( int& natoms ){
       }
   } 
 
-  if( !usingDistributionFunctions() && keywords.exists("DISTRIBUTION") ) addField("DISTRIBUTION", new Field( "identity",1 ) );
+  if( keywords.exists("DISTRIBUTION") ) addField("DISTRIBUTION", new Field( "identity",1 ) );
   requestAtoms();         // Request the atoms in ActionAtomistic and set up the value sizes
 }
 
@@ -325,34 +325,6 @@ void MultiColvar::readSpeciesKeyword( int& natoms ){
   } 
 }
 
-void MultiColvar::checkRead(){
-  plumed_massert(setperiods, "You must set the periodicity of the various component functions");
-  Action::checkRead();
-}
-
-void MultiColvar::setNotPeriodic(){
-  setperiods=true;
-  if( !usingDistributionFunctions() ){
-      std::string num;
-      for(unsigned i=0;i<getNumberOfComponents();++i){
-          Tools::convert(i+1,num);
-          componentIsNotPeriodic( "fval"+num );
-      }
-  }
-}
-
-void MultiColvar::setPeriodicDomain( const double& min, const double max ){
-  setperiods=true;
-  if( !usingDistributionFunctions() ){
-     std::string num;
-     for(unsigned i=0;i<getNumberOfComponents();++i){ 
-         Tools::convert(i+1,num);
-         componentIsPeriodic( "fval"+num , min, max );
-     }
-  }     
-}
-
-
 void MultiColvar::prepareForNeighborListUpdate(){
    for(unsigned i=0;i<colvar_atoms.size();++i){
       colvar_atoms[i].activateAll(); colvar_atoms[i].updateActiveMembers();
@@ -372,9 +344,7 @@ void MultiColvar::completeNeighborListUpdate(){
 
 void MultiColvar::requestAtoms(){
    ActionAtomistic::requestAtoms( all_atoms.retrieveActiveList() );
-   if( usingDistributionFunctions() ){
-       for(unsigned i=0;i<getNumberOfComponents();++i) getPntrToComponent(i)->resizeDerivatives(3*getNumberOfAtoms()+9);
-   } 
+   for(unsigned i=0;i<getNumberOfComponents();++i) getPntrToComponent(i)->resizeDerivatives(3*getNumberOfAtoms()+9);
 }
 
 void MultiColvar::getCentralAtom( const std::vector<Vector>& pos, Vector& cpos, std::vector<Tensor>& deriv ){
diff --git a/src/MultiColvar.h b/src/MultiColvar.h
index 41f830043864af787066408bf93760694ddcc365..1b74409501cc4510545b4af6f06344f65b4d89bb 100644
--- a/src/MultiColvar.h
+++ b/src/MultiColvar.h
@@ -161,7 +161,6 @@ class MultiColvar :
 private:
   bool usepbc;
   bool readatoms;
-  bool setperiods;
   bool needsCentralAtomPosition;
 /// Constants for fields
   double fsigma2, fnorm;
@@ -193,11 +192,6 @@ protected:
   void removeAtomRequest( const unsigned& aa );
 /// Do we use pbc to calculate this quantity
   bool usesPbc() const ;
-/// Check the readin
-  void checkRead();
-/// Set the periodicities of the base quantities for the fields
-  void setNotPeriodic();
-  void setPeriodicDomain( const double& min, const double max );
 public:
   MultiColvar(const ActionOptions&);
   ~MultiColvar(){};
diff --git a/src/MultiColvarCoordination.cpp b/src/MultiColvarCoordination.cpp
index e69aa2e2b672753baa2c15a6efaedabcc4b08fad..31f02678c299830d06fa9ca9bfa93f619660c640 100644
--- a/src/MultiColvarCoordination.cpp
+++ b/src/MultiColvarCoordination.cpp
@@ -100,8 +100,6 @@ PLUMED_MULTICOLVAR_INIT(ao)
   if(nl_cut>0.0){
      log.printf("  ignoring distances greater than %lf in neighbor list\n",nl_cut); 
   } 
-  // Functon is not periodic
-  setNotPeriodic();
   // And check everything has been read in correctly
   checkRead();
 }
diff --git a/src/MultiColvarDensity.cpp b/src/MultiColvarDensity.cpp
index 5073ffc22f67d525497b12265977978214c616c9..4ac94357d41baeb784351894be08fecc68aea783 100644
--- a/src/MultiColvarDensity.cpp
+++ b/src/MultiColvarDensity.cpp
@@ -50,8 +50,6 @@ PLUMED_MULTICOLVAR_INIT(ao)
 {
   int nat; readAtoms( nat ); 
   requestDistribution();
-  // Functon is not periodic
-  setNotPeriodic();
   // And check everything has been read in correctly
   checkRead(); 
 }
diff --git a/src/MultiColvarDistance.cpp b/src/MultiColvarDistance.cpp
index 908b097ea3618d1651230affe8c2cfec7e6eea2b..e6a2713dd88790c8c7e7b80b1b5e78c6a384c475 100644
--- a/src/MultiColvarDistance.cpp
+++ b/src/MultiColvarDistance.cpp
@@ -79,8 +79,6 @@ rcut(-1)
       parse("NL_CUTOFF",rcut);
       if( rcut>0 ) log.printf("  ignoring distances greater than %lf in neighbor list\n",rcut);
   }
-  // Functon is not periodic
-  setNotPeriodic();
   // And check everything has been read in correctly
   checkRead();
 }
diff --git a/src/PDB.cpp b/src/PDB.cpp
index 083654c58784e06de01a3c0ef8428a78e3b06f58..d4d8d8b1046a45c2793d7d75630ff2fc4bd45fcf 100644
--- a/src/PDB.cpp
+++ b/src/PDB.cpp
@@ -7,6 +7,15 @@ using namespace std;
 
 namespace PLMD{
 
+std::string PDB::documentation(){
+  std::ostringstream ostr;
+  ostr<<"In PDB files the atomic coordinates and box lengths should be in Angstroms unless you are working with natural units. ";
+  ostr<<"If you are working with natural units then the coordinates should be in your natural length unit. In the PDB files used ";
+  ostr<<"by plumed the beta column is used to specify the charges on the atoms and the occupancy column is used to specify the atomic masses. ";
+  ostr<<"For more details on the PDB file format visit http://www.wwpdb.org/docs.html";
+  return ostr.str();
+}
+
 const std::vector<Vector> & PDB::getPositions()const{
   return positions;
 }
@@ -27,7 +36,8 @@ unsigned PDB::size()const{
   return positions.size();
 }
 
-void PDB::read(const std::string&file,double scale){
+void PDB::read(const std::string&file,bool naturalUnits,double scale){
+  if(naturalUnits) scale=1.0;
   FILE* fp=fopen(file.c_str(),"r");
   //cerr<<file<<endl;
   string line;
diff --git a/src/PDB.h b/src/PDB.h
index 09c652bf451d1db1cd17e98f80a9279bb730dadc..a038c90224cee546314066fb73ab6717ed1160d8 100644
--- a/src/PDB.h
+++ b/src/PDB.h
@@ -17,8 +17,10 @@ class PDB{
   std::vector<double> beta;
   std::vector<AtomNumber> numbers;
 public:
+/// The documentation for PDB file parser
+  static std::string documentation();
 /// Read the pdb from a file, scaling positions by a factor scale
-  void read(const std::string&file,double scale);
+  void read(const std::string&file,bool naturalUnits,double scale);
 /// Access to the position array
   const std::vector<Vector>     & getPositions()const;
 /// Access to the occupancy array
diff --git a/src/TargetDist.cpp b/src/TargetDist.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cea6f5c37b9e6df91b5f34aa92ab4eb8ddda47c4
--- /dev/null
+++ b/src/TargetDist.cpp
@@ -0,0 +1,47 @@
+#include "TargetDist.h"
+
+namespace PLMD {
+
+void TargetDist::read( const PDB& pdb, std::vector<Value*> ar ){
+  // Clear values in target actions
+  for(unsigned i=0;i<ar.size();++i){
+     (ar[i]->getPntrToAction())->clearInputForces();
+     (ar[i]->getPntrToAction())->clearDerivatives();
+  }
+
+  // Caclulate target actions from input in PDB file
+  std::vector<double> targ( ar.size() );
+  for(unsigned i=0;i<ar.size();++i){
+      if( ar[i]->valueHasBeenSet() ){ 
+         targ[i]=ar[i]->get();
+      } else {
+         ActionWithValue* vv=ar[i]->getPntrToAction();
+         (ar[i]->getPntrToAction())->calculateFromPDB( pdb );
+         targ[i]=ar[i]->get();
+      }
+  }
+  read( targ, ar );
+}
+
+void TargetDist::read( const std::vector<double>& targ, std::vector<Value*> ar ){
+  plumed_assert( targ.size()==ar.size() );
+
+  target.resize( ar.size() ); args.resize( ar.size() );
+  log.printf("  distance from this point in cv space : ");
+  for(unsigned i=0;i<target.size();++i){ log.printf("%f ", targ[i]); target[i]=targ[i]; args[i]=ar[i]; }
+  log.printf("\n");
+}
+
+double TargetDist::calculate( std::vector<double>& derivs ){
+  plumed_assert( derivs.size()==args.size() );
+  double dist=0, tmp;
+  for(unsigned i=0;i<args.size();++i){
+      tmp=args[i]->difference( target[i], args[i]->get() );
+      derivs[i]=tmp; dist+=tmp*tmp; 
+  }
+  dist=sqrt(dist);
+  for(unsigned i=0;i<args.size();++i) derivs[i]/=dist;
+  return dist;
+}
+
+}
diff --git a/src/TargetDist.h b/src/TargetDist.h
new file mode 100644
index 0000000000000000000000000000000000000000..23a931f0398cdf87591d0161d1f21dd9e11d5399
--- /dev/null
+++ b/src/TargetDist.h
@@ -0,0 +1,29 @@
+#ifndef __PLUMED_TargetDist_h
+#define __PLUMED_TargetDist_h
+
+#include "Value.h"
+#include "ActionWithValue.h"
+#include "PDB.h"
+#include <vector>
+#include <string>
+
+namespace PLMD{
+
+class Log;
+class PDB;
+
+class TargetDist {
+private:
+  std::vector<Value*> args; 
+  std::vector<double> target;
+  Log &log;
+public:
+  TargetDist(Log& log) : log(log) {};
+  void read( const PDB& pdb, std::vector<Value*> args ); 
+  void read( const std::vector<double>& targ, std::vector<Value*> ar );
+  double calculate( std::vector<double>& derivs );
+};
+
+}
+
+#endif
diff --git a/src/Value.cpp b/src/Value.cpp
index 9c12ce19e6c1f61541694f06a749d0ee689f8ba4..09fd5fb0c3e30a15a5f12b7b569bfd5bf1b8e8a8 100644
--- a/src/Value.cpp
+++ b/src/Value.cpp
@@ -103,6 +103,12 @@ double Value::projection(const Value& v1,const Value&v2){
   return proj;
 }
 
+ActionWithValue* Value::getPntrToAction(){
+  plumed_assert( action!=NULL );
+  return action;
+}
+
+
 
 
 
diff --git a/src/Value.h b/src/Value.h
index 6dafffffcab8d162c105ad13796d468307599a78..55b5e4b46fd876ab00ea93d7783df35a8c15fcbd 100644
--- a/src/Value.h
+++ b/src/Value.h
@@ -109,6 +109,8 @@ public:
   double difference(double)const;
 /// Calculate the difference between two values of this function
   double difference(double,double)const;
+/// This returns the pointer to the action where this value is calculated
+  ActionWithValue* getPntrToAction();
 /// This sets up the gradients
   void setGradients();
   static double projection(const Value&,const Value&);