Subgraph but with clean changes

build/makefile

@@ -106,7 +106,7 @@ test:	./test/RBT_HOME/1koc.as ./test/RBT_HOME/1YET.as ./test/RBT_HOME/1YET_test.
 	@echo ""
 	@echo "Running rDock unit tests..."
 	@echo ""
-	@$(RBT_ROOT)/bin/rbdock -r1YET_test.prm -i ./test/RBT_HOME/1YET_c.sd -p dock.prm -n 1 -s 48151623 -o 1YET_test_out > 1YET_test_out.log
+	@$(RBT_ROOT)/bin/rbdock -r1YETmine.prm -i ./test/RBT_HOME/1YET_c.sd -p dock.prm -n 20 -s 48151623 -o 1YET_test_out > 1YET_test_out.log
 	@mv 1YET_test_out.log 1YET_test_out.sd ./test/RBT_HOME/
 	@python ./test/RBT_HOME/check_test.py ./test/RBT_HOME/1YET_reference_out.sd ./test/RBT_HOME/1YET_test_out.sd
 
@@ -129,8 +129,8 @@ test:	./test/RBT_HOME/1koc.as ./test/RBT_HOME/1YET.as ./test/RBT_HOME/1YET_test.
 	@echo ""
 	@echo "Cavity mapping of 1YET.prm"
 	@echo ""
-	@$(RBT_ROOT)/bin/rbcavity -r1YET.prm -was
-	@mv 1YET.as ./test/RBT_HOME/
+	@$(RBT_ROOT)/bin/rbcavity -r1YETmine.prm -was
+	@mv 1YETmine.as ./test/RBT_HOME/
 
 # Generate the .as file for 1koc
 ./test/RBT_HOME/1YET_test.as:	./test/RBT_HOME/1YET_test.prm

build/test/RBT_HOME/1YET_c.sd

@@ -92,4 +92,8 @@ L_1YET
  11 43  1  0  0  0
  11 44  1  0  0  0
 M  END
+
+>  <SUBGRAPH_CONSTR_0>
+5
+
 $$$$

include/RbtAtom.h

@@ -732,6 +732,17 @@ namespace Rbt
     RbtBool operator()(const RbtAtom *pAtom) const { return (pAtom->GetGroupCharge() > 0.001); }
   };
 
+  // TODO
+  // Is this atom part of the subgraph specified by the user?
+  class isAtomSubgraph : public std::unary_function<RbtAtom *, RbtBool>
+  {
+    RbtInt n;
+
+  public:
+    explicit isAtomSubgraph(RbtInt nn) : n(nn) {}
+    RbtBool operator()(const RbtAtom *pAtom) const { return pAtom->GetAtomId() == n; }
+  };
+
   // Is atom defined as an anionic interaction center ?
   // Checks if group charge is < zero
   class isAtomAnionic : public std::unary_function<RbtAtom *, RbtBool>

include/RbtConstraint.h

@@ -43,7 +43,7 @@ class RbtConstraint
 
   RbtCoord GetCoords() const { return coord; };
   RbtDouble GetTolerance() const { return tolerance; };
-  virtual void AddAtomList(RbtModelPtr, RbtBool bCheck = true) = 0;
+  virtual void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt atomId = -1) = 0;
   RbtDouble Score() const;
 
   ////////////////////
