double
  OPMSD::calcStructMSD(const Topology& Itop) const
  {
    //Required to get the correct results
    Sim->dynamics->updateAllParticles();

    double acc = 0.0;
    for (const shared_ptr<IDRange>& molRange : Itop.getMolecules())
      {
	Vector origPos{0,0,0}, currPos{0,0,0};
	double totmass = 0.0;
	for (const unsigned long& ID : *molRange)
	  {
	    double pmass = Sim->species[Sim->particles[ID]]->getMass(ID);

	    totmass += pmass;
	    currPos += Sim->particles[ID].getPosition() * pmass;
	    origPos += initPos[ID] * pmass;
	  }
      
	currPos /= totmass;
	origPos /= totmass;

	acc += (currPos - origPos).nrm2();
      }

    acc /= Itop.getMoleculeCount();
  
    return acc;
  }

Action::RetType Action_CreateCrd::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  // Keywords
  Topology* parm = init.DSL().GetTopology( actionArgs );
  if (parm == 0) {
    mprinterr("Error: createcrd: No parm files loaded.\n");
    return Action::ERR;
  }
  pindex_ = parm->Pindex();
  check_ = !actionArgs.hasKey("nocheck");
  // DataSet
  std::string setname = actionArgs.GetStringNext();
  if (setname == "_DEFAULTCRD_") {
    // Special case: Creation of COORDS DataSet has been requested by an
    //               analysis and should already be present.
    coords_ = (DataSet_Coords_CRD*)init.DSL().FindSetOfType(setname, DataSet::COORDS);
  } else 
    coords_ = (DataSet_Coords_CRD*)init.DSL().AddSet(DataSet::COORDS, setname, "CRD");
  if (coords_ == 0) return Action::ERR;
  // Do not set topology here since it may be modified later.

  mprintf("    CREATECRD: Saving coordinates from Top %s to \"%s\"\n",
          parm->c_str(), coords_->legend());
  if (!check_)
    mprintf("\tNot strictly enforcing that all frames have same # atoms.\n");
# ifdef MPI
  if (init.TrajComm().Size() > 1)
    mprintf("Warning: Synchronization of COORDS data sets over multiple threads is\n"
            "Warning:   experimental and may be slower than reading in via a single\n"
            "Warning:   thread. Users are encouraged to run benchmarks before\n"
            "Warning:   extensive usage.\n");
# endif
  return Action::OK;
}

/** Write file containing only cut atoms and energies as charges. */
int Action_Pairwise::WriteCutFrame(int frameNum, Topology const& Parm, AtomMask const& CutMask, 
                                   Darray const& CutCharges,
                                   Frame const& frame, std::string const& outfilename) 
{
  if (CutMask.Nselected() != (int)CutCharges.size()) {
    mprinterr("Error: WriteCutFrame: # of charges (%u) != # mask atoms (%i)\n",
              CutCharges.size(), CutMask.Nselected());
    return 1;
  }
  Frame CutFrame(frame, CutMask);
  Topology* CutParm = Parm.modifyStateByMask( CutMask );
  if (CutParm == 0) return 1;
  // Set new charges
  for (int i = 0; i != CutParm->Natom(); i++)
    CutParm->SetAtom(i).SetCharge( CutCharges[i] );
  int err = 0;
  Trajout_Single tout;
  if (tout.PrepareTrajWrite(outfilename, "multi", CutParm, CoordinateInfo(), 1,
                            TrajectoryFile::MOL2FILE))
  {
    mprinterr("Error: Could not set up cut mol2 file %s\n", outfilename.c_str());
    err = 1;
  } else {
    tout.WriteSingle(frameNum, CutFrame);
    tout.EndTraj();
  }
  delete CutParm;
  return err;
}

/** Based on the given atom mask expression determine what molecules are
  * selected by the mask.
  * \return A list of atom pairs that mark the beginning and end of each
  *         selected molecule.
  */
Action_AutoImage::pairList
  Action_AutoImage::SetupAtomRanges( Topology const& currentParm, std::string const& maskexpr )
{
  pairList imageList;
  CharMask Mask1( maskexpr.c_str() );

  if (currentParm.SetupCharMask( Mask1 )) return imageList;
  if (Mask1.None()) return imageList;
  for (Topology::mol_iterator mol = currentParm.MolStart(); mol != currentParm.MolEnd(); mol++)
  {
    int firstAtom = mol->BeginAtom();
    int lastAtom = mol->EndAtom();
    // Check that each atom in the range is in Mask1
    bool rangeIsValid = true;
    for (int atom = firstAtom; atom < lastAtom; ++atom) {
      if (!Mask1.AtomInCharMask(atom)) {
        rangeIsValid = false;
        break;
      }
    }
    if (rangeIsValid) {
      imageList.push_back( firstAtom );
      imageList.push_back( lastAtom );
    }
  }
  mprintf("\tMask [%s] corresponds to %zu molecules\n", Mask1.MaskString(), imageList.size()/2);
  return imageList;
}

