From 879363a16c09eeaab3f0151f74c50f3587198d5b Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Sat, 31 Jan 2015 16:50:07 +0100
Subject: [PATCH] make it possible the merge topologies with corresponding
 entities, but also the raw topologies

---
 modules/mol/mm/pymod/export_topology.cc |  17 +-
 modules/mol/mm/src/topology.cc          | 201 ++++++++++++++----------
 modules/mol/mm/src/topology.hh          |  11 ++
 3 files changed, 142 insertions(+), 87 deletions(-)

diff --git a/modules/mol/mm/pymod/export_topology.cc b/modules/mol/mm/pymod/export_topology.cc
index eb7a16fcb..05198b22d 100644
--- a/modules/mol/mm/pymod/export_topology.cc
+++ b/modules/mol/mm/pymod/export_topology.cc
@@ -299,6 +299,14 @@ namespace{
     return VecToList<uint>(return_vec);
   }
 
+  void MergeTop(ost::mol::mm::TopologyPtr top, ost::mol::mm::TopologyPtr other){
+    top->Merge(other);
+  }
+
+  void MergeTopEnt(ost::mol::mm::TopologyPtr top, ost::mol::EntityHandle& ent, 
+                   ost::mol::mm::TopologyPtr other, ost::mol::EntityHandle& other_ent){
+    top->Merge(ent,other,other_ent);
+  }
 }
 
 void export_Topology()
@@ -313,6 +321,7 @@ void export_Topology()
     .def("__init__",make_constructor(&WrapTopologyConstructor))
     .def("Save",&ost::mol::mm::Topology::Save)
     .def("Load",&ost::mol::mm::Topology::Load).staticmethod("Load")
+ 
     //interaction adding functions
     .def("AddHarmonicBond",&ost::mol::mm::Topology::AddHarmonicBond)
     .def("AddHarmonicAngle",&ost::mol::mm::Topology::AddHarmonicAngle)
@@ -329,6 +338,7 @@ void export_Topology()
     .def("ResetExclusions",&ost::mol::mm::Topology::ResetExclusions)
     .def("AddHarmonicPositionRestraint",&ost::mol::mm::Topology::AddHarmonicPositionRestraint,(arg("index"),arg("ref_position"),arg("k"),arg("x_scale")=1.0,arg("y_scale")=1.0,arg("z_scale")=1.0))
     .def("AddHarmonicDistanceRestraint",&ost::mol::mm::Topology::AddHarmonicDistanceRestraint)
+ 
     //single atom parameter getter and setter functions
     .def("SetSigmas",&WrapSetSigmas)
     .def("SetSigma",&ost::mol::mm::Topology::SetSigma)
@@ -352,11 +362,13 @@ void export_Topology()
     .def("GetMass",&ost::mol::mm::Topology::GetMass)
     .def("GetOBCScaling",&ost::mol::mm::Topology::GetOBCScaling)
     .def("GetGBSARadius",&ost::mol::mm::Topology::GetGBSARadius)
+
     //getter and setter functions for nonbonded fudge parameters
     .def("SetFudgeLJ",&ost::mol::mm::Topology::SetFudgeLJ)
     .def("SetFudgeQQ",&ost::mol::mm::Topology::SetFudgeQQ)
     .def("GetFudgeLJ",&ost::mol::mm::Topology::GetFudgeLJ)
     .def("GetFudgeQQ",&ost::mol::mm::Topology::GetFudgeQQ)
+ 
     //getter functions for interaction parameters
     .def("GetHarmonicBondParameters",&WrapGetHarmonicBondParam)
     .def("GetHarmonicAngleParameters",&WrapGetHarmonicAngleParam)
@@ -369,6 +381,7 @@ void export_Topology()
     .def("GetDistanceConstraintParameters",&WrapGetDistanceConstraintParam)
     .def("GetHarmonicPositionRestraintParameters",&WrapGetHarmonicPositionRestraintParam)
     .def("GetHarmonicDistanceRestraintParameters",&WrapGetHarmonicDistanceRestraintParam)
+
     //setter functions for interaction parameters
     .def("SetHarmonicBondParameters",&ost::mol::mm::Topology::SetHarmonicBondParameters)
     .def("SetHarmonicAngleParameters",&ost::mol::mm::Topology::SetHarmonicAngleParameters)