@@ -84,71 +84,79 @@ class RbtHeavyConstraint : public RbtConstraint
 {
 public:
   RbtHeavyConstraint(RbtCoord c, RbtDouble t) : RbtConstraint(c, t){};
-  void AddAtomList(RbtModelPtr, RbtBool bCheck = true) throw(RbtError);
+  void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt = -1) throw(RbtError);
   static RbtInt counter;
 };
 
 class RbtHBAConstraint : public RbtConstraint
 {
 public:
   RbtHBAConstraint(RbtCoord c, RbtDouble t) : RbtConstraint(c, t){};
-  void AddAtomList(RbtModelPtr, RbtBool bCheck = true) throw(RbtError);
+  void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt = -1) throw(RbtError);
   static RbtInt counter;
 };
 
 class RbtHBDConstraint : public RbtConstraint
 {
 public:
   RbtHBDConstraint(RbtCoord c, RbtDouble t) : RbtConstraint(c, t){};
-  void AddAtomList(RbtModelPtr, RbtBool bCheck = true) throw(RbtError);
+  void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt = -1) throw(RbtError);
   static RbtInt counter;
 };
 
 class RbtHydroConstraint : public RbtConstraint
 {
 public:
   RbtHydroConstraint(RbtCoord c, RbtDouble t) : RbtConstraint(c, t){};
-  void AddAtomList(RbtModelPtr, RbtBool bCheck = true) throw(RbtError);
+  void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt = -1) throw(RbtError);
   static RbtInt counter;
 };
 
 class RbtHydroAliphaticConstraint : public RbtConstraint
 {
 public:
   RbtHydroAliphaticConstraint(RbtCoord c, RbtDouble t) : RbtConstraint(c, t){};
-  void AddAtomList(RbtModelPtr, RbtBool bCheck = true) throw(RbtError);
+  void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt = -1) throw(RbtError);
   static RbtInt counter;
 };
 
 class RbtHydroAromaticConstraint : public RbtConstraint
 {
 public:
   RbtHydroAromaticConstraint(RbtCoord c, RbtDouble t) : RbtConstraint(c, t){};
-  void AddAtomList(RbtModelPtr, RbtBool bCheck = true) throw(RbtError);
+  void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt = -1) throw(RbtError);
   static RbtInt counter;
 };
 
 class RbtNegChargeConstraint : public RbtConstraint
 {
 public:
   RbtNegChargeConstraint(RbtCoord c, RbtDouble t) : RbtConstraint(c, t){};
-  void AddAtomList(RbtModelPtr, RbtBool bCheck = true) throw(RbtError);
+  void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt = -1) throw(RbtError);
   static RbtInt counter;
 };
 
 class RbtPosChargeConstraint : public RbtConstraint
 {
 public:
   RbtPosChargeConstraint(RbtCoord c, RbtDouble t) : RbtConstraint(c, t){};
-  void AddAtomList(RbtModelPtr, RbtBool bCheck = true) throw(RbtError);
+  void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt = -1) throw(RbtError);
+  static RbtInt counter;
+};
+
+class RbtSubgraphConstraint : public RbtConstraint
+{
+public:
+  RbtSubgraphConstraint(RbtCoord c, RbtDouble t) : RbtConstraint(c, t){};
+  void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt = -1) throw(RbtError);
   static RbtInt counter;
 };
 
 class RbtRingAromaticConstraint : public RbtConstraint
 {
 public:
   RbtRingAromaticConstraint(RbtCoord c, RbtDouble t) : RbtConstraint(c, t){};
-  void AddAtomList(RbtModelPtr, RbtBool bCheck = true) throw(RbtError);
+  void AddAtomList(RbtModelPtr, RbtBool bCheck = true, RbtInt = -1) throw(RbtError);
   static RbtInt counter;
 };
 

src/lib/RbtCavityGridSF.cxx

@@ -54,10 +54,15 @@ void RbtCavityGridSF::SetupReceptor()
   m_spGrid = spDS->GetGrid();
   if (m_spGrid.Null())
     return;
+
+  // RbtStringList strTetherAtomsL = GetLigand()->GetDataValue("SYMMETRIC_BONDS");
   m_maxDist = m_spGrid->MaxValue();
 }
 
-void RbtCavityGridSF::SetupLigand() {}
+void RbtCavityGridSF::SetupLigand()
+{
+  RbtStringList strTetherAtomsL = GetLigand()->GetDataValue("SYMMETRIC_BONDS");
+}
 
 void RbtCavityGridSF::SetupSolvent() {}
 

src/lib/RbtConstraint.cxx

@@ -27,6 +27,7 @@ RbtInt RbtHydroAliphaticConstraint::counter = 0;
 RbtInt RbtHydroAromaticConstraint::counter = 0;
 RbtInt RbtNegChargeConstraint::counter = 0;
 RbtInt RbtPosChargeConstraint::counter = 0;
