C++ OBBond类代码示例

  vector<OBBond*> OBResidue::GetBonds(bool exterior) const
  {
    OBAtom         *atom;
    vector<OBBond*> bonds;
    OBBitVec        idxs;
    unsigned int    sz;

    sz = (unsigned int) _atoms.size();
    for ( unsigned int i = 0 ; i < sz ; ++i )
      {
        atom = _atoms[i];
        OBBond *bond;
        vector<OBBond*>::iterator b;
        for (bond = atom->BeginBond(b) ; bond ; bond = atom->NextBond(b))
          {
            if (!idxs.BitIsOn(bond->GetIdx()))
              {
                if (!exterior)
                  {
                    if (bond->GetNbrAtom(atom)->GetResidue() == this)
                      bonds.push_back(&(*bond));
                  }
                else
                  bonds.push_back(&(*bond));

                idxs.SetBitOn(bond->GetIdx());
              }
          }
      }

    return bonds;
  }

   bool OBBond::IsTertiaryAmide()
   {
      OBAtom *c,*n;
      c = n = NULL;

      // Look for C-N bond
      if (_bgn->GetAtomicNum() == 6 && _end->GetAtomicNum() == 7)
      {
         c = (OBAtom*)_bgn;
         n = (OBAtom*)_end;
      }
      if (_bgn->GetAtomicNum() == 7 && _end->GetAtomicNum() == 6)
      {
         c = (OBAtom*)_end;
         n = (OBAtom*)_bgn;
      }
      if (!c || !n) return(false);
      if (GetBondOrder() != 1) return(false);
      if (n->GetImplicitValence() != 3) return(false);

      // Make sure that N is connected to three non-H atoms
      if (n->GetHvyValence() != 3) return(false);

      // Make sure C is attached to =O
      OBBond *bond;
      vector<OBBond*>::iterator i;
      for (bond = c->BeginBond(i); bond; bond = c->NextBond(i))
      {
         if (bond->IsCarbonyl()) return(true);
      }

      return(false);
   }

  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);
  }

//preprocess molecule into a standardized state for heavy atom rmsd computation
static void processMol(OBMol& mol)
{
	//isomorphismmapper wants isomorphic atoms to have the same aromatic and ring state,
	//but these proporties aren't reliable enough to be trusted in evaluating molecules
	//should be considered the same based solely on connectivity
	
	mol.DeleteHydrogens(); //heavy atom rmsd
	for(OBAtomIterator aitr = mol.BeginAtoms(); aitr != mol.EndAtoms(); aitr++)
	{
		OBAtom *a = *aitr;
		a->UnsetAromatic();
		a->SetInRing();
	}
	for(OBBondIterator bitr = mol.BeginBonds(); bitr != mol.EndBonds(); bitr++)
	{
		OBBond *b = *bitr;
		b->UnsetAromatic();
		b->SetBondOrder(1);
		b->SetInRing();
	}
	//avoid recomputations
	mol.SetHybridizationPerceived();
	mol.SetRingAtomsAndBondsPerceived();
	mol.SetAromaticPerceived();
}

