void 
  DynNewtonianMC::swapSystem(Dynamics& oDynamics)
  {
#ifdef DYNAMO_DEBUG
    if (dynamic_cast<const DynNewtonianMC*>(&oDynamics) == NULL)
      M_throw() << "Trying to swap Dynamicss with different derived types!";
#endif

    DynNewtonianMC& ol(static_cast<DynNewtonianMC&>(oDynamics));

    std::swap(EnergyPotentialStep, ol.EnergyPotentialStep);
    std::swap(_W, ol._W);
  }

  void 
  GCellsShearing::initialise(size_t nID)
  {
    ID=nID;
   
    if (!(Sim->dynamics.BCTypeTest<BCLeesEdwards>()))
      derr << "You should not use the shearing neighbour list"
	   << " in a system without Lees Edwards BC's" << std::endl;

    if (overlink != 1) M_throw() << "Cannot shear with overlinking yet";

    reinitialise();
  }

  Event 
  IHardSphere::getEvent(const Particle &p1, const Particle &p2) const 
  { 
#ifdef DYNAMO_DEBUG
    if (!Sim->dynamics->isUpToDate(p1))
      M_throw() << "Particle 1 is not up to date: ID1=" << p1.getID() << ", ID2=" << p2.getID() << ", delay1=" << Sim->dynamics->getParticleDelay(p1);
  
    if (!Sim->dynamics->isUpToDate(p2))
      M_throw() << "Particle 2 is not up to date: ID1=" << p1.getID() << ", ID2=" << p2.getID() << ", delay2=" << Sim->dynamics->getParticleDelay(p2);

    if (p1 == p2)
      M_throw() << "You shouldn't pass p1==p2 events to the interactions!";
#endif 

    const double d = _diameter->getProperty(p1, p2);
    const double dt = Sim->dynamics->SphereSphereInRoot(p1, p2, d);

    if (dt != HUGE_VAL)
      return Event(p1, dt, INTERACTION, CORE, ID, p2);
  
    return Event(p1, HUGE_VAL, INTERACTION, NONE, ID, p2);
  }

      /*! \brief Destroys any OpenGL resources associated with this
       * object.
       */
      inline void deinit() 
      {
#ifdef MAGNET_DEBUG
	if (_cl_buffer_acquired)
	  M_throw() << "Deinitialising a buffer which is acquired by the OpenCL system!";
#endif
	_cl_handle = ::cl::BufferGL();
	_cl_handle_init = false;
	if (_size)
	  glDeleteBuffersARB(1, &_buffer);
	_context = NULL;
	_size = 0;
      }

  void
  IThinThread::runEvent(Particle& p1, Particle& p2, const IntEvent& iEvent)
  {
    ++Sim->eventCount;

    double d = (_diameter->getProperty(p1.getID())
		+ _diameter->getProperty(p2.getID())) * 0.5;
    double d2 = d * d;

    double e = (_e->getProperty(p1.getID())
		+ _e->getProperty(p2.getID())) * 0.5;

    double l = (_lambda->getProperty(p1.getID())
		+ _lambda->getProperty(p2.getID())) * 0.5;
    double ld2 = d * l * d * l;

    double wd = (_wellDepth->getProperty(p1.getID())
		 + _wellDepth->getProperty(p2.getID())) * 0.5;
    switch (iEvent.getType())
      {
      case CORE:
	{
	  PairEventData retVal(Sim->dynamics->SmoothSpheresColl(iEvent, e, d2, CORE));
	  IntEvent event(iEvent);
	  if (!isCaptured(p1, p2))
	    {
	      event.setType(STEP_IN);
	      retVal.setType(STEP_IN);
	      ICapture::add(p1, p2);
	    }
	  
	  (*Sim->_sigParticleUpdate)(retVal);
	  Sim->ptrScheduler->fullUpdate(p1, p2);
	  for (shared_ptr<OutputPlugin> & Ptr : Sim->outputPlugins)
	    Ptr->eventUpdate(event, retVal);
	  break;
	}
      case STEP_OUT:
	{
	  PairEventData retVal(Sim->dynamics->SphereWellEvent(iEvent, -wd, ld2, 0));
	  if (retVal.getType() != BOUNCE) ICapture::remove(p1, p2);
	  (*Sim->_sigParticleUpdate)(retVal);
	  Sim->ptrScheduler->fullUpdate(p1, p2);
	  for (shared_ptr<OutputPlugin> & Ptr : Sim->outputPlugins)
	    Ptr->eventUpdate(iEvent, retVal);
	  break;
	}
      default:
	M_throw() << "Unknown collision type";
      } 
  }