@@ -381,6 +394,7 @@ void export_Topology()
     .def("SetDistanceConstraintParameters",&ost::mol::mm::Topology::SetDistanceConstraintParameters)
     .def("SetHarmonicPositionRestraintParameters",&ost::mol::mm::Topology::SetHarmonicPositionRestraintParameters)
     .def("SetHarmonicDistanceRestraintParameters",&ost::mol::mm::Topology::SetHarmonicDistanceRestraintParameters)
+
     //functions to find interactions certain atoms are involved in
     .def("GetHarmonicBondIndices",&GetHarmonicBondIndices)
     .def("GetHarmonicBondIndices",&GetHarmonicBondIndicesSingleIndex)
@@ -418,7 +432,8 @@ void export_Topology()
     .def("GetNumHarmonicPositionRestraints",&ost::mol::mm::Topology::GetNumHarmonicPositionRestraints)
     .def("GetNumHarmonicDistanceRestraints",&ost::mol::mm::Topology::GetNumHarmonicDistanceRestraints)  
 
-    .def("Merge",&ost::mol::mm::Topology::Merge)                                                           
+    .def("Merge",&MergeTop)
+    .def("Merge",&MergeTopEnt)                                                           
   ;
 
   boost::python::register_ptr_to_python<ost::mol::mm::TopologyPtr>();
diff --git a/modules/mol/mm/src/topology.cc b/modules/mol/mm/src/topology.cc
index d5b238b84..99dc19c47 100644
--- a/modules/mol/mm/src/topology.cc
+++ b/modules/mol/mm/src/topology.cc
@@ -911,99 +911,23 @@ std::vector<uint> Topology::GetHarmonicDistanceRestraintIndices(uint atom_index)
 void Topology::Merge(ost::mol::EntityHandle& ent, TopologyPtr other, 
                      const ost::mol::EntityHandle& other_ent){
 
-  //check whether the particle numbers match
-  if(num_particles_ != static_cast<uint>(ent.GetAtomCount())){
-    throw ost::Error("Num Particles is not consistent with num atoms of provided ent!");
-  }
-
-  if(other->GetNumParticles() != static_cast<uint>(other_ent.GetAtomCount())){
-    throw ost::Error("Num Particles is not consistent with num atoms of provided other_ent!");
-  } 
-
-  //check whether there is chain from the new entity is already present in the actual entity 
-  ost::mol::ChainHandleList other_chains = other_ent.GetChainList();
-  for(ost::mol::ChainHandleList::iterator i = other_chains.begin(), 
-      e = other_chains.end(); i != e; ++i){
-    if(ent.FindChain(i->GetName()).IsValid()){
-      std::stringstream ss;
-      ss << "Chain with name \"" << i->GetName() << "\" from other topology";
-      ss << "is already present in destination topology!";
-      throw ost::Error(ss.str());
-    }
-  }
-  //check whether the fudge parameters are consistent
-  if(other->fudge_lj_ != fudge_lj_ || other->fudge_qq_ != fudge_qq_){
-    throw ost::Error("Expect the fudge parameters of source topology to consistent with the fudge parameters from the destination topology!");
-  }
-
-  if(!charges_.empty()){
-    if(other->charges_.empty()){
-      throw ost::Error("Cannot merge topology without charges into a topology with defined charges!");
-    }
-  }
+  this->CheckEntToAdd(ent,other,other_ent);
+  this->CheckTopToAdd(other);
+  this->MergeEnt(ent, other_ent);
+  this->MergeTop(other);
+}
 
-  if(!sigmas_.empty()){
-    if(other->sigmas_.empty()){
-      throw ost::Error("Cannot merge topology without lj sigmas into a topology with defined sigmas!");
-    }
-  }
+void Topology::Merge(TopologyPtr other){
 
-  if(!epsilons_.empty()){
-    if(other->epsilons_.empty()){
-      throw ost::Error("Cannot merge topology without lj epsilons into a topology with defined epsilons!");
-    }
-  }
+  this->CheckTopToAdd(other);
+  this->MergeTop(other);
+}
 
