system_creator.cc

#include <ost/mol/mm/system_creator.hh>

namespace ost{ namespace mol{ namespace mm{

SystemPtr SystemCreator::Create(const TopologyPtr top, 
                                const MMSettingsPtr settings,
                                std::map<FuncType,uint>& mapper){

  uint mapper_index = 0;

  SystemPtr sys(new OpenMM::System);

  //add particles to the system
  std::vector<Real> masses = top->GetMasses();

  //set masses of position constraints to zero
  const std::vector<uint>& position_constraints = top->GetPositionConstraints();
  if(!position_constraints.empty()){
    for(std::vector<uint>::const_iterator i = position_constraints.begin();
        i != position_constraints.end(); ++i){
      masses[*i] = 0.0;
    }
  }  

  for(std::vector<Real>::iterator i = masses.begin(); i != masses.end(); ++i){
    sys->addParticle(*i);
  }

  //parameters of all added interactions will temporarily be stored in here
  std::vector<Real> parameters;

  const std::vector<std::pair<Index<2>,std::vector<Real> > >& harmonic_bonds = top->GetHarmonicBonds();
  if(!harmonic_bonds.empty()){
    OpenMM::HarmonicBondForce& harmonic_bond_force = *new OpenMM::HarmonicBondForce();
    sys->addForce(&harmonic_bond_force);
    mapper[HarmonicBond] = mapper_index++;
    for(std::vector<std::pair<Index<2>,std::vector<Real> > >::const_iterator i = harmonic_bonds.begin();
        i != harmonic_bonds.end(); ++i){
      harmonic_bond_force.addBond(i->first[0],i->first[1],i->second[0],i->second[1]);
    }
  }

  const std::vector<std::pair<Index<3>, std::vector<Real> > >& harmonic_angles = top->GetHarmonicAngles();
  if(!harmonic_angles.empty()){
    OpenMM::HarmonicAngleForce& harmonic_angle_force = *new OpenMM::HarmonicAngleForce();
    sys->addForce(&harmonic_angle_force);
    mapper[HarmonicAngle] = mapper_index++;
    for(std::vector<std::pair<Index<3>,std::vector<Real> > >::const_iterator i = harmonic_angles.begin();
        i != harmonic_angles.end(); ++i){
      harmonic_angle_force.addAngle(i->first[0],i->first[1],i->first[2],
                                    i->second[0],i->second[1]);
    }
  }
  
  const std::vector<std::pair<Index<3>,std::vector<Real> > >& urey_bradley_angles = top->GetUreyBradleyAngles();
  if(!urey_bradley_angles.empty()){
    //urey bradley is a mixture between a harmonic angle and harmonic bond term,
    //that gets split up
    OpenMM::HarmonicAngleForce& urey_angle_force = *new OpenMM::HarmonicAngleForce();
    OpenMM::HarmonicBondForce& urey_bond_force = *new OpenMM::HarmonicBondForce();
    sys->addForce(&urey_angle_force);
    sys->addForce(&urey_bond_force);
    mapper[UreyBradleyAngle] = mapper_index++;
    ++mapper_index;
    for(std::vector<std::pair<Index<3>,std::vector<Real> > >::const_iterator i = urey_bradley_angles.begin();
        i != urey_bradley_angles.end(); ++i){
      urey_angle_force.addAngle(i->first[0],i->first[1],i->first[2],
                                i->second[0],i->second[1]);
      urey_bond_force.addBond(i->first[0],i->first[2],i->second[2],i->second[3]);
    }