int
DynNewtonian::getSquareCellCollision3(const Particle& part, const Vector & origin, const Vector & width) const
{
    Vector  rpos(part.getPosition() - origin);
    Vector  vel(part.getVelocity());

    Sim->BCs->applyBC(rpos, vel);

    int retVal(0);
    double time(HUGE_VAL);

#ifdef DYNAMO_DEBUG
    for (size_t iDim = 0; iDim < NDIM; ++iDim)
        if ((vel[iDim] == 0) && (std::signbit(vel[iDim])))
            M_throw() << "You have negative zero velocities, dont use them."
                      << "\nPlease think of the neighbour lists.";
#endif

    for (size_t iDim = 0; iDim < NDIM; ++iDim)
    {
        double tmpdt = ((vel[iDim] < 0)
                        ? -rpos[iDim]/vel[iDim]
                        : (width[iDim]-rpos[iDim]) / vel[iDim]);

        if (tmpdt < time)
        {
            time = tmpdt;
            retVal = (vel[iDim] < 0) ? -(iDim+1) : (iDim+1);
        }
    }

    if (((retVal < 0) && (vel[abs(retVal)-1] > 0))
            || ((retVal > 0) && (vel[abs(retVal)-1] < 0)))
        M_throw() << "Found an error! retVal " << retVal
                  << " vel is " << vel[abs(retVal)-1];

    return retVal;
}

