C++ OBMol::NextBond方法代码示例

  static int DetermineFRJ(OBMol &mol)
  {
    vector<vector<int> >::iterator i;
    vector<vector<int> > cfl;
    //find all continuous graphs in the mol area
    mol.ContigFragList(cfl);

    if (cfl.empty())
      return(0);
    if (cfl.size() == 1)
      return(mol.NumBonds() - mol.NumAtoms() + 1);

    //count up the atoms and bonds belonging to each graph
    OBBond *bond;
    vector<OBBond*>::iterator j;
    int numatoms,numbonds,frj=0;
    OBBitVec frag;
    for (i = cfl.begin();i != cfl.end();++i)
      {
        frag.Clear();
        frag.FromVecInt(*i);
        numatoms = (*i).size();
        numbonds=0;
        for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
          if (frag.BitIsOn(bond->GetBeginAtomIdx()) &&
              frag.BitIsOn(bond->GetEndAtomIdx()))
            numbonds++;
        frj += numbonds - numatoms + 1;
      }

    return(frj);
  }

  void ReportFormat::WriteTorsions(ostream &ofs,OBMol &mol)
  {
    vector<OBBond*>::iterator bi1,bi2,bi3;
    OBBond* bond;
    OBAtom *a,*b,*c,*d;
    char buffer[BUFF_SIZE];

    //loop through all bonds generating torsions
    for(bond = mol.BeginBond(bi1); bond; bond = mol.NextBond(bi1))
      {
        b = bond->GetBeginAtom();
        c = bond->GetEndAtom();

        for(a = b->BeginNbrAtom(bi2);a;a = b->NextNbrAtom(bi2))
          {
            if(a == c)
              continue;

            for(d = c->BeginNbrAtom(bi3);d;d = c->NextNbrAtom(bi3))
              {
                if(d == b)
                  continue;

                snprintf(buffer, BUFF_SIZE, "%4d %4d %4d %4d %10.3f",
                        a->GetIdx(), b->GetIdx(),c->GetIdx(),d->GetIdx(),
                        CalcTorsionAngle(a->GetVector(), b->GetVector(),
                                         c->GetVector(), d->GetVector()));
                ofs << buffer << "\n";
              }
          }
      }
  }

  bool OBBond::IsClosure()
  {
    OBMol *mol = (OBMol*)GetParent();
    if (!mol)
      return(false);
    if (mol->HasClosureBondsPerceived())
      return(HasFlag(OB_CLOSURE_BOND));

    mol->SetClosureBondsPerceived();

    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::PerceiveClosureBonds", obAuditMsg);

    OBBond *bond;
    OBAtom *atom,*nbr;
    OBBitVec uatoms,ubonds;
    vector<OBAtom*> curr,next;
    vector<OBAtom*>::iterator i;
    vector<OBBond*>::iterator j;

    uatoms.Resize(mol->NumAtoms()+1);
    ubonds.Resize(mol->NumAtoms()+1);

    for (;static_cast<unsigned int>(uatoms.CountBits()) < mol->NumAtoms();)
      {
        if (curr.empty())
          for (atom = mol->BeginAtom(i);atom;atom = mol->NextAtom(i))
            if (!uatoms[atom->GetIdx()])
              {
                uatoms |= atom->GetIdx();
                curr.push_back(atom);
                break;
              }

        for (;!curr.empty();)
          {
            for (i = curr.begin();i != curr.end();++i)
              for (nbr = ((OBAtom*)*i)->BeginNbrAtom(j);nbr;nbr = ((OBAtom*)*i)->NextNbrAtom(j))
                if (!uatoms[nbr->GetIdx()])
                  {
                    uatoms |= nbr->GetIdx();
                    ubonds |= (*j)->GetIdx();
                    next.push_back(nbr);
                  }

            curr = next;
            next.clear();
          }
      }

    for (bond = mol->BeginBond(j);bond;bond = mol->NextBond(j))
      if (!ubonds[bond->GetIdx()])
        bond->SetClosure();

    return(HasFlag(OB_CLOSURE_BOND));
  }

  bool OBRotorList::FindRotors(OBMol &mol, bool sampleRingBonds)
  {
    // Find ring atoms & bonds
    // This function will set OBBond::IsRotor().
    mol.FindRingAtomsAndBonds();

    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::FindRotors", obAuditMsg);

    //
    // Score the bonds using the graph theoretical distance (GTD).
    // The GTD is the distance from atom i to every other atom j.
    // Atoms on the "inside" of the molecule will have a lower GTD
    // value than atoms on the "outside"
    //
    // The scoring will rank "inside" bonds first.
    //
    vector<int> gtd;
    mol.GetGTDVector(gtd);
    // compute the scores
    vector<OBBond*>::iterator i;
    vector<pair<OBBond*,int> > vtmp;
    for (OBBond *bond = mol.BeginBond(i);bond;bond = mol.NextBond(i)) {
      // check if the bond is "rotatable"
      if (bond->IsRotor(sampleRingBonds)) {
        // check if the bond is fixed (using deprecated fixed atoms or new fixed bonds)
        if ((HasFixedAtoms() || HasFixedBonds()) && IsFixedBond(bond))
          continue;

        if (bond->IsInRing()) {
          //otherwise mark that we have them and add it to the pile
          _ringRotors = true;
        }

        int score = gtd[bond->GetBeginAtomIdx()-1] + gtd[bond->GetEndAtomIdx()-1];
        // compute the GTD bond score as sum of atom GTD scores
        vtmp.push_back(pair<OBBond*,int> (bond,score));
      }
    }

    // sort the rotatable bonds by GTD score
    sort(vtmp.begin(),vtmp.end(),CompareRotor);

    // create rotors for the bonds
    int count = 0;
    vector<pair<OBBond*,int> >::iterator j;
    for (j = vtmp.begin(); j != vtmp.end(); ++j, ++count) {
      OBRotor *rotor = new OBRotor;
      rotor->SetBond((*j).first);
      rotor->SetIdx(count);
      rotor->SetNumCoords(mol.NumAtoms()*3);
      _rotor.push_back(rotor);
    }

    return true;
  }

  static int DetermineFRJ(OBMol &mol)
  {
    if (!mol.HasClosureBondsPerceived())
      return (int)FindRingAtomsAndBonds2(mol);

    int frj = 0;
    OBBond *bond;
    vector<OBBond*>::iterator j;
    for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
      if (bond->IsClosure()) // bond->HasFlag(OB_CLOSURE_BOND)?
        frj++;
    return frj;
  }