-  if(!gbsa_radii_.empty()){
-    if(other->gbsa_radii_.empty()){
-      throw ost::Error("Cannot merge topology without gbsa radii into a topology with defined radii!");
-    }
-  }
-
-  if(!obc_scaling_.empty()){
-    if(other->obc_scaling_.empty()){
-      throw ost::Error("Cannot merge topology without obc scaling into a topology with defined scaling!");
-    }
-  }
+void Topology::MergeTop(TopologyPtr other){
 
   uint old_num_particles = num_particles_;
   num_particles_ = num_particles_ + other->GetNumParticles();
 
-  //mapper of hashcode from source atom to index in added_atom list
-  std::map<long,int> index_mapper;
-  ost::mol::AtomHandleList added_atoms;
-
-
-  //let's create an editor to copy over all chains, residues and atoms
-  //from the other topologies entity
-  ost::mol::XCSEditor ed = ent.EditXCS(ost::mol::BUFFERED_EDIT);
-  for(ost::mol::ChainHandleList::iterator i = other_chains.begin(),
-      e = other_chains.end(); i != e; ++i){
-    ost::mol::ResidueHandleList res_list = i->GetResidueList();
-    ost::mol::ChainHandle added_chain = ed.InsertChain(i->GetName());
-
-    for(ost::mol::ResidueHandleList::iterator j = res_list.begin();
-        j != res_list.end(); ++j){
-      ost::mol::AtomHandleList atom_list = j->GetAtomList();
-      ost::mol::ResidueHandle added_residue = ed.AppendResidue(added_chain,*j);
-
-      for(ost::mol::AtomHandleList::iterator k = atom_list.begin();
-          k != atom_list.end(); ++k){
-        index_mapper[k->GetHashCode()] = added_atoms.size();
-        added_atoms.push_back(ed.InsertAtom(added_residue,*k));
-      }
-    }
-  }
-
-  //let's rebuild the connectivity
-  ost::mol::BondHandleList bond_list = other_ent.GetBondList();
-  for(ost::mol::BondHandleList::iterator i = bond_list.begin();
-      i != bond_list.end(); ++i){
-    ed.Connect(added_atoms[index_mapper[i->GetFirst().GetHashCode()]],
-               added_atoms[index_mapper[i->GetSecond().GetHashCode()]]);
-  }
-  ed.UpdateICS();
-
   //let's map over masses
   atom_masses_.resize(old_num_particles + other->GetNumParticles());
   memcpy(&atom_masses_[old_num_particles],&other->atom_masses_[0],other->GetNumParticles() * sizeof(Real));
@@ -1200,6 +1124,111 @@ void Topology::Merge(ost::mol::EntityHandle& ent, TopologyPtr other,
                                       old_num_particles + (*i)[1]));
   }
 
