diff --git a/src/gridtools/ActionWithInputGrid.cpp b/src/gridtools/ActionWithInputGrid.cpp
index 5cb82f3e7e05dc9f714189e15e8f33eb213d3026..e409d9662bec7b510cf9d5b332a06569c9ea8e4e 100644
--- a/src/gridtools/ActionWithInputGrid.cpp
+++ b/src/gridtools/ActionWithInputGrid.cpp
@@ -59,7 +59,7 @@ ingrid(NULL)
 }
 
 void ActionWithInputGrid::clearAverage(){
-  mygrid->setBounds( ingrid->getMin(), ingrid->getMax(), mygrid->getNbin(), mygrid->getGridSpacing() );
+  if( mygrid->getType()=="flat" ) mygrid->setBounds( ingrid->getMin(), ingrid->getMax(), mygrid->getNbin(), mygrid->getGridSpacing() );
   ActionWithAveraging::clearAverage();
 }
 
diff --git a/src/gridtools/ContourFindingBase.cpp b/src/gridtools/ContourFindingBase.cpp
index 0ce55ccc1172cbe8758531980896c4ff5776709c..dd9cfc50f8eab3db1d15e10205bd7424ee9696dc 100644
--- a/src/gridtools/ContourFindingBase.cpp
+++ b/src/gridtools/ContourFindingBase.cpp
@@ -20,7 +20,6 @@
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
 #include "ContourFindingBase.h"
-#include "vesselbase/StoreDataVessel.h"
 #include "core/PlumedMain.h"
 #include "core/Atoms.h"
 
@@ -39,8 +38,7 @@ void ContourFindingBase::registerKeywords( Keywords& keys ){
 ContourFindingBase::ContourFindingBase(const ActionOptions&ao):
 Action(ao),
 ActionWithInputGrid(ao),
-mymin(this),
-mydata(NULL)
+mymin(this)
 {
   if( ingrid->noDerivatives() ) error("cannot find contours if input grid has no derivatives");
   parse("CONTOUR",contour); 
@@ -73,26 +71,9 @@ mydata(NULL)
          else lenunit=1.0; 
 
          of.link(*this); of.open(file); 
-
-         // Now create a store data vessel to hold the points
-         mydata=buildDataStashes( NULL );
       }
   }
 }
 
-void ContourFindingBase::finishAveraging(){
-  // And this looks after the output
-  if( mydata ){
-      std::vector<double> point( 1 + ingrid->getDimension() );
-      of.printf("%u\n",mydata->getNumberOfStoredValues());  
-      of.printf("Points found on isocontour\n");
-      for(unsigned i=0;i<mydata->getNumberOfStoredValues();++i){
-          mydata->retrieveSequentialValue( i, false, point ); of.printf("X");
-          for(unsigned j=0;j<ingrid->getDimension();++j) of.printf( (" " + fmt_xyz).c_str(), lenunit*point[1+j] );
-          of.printf("\n");
-      }   
-  }
-}
-
 }
 }
