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

void testIdsDeleteAtom()
{
  OBMol mol;
  for (int i = 0; i < 10; ++i)
    mol.NewAtom();

  OB_REQUIRE( mol.GetAtomById(3) );
  OB_REQUIRE( mol.GetAtomById(4) );
  OB_REQUIRE( mol.GetAtomById(5) );

  mol.DeleteAtom(mol.GetAtomById(4));

  OB_REQUIRE( mol.GetAtomById(3) );
  OB_REQUIRE( mol.GetAtomById(3)->GetId() == 3 );
  OB_REQUIRE( !mol.GetAtomById(4) );
  OB_REQUIRE( mol.GetAtomById(5) );
  OB_REQUIRE( mol.GetAtomById(5)->GetId() == 5 );
}

void AliasData::DeleteExpandedAtoms(OBMol& mol)
{
  //The atom that carries the AliasData object remains as an Xx atom with no charge;
  //the others are deleted. All the attached hydrogens are also deleted.
  for(unsigned i=0;i<_expandedatoms.size();++i)
  {
    OBAtom* at = mol.GetAtomById(_expandedatoms[i]);
    if(!at)
      continue;
    mol.DeleteHydrogens(at);
    if(at->HasData(AliasDataType))
    {
      at->SetAtomicNum(0);
      at->SetFormalCharge(0);
      at->SetSpinMultiplicity(0);
    }
    else
      mol.DeleteAtom(at);
  }
  _expandedatoms.clear();
}

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

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

//.........这里部分代码省略.........
  //! beginning).
  //! Implements <a href="http://qsar.sourceforge.net/dicts/blue-obelisk/index.xhtml#zmatrixCoordinatesIntoCartesianCoordinates">blue-obelisk:zmatrixCoordinatesIntoCartesianCoordinates</a>
  void InternalToCartesian(std::vector<OBInternalCoord*> &vic,OBMol &mol)
  {
    vector3 n,nn,v1,v2,v3,avec,bvec,cvec;
    double dst = 0.0, ang = 0.0, tor = 0.0;
    OBAtom *atom;
    vector<OBAtom*>::iterator i;
    unsigned int index;

    if (vic.empty())
      return;

    if (vic[0] != NULL) {
      std::vector<OBInternalCoord*>::iterator it = vic.begin();
      vic.insert(it, static_cast<OBInternalCoord*>(NULL));
    }

    if (vic.size() != mol.NumAtoms() + 1) {
      string error = "Number of internal coordinates is not the same as";
      error += " the number of atoms in molecule";
      obErrorLog.ThrowError(__FUNCTION__, error, obError);
      return;
    }

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

    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      {
        index = atom->GetIdx();

        // make sure we always have valid pointers
        if (index >= vic.size() || !vic[index])
          return;

        if (vic[index]->_a) // make sure we have a valid ptr
          {
            avec = vic[index]->_a->GetVector();
            dst = vic[index]->_dst;
          }
        else
          {
            // atom 1
            atom->SetVector(0.0, 0.0, 0.0);
            continue;
          }

        if (vic[index]->_b)
          {
            bvec = vic[index]->_b->GetVector();
            ang = vic[index]->_ang * DEG_TO_RAD;
          }
        else
          {
            // atom 2
            atom->SetVector(dst, 0.0, 0.0);
            continue;
          }

        if (vic[index]->_c)
          {
            cvec = vic[index]->_c->GetVector();
            tor = vic[index]->_tor * DEG_TO_RAD;
          }
        else
          {
            // atom 3
            cvec = VY;
            tor = 90. * DEG_TO_RAD;
          }

        v1 = avec - bvec;
        v2 = avec - cvec;
        n = cross(v1,v2);
        nn = cross(v1,n);
        n.normalize();
        nn.normalize();

        n  *= -sin(tor);
        nn *= cos(tor);
        v3 = n + nn;
        v3.normalize();
        v3 *= dst * sin(ang);
        v1.normalize();
        v1 *= dst * cos(ang);
        v2 = avec + v3 - v1;

        atom->SetVector(v2);
      }

    // Delete dummy atoms
    vector<OBAtom*> for_deletion;
    FOR_ATOMS_OF_MOL(a, mol)
      if (a->GetAtomicNum() == 0)
        for_deletion.push_back(&(*a));
    for(vector<OBAtom*>::iterator a_it=for_deletion.begin(); a_it!=for_deletion.end(); ++a_it)
      mol.DeleteAtom(*a_it);

  }

//.........这里部分代码省略.........
          v31, v32, v33;

        ifs.getline(buffer,BUFF_SIZE); // v1
        tokenize(vs, buffer);
        v11 = atof(vs.at(3).c_str()) * conv;
        v12 = atof(vs.at(4).c_str()) * conv;
        v13 = atof(vs.at(5).c_str()) * conv;

        ifs.getline(buffer,BUFF_SIZE); // v2
        tokenize(vs, buffer);
        v21 = atof(vs.at(3).c_str()) * conv;
        v22 = atof(vs.at(4).c_str()) * conv;
        v23 = atof(vs.at(5).c_str()) * conv;

        ifs.getline(buffer,BUFF_SIZE); // v3
        tokenize(vs, buffer);
        v31 = atof(vs.at(3).c_str()) * conv;
        v32 = atof(vs.at(4).c_str()) * conv;
        v33 = atof(vs.at(5).c_str()) * conv;

        // Build unit cell
        cell->SetData(vector3(v11,v12,v13),
                      vector3(v21,v22,v23),
                      vector3(v31,v32,v33));
      }

      // Atoms info
      if (strstr(buffer, "ATOMIC_POSITIONS")) {
        // Clear old atoms from pmol
        vector<OBAtom*> toDelete;
        FOR_ATOMS_OF_MOL(a, *pmol)
          toDelete.push_back(&*a);
        for (size_t i = 0; i < toDelete.size(); i++)
          pmol->DeleteAtom(toDelete.at(i));


        // Discover units
        matrix3x3 conv (1);
        tokenize(vs, buffer);

        if (strstr(vs[1].c_str(), "alat")) {
          conv *= (alat * BOHR_TO_ANGSTROM);
        }
        else if (strstr(vs[1].c_str(), "crystal")) {
          // Set to the zero matrix and test below.
          conv = matrix3x3 (0.0);
        }
        // Add others if needed

        // Load new atoms from molecule
        ifs.getline(buffer,BUFF_SIZE); // First entry
        tokenize(vs, buffer);
        int size = vs.size();
        while (size == 4) {
          atomicNum = OBElements::GetAtomicNum(vs[0].c_str());
          x = atof((char*)vs[1].c_str());
          y = atof((char*)vs[2].c_str());
          z = atof((char*)vs[3].c_str());
          // Add atom
          OBAtom *atom = pmol->NewAtom();
          atom->SetAtomicNum(atomicNum);
          vector3 coords (x,y,z);
          if (conv.determinant() == 0.0) { // Fractional coords
            atom->SetVector(cell->FractionalToCartesian(coords));
          }
          else {