// ParmFile::WriteTopology()
int ParmFile::WriteTopology(Topology const& Top, FileName const& fnameIn, 
                            ArgList const& argListIn, ParmFormatType fmtIn, int debugIn)
{
  parmName_ = fnameIn;
  ArgList argIn = argListIn;
  ParmFormatType fmt = fmtIn;
  if (fmt == UNKNOWN_PARM) {
    // Check arg list to see if format specified.
    fmt = (ParmFormatType)FileTypes::GetFormatFromArg(PF_KeyArray, argIn, UNKNOWN_PARM);
    // If still UNKNOWN check file extension. Default to AMBERPARM
    if (fmt == UNKNOWN_PARM)
      fmt = (ParmFormatType)FileTypes::GetTypeFromExtension(PF_KeyArray, parmName_.Ext(),
                                                            AMBERPARM);
  }
  ParmIO* parmio = (ParmIO*)FileTypes::AllocIO(PF_AllocArray, fmt, true);
  if (parmio == 0) return 1;
  parmio->SetDebug( debugIn );
  parmio->processWriteArgs( argIn );
  mprintf("\tWriting topology %i (%s) to '%s' with format %s\n", Top.Pindex(),
          Top.c_str(), parmName_.full(), FileTypes::FormatDescription(PF_AllocArray, fmt));
  int err = parmio->WriteParm( parmName_.Full(), Top );
  delete parmio;
  if (err != 0 ) {
    mprinterr("Error: writing topology file '%s'\n", parmName_.full());
    return 1;
  }
  return 0;
}

int Parm_PDB::ReadParm(std::string const& fname, Topology &TopIn) {
  PDBfile infile;
  double XYZ[3];
  int current_res = 0;
  int last_res = -1;
  if (infile.OpenRead(fname)) return 1;
  // Loop over PDB records 
  while ( infile.NextLine() != 0 ) {
    if (infile.IsPDBatomKeyword()) {
      // If this is an ATOM / HETATM keyword, add to topology
      infile.pdb_XYZ(XYZ);
      NameType pdbresname = infile.pdb_Residue( current_res );
      TopIn.AddTopAtom(infile.pdb_Atom(), pdbresname, current_res, last_res, XYZ);
    } else if (infile.IsPDB_TER() || infile.IsPDB_END()) {
      // Indicate end of molecule for TER/END. Finish if END.
      TopIn.StartNewMol();
      if (infile.IsPDB_END()) break;
    }
  }
  // If Topology name not set with TITLE etc, use base filename.
  // TODO: Read in title.
  std::string pdbtitle;
  TopIn.SetParmName( pdbtitle, infile.Filename().Base() );

  infile.CloseFile();
  return 0;
}

void Node::SetDistances(const Topology &top)
{
	distances.reserve(top.GetSize());
	for(int i = 0; i < top.GetSize(); i++)
		distances[i] = top.GetNumHops(_logical, i);

}

/** Determine VDW long range correction prefactor. */
void Ewald::Setup_VDW_Correction(Topology const& topIn, AtomMask const& maskIn) {
  Vdw_Recip_term_ = 0.0;
  NB_ = static_cast<NonbondParmType const*>( &(topIn.Nonbond()) );
  if (!NB_->HasNonbond()) {
    mprintf("Warning: '%s' has no nonbonded parameters. Cannot calculate VDW correction.\n",
            topIn.c_str());
    return;
  }
  // Count the number of each unique nonbonded type.
  Iarray N_vdw_type( NB_->Ntypes(), 0 );
  for (AtomMask::const_iterator atm = maskIn.begin(); atm != maskIn.end(); ++atm)
    N_vdw_type[ topIn[*atm].TypeIndex() ]++;
  if (debug_ > 0) {
    mprintf("DEBUG: %zu VDW types.\n", N_vdw_type.size());
    for (Iarray::const_iterator it = N_vdw_type.begin(); it != N_vdw_type.end(); ++it)
      mprintf("\tType %li = %i\n", it-N_vdw_type.begin(), *it);
  }
  // Determine correction term from types and LJ B parameters
  for (unsigned int itype = 0; itype != N_vdw_type.size(); itype++)
  {
    unsigned int offset = N_vdw_type.size() * itype;
    for (unsigned int jtype = 0; jtype != N_vdw_type.size(); jtype++)
    {
      unsigned int idx = offset + jtype;
      int nbidx = NB_->NBindex()[ idx ];
      if (nbidx > -1)
        Vdw_Recip_term_ += N_vdw_type[itype] * N_vdw_type[jtype] * NB_->NBarray()[ nbidx ].B();
    }
  }
}