diff --git a/src/gridtools/ContourFindingBase.h b/src/gridtools/ContourFindingBase.h
index 57611ba2ef1913dfd8bad496b42b4957af8c3b13..958e31f6f383a5cce29d0ee69423e7f8ea034833 100644
--- a/src/gridtools/ContourFindingBase.h
+++ b/src/gridtools/ContourFindingBase.h
@@ -30,15 +30,13 @@ namespace gridtools {
 
 class ContourFindingBase : public ActionWithInputGrid {
 private:
+/// This is the object that does the root finding
+  RootFindingBase<ContourFindingBase> mymin;
+protected:
 /// Stuff for output 
   OFile of;
   double lenunit;
   std::string fmt_xyz;
-/// This is the object that does the root finding
-  RootFindingBase<ContourFindingBase> mymin;
-/// The data is stored in a grid
-  vesselbase::StoreDataVessel* mydata;
-protected:
 /// Where you would like to find the contour
   double contour;
 /// Find a contour along line specified by direction
@@ -52,10 +50,6 @@ public:
   unsigned getNumberOfDerivatives(){ return 0; }
 /// This is not periodic
   bool isPeriodic(){ return false; }
-/// Number of quantities is the number of points in each point on the grid
-  virtual unsigned getNumberOfQuantities() const { return 1 + ingrid->getDimension(); }
-/// This does output if needs be
-  void finishAveraging();
 };
 
 inline
diff --git a/src/gridtools/FindContour.cpp b/src/gridtools/FindContour.cpp
index 45d3d7416a6b45c6dbd9a7746f73354a4697cc2c..097645f50d251125cb99b991c5c08c07445c07b9 100644
--- a/src/gridtools/FindContour.cpp
+++ b/src/gridtools/FindContour.cpp
@@ -20,6 +20,7 @@
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
 #include "core/ActionRegister.h"
+#include "vesselbase/StoreDataVessel.h"
 #include "ContourFindingBase.h"
 
 //+PLUMEDOC GRIDANALYSIS FIND_CONTOUR
@@ -38,12 +39,15 @@ class FindContour : public ContourFindingBase {
 private:
   bool firsttime;
   unsigned gbuffer;
+/// The data is stored in a grid
+  vesselbase::StoreDataVessel* mydata;
 public:
   static void registerKeywords( Keywords& keys );
   explicit FindContour(const ActionOptions&ao);
   bool checkAllActive() const { return gbuffer==0; }
   void prepareForAveraging();
   bool isPeriodic(){ return false; }
+  unsigned getNumberOfQuantities() const { return 1 + ingrid->getDimension(); }
   void compute( const unsigned& current, MultiValue& myvals ) const ;
   void finishAveraging();
 };
@@ -64,7 +68,7 @@ firsttime(true)
 
   parse("BUFFER",gbuffer);
   if( gbuffer>0 ) log.printf("  after first step a subset of only %u grid points around where the countour was found will be checked\n",gbuffer);
-  checkRead();
+  checkRead(); mydata=buildDataStashes( NULL );
 }
 
 void FindContour::prepareForAveraging(){
@@ -127,7 +131,6 @@ void FindContour::compute( const unsigned& current, MultiValue& myvals ) const {
 }
 
 void FindContour::finishAveraging(){
-  ContourFindingBase::finishAveraging();
   // And update the list of active grid points
   if( gbuffer>0 ){
       std::vector<unsigned> neighbours; unsigned num_neighbours;
@@ -145,6 +148,14 @@ void FindContour::finishAveraging(){
       }
       ingrid->activateThesePoints( active );
   }
+  std::vector<double> point( 1 + ingrid->getDimension() );
+  of.printf("%u\n",mydata->getNumberOfStoredValues());
+  of.printf("Points found on isocontour\n");
+  for(unsigned i=0;i<mydata->getNumberOfStoredValues();++i){
+      mydata->retrieveSequentialValue( i, false, point ); of.printf("X");
+      for(unsigned j=0;j<ingrid->getDimension();++j) of.printf( (" " + fmt_xyz).c_str(), lenunit*point[1+j] );
+      of.printf("\n");
+  }
 }
 
 }
diff --git a/src/gridtools/FindSphericalContour.cpp b/src/gridtools/FindSphericalContour.cpp
index ab494a748c0217ad2d8e5ae08d595a3cb563dac8..fb8939ce4c41c83706f618f38778af56ee3238e9 100644
--- a/src/gridtools/FindSphericalContour.cpp
+++ b/src/gridtools/FindSphericalContour.cpp
@@ -39,12 +39,13 @@ class FindSphericalContour : public ContourFindingBase {
 private:
   int rnd;
   unsigned nbins;
-  double offset, increment;
   double min, max;
 public:
   static void registerKeywords( Keywords& keys );
   explicit FindSphericalContour(const ActionOptions&ao);
+  unsigned getNumberOfQuantities() const { return 2; }
   void compute( const unsigned& current, MultiValue& myvals ) const ;
+  void finishAveraging();
 };
 
 PLUMED_REGISTER_ACTION(FindSphericalContour,"FIND_SPHERICAL_CONTOUR")