void fingerprint2::getFragments(vector<int> levels, vector<int> curfrag,
					int level, OBAtom* patom, OBBond* pbond)
{
	//Recursive routine to analyse schemical structure and populate fragset and ringset
	//Hydrogens,charges(except dative bonds), spinMultiplicity ignored
	const int Max_Fragment_Size = 7;
	int bo=0;
	if(pbond)
	{
		bo = pbond->IsAromatic() ? 5 : pbond->GetBO();

//		OBAtom* pprevat = pbond->GetNbrAtom(patom);
//		if(patom->GetFormalCharge() && (patom->GetFormalCharge() == -pprevat->GetFormalCharge()))
//			++bo; //coordinate (dative) bond eg C[N+]([O-])=O is seen as CN(=O)=O
	}
	curfrag.push_back(bo);
	curfrag.push_back(patom->GetAtomicNum());
	levels[patom->GetIdx()-1] = level;

	vector<OBEdgeBase*>::iterator itr;
	OBBond *pnewbond;
//	PrintFpt(curfrag,(int)patom);
	for (pnewbond = patom->BeginBond(itr);pnewbond;pnewbond = patom->NextBond(itr))
	{
		if(pnewbond==pbond) continue; //don't retrace steps
		OBAtom* pnxtat = pnewbond->GetNbrAtom(patom);
		if(pnxtat->GetAtomicNum() == OBElements::Hydrogen) continue;

		int atlevel = levels[pnxtat->GetIdx()-1];
		if(atlevel) //ring
		{
			if(atlevel==1)
			{
				//If complete ring (last bond is back to starting atom) add bond at front
				//and save in ringset
				curfrag[0] = pnewbond->IsAromatic() ? 5 : pnewbond->GetBO();
				ringset.insert(curfrag);
 				curfrag[0] = 0;
			}
		}
		else //no ring
		{
			if(level<Max_Fragment_Size)
			{
//				TRACE("level=%d size=%d %p frag[0]=%p\n",level, curfrag.size(),&curfrag, &(curfrag[0]));
				//Do the next atom; levels, curfrag are passed by value and hence copied
				getFragments(levels, curfrag, level+1, pnxtat, pnewbond);
			}
		}
	}

	//do not save C,N,O single atom fragments
	if(curfrag[0]==0 &&
		(level>1 || patom->GetAtomicNum()>8  || patom->GetAtomicNum()<6))
	{
		fragset.insert(curfrag); //curfrag ignored if an identical fragment already present
//		PrintFpt(curfrag,level);
	}
}

  //! Sets atom->IsChiral() to true for chiral atoms
  void OBMol::FindChiralCenters()
  {
    if (HasChiralityPerceived())
      return;
    SetChiralityPerceived();

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

    //do quick test to see if there are any possible chiral centers
    bool mayHaveChiralCenter=false;
    OBAtom *atom,*nbr;
    vector<OBAtom*>::iterator i;
    for (atom = BeginAtom(i);atom;atom = NextAtom(i))
      if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3)
        {
          mayHaveChiralCenter=true;
          break;
        }

    if (!mayHaveChiralCenter)
      return;

    OBBond *bond;
    vector<OBBond*>::iterator j;
    for (bond = BeginBond(j);bond;bond = NextBond(j))
      if (bond->IsWedge() || bond->IsHash())
        (bond->GetBeginAtom())->SetChiral();

    vector<unsigned int> vgid;
    GetGIDVector(vgid);
    vector<unsigned int> tlist;
    vector<unsigned int>::iterator k;

    bool ischiral;
    for (atom = BeginAtom(i);atom;atom = NextAtom(i))
      if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3 && !atom->IsChiral())
        {
          tlist.clear();
          ischiral = true;

          for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
            {
              for (k = tlist.begin();k != tlist.end();++k)
                if (vgid[nbr->GetIdx()-1] == *k)
                  ischiral = false;

              if (ischiral)
                tlist.push_back(vgid[nbr->GetIdx()-1]);

              else
                break;
            }

          if (ischiral)
            atom->SetChiral();
        }
  }

  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;
  }