bool WriteAlchemy(ostream &ofs,OBMol &mol)
{
  unsigned int i;
  char buffer[BUFF_SIZE];
  char bond_string[10];
  
  snprintf(buffer, BUFF_SIZE, "%5d ATOMS, %5d BONDS,     0 CHARGES",
	  mol.NumAtoms(),
	  mol.NumBonds());
  ofs << buffer << endl;
  ttab.SetFromType("INT"); ttab.SetToType("ALC");

  OBAtom *atom;
  string str,str1;
  for(i = 1;i <= mol.NumAtoms(); i++)
  {
    atom = mol.GetAtom(i);
    str = atom->GetType();
    ttab.Translate(str1,str);
    snprintf(buffer, BUFF_SIZE, "%5d %-6s%8.4f %8.4f %8.4f     0.0000",
	    i,
	    (char*)str1.c_str(),
	    atom->GetX(),
	    atom->GetY(),
	    atom->GetZ());
    ofs << buffer << endl;
  }

  OBBond *bond;
  vector<OBEdgeBase*>::iterator j;

  for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
  {
    switch(bond->GetBO())
    {
    case 1 :  strcpy(bond_string,"SINGLE");  break;
    case 2 :  strcpy(bond_string,"DOUBLE");  break;
    case 3 :  strcpy(bond_string,"TRIPLE");  break;
    case 5 :  strcpy(bond_string,"AROMATIC"); break;
    default : strcpy(bond_string,"SINGLE");
    }
    snprintf(buffer, BUFF_SIZE, "%5d  %4d  %4d  %s",
	    bond->GetIdx()+1,
	    bond->GetBeginAtomIdx(),
	    bond->GetEndAtomIdx(),
	    bond_string);
    ofs << buffer << endl;
  }
  return(true);
}

  void WriteAngles(ostream &ofs,OBMol &mol)
  {
    // Alas, we still need to sort these to only list unique entries...
    vector3 v1, v2;
    OBAtom *a, *b, *c, *d;
    OBBond *bond1, *bond2, *bond3;
    vector<OBEdgeBase*>::iterator i, j, k;
    char buffer[BUFF_SIZE];

    for (bond1 = mol.BeginBond(i); bond1; bond1 = mol.NextBond(i))
      {
	b = bond1->GetBeginAtom();
	c = bond1->GetEndAtom();

	for (bond2 = b->BeginBond(j); bond2; bond2 = b->NextBond(j))
	  {
	    if (bond2->GetEndAtomIdx() != c->GetIdx() 
		&& bond2->GetEndAtomIdx() != b->GetIdx())
	      {
		a = bond2->GetEndAtom();

		v1 = a->GetVector() - b->GetVector();
		v2 = c->GetVector() - b->GetVector();

		sprintf(buffer,"%4d %4d %4d %4s %4s %4s %10.3f",
			a->GetIdx(),b->GetIdx(),c->GetIdx(),
			a->GetType(),b->GetType(),c->GetType(),
			vectorAngle(v1, v2));
		ofs << buffer << endl;

		for (bond3 = c->BeginBond(k); bond3; bond3 = c->NextBond(k))
		  if (bond3->GetEndAtomIdx() != b->GetIdx()
		      && bond3->GetEndAtomIdx() != c->GetIdx())
		    {
		      d = bond3->GetEndAtom();

		      v1 = b->GetVector() - c->GetVector();
		      v2 = d->GetVector() - c->GetVector();

		      sprintf(buffer,"%4d %4d %4d %4s %4s %4s %10.3f",
			      b->GetIdx(),c->GetIdx(),d->GetIdx(),
			      b->GetType(),c->GetType(),d->GetType(),
			      vectorAngle(v1, v2));
		      ofs << buffer << endl;
		    }
	      }
	  }
      }
  }