@@ -55,7 +56,6 @@ void FindSphericalContour::registerKeywords( Keywords& keys ){
   keys.add("compulsory","INNER_RADIUS","the minimum radius on which to look for the contour");
   keys.add("compulsory","OUTER_RADIUS","the outer radius on which to look for the contour");
   keys.add("compulsory","NBINS","1","the number of discrete sections in which to divide the distance between the inner and outer radius when searching for a contour");
-  keys.remove("CLEAR");
 }
 
 FindSphericalContour::FindSphericalContour(const ActionOptions&ao):
@@ -70,9 +70,9 @@ ContourFindingBase(ao)
   log.printf("  expecting to find dividing surface at radii between %f and %f \n",min,max);
   log.printf("  looking for contour in windows of length %f \n", (max-min)/nbins);
   // Set this here so the same set of grid points are used on every turn
-  Random random; rnd = std::floor( npoints*random.RandU01() );
-  offset = 2 / static_cast<double>( npoints );
-  increment = pi*( 3 - sqrt(5) );
+  std::string vstring = "TYPE=fibonacci COMPONENTS=" + getLabel() + " COORDINATES=x,y,z";
+  createGrid( "grid", vstring ); mygrid->setNoDerivatives();
+  setAveragingAction( mygrid, true ); mygrid->setupFibonacciGrid( npoints );
 
   checkRead();
   // Create a task list
@@ -85,19 +85,9 @@ ContourFindingBase(ao)
 void FindSphericalContour::compute( const unsigned& current, MultiValue& myvals ) const {
   // Generate contour point on inner sphere
   std::vector<double> contour_point(3), direction(3), der(3), tmp(3);
-  contour_point[1] = ((current*offset) - 1) + (offset/2);
-  double r = sqrt( 1 - pow(contour_point[1],2) );
-  double phi = ((current + rnd)%getFullNumberOfTasks())*increment;
-  contour_point[0] = r*cos(phi);
-  contour_point[2] = r*sin(phi);
-
-  // normalize direction vector 
-  double norm=0;
-  for(unsigned j=0;j<3;++j) norm+=contour_point[j]*contour_point[j];
-  norm = sqrt( norm );
-  for(unsigned j=0;j<3;++j) direction[j] = contour_point[j] / norm;
- 
-  // Now set up direction as vector from inner sphere to outer sphere
+  // Retrieve this contour point from grid
+  mygrid->getGridPointCoordinates( current, direction );
+  // Now setup contour point on inner sphere
   for(unsigned j=0;j<3;++j){
      contour_point[j] = min*direction[j];
      direction[j] = (max-min)*direction[j] / static_cast<double>(nbins);
@@ -109,13 +99,25 @@ void FindSphericalContour::compute( const unsigned& current, MultiValue& myvals
      double val2 = getDifferenceFromContour( tmp, der ); 
      if( val1*val2<0 ){
          findContour( direction, contour_point );
-         for(unsigned j=0;j<3;++j) myvals.setValue( 1+j, contour_point[j] );
-         found=true; break;
+         double norm=0; for(unsigned j=0;j<3;++j) norm += contour_point[j]*contour_point[j];
+         myvals.setValue( 1, sqrt(norm) ); found=true; break;
      }   
      for(unsigned j=0;j<3;++j) contour_point[j] = tmp[j]; 
   } 
   if( !found ) error("range does not bracket the dividing surface");
 }
 
+void FindSphericalContour::finishAveraging(){
+   std::vector<double> point( 3 );
+   of.printf("%u\n",mygrid->getNumberOfPoints());  
+   of.printf("Points found on isocontour\n");
+   for(unsigned i=0;i<mygrid->getNumberOfPoints();++i){
+       mygrid->getGridPointCoordinates( i, point );
+       of.printf("X"); double val=mygrid->getGridElement( i, 0 );
+       for(unsigned j=0;j<3;++j) of.printf( (" " + fmt_xyz).c_str(), val*lenunit*point[j] );
+       of.printf("\n");
+   } 
+}
+
 }
 }
diff --git a/src/gridtools/GridVessel.cpp b/src/gridtools/GridVessel.cpp
index 08b9bb19d29ada0f93dcc5fd51ccfd321a6a1272..4ed72bf1fb3c8d03da2fb417204f585a57f60956 100644
--- a/src/gridtools/GridVessel.cpp
+++ b/src/gridtools/GridVessel.cpp
@@ -21,6 +21,7 @@
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
 #include "GridVessel.h"
 #include "vesselbase/ActionWithVessel.h"
