make sure, that disulfid bridge is only built to ONE other sulfur

modules/mol/mm/src/modeller.cc

@@ -19,6 +19,32 @@
 
 #include <ost/mol/mm/modeller.hh>
 
+namespace{
+
+bool IsCysteine(const ost::mol::ResidueHandle& res){
+  String res_name = res.GetName();
+  if(res_name[0] == 'C'){
+    if(res_name == "CYS" || res_name == "CYX" || 
+      res_name == "CYS2" || res_name == "CYM"){
+      return true;
+    }
+  }
+  return false;
+}
+
+bool BoundToSG(const ost::mol::AtomHandle& atom){
+  ost::mol::AtomHandleList bound_atoms = atom.GetBondPartners();
+  for(ost::mol::AtomHandleList::iterator i = bound_atoms.begin();
+      i != bound_atoms.end(); ++i){
+    if(i->GetResidue() != atom.GetResidue()) continue;
+    if(i->GetName() != "SG") continue;
+    if(IsCysteine(i->GetResidue())) return true;
+  }
+  return false;
+}
+
+}
+
 
 namespace ost{ namespace mol{ namespace mm{
 
@@ -29,39 +55,29 @@ void Modeller::GenerateDisulfidBonds(ost::mol::EntityHandle& handle){
 
   for(ost::mol::ResidueHandleList::iterator i = res_list.begin();
       i != res_list.end(); ++i){
-    if(i->GetName() == "CYS" || i->GetName() == "CYX" || i->GetName() == "CYS2" || i->GetName() == "CYM"){
-      ost::mol::AtomHandle s = i->FindAtom("SG");
-      if(s.IsValid()){
-        //quick check, whether max bond partners is already reached
-        if(s.GetBondCount() >= 2) continue;
-        //check, whether SG is already bound to a SG
-        ost::mol::AtomHandleList bonded_atoms = s.GetBondPartners();
-        bool already_bound = false;
-        for(ost::mol::AtomHandleList::iterator bound_it = bonded_atoms.begin();
-            bound_it != bonded_atoms.end(); ++bound_it){
-          if(bound_it->GetName() == "SG"){
-            if(bound_it->GetResidue().GetName() == "CYS" || bound_it->GetResidue().GetName() == "CYX" ||
-               bound_it->GetResidue().GetName() == "CYS2" || bound_it->GetResidue().GetName() == "CYM"){
-              already_bound = true;
-              break;
-            }
-          }
-        }
-        if(already_bound) continue;
-        ost::mol::AtomHandleList in_reach = handle.FindWithin(s.GetPos(),2.5);
-        for(ost::mol::AtomHandleList::iterator j = in_reach.begin();
-            j != in_reach.end(); ++j){
-          if(j->GetName() == "SG" && j->GetResidue() != *i){
-            if(j->GetBondCount() >= 2) continue;
-            if(!ost::mol::BondExists(s,*j)){
-              ed.Connect(s,*j);
-            }
-            ed.RenameResidue(*i,"CYS2");
-            ed.RenameResidue(j->GetResidue(),"CYS2");
-          }
+
+    if(!IsCysteine(*i)) continue;
+    ost::mol::AtomHandle s = i->FindAtom("SG");
+    if(!s.IsValid()) continue;
+    if(BoundToSG(s)) continue; //it's already connected to another sulfur
+
+    //we have now a cysteine gamma-sulfur and can look around for potential
+    //disulfid bridges...
+    ost::mol::AtomHandleList in_reach = handle.FindWithin(s.GetPos(),2.5);
+    for(ost::mol::AtomHandleList::iterator j = in_reach.begin();
+        j != in_reach.end(); ++j){
+      if(j->GetResidue() != *i && j->GetName() == "SG" &&
+         IsCysteine(j->GetResidue())){
+        //the other atom is a gamma-sulfur as well!
+        if(BoundToSG(*j)) continue; //it's already connected to another sulfur
+        if(!ost::mol::BondExists(s,*j)){
+          ed.Connect(s,*j);
         }
+        ed.RenameResidue(*i,"CYS2");
+        ed.RenameResidue(j->GetResidue(),"CYS2");
+        break; //sg can only be bound to one other sg
       }
-    }
+    } 
   }
 }