void GenerateRingReference()
{
  std::ifstream ifs;
  if (!SafeOpen(ifs,"attype.00.smi")) return;

  std::ofstream ofs;
  if (!SafeOpen(ofs,"ringresults.txt")) return;

  int count;
  OBAtom *atom;
  OBBond *bond;
  char buffer[BUFF_SIZE];
  vector<OBRing*> vr;
  vector<OBEdgeBase*>::iterator i;
  vector<OBNodeBase*>::iterator j;
  vector<OBRing*>::iterator k;
  OBMol mol(SMI,SMI);
  OBFileFormat ff;

  for (;ifs;)
    {
      mol.Clear();
      ff.ReadMolecule(ifs, mol);
      if (mol.Empty()) continue;
      
      //write out ring bonds
      for (bond = mol.BeginBond(i);bond;bond = mol.NextBond(i))
	if (bond->IsInRing())
	  {
	    sprintf(buffer,"%3d",bond->GetIdx());
	    ofs << buffer;
	  }
      ofs << endl;

      vr = mol.GetSSSR();
      //write the total number of rings
      ofs << vr.size() << endl;

      //write the number of rings that each atom is a member of
      for (atom = mol.BeginAtom(j);atom;atom = mol.NextAtom(j))
	{
	  count = 0;
	  for (k = vr.begin();k != vr.end();k++)
	    if ((*k)->_pathset[atom->GetIdx()])
	      count++;

	  sprintf(buffer,"%3d",count);
	  ofs << buffer;
	}
      ofs << endl;

    }

  ThrowError("Ring perception test results written successfully");
}

  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 GetDFFVector(OBMol &mol,vector<int> &dffv,OBBitVec &bv)
  {
    dffv.clear();
    dffv.resize(mol.NumAtoms());

    int dffcount,natom;
    OBBitVec used,curr,next;
    OBAtom *atom,*atom1;
    OBBond *bond;
    vector<OBAtom*>::iterator i;
    vector<OBBond*>::iterator j;

    next.Clear();

    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      {
        if (bv[atom->GetIdx()])
          {
            dffv[atom->GetIdx()-1] = 0;
            continue;
          }

        dffcount = 0;
        used.Clear();
        curr.Clear();
        used.SetBitOn(atom->GetIdx());
        curr.SetBitOn(atom->GetIdx());

        while (!curr.IsEmpty() && (bv&curr).Empty())
          {
            next.Clear();
            for (natom = curr.NextBit(-1);natom != curr.EndBit();natom = curr.NextBit(natom))
              {
                atom1 = mol.GetAtom(natom);
                for (bond = atom1->BeginBond(j);bond;bond = atom1->NextBond(j))
                  if (!used.BitIsOn(bond->GetNbrAtomIdx(atom1)) &&
                      !curr.BitIsOn(bond->GetNbrAtomIdx(atom1)))
                    if (!(bond->GetNbrAtom(atom1))->IsHydrogen())
                      next.SetBitOn(bond->GetNbrAtomIdx(atom1));
              }

            used |= next;
            curr = next;
            dffcount++;
          }

        dffv[atom->GetIdx()-1] = dffcount;
      }

    return(true);
  }

  void OBDepictPrivate::SetWedgeAndHash(OBMol* mol)  {
    std::map<OBBond*, enum OBStereo::BondDirection> updown;
    std::map<OBBond*, OBStereo::Ref> from;
    std::map<OBBond*, OBStereo::Ref>::const_iterator from_cit;
    TetStereoToWedgeHash(*mol, updown, from);

    for(from_cit=from.begin();from_cit!=from.end();++from_cit) {
      OBBond* pbond = from_cit->first;
      if(updown[pbond]==OBStereo::UpBond)
        pbond->SetHash();
      else if(updown[pbond]==OBStereo::DownBond)
        pbond->SetWedge();
    }
  }

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 OBRotorList::RemoveSymVals(OBMol &mol)
  {
    OBGraphSym gs(&mol);
    vector<unsigned int> sym_classes;
    gs.GetSymmetry(sym_classes);

    OBRotor *rotor;
    vector<OBRotor*>::iterator i;
    std::set<unsigned int> syms;
    for (rotor = BeginRotor(i);rotor;rotor = NextRotor(i)) {
      OBBond* bond = rotor->GetBond();
      OBAtom* end = bond->GetEndAtom();
      OBAtom* begin = bond->GetBeginAtom();
      int N_fold_symmetry = 1;
      for (int here=0; here <= 1; ++here) { // Try each side of the bond in turn
        
        OBAtom *this_side, *other_side;
        if (here == 0) {
          this_side = begin; other_side = end;
        }
        else {
          this_side = end; other_side = begin;
        }

        for (int hyb=2; hyb<=3; ++hyb) { // sp2 and sp3 carbons, with explicit Hs
          if (this_side->GetAtomicNum() == 6 && this_side->GetHyb() == hyb && this_side->GetValence() == (hyb + 1) ) {
            syms.clear();
            FOR_NBORS_OF_ATOM(nbr, this_side) {
              if ( &(*nbr) == other_side ) continue;
              syms.insert(sym_classes[nbr->GetIdx() - 1]);
            }
            if (syms.size() == 1) // All of the rotated atoms have the same symmetry class
              N_fold_symmetry *= hyb;
          }
        }
      }

      if (N_fold_symmetry  > 1) {
        size_t old_size = rotor->Size();
        rotor->RemoveSymTorsionValues(N_fold_symmetry);
        if (!_quiet) {
          cout << "...." << N_fold_symmetry << "-fold symmetry at rotor between " << 
                 begin->GetIdx() << " and " << end->GetIdx();
          cout << " - reduced from " << old_size << " to " << rotor->Size() << endl;
        }
      }
    }

// OBMol::NewBond()
void testIdsNewBond1()
{
  OBMol mol;
  for (int i = 0; i < 10; ++i) {
    OBBond *bond = mol.NewBond();
    OB_REQUIRE(bond->GetId() == i);
  }

  OB_REQUIRE( mol.GetBondById(0) );
  OB_REQUIRE( mol.GetBondById(4) );
  OB_REQUIRE( mol.GetBondById(9) );
  OB_REQUIRE( !mol.GetBondById(10) );

  OB_REQUIRE( mol.GetBondById(0)->GetId() == 0 );
  OB_REQUIRE( mol.GetBondById(4)->GetId() == 4 );
  OB_REQUIRE( mol.GetBondById(9)->GetId() == 9 );
}