void ring_test()
{
  ostringstream os;

#ifdef TESTDATADIR
  string testdatadir = TESTDATADIR;
  string results_file = testdatadir + "ringresults.txt";
  string smilestypes_file = testdatadir + "attype.00.smi";
#else
  string results_file = "files/ringresults.txt";
  string smilestypes_file = "files/attype.00.smi";
#endif

  cout << "# Testing ring perception..." << endl;

  std::ifstream mifs;
  os << "Bail out! Cannot read test data " << smilestypes_file.c_str();
  BOOST_REQUIRE_MESSAGE( SafeOpen(mifs, smilestypes_file.c_str()), os.str().c_str() );

  std::ifstream rifs;
  os.str("");
  os << "Bail out! Cannot read test data " << results_file.c_str();
  BOOST_REQUIRE_MESSAGE( SafeOpen(rifs, results_file.c_str()), os.str().c_str() );

  unsigned int size;
  OBBond *bond;
  OBAtom *atom;
  int count;
  char buffer[BUFF_SIZE];
  vector<string> vs;
  vector<OBRing*> vr;
  vector<bool> vb;
  vector<int> vi;
  OBMol mol;
  vector<string>::iterator i;
  vector<OBBond*>::iterator j;
  vector<OBAtom*>::iterator k;
  vector<OBRing*>::iterator m;
  OBConversion conv(&mifs, &cout);
  unsigned int currentTest = 0;

  BOOST_REQUIRE_MESSAGE( conv.SetInAndOutFormats("SMI","SMI"), "Bail out! SMILES format is not loaded" );

  for (;mifs;)
    {
      mol.Clear();
      conv.Read(&mol);
      if (mol.Empty())
        continue;
      
      BOOST_REQUIRE_MESSAGE( rifs.getline(buffer,BUFF_SIZE), "Bail out! error reading reference data" );

      vb.clear();
      vb.resize(mol.NumBonds(),false);
      //check ring bonds
      tokenize(vs,buffer);
      for (i = vs.begin();i != vs.end();i++)
        vb[atoi(i->c_str())] = true;

      for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
        {
          os.str("");
	  os << "ring bond data different than reference # Molecule: " << mol.GetTitle();
          BOOST_CHECK_MESSAGE( vb[bond->GetIdx()] == bond->IsInRing(), os.str().c_str() );
        }

      vr = mol.GetSSSR();
      BOOST_REQUIRE_MESSAGE( rifs.getline(buffer,BUFF_SIZE), "Bail out! error reading reference data" );
      sscanf(buffer,"%d",&size);

      os.str("");
      os << "SSSR size different than reference # Molecule: " << mol.GetTitle();
      BOOST_CHECK_MESSAGE( vr.size() == size, os.str().c_str() ); //check SSSR size

      BOOST_REQUIRE_MESSAGE( rifs.getline(buffer,BUFF_SIZE), "Bail out! error reading reference data" );

      tokenize(vs,buffer);
      i = vs.begin();
      for (atom = mol.BeginAtom(k);atom;atom = mol.NextAtom(k))
        {
          os.str("");
          os << "error in SSSR count # Molecule: " << mol.GetTitle();
          BOOST_CHECK_MESSAGE( i != vs.end(), os.str().c_str() ); //check SSSR size

          count = 0;
          for (m = vr.begin();m != vr.end();m++)
            if ((*m)->_pathset[atom->GetIdx()])
              count++;

	  os.str("");
          os << "ring membership test failed # Molecule: " << mol.GetTitle();
          BOOST_CHECK_MESSAGE( atoi(i->c_str()) == count, os.str().c_str() ); 

          i++;
        }
    }

}

  /**
   * \brief Select the root atoms for traversing atoms in rings.
   *
   * Picking only the begin atom of a closure bond can cause
   * difficulties when the selected atom is an inner atom
   * with three neighbour ring atoms. Why ? Because this atom
   * can get trapped by the other atoms when determining aromaticity,
   * because a simple visited flag is used in the
   * OBAromaticTyper::TraverseCycle() method.
   *
   * Ported from JOELib, copyright Joerg Wegner, 2003 under the GPL version 2
   * Improved by Fabian (fab5) in 2009 -- PR#2889708
   *
   * @param mol the molecule
   * @param avoidInnerRingAtoms inner closure ring atoms with more than 2 neighbours will be avoided
   *
   */
  void OBAromaticTyper::SelectRootAtoms(OBMol &mol, bool avoidInnerRingAtoms)
  {
    OBBond *bond;
    OBAtom *atom, *nbr, *nbr2;
    OBRing *ring;
    //    vector<OBAtom*>::iterator i;
    vector<OBBond*>::iterator j, l, nbr2Iter;
    vector<OBRing*> sssRings = mol.GetSSSR();
    vector<OBRing*>::iterator k;
    
    int rootAtom;
    int ringNbrs;
    int heavyNbrs;
    int newRoot = -1;
    vector<int> tmpRootAtoms;
    vector<int> tmp;

    vector<OBBond*> cbonds;
    vector< vector<OBRing*> > ringAtoms; // store ring pointers on an atom basis
    
    //generate list of closure bonds 
    for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
      {
        if( bond->IsClosure() )
          {
            cbonds.push_back(bond);
            if(avoidInnerRingAtoms)
              {
                tmpRootAtoms.push_back(bond->GetBeginAtomIdx());
              }
          }
      }
    
    if(avoidInnerRingAtoms)
      {
        //for every atom fill vector with ring pointer it's associated with
        ringAtoms.resize(mol.NumAtoms()+1);
        for (k = sssRings.begin();k != sssRings.end();++k)
          {
            tmp = (*k)->_path;
            for (unsigned int j (0),j_end(tmp.size()); j < j_end; ++j)
              {
                ringAtoms[tmp[j]].push_back(*k);
              }
          }
      }
    
    
    //loop over closure bonds
    for(OBBondIterator bd(cbonds.begin()),bd_end(cbonds.end());bd!=bd_end;++bd)
      {
        bond = *bd;

        // BASIC APPROACH
        // pick beginning atom at closure bond
        // this is really ready, isn't it ! ;-)
        rootAtom = bond->GetBeginAtomIdx();
        _root[rootAtom] = true;
	
        // EXTENDED APPROACH
        if (avoidInnerRingAtoms)
          {
            // count the number of neighbor ring atoms
            atom = mol.GetAtom(rootAtom);
            ringNbrs = heavyNbrs = 0;
	    
            for (nbr = atom->BeginNbrAtom(l);nbr;nbr = atom->NextNbrAtom(l))
              {
                // we can get this from atom->GetHvyValence()
                // but we need to find neighbors in rings too
                // so let's save some time
                if (!nbr->IsHydrogen())
                  {
                    heavyNbrs++;
                    if (nbr->IsInRing())
                      ringNbrs++;
                  }
		
                // if this atom has more than 2 neighbor ring atoms
                // we could get trapped later when traversing cycles
                // which can cause aromaticity false detection
                newRoot = -1;
		
//.........这里部分代码省略.........

  void OBAromaticTyper::AssignAromaticFlags(OBMol &mol)
  {
    if (!_init)
      Init();

    if (mol.HasAromaticPerceived())
      return;
    mol.SetAromaticPerceived();
    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::AssignAromaticFlags", obAuditMsg);

    _vpa.clear();
    _vpa.resize(mol.NumAtoms()+1);
    _velec.clear();
    _velec.resize(mol.NumAtoms()+1);
    _root.clear();
    _root.resize(mol.NumAtoms()+1);

    OBBond *bond;
    OBAtom *atom;
    vector<OBAtom*>::iterator i;
    vector<OBBond*>::iterator j;

    //unset all aromatic flags
    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      atom->UnsetAromatic();
    for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
      bond->UnsetAromatic();

    int idx;
    vector<vector<int> >::iterator m;
    vector<OBSmartsPattern*>::iterator k;

    //mark atoms as potentially aromatic
    for (idx=0,k = _vsp.begin();k != _vsp.end();++k,++idx)
      if ((*k)->Match(mol))
        {
          _mlist = (*k)->GetMapList();
          for (m = _mlist.begin();m != _mlist.end();++m)
            {
              _vpa[(*m)[0]] = true;
              _velec[(*m)[0]] = _verange[idx];
            }
        }

    //sanity check - exclude all 4 substituted atoms and sp centers
    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      {
        if (atom->GetImplicitValence() > 3)
          {
            _vpa[atom->GetIdx()] = false;
            continue;
          }

        switch(atom->GetAtomicNum())
          {
            //phosphorus and sulfur may be initially typed as sp3
          case 6:
          case 7:
          case 8:
            if (atom->GetHyb() != 2)
              _vpa[atom->GetIdx()] = false;
            break;
          }
      }

    //propagate potentially aromatic atoms
    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      if (_vpa[atom->GetIdx()])
        PropagatePotentialAromatic(atom);

    //select root atoms
    SelectRootAtoms(mol);

    ExcludeSmallRing(mol); //remove 3 membered rings from consideration

    //loop over root atoms and look for aromatic rings
    _visit.clear();
    _visit.resize(mol.NumAtoms()+1);
    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      if (_root[atom->GetIdx()])
        CheckAromaticity(atom,14);

    //for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
    //	  if (atom->IsAromatic())
    //		  cerr << "aro = " <<atom->GetIdx()  << endl;

    //for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
    //if (bond->IsAromatic())
    //cerr << bond->GetIdx() << ' ' << bond->IsAromatic() << endl;
  }

//.........这里部分代码省略.........
//            cerr << nodeid1 << inputedgelabel << "(" << (*bond)->IsAromatic() << ")" << nodeid2 << " ";
            NodeLabel node1label = tree->nodes[nodeid1].nodelabel;
            NodeLabel node2label = tree->nodes[nodeid2].nodelabel;
            
            //cerr << "(" << (*bond)->GetBeginAtomIdx() << "[" << (int) node1label << "] " << (*bond)->GetBondOrder() << " " << (*bond)->GetEndAtomIdx() << "[" << (int) node2label << "]"<< ")" << endl;

            
            // Direction of edge always from low to high
            if ( node1label > node2label ) {
                NodeLabel temp = node1label;
                node1label = node2label;
                node2label = temp;
            }

            // Create combined input label for edge
            CombinedInputLabel combinedinputlabel;
            combinedinputlabel = combineInputLabels ( inputedgelabel, node1label, node2label );

            // Insert into map, analoguous to nodes, using subsequent numbering for internal labels:
            // edge nr. 1, edge nr. 2, ...
            
            // Direction inputlabel -> internal label
            //cerr << "EdgeLabelMap: insert " << combinedinputlabel << "-->" << edgelabels.size() << endl;
            map_insert_pair ( edgelabelmap ) p = edgelabelmap.insert ( make_pair ( combinedinputlabel, edgelabels.size () ) );
     
            // Direction internal label -> edge
            if ( p.second ) {
              vector_push_back ( DatabaseEdgeLabel, edgelabels, edgelabel );
              edgelabel.fromnodelabel = node1label;	// directed edges
              edgelabel.tonodelabel = node2label;
              edgelabel.inputedgelabel = inputedgelabel;
              edgelabel.lasttid = tid;
              //cerr << "Created edge " << edgelabel.inputedgelabel << " (" << (int) edgelabel.fromnodelabel << "-->" << (int) edgelabel.tonodelabel << ")" << ":" << edgelabel.frequency << endl;
            }
            else {
              DatabaseEdgeLabel &edgelabel = edgelabels[p.first->second];
              if ( edgelabel.lasttid != tid )
                edgelabel.frequency++;
                edgelabel.lasttid = tid;
                //cerr << "Updated edge " << edgelabel.inputedgelabel << " (" << (int) edgelabel.fromnodelabel << "-->" << (int) edgelabel.tonodelabel << ")" << ":" << edgelabel.frequency << endl;
            }
        
            // Tree edge (2 versions)
            vector_push_back ( DatabaseTreeEdge, edges[nodeid1], edge );
            edge.edgelabel = p.first->second;
            edge.tonode = nodeid2;
    //        edge.bond = (*bond);
                
            vector_push_back ( DatabaseTreeEdge, edges[nodeid2], edge2 );
            edge2.edgelabel = p.first->second;
            edge2.tonode = nodeid1;
    //        edge2.bond = (*bond);
           
            edgessize++;

        } while (mol.NextBond(bond));
    }

    // copy edges to tree
    tree->edges = new DatabaseTreeEdge[edgessize * 2];	// allocate space
                                                        // two instances per edge (both directions)
    int pos = 0;
    for ( int i = 0; i < nodessize; i++ ) {				// for all nodes...
        int s = edges[i].size ();
        tree->nodes[i].edges._size = s;						// edges stored in s steps from ...
        tree->nodes[i].edges.array = tree->edges + pos;		// ... offset for that node's edges
        for ( int j = 0; j < s; j++, pos++ ) {
            //cerr << i << " " << j << " L: " << edges[i][j].tonode << endl;
            tree->edges[pos] = edges[i][j];					// flattened storage
        }
    }
  
    ////////////
    // CYCLES //
    ////////////

    static vector<int> nodestack;
    static vector<bool> visited1, visited2;
    nodestack.resize ( 0 );
    visited1.resize ( 0 );
    visited1.resize ( nodessize, false );
    visited2.resize ( 0 );
    visited2.resize ( nodessize, false );
  


    // for every node...
    for ( int i = 0; i < nodessize; i++ ) {
        if ( !visited1[i] ) {
            nodestack.push_back ( i );
            visited1[i] = visited2[i] = true;		// visited1: has been reached by DFS
											        // visited2: is on the current path (indicator for cycle)
	        // ...perform DFS
            determineCycledNodes ( tree, nodestack, visited1, visited2 );
            visited2[i] = false;
            nodestack.pop_back ();
        }
    }
    return(1);
}