+RbtInt RbtSubgraphConstraint::counter = 0;
 RbtInt RbtRingAromaticConstraint::counter = 0;
 
 RbtConstraint::RbtConstraint(RbtCoord c, RbtDouble t)
@@ -99,6 +100,7 @@ RbtDouble RbtConstraint::Score() const
 
 RbtConstraintPtr Rbt::CreateConstraint(RbtCoord &c, RbtDouble &t, RbtString &n, RbtBool bCount)
 {
+  // throw RbtError(_WHERE_, "Problems reading constraint ABC " + n);
   if (n == "Any")
   {
     if (bCount)
@@ -153,6 +155,12 @@ RbtConstraintPtr Rbt::CreateConstraint(RbtCoord &c, RbtDouble &t, RbtString &n,
       RbtHydroAromaticConstraint::counter++;
     return new RbtHydroAromaticConstraint(c, t);
   }
+  else if (n == "Sub")
+  {
+    if (bCount)
+      RbtSubgraphConstraint::counter++;
+    return new RbtSubgraphConstraint(c, t);
+  }
   else
     throw RbtError(_WHERE_, "Constraint " + n + " not recognized");
 }
@@ -168,6 +176,7 @@ void Rbt::ZeroCounters()
   RbtNegChargeConstraint::counter = 0;
   RbtPosChargeConstraint::counter = 0;
   RbtRingAromaticConstraint::counter = 0;
+  RbtSubgraphConstraint::counter = 0;
 }
 
 void Rbt::ReadConstraint(istream &ifile, RbtConstraintPtr &cnt, RbtBool bCount)
@@ -256,7 +265,7 @@ void Rbt::ReadConstraints(istream &ifile, RbtConstraintList &cl, RbtBool bCount)
 }
 
 // 07 Feb 2005 (DM) - new constraint type, any heavy atom
-void RbtHeavyConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtError)
+void RbtHeavyConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck, RbtInt atomID) throw(RbtError)
 {
   m_atomList = Rbt::GetAtomList(lig->GetAtomList(), std::not1(Rbt::isAtomicNo_eq(1)));
   if (bCheck && (m_atomList.size() < counter))
@@ -269,7 +278,7 @@ void RbtHeavyConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtE
 }
 
 // 16 May 2003 (DM) - limit to neutral H-bond donor hydrogens
-void RbtHBDConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtError)
+void RbtHBDConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck, RbtInt atomID) throw(RbtError)
 {
   m_atomList = Rbt::GetAtomList(lig->GetAtomList(), Rbt::isAtomHBondDonor());
   m_atomList = Rbt::GetAtomList(m_atomList, std::not1(Rbt::isAtomCationic()));
@@ -283,7 +292,7 @@ void RbtHBDConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtErr
 }
 
 // 16 May 2003 (DM) - limit to neutral H-bond acceptors
-void RbtHBAConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtError)
+void RbtHBAConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck, RbtInt atomID) throw(RbtError)
 {
   m_atomList = Rbt::GetAtomList(lig->GetAtomList(), Rbt::isAtomHBondAcceptor());
   m_atomList = Rbt::GetAtomList(m_atomList, std::not1(Rbt::isAtomAnionic()));
@@ -296,7 +305,7 @@ void RbtHBAConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtErr
   }
 }
 
-void RbtHydroConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtError)
+void RbtHydroConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck, RbtInt atomID) throw(RbtError)
 {
   m_atomList = Rbt::GetAtomList(lig->GetAtomList(), Rbt::isAtomLipophilic());
   if (bCheck && (m_atomList.size() < counter))
@@ -308,7 +317,7 @@ void RbtHydroConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtE
   }
 }
 