IntEvent 
IRoughHardSphere::getEvent(const Particle& p1, const Particle& p2) const 
{ 
#ifdef DYNAMO_DEBUG
  if (!Sim->dynamics.getLiouvillean().isUpToDate(p1))
    M_throw() << "Particle 1 is not up to date";
  
  if (!Sim->dynamics.getLiouvillean().isUpToDate(p2))
    M_throw() << "Particle 2 is not up to date";
#endif

#ifdef DYNAMO_DEBUG
  if (p1 == p2)
    M_throw() << "You shouldn't pass p1==p2 events to the interactions!";
#endif 

  CPDData colldat(*Sim, p1, p2);

  double d2 = (_diameter->getProperty(p1.getID())
	       + _diameter->getProperty(p2.getID())) * 0.5;
  d2 *= d2;

  if (Sim->dynamics.getLiouvillean()
      .SphereSphereInRoot(colldat, d2,
			  p1.testState(Particle::DYNAMIC), p2.testState(Particle::DYNAMIC)))
    {
#ifdef DYNAMO_OverlapTesting
      if (Sim->dynamics.getLiouvillean().sphereOverlap(colldat, d2))
	M_throw() << "Overlapping particles found" 
		  << ", particle1 " << p1.getID() << ", particle2 " 
		  << p2.getID() << "\nOverlap = " << (sqrt(colldat.r2) - sqrt(d2))/Sim->dynamics.units().unitLength();
#endif

      return IntEvent(p1, p2, colldat.dt, CORE, *this);
    }
  
  return IntEvent(p1,p2,HUGE_VAL, NONE, *this);  
}

  void 
  SpSphericalTop::operator<<(const magnet::xml::Node& XML)
  {
    SpPoint::operator<<(XML);

    try {
      inertiaConstant 
	= XML.getAttribute("InertiaConstant").as<double>() * Sim->dynamics.units().unitArea();
    } 
    catch (boost::bad_lexical_cast &)
      {
	M_throw() << "Failed a lexical cast in CSSphericalTop";
      }
  }

  void 
  SpFixedCollider::operator<<(const magnet::xml::Node& XML)
  {
    range = std::tr1::shared_ptr<CRange>(CRange::getClass(XML, Sim));
  
    try {
      spName = XML.getAttribute("Name");
      intName = XML.getAttribute("IntName");
    } 
    catch (boost::bad_lexical_cast &)
      {
	M_throw() << "Failed a lexical cast in SpFixedCollider";
      }
  }

  void 
  IRotatedParallelCubes::operator<<(const magnet::xml::Node& XML)
  { 
    if (strcmp(XML.getAttribute("Type"),"RotatedParallelCubes"))
      M_throw() << "Attempting to load RotatedParallelCubes from " 
		<< XML.getAttribute("Type") << " entry";
  
    Interaction::operator<<(XML);
  
    try 
      {
	_diameter = Sim->_properties.getProperty(XML.getAttribute("Diameter"),
						 Property::Units::Length());
	_e = Sim->_properties.getProperty(XML.getAttribute("Elasticity"),
					  Property::Units::Dimensionless());
	intName = XML.getAttribute("Name");
	magnet::math::operator<<(Rotation, XML.getNode("Rotation"));
      }
    catch (boost::bad_lexical_cast &)
      {
	M_throw() << "Failed a lexical cast in CIRotatedParallelCubes";
      }
  }

    double eval() const
    {
      switch (deriv)
	{
	case 0:
	  return (rp | nhat) - ( Sigma + wallnHatPosition());
	case 1:
	  return (vp | nhat) - velnHatWall();
	case 2:
	  return Delta * Omega * Omega * std::cos(Omega * t); 
	default:
	  M_throw() << "Invalid access";
	}
    }

      size_t padding()
      {
	switch(_mode)
	  {
	  case CPU:
	    return 1;
	  case NVIDIA:
	    return 1024;
	  case AMD:
	    return 64 * 256;
	  default:
	    M_throw() << "Functor has not yet been built";
	  }
      }

  void 
  TChain::operator<<(const magnet::xml::Node& XML) 
  {
    Topology::operator<<(XML);

    for (magnet::xml::Node node = XML.findNode("Molecule"); node.valid(); ++node)
      ranges.push_back(shared_ptr<IDRange>(IDRange::getClass(node.getNode("IDRange"), Sim)));
  
    size_t Clength = (*ranges.begin())->size();
    for (const shared_ptr<IDRange>& nRange : ranges)
      if (nRange->size() != Clength)
	M_throw() << "Size mismatch in loading one of the ranges in Chain topology \"" 
		  << _name << "\"";
  }

    IDPairRangeChainEnds(const magnet::xml::Node& XML, const dynamo::Simulation*):
      rangeStart(0),rangeEnd(0), interval(0) 
    { 
      rangeStart = XML.getAttribute("Start").as<size_t>();
      rangeEnd = XML.getAttribute("End").as<size_t>();
      interval = XML.getAttribute("Interval").as<size_t>();

      //Guarrantee that they are ordered
      if (rangeStart > rangeEnd) std::swap(rangeStart, rangeEnd);
  
      if ((rangeEnd - rangeStart + 1) % interval)
	M_throw() << "Length of range does not split into an integer"
		  << " number of intervals";
    }

  void 
  Liouvillean::loadParticleXMLData(const magnet::xml::Node& XML)
  {
    dout << "Loading Particle Data" << std::endl;

    bool outofsequence = false;  
  
    for (magnet::xml::Node node = XML.getNode("ParticleData").fastGetNode("Pt"); 
	 node.valid(); ++node)
      {
	if (!node.hasAttribute("ID")
	    || node.getAttribute("ID").as<size_t>() != Sim->particleList.size())
	  outofsequence = true;
      
	Particle part(node, Sim->particleList.size());
	part.getVelocity() *= Sim->dynamics.units().unitVelocity();
	part.getPosition() *= Sim->dynamics.units().unitLength();
	Sim->particleList.push_back(part);
      }

    if (outofsequence)
      dout << "Particle ID's out of sequence!\n"
	   << "This can result in incorrect capture map loads etc.\n"
	   << "Erase any capture maps in the configuration file so they are regenerated." << std::endl;

    Sim->N = Sim->particleList.size();

    dout << "Particle count " << Sim->N << std::endl;

    if (XML.getNode("ParticleData").hasAttribute("OrientationData"))
      {
	orientationData.resize(Sim->N);
	size_t i(0);
	for (magnet::xml::Node node = XML.getNode("ParticleData").fastGetNode("Pt"); 
	     node.valid(); ++node, ++i)
	  {
	    orientationData[i].orientation << node.getNode("U");
	    orientationData[i].angularVelocity << node.getNode("O");
      
	    double oL = orientationData[i].orientation.nrm();
      
	    if (!(oL > 0.0))
	      M_throw() << "Particle ID " << i 
			<< " orientation vector is zero!";
      
	    //Makes the vector a unit vector
	    orientationData[i].orientation /= oL;
	  }
      }
  }

  std::unique_ptr<IDRange>
  SNeighbourList::getParticleNeighbours(const Particle& part) const
  {
#ifdef DYNAMO_DEBUG
    if (!std::dynamic_pointer_cast<GNeighbourList>(Sim->globals[NBListID]))
      M_throw() << "Not a GNeighbourList!";
#endif

    //Grab a reference to the neighbour list
    const GNeighbourList& nblist(*static_cast<const GNeighbourList*>(Sim->globals[NBListID].get()));
    IDRangeList* range_ptr = new IDRangeList();
    nblist.getParticleNeighbours(part, range_ptr->getContainer());
    return std::unique_ptr<IDRange>(range_ptr);
  }

    double max() const
    {
      switch (deriv)
	{
	case 1:
	  return _f1max;
	case 2:
	  return _f2max;
	case 3:
	  return _f3max;
	default:
	  M_throw() << "Invalid access";
	}
    }

  void 
  ISWSequence::operator<<(const magnet::xml::Node& XML)
  {
    Interaction::operator<<(XML);

    _diameter = Sim->_properties.getProperty(XML.getAttribute("Diameter"), Property::Units::Length());
    _lambda = Sim->_properties.getProperty(XML.getAttribute("Lambda"), Property::Units::Dimensionless());
    
    if (XML.hasAttribute("Elasticity"))
      _e = Sim->_properties.getProperty(XML.getAttribute("Elasticity"), Property::Units::Dimensionless());
    else
      _e = Sim->_properties.getProperty(1.0, Property::Units::Dimensionless());
    
    ICapture::loadCaptureMap(XML);
    
    //Load the sequence
    sequence.clear();
    std::set<size_t> letters;
    
    for (magnet::xml::Node node = XML.getNode("Sequence").findNode("Element");
	 node.valid(); ++node)
      {
	if (node.getAttribute("seqID").as<size_t>() != sequence.size())
	  M_throw() << "Sequence of letters not in order, missing element " << sequence.size();
	
	size_t letter = node.getAttribute("Letter").as<size_t>();
	letters.insert(letter);
	sequence.push_back(letter);
      }
    
    //Initialise all the well depths to 1.0
    alphabet.resize(letters.size());
    
    for (std::vector<double>& vec : alphabet)
      vec.resize(letters.size(), 0.0);
    
    for (magnet::xml::Node node = XML.getNode("Alphabet").findNode("Word");
	 node.valid(); ++node)
      {
	alphabet
	  .at(node.getAttribute("Letter1").as<size_t>())
	  .at(node.getAttribute("Letter2").as<size_t>())
	  = node.getAttribute("Depth").as<double>();
	
	alphabet
	  .at(node.getAttribute("Letter2").as<size_t>())
	  .at(node.getAttribute("Letter1").as<size_t>())
	  = node.getAttribute("Depth").as<double>();
      }
  }

    /*! \brief Method which loads the properties from the XML configuration file.
      \param node A xml Node at the root dynamoconfig Node of the config file.
    */
    inline PropertyStore& operator<<(const magnet::xml::Node& node)
    {
      if (node.hasNode("Properties"))
	for (magnet::xml::Node propNode = node.getNode("Properties").fastGetNode("Property");
	     propNode.valid(); ++propNode)
	  {
	    if (!std::string("PerParticle").compare(propNode.getAttribute("Type")))
	      _namedProperties.push_back(Value(new ParticleProperty(propNode)));
	    else
	      M_throw() << "Unsupported Property type, " << propNode.getAttribute("Type").getValue();
	  }
    
      return *this;
    }

  void 
  SCENBList::getParticleNeighbourhood(const Particle& part, 
				       const GNeighbourList::nbHoodFunc& func) const
  {
#ifdef DYNAMO_DEBUG
    if (!isApplicable(part))
      M_throw() << "This complexNBlist entry ("
		<< name << ") is not valid for this particle (" 
		<< part.getID() << ") yet it is being used anyway!";
#endif

//    static_cast<const GNeighbourList&>(*Sim->globals[nblistID])
//      .getParticleNeighbourhood(part, func);
  }