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

  void FastSearchFormat::AddPattern(vector<OBMol>& patternMols, OBMol patternMol, int idx)
  {
    //Recursive function to generate all combinations of aromatic/single bonds for each tilde bond
    //Copying an OBMol, which happens when adding it to a vector, kekulizes it,
    // changing aromatic (bo=5) bonds. So set order after adding. Should work here,
    // but is dangerous if the vector needs to be reallocated.

    if(idx>=patternMol.NumBonds())
      return;
    if(patternMol.GetBond(idx)->GetBO()==4)
    {
      patternMol.GetBond(idx)->SetBO(1);
      patternMols.push_back(patternMol);
      AddPattern(patternMols, patternMol,idx+1);

      patternMols.push_back(patternMol);
      patternMols.back().GetBond(idx)->SetBO(5);
    }
    AddPattern(patternMols, patternMol,idx+1);
  }

  static void FindRings(OBMol &mol,vector<int> &path,OBBitVec &avisit,
                        OBBitVec &bvisit, int natom,int depth )
  {
    OBAtom *atom;
    OBBond *bond;
    vector<OBBond*>::iterator k;

    // don't return if all atoms are visited
    // (For example, some atoms are in multiple rings!) -GRH
      
    if (avisit[natom])
      {
        int j = depth-1;
        bond=mol.GetBond(path[j--]);
        bond->SetInRing();
        while( j >= 0 )
          {
            bond=mol.GetBond(path[j--]);
            bond->SetInRing();
            (bond->GetBeginAtom())->SetInRing();
            (bond->GetEndAtom())->SetInRing();
            if(bond->GetBeginAtomIdx()==static_cast<unsigned int>(natom) || bond->
               GetEndAtomIdx()==static_cast<unsigned int>(natom))
              break;
          }
      }
    else
      {
        avisit.SetBitOn(natom);
        atom = mol.GetAtom(natom);
        for(bond = atom->BeginBond(k);bond;bond=atom->NextBond(k))
          if( !bvisit[bond->GetIdx()])
            {
              path[depth] = bond->GetIdx();
              bvisit.SetBitOn(bond->GetIdx());
              FindRings(mol,path,avisit,bvisit,bond->GetNbrAtomIdx(atom),
                        depth+1);
            }
      }
  }

  void OutputCSTBonds(ostream &ofs, OBMol &mol, string prefix)
  {
    /* ---- Write povray-description of all bonds---- */
    for(unsigned int i = 0; i < mol.NumBonds(); ++i)
      {

        double x1,y1,z1,x2,y2,z2; /* Start and stop coordinates of a bond       */
        double dist;              /* Distance between (x1|y1|z1) and (x2|y2|z2) */
        double phi,theta;         /* Angles between (x1|y1|z1) and (x2|y2|z2)   */
        double dy;                /* Distance between (x1|0|z1) and (x2|0|z2)   */

        /* ---- Get a pointer to ith atom ---- */
        OBBond *bond = mol.GetBond(i);

        /* ---- Assign start of bond i ---- */
        x1 = (bond -> GetBeginAtom()) -> GetX();
        y1 = (bond -> GetBeginAtom()) -> GetY();
        z1 = (bond -> GetBeginAtom()) -> GetZ();

        /* ---- Assign end of bond i ---- */
        x2 = (bond -> GetEndAtom()) -> GetX();
        y2 = (bond -> GetEndAtom()) -> GetY();
        z2 = (bond -> GetEndAtom()) -> GetZ();

        /* ---- Calculate length of bond and (x1|0|z1) - (x2|0|z2) ---- */
        dist = sqrt(SQUARE(x2-x1) + SQUARE(y2-y1) + SQUARE(z2-z1));
        dy = sqrt(SQUARE(x2-x1) + SQUARE(z2-z1));

        /* ---- Calculate Phi and Theta ---- */
        phi = (double) 0.0;
        theta = (double) 0.0;
        if (fabs(dist) >= EPSILON)
          phi = acos((y2-y1)/dist);
        if (fabs(dy) >= EPSILON)
          theta = acos((x2-x1)/dy);

        /* ---- Begin of description of bond i (for a capped sticks model) ---- */
        ofs << "#declare " << prefix << "_bond" << i << " = object {" << endl;
        ofs << "\t  union {" << endl;

        /* ---- Begin of Start-Half of Bond (i) ---- */
        ofs << "\t   object {" << endl << "\t    bond_" << bond -> GetBondOrder()  << "\n";

        /* ---- Add a pigment - statement for start-atom of bond ---- */
        ofs << "\t    pigment{color Color_"
            << bond -> GetBeginAtom() -> GetType()
            << "}" << endl;

        /* ---- Scale bond if needed ---- */
        if (fabs((double) 2.0 * dist) >= EPSILON)
          {

            /* ---- Print povray scale-statement (x-Axis) ---- */
            ofs << "\t    scale <" << (double) 0.5 * dist << ",1.0000,1.0000>" << endl;

          }

        /* ---- Rotate (Phi) bond if needed ---- */
        if (fabs(RAD2DEG(-phi) + (double) 90.0) >= EPSILON)
          {

            /* ---- Rotate along z-axis ---- */
            ofs << "\t    rotate <0.0000,0.0000,"
                << RAD2DEG(-phi) + (double) 90.0
                << ">" << endl;

          }

        /* ---- Check angle between (x1|0|z1) and (x2|0|z2) ---- */
        if (theta >= EPSILON)
          {

            /* ---- Check direction ---- */
            if ((z2 - z1) >= (double) 0.0)
              {

                /* ---- Rotate along y-Axis (negative) ---- */
                ofs << "\t    rotate <0.0000,"
                    << RAD2DEG((double) -1.0 *theta) << ",0.0000>"
                    << endl;

              }
            else
              {

                /* ---- Rotate along y-Axis (positive) ---- */
                ofs << "\t    rotate <0.0000," << RAD2DEG(theta) << ",0.0000>" << endl;

              }

          }

        /* ---- Translate bond to start ---- */

        ofs << "\t    translate " << prefix << "_pos_" << bond -> GetBeginAtomIdx() << endl;

        /* ---- End of description of Start-Bond ---- */
        ofs << "\t   }" << endl;

        /* ---- Begin of End-Half of Bond i ---- */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                                GetSerialNum(a1)) == startAtomSerialNumber)
        {
          firstAtom = a1;
          break;
        }
    }

    if (firstAtom == NULL)
      {
        errorMsg << "WARNING: Problems reading a PDB file:\n"
                 << "  Problems reading a CONECT record.\n"
                 << "  According to the PDB specification,\n"
                 << "  columns 7-11 should contain the serial number of an atom.\n"
                 << "  No atom was found with this serial number.\n"
                 << "  THIS CONECT RECORD WILL BE IGNORED." << endl;
        obErrorLog.ThrowError(__FUNCTION__, errorMsg.str() , obWarning);
        return(false);
      }

    if (mol.NumAtoms() < 9999)
      {
        if (vs.size() > 1) boundedAtomsSerialNumbers[0] = atoi(vs[1].c_str());
        if (vs.size() > 2) boundedAtomsSerialNumbers[1] = atoi(vs[2].c_str());
        if (vs.size() > 3) boundedAtomsSerialNumbers[2] = atoi(vs[3].c_str());
        if (vs.size() > 4) boundedAtomsSerialNumbers[3] = atoi(vs[4].c_str());

        unsigned int limit = 4;
        if (vs.size() <= 4)
          limit = vs.size() - 1;

        for (unsigned int s = 0; s < limit; ++s)
          boundedAtomsSerialNumbersValid[s] = true;
      }
    else
      {
        // Now read the serial numbers. If the first serial number is not
        // present, this connect record probably contains only hydrogen
        // bonds and salt bridges, which we ignore. In that case, we just
        // exit gracefully.
        boundedAtomsSerialNumbersValid[0] = readIntegerFromRecord(buffer, 12, boundedAtomsSerialNumbers+0);
        if (boundedAtomsSerialNumbersValid[0] == false)
          return(true);
        boundedAtomsSerialNumbersValid[1] = readIntegerFromRecord(buffer, 17, boundedAtomsSerialNumbers+1);
        boundedAtomsSerialNumbersValid[2] = readIntegerFromRecord(buffer, 22, boundedAtomsSerialNumbers+2);
        boundedAtomsSerialNumbersValid[3] = readIntegerFromRecord(buffer, 27, boundedAtomsSerialNumbers+3);
      }

    // Now iterate over the VALID boundedAtomsSerialNumbers and connect
    // the atoms.
    for(unsigned int k=0; boundedAtomsSerialNumbersValid[k]; k++)
      {
        // Find atom that is connected to, write an error message
        OBAtom *connectedAtom = 0L;
        for (OBAtom *a1 = mol.BeginAtom(i);a1;a1 = mol.NextAtom(i)) {
          // again, atoms may not have residues, but if they do, check serials
          if (a1->GetResidue() != NULL &&
              static_cast<long int>(a1->GetResidue()->
                                    GetSerialNum(a1)) == boundedAtomsSerialNumbers[k])
            {
              connectedAtom = a1;
              break;
            }
        }
        if (connectedAtom == 0L)
          {
            errorMsg << "WARNING: Problems reading a PDB file:\n"
                     << "  Problems reading a CONECT record.\n"
                     << "  According to the PDB specification,\n"
                     << "  Atoms with serial #" << startAtomSerialNumber
                     << " and #" << boundedAtomsSerialNumbers[k]
                     << " should be connected\n"
                     << "  However, an atom with serial #" << boundedAtomsSerialNumbers[k] << " was not found.\n"
                     << "  THIS CONECT RECORD WILL BE IGNORED." << endl;
            obErrorLog.ThrowError(__FUNCTION__, errorMsg.str() , obWarning);
            return(false);
          }

        // Figure the bond order
        unsigned char order = 0;
        while(boundedAtomsSerialNumbersValid[k+order+1] && (boundedAtomsSerialNumbers[k+order]
                                                            == boundedAtomsSerialNumbers[k+order+1]))
          order++;
        k += order;

        // Generate the bond
        if (firstAtom->GetIdx() < connectedAtom->GetIdx()) { // record the bond 'in one direction' only
          OBBond *bond = mol.GetBond(firstAtom, connectedAtom);
          if (!bond)
            mol.AddBond(firstAtom->GetIdx(), connectedAtom->GetIdx(), order+1);
          else // An additional CONECT record with the same firstAtom that references
               // a bond created in the previous CONECT record.
               // For example, the 1136->1138 double bond in the following:
               //   CONECT 1136 1128 1137 1137 1138
               //   CONECT 1136 1138 1139
            bond->SetBondOrder(bond->GetBondOrder() + order+1);
        }

      }
    return(true);
  }