// Action_CheckStructure::SeparateSetup()
int Action_CheckStructure::SeparateSetup(Topology const& top, Box::BoxType btype, bool checkBonds)
{
  image_.SetupImaging( btype );
  bondList_.clear();
  // Set up masks
  if ( top.SetupIntegerMask( Mask1_ ) ) return 1;
  Mask1_.MaskInfo();
  if (Mask1_.None()) {
    mprinterr("Error: Mask '%s' has no atoms.\n", Mask1_.MaskString());
    return 1;
  }
  if (checkBonds) SetupBondList(Mask1_, top);
  if ( Mask2_.MaskStringSet() ) {
    if (top.SetupIntegerMask( Mask2_ ) ) return 1;
    Mask2_.MaskInfo();
    if (Mask2_.None()) {
      mprinterr("Error: Mask '%s' has no atoms.\n", Mask2_.MaskString());
      return 1;
    }
    int common = Mask1_.NumAtomsInCommon( Mask2_ );
    if (common > 0)
      mprintf("Warning: '%s' has %i atoms in common with '%s'. Some problems may be reported\n"
              "Warning:   more than once.\n", Mask1_.MaskString(), common, Mask2_.MaskString());
    // Outer mask should be the one with the most atoms.
    if ( Mask2_.Nselected() > Mask1_.Nselected() ) {
      OuterMask_ = Mask2_;
      InnerMask_ = Mask1_;
    } else {
      OuterMask_ = Mask1_;
      InnerMask_ = Mask2_;
    }
    if (checkBonds) SetupBondList(Mask2_, top);
  }
  return 0;
}

Channel_Generator_Base::Channel_Generator_Base(int _nin,int _nout,
					       Point * _plist) 
  : nin(_nin), nout(_nout), m_valid(1)
{
  Topology top;
  plist  = new Point[2*(nout+1)];
  int ll = 0;
  top.Copy(_plist,plist,ll);
}

	FilesystemID register_filesystemID( const string & global_unique_identifier ) {
		TopologyObject d = topology->registerObject( pluginTopoObjectID, fsID, global_unique_identifier, fsID);
		if ( ! d ){
			throw IllegalStateError("Could not register filesystem: " + global_unique_identifier);
		}
		topology->setAttribute(d.id(), attrFSNameID, global_unique_identifier);		

		return d.id();
	}