+
+
+}
+
+void Topology::MergeEnt(ost::mol::EntityHandle& ent, const ost::mol::EntityHandle& other_ent){
+
+  //mapper of hashcode from source atom to index in added_atom list
+  std::map<long,int> index_mapper;
+  ost::mol::AtomHandleList added_atoms;
+
+  ost::mol::ChainHandleList other_chains = other_ent.GetChainList();
+
+  //let's create an editor to copy over all chains, residues and atoms
+  //from the other topologies entity
+  ost::mol::XCSEditor ed = ent.EditXCS(ost::mol::BUFFERED_EDIT);
+  for(ost::mol::ChainHandleList::iterator i = other_chains.begin(),
+      e = other_chains.end(); i != e; ++i){
+    ost::mol::ResidueHandleList res_list = i->GetResidueList();
+    ost::mol::ChainHandle added_chain = ed.InsertChain(i->GetName());
+
+    for(ost::mol::ResidueHandleList::iterator j = res_list.begin();
+        j != res_list.end(); ++j){
+      ost::mol::AtomHandleList atom_list = j->GetAtomList();
+      ost::mol::ResidueHandle added_residue = ed.AppendResidue(added_chain,*j);
+
+      for(ost::mol::AtomHandleList::iterator k = atom_list.begin();
+          k != atom_list.end(); ++k){
+        index_mapper[k->GetHashCode()] = added_atoms.size();
+        added_atoms.push_back(ed.InsertAtom(added_residue,*k));
+      }
+    }
+  }
+
+  //let's rebuild the connectivity
+  ost::mol::BondHandleList bond_list = other_ent.GetBondList();
+  for(ost::mol::BondHandleList::iterator i = bond_list.begin();
+      i != bond_list.end(); ++i){
+    ed.Connect(added_atoms[index_mapper[i->GetFirst().GetHashCode()]],
+               added_atoms[index_mapper[i->GetSecond().GetHashCode()]]);
+  }
+  ed.UpdateICS();
+}
+
+
+void Topology::CheckTopToAdd(TopologyPtr other){
+
+  if(other->fudge_lj_ != fudge_lj_ || other->fudge_qq_ != fudge_qq_){
+    throw ost::Error("Expect the fudge parameters of source topology to consistent with the fudge parameters from the destination topology!");
+  }
+
+  if(!charges_.empty()){
+    if(other->charges_.empty()){
+      throw ost::Error("Cannot merge topology without charges into a topology with defined charges!");
+    }
+  }
+
+  if(!sigmas_.empty()){
+    if(other->sigmas_.empty()){
+      throw ost::Error("Cannot merge topology without lj sigmas into a topology with defined sigmas!");
+    }
+  }
+
+  if(!epsilons_.empty()){
+    if(other->epsilons_.empty()){
+      throw ost::Error("Cannot merge topology without lj epsilons into a topology with defined epsilons!");
+    }
+  }
+
+  if(!gbsa_radii_.empty()){
+    if(other->gbsa_radii_.empty()){
+      throw ost::Error("Cannot merge topology without gbsa radii into a topology with defined radii!");
+    }
+  }
+
+  if(!obc_scaling_.empty()){
+    if(other->obc_scaling_.empty()){
+      throw ost::Error("Cannot merge topology without obc scaling into a topology with defined scaling!");
+    }
+  }
+}
+
+void Topology::CheckEntToAdd(ost::mol::EntityHandle& ent, TopologyPtr other,
+                             const ost::mol::EntityHandle& other_ent){
+
+  //check whether the particle numbers match
+  if(num_particles_ != static_cast<uint>(ent.GetAtomCount())){
+    throw ost::Error("Num Particles is not consistent with num atoms of provided ent!");
+  }
+
+  if(other->GetNumParticles() != static_cast<uint>(other_ent.GetAtomCount())){
+    throw ost::Error("Num Particles is not consistent with num atoms of provided other_ent!");
+  } 
+
+  //check whether there is a chain from the new entity is already 
+  //present in the actual entity 
+  ost::mol::ChainHandleList other_chains = other_ent.GetChainList();
+  for(ost::mol::ChainHandleList::iterator i = other_chains.begin(), 
+      e = other_chains.end(); i != e; ++i){
+    if(ent.FindChain(i->GetName()).IsValid()){
+      std::stringstream ss;
+      ss << "Chain with name \"" << i->GetName() << "\" from other topology";
+      ss << "is already present in destination topology!";
+      throw ost::Error(ss.str());
+    }
+  }
 }
 
 }}}//ns
diff --git a/modules/mol/mm/src/topology.hh b/modules/mol/mm/src/topology.hh
index 5298a6007..0c9b800d9 100644
--- a/modules/mol/mm/src/topology.hh
+++ b/modules/mol/mm/src/topology.hh
@@ -356,6 +356,8 @@ public:
 
   void Merge(ost::mol::EntityHandle& ent, TopologyPtr other, const ost::mol::EntityHandle& other_ent);
 
+  void Merge(TopologyPtr other);
+
   template <typename DS>
   void Serialize(DS& ds){
 
@@ -665,6 +667,15 @@ public:
 
 private:
 
+  void MergeTop(TopologyPtr other);
+
+  void MergeEnt(ost::mol::EntityHandle& ent, const ost::mol::EntityHandle& other_ent);
+
+  void CheckTopToAdd(TopologyPtr other);
+
+  void CheckEntToAdd(ost::mol::EntityHandle& ent, TopologyPtr other, 
+                     const ost::mol::EntityHandle& other_ent);
+
   uint num_particles_;
 
   //fudge parameters for lj 1,4 pairs
-- 
GitLab