+#include "tools/Random.h"
 #include "tools/Tools.h"
 
 namespace PLMD {
@@ -28,6 +29,7 @@ namespace gridtools {
 
 void GridVessel::registerKeywords( Keywords& keys ){
   AveragingVessel::registerKeywords( keys );
+  keys.add("compulsory","TYPE","flat","how the grid points are being generated");
   keys.add("compulsory","COMPONENTS","the names of the components in the vector");
   keys.add("compulsory","COORDINATES","the names of the coordinates of the grid");
   keys.add("compulsory","PBC","is the grid periodic in each direction or not");
@@ -41,13 +43,27 @@ noderiv(false),
 npoints(0),
 wasforced(false)
 {
+  std::string geom; parse("TYPE",geom);
+  if( geom=="flat" ) gtype=flat;
+  else if( geom=="fibonacci" ) gtype=fibonacci;
+  else plumed_merror( geom + " is invalid geometry type");
   std::vector<std::string> compnames; parseVector("COMPONENTS",compnames);
   std::vector<std::string> coordnames; parseVector("COORDINATES",coordnames);
-  dimension=coordnames.size();
-  std::vector<std::string> spbc( dimension ); parseVector("PBC",spbc); 
-  str_min.resize( dimension);  str_max.resize( dimension ); stride.resize( dimension ); 
-  max.resize( dimension ); dx.resize( dimension ); nbin.resize( dimension ); min.resize( dimension );
- 
+  if( gtype==flat ){
+      dimension=coordnames.size();
+      std::vector<std::string> spbc( dimension ); parseVector("PBC",spbc);
+      str_min.resize( dimension);  str_max.resize( dimension ); stride.resize( dimension ); 
+      max.resize( dimension ); dx.resize( dimension ); nbin.resize( dimension ); min.resize( dimension ); 
+       pbc.resize( dimension );
+       for(unsigned i=0;i<dimension;++i){
+           if( spbc[i]=="F" ) pbc[i]=false;
+           else if( spbc[i]=="T" ) pbc[i]=true;
+           else plumed_error();
+       }
+  } else if( gtype==fibonacci ){
+      if( coordnames.size()!=3 ) error("cannot generate fibonacci grid points on surface of sphere if not 3 input coordinates");
+      dimension=3; 
+  }
 
   unsigned n=0; nper=compnames.size()*( 1 + coordnames.size() );
   arg_names.resize( coordnames.size() + compnames.size()*( 1 + coordnames.size() ) );
@@ -56,13 +72,6 @@ wasforced(false)
       arg_names[n]=compnames[i]; n++;
       for(unsigned j=0;j<coordnames.size();++j){ arg_names[n] = "d" + compnames[i] + "_" + coordnames[j]; n++; }
   }
-
-  pbc.resize( dimension ); 
-  for(unsigned i=0;i<dimension;++i){
-      if( spbc[i]=="F" ) pbc[i]=false;   
-      else if( spbc[i]=="T" ) pbc[i]=true;
-      else plumed_error();
-  }
 }
 
 void GridVessel::setNoDerivatives(){
@@ -78,7 +87,7 @@ void GridVessel::setNoDerivatives(){
 void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
                             const std::vector<unsigned>& binsin, const std::vector<double>& spacing ){
   plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
-  plumed_assert( (spacing.size()==dimension || binsin.size()==dimension) );
+  plumed_assert( gtype==flat && (spacing.size()==dimension || binsin.size()==dimension) );
 
   npoints=1; bounds_set=true;
   for(unsigned i=0;i<dimension;++i){
@@ -101,21 +110,35 @@ void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vec
   resize();  // Always resize after setting new bounds as grid size may have have changed 
 } 
 
+void GridVessel::setupFibonacciGrid( const unsigned& np ){
+  bounds_set=true;
+  npoints = np; fib_increment = pi*( 3 - sqrt(5) );
+  fib_offset = 2 / static_cast<double>( npoints ); 
+  Random random; fib_rnd = std::floor( npoints*random.RandU01() );
+  resize();
+}
+
 std::string GridVessel::description(){
   if( !bounds_set ) return "";
  
-  std::string des="grid of "; std::string num;
-  for(unsigned i=0;i<dimension-1;++i){
-      Tools::convert( nbin[i], num );
-      des += num + " X ";
+  std::string des;
+  if( gtype==flat ){
+      des="grid of "; std::string num;
+      for(unsigned i=0;i<dimension-1;++i){
+          Tools::convert( nbin[i], num );
+          des += num + " X ";
+      }
+      Tools::convert( nbin[dimension-1], num );
+      des += num + " equally spaced points between (";
+      for(unsigned i=0;i<dimension-1;++i) des += str_min[i] + ",";
+      Tools::convert( nbin[dimension-1], num );
+      des += str_min[dimension-1] + ") and (";
+      for(unsigned i=0;i<dimension-1;++i) des += str_max[i] + ",";
+      des += str_max[dimension-1] + ")";
+  } else if( gtype==fibonacci ){
+     std::string num; Tools::convert( npoints, num );
+     des += "fibonacci grid of " + num + " points on spherical surface"; 
   }
-  Tools::convert( nbin[dimension-1], num );
-  des += num + " equally spaced points between (";
-  for(unsigned i=0;i<dimension-1;++i) des += str_min[i] + ",";
-  Tools::convert( nbin[dimension-1], num );
-  des += str_min[dimension-1] + ") and (";
-  for(unsigned i=0;i<dimension-1;++i) des += str_max[i] + ",";
-  des += str_max[dimension-1] + ")";
   return des;
 }
 
@@ -127,7 +150,7 @@ void GridVessel::resize(){
 }
 
 unsigned GridVessel::getIndex( const std::vector<unsigned>& indices ) const {
-  plumed_dbg_assert( bounds_set && indices.size()==dimension );
+  plumed_dbg_assert( gtype==flat && bounds_set && indices.size()==dimension );
   // indices are flattended using a column-major order
   unsigned index=indices[dimension-1];
   for(unsigned i=dimension-1;i>0;--i){
@@ -137,7 +160,7 @@ unsigned GridVessel::getIndex( const std::vector<unsigned>& indices ) const {
 }
 
 void GridVessel::getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const { 
-  plumed_dbg_assert( bounds_set && point.size()==dimension && indices.size()==dimension );
+  plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension && indices.size()==dimension );
   for(unsigned i=0;i<dimension;++i){
       indices[i]=std::floor( (point[i] - min[i])/dx[i] );
       if( pbc[i] ) indices[i]=indices[i]%nbin[i];
@@ -145,14 +168,13 @@ void GridVessel::getIndices( const std::vector<double>& point, std::vector<unsig
 }
 
 unsigned GridVessel::getIndex( const std::vector<double>& point ) const {
-  plumed_dbg_assert( bounds_set && point.size()==dimension );
+  plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension );
   std::vector<unsigned> indices(dimension); getIndices( point, indices );
   return getIndex( indices );
 }
 
 void GridVessel::convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const {
- unsigned kk=index;
- indices[0]=index%nnbin[0];
+ plumed_dbg_assert( gtype==flat ); unsigned kk=index; indices[0]=index%nnbin[0];
  for(unsigned i=1;i<dimension-1;++i){
     kk=(kk-indices[i-1])/nnbin[i-1];
     indices[i]=kk%nnbin[i];
@@ -164,17 +186,30 @@ void GridVessel::convertIndexToIndices( const unsigned& index, const std::vector
 }
 
 void GridVessel::getIndices( const unsigned& index, std::vector<unsigned>& indices ) const {
- convertIndexToIndices( index, nbin, indices );
+ plumed_dbg_assert( gtype==flat ); convertIndexToIndices( index, nbin, indices );
 }
 
 void GridVessel::getGridPointCoordinates( const unsigned& ipoint , std::vector<double>& x ) const {
-  plumed_dbg_assert( x.size()==dimension && ipoint<npoints );
-  std::vector<unsigned> tindices( dimension ); getIndices( ipoint, tindices ); 
-  for(unsigned i=0;i<dimension;++i) x[i] = min[i] + dx[i]*tindices[i];
+  plumed_dbg_assert( bounds_set && x.size()==dimension && ipoint<npoints );
+  if( gtype==flat ){
+      std::vector<unsigned> tindices( dimension ); getIndices( ipoint, tindices ); 
+      for(unsigned i=0;i<dimension;++i) x[i] = min[i] + dx[i]*tindices[i];
+  } else if( gtype==fibonacci ){
+      x[1] = ((ipoint*fib_offset) - 1) + (fib_offset/2);
+      double r = sqrt( 1 -x[1]*x[1] );
+      double phi = ((ipoint+fib_rnd)%npoints)*fib_increment;
+      x[0] = r*cos(phi); x[2] = r*sin(phi);
+      double norm=0; 
+      for(unsigned j=0;j<3;++j) norm+=x[j]*x[j]; 
+      norm = sqrt(norm);
+      for(unsigned j=0;j<3;++j) x[j] = x[j] / norm;
+  } else {
+      plumed_error();
+  }
 }
 
 void GridVessel::getSplineNeighbors( const unsigned& mybox, std::vector<unsigned>& mysneigh ) const {
-  mysneigh.resize( static_cast<unsigned>(pow(2.,dimension)) );
+  plumed_dbg_assert( gtype==flat ); mysneigh.resize( static_cast<unsigned>(pow(2.,dimension)) );
 
   std::vector<unsigned> tmp_indices( dimension );
   std::vector<unsigned> my_indices( dimension );
@@ -221,15 +256,15 @@ void GridVessel::setGridElement( const std::vector<unsigned>& indices, const uns
 }
 
 std::vector<std::string> GridVessel::getMin() const {
-  return str_min;
+  plumed_dbg_assert( gtype==flat ); return str_min;
 }
   
 std::vector<std::string> GridVessel::getMax() const {
-  return str_max;
+  plumed_dbg_assert( gtype==flat ); return str_max;
 }
 
 std::vector<unsigned> GridVessel::getNbin() const {
-  plumed_dbg_assert( bounds_set );
+  plumed_dbg_assert( gtype==flat && bounds_set );
   std::vector<unsigned> ngrid( dimension );
   for(unsigned i=0;i<dimension;++i){
       if( !pbc[i] ) ngrid[i]=nbin[i] - 1;
@@ -240,7 +275,7 @@ std::vector<unsigned> GridVessel::getNbin() const {
 
 void GridVessel::getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
                                unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
-  plumed_dbg_assert( bounds_set && nneigh.size()==dimension );
+  plumed_dbg_assert( gtype==flat && bounds_set && nneigh.size()==dimension );
 
   std::vector<unsigned> indices( dimension );
   for(unsigned i=0;i<dimension;++i) indices[i] = std::floor( (pp[i]-min[i])/dx[i] );
@@ -249,7 +284,7 @@ void GridVessel::getNeighbors( const std::vector<double>& pp, const std::vector<
 
 void GridVessel::getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh, 
                                unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
-  plumed_dbg_assert( bounds_set && nneigh.size()==dimension );
+  plumed_dbg_assert( gtype==flat && bounds_set && nneigh.size()==dimension );
 
   unsigned num_neigh=1; std::vector<unsigned> small_bin( dimension );
   for(unsigned i=0;i<dimension;++i){
@@ -279,28 +314,32 @@ void GridVessel::getNeighbors( const std::vector<unsigned>& indices, const std::
 }
 
 void GridVessel::setCubeUnits( const double& units ){
-  cube_units=units;
+  plumed_dbg_assert( gtype==flat ); cube_units=units;
 }
 
 double GridVessel::getCubeUnits() const {
-  return cube_units;
+  plumed_dbg_assert( gtype==flat ); return cube_units;
 }
 
 std::string GridVessel::getInputString() const {
   std::string mstring="COORDINATES="+arg_names[0];
   for(unsigned i=1;i<dimension;++i) mstring+="," + arg_names[i];
-  mstring += " PBC=";
-  if( pbc[0] ) mstring +="T";
-  else mstring +="F";
-  for(unsigned i=1;i<dimension;++i){
-     if( pbc[i] ) mstring +=",T";
-     else mstring +=",F";
+  if( gtype==flat ){
+      mstring += " TYPE=flat PBC=";
+      if( pbc[0] ) mstring +="T";
+      else mstring +="F";
+      for(unsigned i=1;i<dimension;++i){
+         if( pbc[i] ) mstring +=",T";
+         else mstring +=",F";
+      }
+  } else if( gtype==fibonacci ){
+      mstring += " TYPE=fibonacci";
   }
   return mstring;
 }
 
 double GridVessel::getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const {
-  plumed_dbg_assert( der.size()==dimension && !noderiv && ind<getNumberOfComponents() );
+  plumed_dbg_assert( gtype==flat && der.size()==dimension && !noderiv && ind<getNumberOfComponents() );
 
   double X,X2,X3,value=0; der.assign(der.size(),0.0);
   std::vector<double> fd(dimension);
diff --git a/src/gridtools/GridVessel.h b/src/gridtools/GridVessel.h
index 7fe9953b10cb52dcab00169ed7d511edb7d0021e..2e2c08a2c073b5d02726681553ed9cf73005c2f4 100644
--- a/src/gridtools/GridVessel.h
+++ b/src/gridtools/GridVessel.h
@@ -34,10 +34,15 @@ class GridVessel : public vesselbase::AveragingVessel {
 friend class ActionWithInputGrid;
 friend class DumpGrid;
 private:
+/// The way that grid points are constructed
+ enum {flat,fibonacci} gtype;
 /// Have the minimum and maximum for the grid been set
  bool bounds_set;
 /// The number of points in the grid
  unsigned npoints;
+/// Stuff for fibonacci grids
+ int fib_rnd;
+ double fib_offset, fib_increment;
 /// Units for Gaussian Cube file
  double cube_units;
 /// This flag is used to check if the user has created a valid input 
@@ -83,8 +88,12 @@ public:
  explicit GridVessel( const vesselbase::VesselOptions& );
 /// Remove the derivatives 
  void setNoDerivatives();
+/// Get the type of grid we are using
+ std::string getType() const ;
 /// Set the minimum and maximum of the grid
  virtual void setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax, const std::vector<unsigned>& nbins, const std::vector<double>& spacing );
+/// Setup the grid if it is a fibonacci grid on the surface of a sphere
+ void setupFibonacciGrid( const unsigned& np );
 /// Get a description of the grid to output to the log
  std::string description();
 /// Convert an index into indices
@@ -184,13 +193,20 @@ unsigned GridVessel::getNumberOfPoints() const {
 
 inline
 const std::vector<double>& GridVessel::getGridSpacing() const {
+  if( gtype==flat ) return dx;
+  plumed_merror("dont understand what spacing means for spherical grids");
   return dx;
 }
 
 inline
 double GridVessel::getCellVolume() const {
-  double myvol=1.0; for(unsigned i=0;i<dimension;++i) myvol *= dx[i];
-  return myvol;
+  if( gtype==flat ){
+      double myvol=1.0; for(unsigned i=0;i<dimension;++i) myvol *= dx[i];
+      return myvol;
+  } else {
+      plumed_merror("cell volume for surface grid not worked out yet"); 
+      return 0.0; 
+  }
 }
 
 inline
@@ -200,6 +216,7 @@ unsigned GridVessel::getDimension() const {
 
 inline
 bool GridVessel::isPeriodic( const unsigned& i ) const {
+  plumed_dbg_assert( gtype==flat );
   return pbc[i];
 }
 
@@ -216,6 +233,7 @@ unsigned GridVessel::getNumberOfComponents() const {
 
 inline
 double GridVessel::getGridExtent( const unsigned& i ) const {
+  plumed_dbg_assert( gtype==flat );
   return max[i] - min[i];
 }
 
@@ -232,6 +250,7 @@ bool GridVessel::inactive( const unsigned& ip ) const {
 
 inline
 const std::vector<unsigned>& GridVessel::getStride() const {
+  plumed_dbg_assert( gtype==flat );
   return stride;
 }
 
@@ -240,6 +259,13 @@ unsigned GridVessel::getNumberOfBufferPoints() const {
   return npoints;
 }
 
+inline
+std::string GridVessel::getType() const {
+  if( gtype==flat ) return "flat";
+  else if( gtype==fibonacci ) return "fibonacci";
+  plumed_error();
+}
+
 }
 }
 #endif