void print_graph_layout(const Graph& g, PositionMap position, const Topology& topology)
{
  typedef typename Topology::point_type Point;
  // Find min/max ranges
  Point min_point = position[*vertices(g).first], max_point = min_point;
  BGL_FORALL_VERTICES_T(v, g, Graph) {
    min_point = topology.pointwise_min(min_point, position[v]);
    max_point = topology.pointwise_max(max_point, position[v]);
  }

 double operator()(const PointDiff& a, const Topology& s) const {
   try {
     detail::generic_interpolator_impl<svp_interpolator,Topology,TimeSpaceType> interp;
     interp.initialize(s.origin(), s.adjust(s.origin(),a), 0.0, s, *t_space, *this);
     return interp.get_minimum_travel_time();
   } catch(optim::infeasible_problem& e) { RK_UNUSED(e);
     return std::numeric_limits<double>::infinity();
   };
 };

	NodeID register_nodeID( const string & hostname ) {
		TopologyObject obj = topology->registerObject( pluginTopoObjectID, hostID, hostname, hostID);
		if ( ! obj ){
			throw IllegalStateError("Could not register node: " + hostname);
		}

		topology->setAttribute( obj.id(), attrNameID, hostname );

		return obj.id();
	}

Topology* Topology::GenLineTopo(int len, int content_num, int content_size, int k, int* cacheSizes)
{
	
	Topology* topo = new Topology();
	topo->router_num = len;
	topo->routers = new Router[topo->router_num];
	topo->content_num = content_num;
	topo->size_of_content = content_size;
	
	topo->server_num = 1;
	topo->servers = new Server[topo->server_num];

	topo->sink_num = 1;
	topo->sinks = new Sink[topo->sink_num];



	//sink-100-0-1- ... - 8-9-server-10
	topo->edges.clear();
	//routers
	for (int router_id = 0; router_id < topo->router_num; router_id++) {
		topo->routers[router_id].Init(router_id, cacheSizes[router_id], k, topo->server_num);
	}
	//server
	for (int server_id = 0; server_id < topo->server_num; server_id++) {
		topo->servers[server_id].Init(SERVER_BASE + server_id);
	}
	//sink
	for (int sink_id = 0; sink_id < topo->sink_num; sink_id++) {
		topo->sinks[sink_id].Init(SINK_BASE + sink_id, &(topo->routers[sink_id]));
	}
	//edges
	for (int router_id = 0; router_id < (topo->router_num - 1); router_id++) {
		Edge* e = new Edge(&topo->routers[router_id], &topo->routers[router_id + 1], router_id, MyRandom::NextDouble());
		topo->routers[router_id].AddEdge(e);
		topo->routers[router_id + 1].AddEdge(e);
		topo->edges.push_back(e);
	}
	for (int server_id = 0; server_id < topo->server_num; server_id++) {
		Edge* e = new Edge(&topo->servers[server_id], &topo->routers[topo->router_num - 1 - server_id],server_id-100,MyRandom::NextDouble());
		topo->servers[server_id].AddEdge(e);
		topo->routers[topo->router_num - 1 - server_id].AddEdge(e);
		topo->edges.push_back(e);
	}
	for (int sink_id = 0; sink_id < topo->sink_num; sink_id++) {
		Edge* e = new Edge(&topo->sinks[sink_id], &topo->routers[sink_id],sink_id + SINK_BASE,MyRandom::NextDouble());
		topo->sinks[sink_id].AddEdge(e);
		topo->routers[sink_id].AddEdge(e);
		topo->edges.push_back(e);
	}

	topo->IssueContentOnServer();
	topo->Announce();
	return topo; 
}

double Dnds::change(void) {

	// update conditional likelihood and transition probability flags
	Topology* t = modelPtr->getActiveTopology();
	t->updateAllCls(true);
	t->updateAllTis(true);
	t->flipAllActiveCls();
	t->flipAllActiveTis();

	// set the maximum distance of the move
	double x = 0.1;
	
	// propose new values for all of the omegas
	bool isValid = true;
	std::vector<double> newVals(numCategories);
	std::vector<double> offsets(numCategories);
	do
		{
		// pick the values for each of the categories
		for (int i=0; i<numCategories; i++)
			offsets[i] = x*(ranPtr->uniformRv()-0.5);
			
		// calculate the new points
		for (int i=0; i<numCategories; i++)
			newVals[i] = dndsVals[i] + offsets[i];
		
		// decide if the new point is valid
		isValid = true;
		for (int i=0; i<numCategories; i++)
			{
			if (newVals[i] < 0.0)
				{
				isValid = false;
				break;
				}
			for (int j=i+1; j<numCategories; j++)
				{
				if (newVals[i] > newVals[j])
					{
					isValid = false;
					break;
					}
				}
			if (isValid == false)
				break;
			}
	
		} while (isValid == false);
	
	// commit the change
	for (int i=0; i<numCategories; i++)
		dndsVals[i] = newVals[i];
	
	return 0.0;
}

	DeviceID register_deviceID( NodeID id, const string & local_unique_identifier ) {
		TopologyObject d = topology->registerObject( id, deviceID, local_unique_identifier, deviceID);
		if ( ! d ){
			throw IllegalStateError("Could not register device: " + local_unique_identifier);
		}

		topology->setAttribute( d.id(), attrNodeID, id );
		topology->setAttribute( d.id(), attrDeviceNameID, local_unique_identifier );

		return d.id();
	}

	UniqueComponentActivityID register_activityID( UniqueInterfaceID id, const string & name ) {
		TopologyObject d = topology->registerObject( id, activityID, name, activityID);
		if ( ! d ){
			throw IllegalStateError("Could not register name: " + name);
		}

		topology->setAttribute(d.id(), attrNameID, name);
		topology->setAttribute(d.id(), attrUNIDID, id);

		return d.id();
	}

	UniqueInterfaceID register_interfaceID( const string & interface, const string & implementation ) {
		TopologyObject d = topology->registerObjectByPath(  {{"Interface", interface, "Interface"}, {"Implementation", implementation, "Implementation"} }, pluginTopoObjectID );
		if ( ! d ){
			throw IllegalStateError("Could not register interface: " + interface + " " + implementation);
		}
		const TopologyObjectId objID = d.id();

		topology->setAttribute(objID, attrImplementationNameID, implementation);
		topology->setAttribute(objID, attrInterfaceNameID, interface);

		return d.id();
	}

/** Sets the temporary charge array and makes sure that we have the necessary
  * parameters in our topology to calculate nonbonded energy terms
  */
int Action_LIE::SetupParms(Topology const& ParmIn) {
  if (!ParmIn.Nonbond().HasNonbond()) {
    mprinterr("Error: Topology does not have LJ information.\n");
    return 1;
  }
  // Store the charges
  atom_charge_.clear();
  atom_charge_.reserve( ParmIn.Natom() );
  for (Topology::atom_iterator atom = ParmIn.begin();
                               atom != ParmIn.end(); ++atom)
    atom_charge_.push_back( atom->Charge() * Constants::ELECTOAMBER / sqrt(dielc_) );
  return 0;
}