-void RbtHydroAliphaticConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtError)
+void RbtHydroAliphaticConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck, RbtInt atomID) throw(RbtError)
 {
   m_atomList = Rbt::GetAtomList(lig->GetAtomList(), Rbt::isAtomLipophilic());
   m_atomList = Rbt::GetAtomList(m_atomList,
@@ -322,7 +331,7 @@ void RbtHydroAliphaticConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) t
   }
 }
 
-void RbtHydroAromaticConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtError)
+void RbtHydroAromaticConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck, RbtInt atomID) throw(RbtError)
 {
   m_atomList = Rbt::GetAtomList(lig->GetAtomList(), Rbt::isAtomLipophilic());
   m_atomList = Rbt::GetAtomList(m_atomList,
@@ -336,7 +345,7 @@ void RbtHydroAromaticConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) th
   }
 }
 
-void RbtNegChargeConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtError)
+void RbtNegChargeConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck, RbtInt atomID) throw(RbtError)
 {
   m_atomList = Rbt::GetAtomList(lig->GetAtomList(), Rbt::isAtomAnionic());
   if (bCheck && (m_atomList.size() < counter))
@@ -348,7 +357,7 @@ void RbtNegChargeConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(
   }
 }
 
-void RbtPosChargeConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtError)
+void RbtPosChargeConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck, RbtInt atomID) throw(RbtError)
 {
   m_atomList = Rbt::GetAtomList(lig->GetAtomList(), Rbt::isAtomCationic());
   if (bCheck && (m_atomList.size() < counter))
@@ -360,7 +369,20 @@ void RbtPosChargeConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(
   }
 }
 