int main(int argc,char *argv[])
{
  // turn off slow sync with C-style output (we don't use it anyway).
  std::ios::sync_with_stdio(false);

  OBConversion conv;
  OBFormat *inFormat, *canFormat;
  OBMol mol;
  ifstream ifs;
  vector<OBMol> fragments;
  unsigned int fragmentCount = 0; // track how many in library -- give a running count
  map<string, int> index; // index of cansmi
  string currentCAN;
  unsigned int size;
  OBAtom *atom;
  OBBond *bond;
  bool nonRingAtoms, nonRingBonds;
  char buffer[BUFF_SIZE];

  canFormat = conv.FindFormat("can");
  conv.SetOutFormat(canFormat);

  if (argc < 2)
    {
      cout << "Usage: obfragment <file>" << endl;
      return(-1);
    }

  for (int i = 1; i < argc; i++) {
    cerr << " Reading file " << argv[i] << endl;

    inFormat = conv.FormatFromExt(argv[i]);
    if(inFormat==NULL || !conv.SetInFormat(inFormat))
      {
        cerr << " Cannot read file format for " << argv[i] << endl;
        continue; // try next file
      }
    
    ifs.open(argv[i]);
    
    if (!ifs)
      {
        cerr << "Cannot read input file: " << argv[i] << endl;
        continue;
      }
    
    
    while(ifs.peek() != EOF && ifs.good())
      {
        conv.Read(&mol, &ifs);
        if (!mol.Has3D()) continue; // invalid coordinates!
        mol.DeleteHydrogens(); // remove these before we do anything else
        
        do {
          nonRingAtoms = false;
          size = mol.NumAtoms();
          for (unsigned int i = 1; i <= size; ++i)
            {
              atom = mol.GetAtom(i);
              if (!atom->IsInRing()) {
                mol.DeleteAtom(atom);
                nonRingAtoms = true;
                break; // don't know how many atoms there are
              } 
              // Previously, we changed atoms to carbon here.
              // Now we perform this alchemy in terms of string-rewriting
              // once the canonical SMILES is generated
            }
        } while (nonRingAtoms);
        
        if (mol.NumAtoms() < 3)
          continue;
        
        if (mol.NumBonds() == 0)
          continue;
        
        do {
          nonRingBonds = false;
          size = mol.NumBonds();
          for (unsigned int i = 0; i < size; ++i)
            {
              bond = mol.GetBond(i);
              if (!bond->IsInRing()) {
                mol.DeleteBond(bond);
                nonRingBonds = true;
                break; // don't know how many bonds there are
              }
            }        
        } while (nonRingBonds);

        fragments = mol.Separate();
        for (unsigned int i = 0; i < fragments.size(); ++i)
          {
            if (fragments[i].NumAtoms() < 3) // too small to care
              continue;
              
            currentCAN = conv.WriteString(&fragments[i], true);
            currentCAN = RewriteSMILES(currentCAN); // change elements to "a/A" for compression
            if (index.find(currentCAN) != index.end()) { // already got this
              index[currentCAN] += 1; // add to the count for bookkeeping
//.........这里部分代码省略.........

  bool OBChemTsfm::Apply(OBMol &mol)
  {
    if (!_bgn.Match(mol))
      return(false);
    mol.BeginModify();
    vector<vector<int> > mlist = _bgn.GetUMapList();

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

    if (!_vchrg.empty()) //modify charges
      {
        vector<vector<int> >::iterator i;
        vector<pair<int,int> >::iterator j;

        for (i = mlist.begin();i != mlist.end();++i)
          for (j = _vchrg.begin();j != _vchrg.end();++j)
            if (j->first < (signed)i->size()) { //goof proofing
              OBAtom *atom = mol.GetAtom((*i)[j->first]);
              int old_charge = atom->GetFormalCharge();
              atom->SetFormalCharge(j->second);
              int new_hcount = atom->GetImplicitHCount() + (j->second - old_charge);
              if (new_hcount < 0)
                new_hcount = 0;
              atom->SetImplicitHCount(new_hcount);
            }
      }

    if (!_vbond.empty()) //modify bond orders
      {
        OBBond *bond;
        vector<vector<int> >::iterator i;
        vector<pair<pair<int,int>,int> >::iterator j;
        for (i = mlist.begin();i != mlist.end();++i)
          for (j = _vbond.begin();j != _vbond.end();++j)
            {
              bond = mol.GetBond((*i)[j->first.first],(*i)[j->first.second]);
              if (!bond)
                {
                  obErrorLog.ThrowError(__FUNCTION__, "unable to find bond", obDebug);
                  continue;
                }
              unsigned int old_bond_order = bond->GetBondOrder();
              bond->SetBondOrder(j->second);
              for (int k = 0; k < 2; ++k) {
                OBAtom* atom = k == 0 ? bond->GetBeginAtom() : bond->GetEndAtom();
                int new_hcount = atom->GetImplicitHCount() - (j->second - old_bond_order);
                if (new_hcount < 0)
                  new_hcount = 0;
                atom->SetImplicitHCount(new_hcount);
              }
            }
      }

    if (!_vadel.empty() || !_vele.empty()) //delete atoms and change elements
      {
        vector<int>::iterator j;
        vector<vector<int> >::iterator i;

        if (!_vele.empty())
          {
            vector<pair<int,int> >::iterator k;
            for (i = mlist.begin();i != mlist.end();++i)
              for (k = _vele.begin();k != _vele.end();++k)
                mol.GetAtom((*i)[k->first])->SetAtomicNum(k->second);
          }

        //make sure same atom isn't deleted twice
        vector<bool> vda;
        vector<OBAtom*> vdel;
        vda.resize(mol.NumAtoms()+1,false);
        for (i = mlist.begin();i != mlist.end();++i)
          for (j = _vadel.begin();j != _vadel.end();++j)
            if (!vda[(*i)[*j]])
              {
                vda[(*i)[*j]] = true;
                vdel.push_back(mol.GetAtom((*i)[*j]));
              }

        vector<OBAtom*>::iterator k;
        for (k = vdel.begin();k != vdel.end();++k)
          mol.DeleteAtom((OBAtom*)*k);
      }

    mol.EndModify();
    return(true);
  }

  void OBRotamerList::Setup(OBMol &mol,OBRotorList &rl)
  {
    //clear the old stuff out if necessary
    _vres.clear();
    vector<unsigned char*>::iterator j;
    for (j = _vrotamer.begin();j != _vrotamer.end();++j)
      delete [] *j;
    _vrotamer.clear();

    vector<pair<OBAtom**,vector<int> > >::iterator k;
    for (k = _vrotor.begin();k != _vrotor.end();++k)
      delete [] k->first;
    _vrotor.clear();
    _vrings.clear();
    _vringTors.clear();

    //create the new list
    OBRotor *rotor;
    vector<OBRotor*>::iterator i;
    vector<int> children;

    int ref[4];
    OBAtom **atomlist;
    for (rotor = rl.BeginRotor(i);rotor;rotor = rl.NextRotor(i))
      {
        atomlist = new OBAtom* [4];
        rotor->GetDihedralAtoms(ref);
        atomlist[0] = mol.GetAtom(ref[0]);
        atomlist[1] = mol.GetAtom(ref[1]);
        atomlist[2] = mol.GetAtom(ref[2]);
        atomlist[3] = mol.GetAtom(ref[3]);
        mol.FindChildren(children,ref[1],ref[2]);
        _vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children));
        _vres.push_back(rotor->GetResolution());
      }

    // if the rotor list has ring bonds, build up an index
    if (rl.HasRingRotors()){
      // go through rings
      // for each step of the path, see if there's a matching rotor
      vector<int> path;
      int pSize;
      vector<double> ringTorsions;
      vector<int> ringRotors;
      FOR_RINGS_OF_MOL(r, mol)
      {
        ringTorsions.clear();
        ringRotors.clear();

        pSize = r->Size();
        if (pSize < 4)
          continue; // not rotatable

        path = r->_path;
        for (int j = 0; j < pSize; ++j) {
          double torsion = mol.GetTorsion(path[(j + pSize - 1) % pSize],
                                       path[(j + pSize) % pSize],
                                       path[(j + pSize + 1) % pSize],
                                       path[(j + pSize + 2) % pSize]);
          ringTorsions.push_back(torsion);

          // now check to see if any of these things are rotors
          int rotorIndex = -1; // not a rotor
          OBBond *bond = mol.GetBond(path[(j + pSize) % pSize], path[(j + pSize + 1) % pSize]);
          for (rotor = rl.BeginRotor(i);rotor;rotor = rl.NextRotor(i))
            {
              if (bond != rotor->GetBond())
                continue; // no match at all

              // Central bond matches, make sure 1..4 atoms are in the path
              rotor->GetDihedralAtoms(ref);
              if ( (ref[0] == path[(j + pSize - 1) % pSize] &&
                    ref[3] == path[(j + pSize + 2) % pSize])
                   ||
                   (ref[3] == path[(j + pSize - 1) % pSize] &&
                    ref[0] == path[(j + pSize + 2) % pSize]) ) {
                rotorIndex = rotor->GetIdx();
              }
            } // end checking all the rotors
          ringRotors.push_back(rotorIndex); // could be -1 if it's not rotatable
        }

        _vringTors.push_back(ringTorsions);
        _vrings.push_back(ringRotors);
      } // finished with the rings

  bool OBChemTsfm::Apply(OBMol &mol)
  {
    if (!_bgn.Match(mol))
      return(false);

    vector<vector<int> > mlist = _bgn.GetUMapList();

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

    if (!_vchrg.empty()) //modify charges
      {
        vector<vector<int> >::iterator i;
        vector<pair<int,int> >::iterator j;

        for (i = mlist.begin();i != mlist.end();++i)
          for (j = _vchrg.begin();j != _vchrg.end();++j)
            if (j->first < (signed)i->size()) //goof proofing
              mol.GetAtom((*i)[j->first])->SetFormalCharge(j->second);

        mol.UnsetImplicitValencePerceived();
      }

    if (!_vbond.empty()) //modify bond orders
      {
        OBBond *bond;
        vector<vector<int> >::iterator i;
        vector<pair<pair<int,int>,int> >::iterator j;
        for (i = mlist.begin();i != mlist.end();++i)
          for (j = _vbond.begin();j != _vbond.end();++j)
            {
              bond = mol.GetBond((*i)[j->first.first],(*i)[j->first.second]);
              if (!bond)
                {
                  obErrorLog.ThrowError(__FUNCTION__, "unable to find bond", obDebug);
                  continue;
                }

              bond->SetBO(j->second);
            }
      }

    if (!_vadel.empty() || !_vele.empty()) //delete atoms and change elements
      {
        vector<int>::iterator j;
        vector<vector<int> >::iterator i;

        if (!_vele.empty())
          {
            vector<pair<int,int> >::iterator k;
            for (i = mlist.begin();i != mlist.end();++i)
              for (k = _vele.begin();k != _vele.end();++k)
                mol.GetAtom((*i)[k->first])->SetAtomicNum(k->second);
          }

        //make sure same atom isn't deleted twice
        vector<bool> vda;
        vector<OBAtom*> vdel;
        vda.resize(mol.NumAtoms()+1,false);
        for (i = mlist.begin();i != mlist.end();++i)
          for (j = _vadel.begin();j != _vadel.end();++j)
            if (!vda[(*i)[*j]])
              {
                vda[(*i)[*j]] = true;
                vdel.push_back(mol.GetAtom((*i)[*j]));
              }

        vector<OBAtom*>::iterator k;
        for (k = vdel.begin();k != vdel.end();++k)
          mol.DeleteAtom((OBAtom*)*k);
      }

    return(true);
  }