-void RbtRingAromaticConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck) throw(RbtError)
+void RbtSubgraphConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck, RbtInt atomId) throw(RbtError)
+{
+  m_atomList = Rbt::GetAtomList(lig->GetAtomList(), Rbt::isAtomSubgraph(atomId));
+
+  if (bCheck && (m_atomList.size() < counter))
+  {
+    ostringstream ostr;
+    ostr << "The ligand has only " << m_atomList.size()
+         << " atom(s) belonging to that subgraph type (" << counter << " required)" << ends;
+    throw RbtLigandError(_WHERE_, ostr.str());
+  }
+}
+
+void RbtRingAromaticConstraint::AddAtomList(RbtModelPtr lig, RbtBool bCheck, RbtInt atomID) throw(RbtError)
 {
   m_atomList.clear();
   RbtAtomList at = lig->GetAtomList();

src/lib/RbtMdlFileSource.cxx

@@ -189,8 +189,10 @@ void RbtMdlFileSource::Parse() throw(RbtError)
       // DM 12 May 1999 - read data records (if any)
       for (; fileIter != fileEnd; fileIter++)
       {
+        // cout << "here"; // Found a data record
         if ((*fileIter).find(">") == 0)
-        {                                                   // Found a data record
+        {
+          // cout << "here record";                            // Found a data record
           RbtString::size_type ob = (*fileIter).find("<");  // First open bracket
           RbtString::size_type cb = (*fileIter).rfind(">"); // Last closed bracket
           if ((ob != RbtString::npos) && (cb != RbtString::npos))
@@ -201,6 +203,7 @@ void RbtMdlFileSource::Parse() throw(RbtError)
             {
               sl.push_back(*fileIter);
             }
+            // cout << "Added record";
             m_dataMap[fieldName] = RbtVariant(sl);
           }
         }
@@ -214,6 +217,16 @@ void RbtMdlFileSource::Parse() throw(RbtError)
         m_dataMap["Name"] = RbtVariant(m_titleList.front());
       }
 
+      // cout << "ALEJANDRO says HELLO WORLD " << (*iter)->GetAtomName() << " to plane of substituents=" << dist << endl;
+      // cout << "ALEJANDRO says \n\n\n HELLO WORLD  \n\n\n ";
+
+      // for (RbtStringVariantMapConstIter iter = m_dataMap.begin(); iter != m_dataMap.end(); iter++)
+      //{
+      //   cout << (*iter).first;
+      //   // cout << (*iter).second;
+      //   cout << '\n\n' << endl;
+      // }
+
       // Setup the atomic params not stored in the file (e.g. hybridisation state etc)
       SetupAtomParams();
 
@@ -553,20 +566,20 @@ void RbtMdlFileSource::SetupPosIonisableGroups()
 #ifdef _DEBUG
         cout << "Possible imidazole/amidine, found 1xN_SP2, 1xN_TRI bonded to " << (*iter)->GetAtomName() << endl;
 #endif //_DEBUG
-        // IMIDAZOLE - Check that the N_SP2 is bonded to 2 x C_SP2 and that
-        // the N_TRI is bonded to 2 x C_SP2 and 1 hydrogen
+       // IMIDAZOLE - Check that the N_SP2 is bonded to 2 x C_SP2 and that
+       // the N_TRI is bonded to 2 x C_SP2 and 1 hydrogen
         if ((Rbt::GetNumAtoms(nsp2BondedAtomList, bIsC_SP2) == 2) &&
             (Rbt::GetNumAtoms(ntriBondedAtomList, bIsC_SP2) == 2) &&
             (Rbt::GetNumAtoms(ntriBondedAtomList, Rbt::isAtomicNo_eq(1)) == 1))
         {
 #ifdef _DEBUG
           cout << "Possible imidazole, bonding requirements for N_SP2 and N_TRI met" << endl;
 #endif //_DEBUG
-          // Now check if any of the atoms bonded to the N_TRI are 1-2 connected to any
-          // of the atoms bonded to the N_SP2. If so, it is a 5-membered ring
-          // DM 25 Jul 2002 - also check whether the atoms are bridgeheads or not
-          // We don't want to protonate imidazoles which are part of larger fused ring systems
-          //(can have very different pKa's)
+       // Now check if any of the atoms bonded to the N_TRI are 1-2 connected to any
+       // of the atoms bonded to the N_SP2. If so, it is a 5-membered ring
+       // DM 25 Jul 2002 - also check whether the atoms are bridgeheads or not
+       // We don't want to protonate imidazoles which are part of larger fused ring systems
+       //(can have very different pKa's)
           RbtAtomList atoms12Conn;
           for (RbtAtomListConstIter iter2 = ntriBondedAtomList.begin(); iter2 != ntriBondedAtomList.end(); iter2++)
           {
@@ -578,7 +591,7 @@ void RbtMdlFileSource::SetupPosIonisableGroups()
 #ifdef _DEBUG
             cout << "5-membered imidazole ring found" << endl;
 #endif //_DEBUG
-            // Remove bridgehead atoms (fused rings)
+       // Remove bridgehead atoms (fused rings)
             atoms12Conn = Rbt::GetAtomList(atoms12Conn, std::not1(Rbt::isAtomBridgehead()));
             if (atoms12Conn.size() == 1)
             { // Imidazole!!

src/lib/RbtPharmaSF.cxx

@@ -129,6 +129,17 @@ void RbtPharmaSF::SetupLigand()
 {
   if (GetLigand().Null())
     return;
+
+  RbtStringList strTetherAtomsL = GetLigand()->GetDataValue("SUBGRAPH_CONSTR_0");
+  // print each item on the stringlist
+  for (RbtStringListIter iter = strTetherAtomsL.begin(); iter != strTetherAtomsL.end(); iter++)
+  {
+    cout << "\n\n\natom to ANCHOR" << *iter << endl;
+  }
+
+  cout << "\n\n\nfirst and only:::" << strTetherAtomsL[0] << endl;
+
+  RbtInt myAtomId = stoi(strTetherAtomsL[0]);
   try
   {
     if (GetTrace() > 0)
@@ -137,7 +148,7 @@ void RbtPharmaSF::SetupLigand()
     }
     for (RbtConstraintListIter iter = m_constrList.begin(); iter != m_constrList.end(); iter++)
     {
-      (*iter)->AddAtomList(GetLigand(), true);
+      (*iter)->AddAtomList(GetLigand(), true, myAtomId); // I would like to pass a stringlist here for one type of constraint (sub)
     }
     if (GetTrace() > 0)
     {