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

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

bool MakeQueriesFromMolInFile(vector<OBQuery*>& queries, const std::string& filename, int* pnAtoms, bool noH)
{
    OBMol patternMol;
    patternMol.SetIsPatternStructure();
    OBConversion patternConv;
    OBFormat* pFormat;
    //Need to distinguish between filename and SMARTS. Not infallable...
    if( filename.empty() ||
        filename.find('.')==string::npos ||
        !(pFormat = patternConv.FormatFromExt(filename.c_str())) ||
        !patternConv.SetInFormat(pFormat) ||
        !patternConv.ReadFile(&patternMol, filename) ||
        patternMol.NumAtoms()==0)
      return false;

    if(noH)
      patternMol.DeleteHydrogens();

    do
    {
      *pnAtoms = patternMol.NumHvyAtoms();
      queries.push_back(CompileMoleculeQuery(&patternMol));
    }while(patternConv.Read(&patternMol));
    return true;
}

// Delete hydrogens should not remove charged or isotopic hydrogens or [H][H] or [Cu][H][Cu]
// or hydrogens with assigned atom classes
void test_Issue178_DeleteHydrogens()
{
  OBConversion conv;
  conv.SetInFormat("smi");
  OBMol mol;
  // Test DeleteHydrogens() and DeleteNonPolarHydrogens()
  static const char *smi[] = { "C[H]", "[H][H]", "C[1H]", "C[H]C", "C[H+]" };
  int numHs[] = { 0, 2, 1, 1, 1 };
  for (int i = 0; i < 5; ++i) {
    for (int j = 0; j < 2; ++j) {
      conv.ReadString(&mol, smi[i]);
      if (j == 0)
        mol.DeleteHydrogens();
      else
        mol.DeleteNonPolarHydrogens();
      int myNumHs = 0;
      FOR_ATOMS_OF_MOL(atom, mol)
        if (atom->IsHydrogen())
          myNumHs++;
      OB_COMPARE(myNumHs, numHs[i]);
    }
  }
  // Test DeletePolarHydrogens()
  static const char *smiB[] = { "N[H]", "[H][H]", "N[1H]", "N[H]C", "N[H+]" };
  int numHsB[] = { 0, 2, 1, 1, 1 };
  for (int i = 0; i < 5; ++i) {
    conv.ReadString(&mol, smiB[i]);
    mol.DeletePolarHydrogens();
    int myNumHs = 0;
    FOR_ATOMS_OF_MOL(atom, mol)
      if (atom->IsHydrogen())
        myNumHs++;
    OB_COMPARE(myNumHs, numHsB[i]);
  }
  // Test atom class
  // Currently, the SMILES parser does not retain atom classes for hydrogens on reading so...
  conv.ReadString(&mol, "C[H]");
  OBAtomClassData *ac = new OBAtomClassData;
  ac->Add(2, 99); // Assign the hydrogen (atom 2) a class of 99
  mol.SetData(ac);
  mol.DeleteHydrogens();
  int myNumHs = 0;
  FOR_ATOMS_OF_MOL(atom, mol)
    if (atom->IsHydrogen())
      myNumHs++;
  OB_COMPARE(myNumHs, 1);
}

  void OBPhModel::CorrectForPH(OBMol &mol, double pH)
  {
    if (!_init)
      Init();
    if (mol.IsCorrectedForPH())
      return;
    if (mol.GetDimension() > 0 && !mol.AutomaticFormalCharge())
      return;

    mol.SetCorrectedForPH();

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

    mol.DeleteHydrogens();

    for (unsigned int i = 0; i < _vtsfm.size(); ++i) {

      if (_vpKa[i] > 1E+9) {
        // always apply when pKa is > 1e+9
        _vtsfm[i]->Apply(mol);
      } else {
        // 10^(pKa - pH) = [HA] / [A-]
        //
        // > 1 : [HA] > [A-]
        // < 1 : [HA] < [A-]
        if (_vtsfm[i]->IsAcid()) {
          //cout << "IsAcid == " << _vtsfm[i]->IsAcid() << endl;
          //cout << "pKa == " << _vpKa[i] << endl;
          //cout << "pow(10, _vpKa[i] - pH) == " << pow(10, _vpKa[i] - pH) << endl;
          if (pow(10, _vpKa[i] - pH) < 1.0) {
            //cout << "APPLY!!" << endl;
            _vtsfm[i]->Apply(mol);
          }
        }

        // 10^(pKa - pH) = [BH+] / [B:]
        //
        // > 1 : [BH+] > [B:]
        // < 1 : [BH+] < [B:]
        if (_vtsfm[i]->IsBase()) {
          //cout << "IsBase == " << _vtsfm[i]->IsBase() << endl;
          //cout << "pKa == " << _vpKa[i] << endl;
          //cout << "pow(10, _vpKa[i] - pH) == " << pow(10, _vpKa[i] - pH) << endl;
          if (pow(10, _vpKa[i] - pH) > 1.0) {
            //cout << "APPLY!!" << endl;
            _vtsfm[i]->Apply(mol);
          }
        }
      }
    }

    atomtyper.CorrectAromaticNitrogens(mol);
  }

extern "C" char *
ob_delete_hydrogens (char *smiles, int nonpolaronly)
{
  OBMol mol;
  OBConversion conv;
  string tmpStr (smiles);
  string outstring;
  istringstream molstream1 (tmpStr);
  ostringstream molstream2;
  char *tmpMolfile;

  conv.SetInAndOutFormats ("SMI", "SMI");

  conv.Read (&mol, &molstream1);

  if(mol.NumHvyAtoms () > 0) {
    if (nonpolaronly != 0)
      {
	mol.DeleteNonPolarHydrogens ();
      }
    else
      {
	mol.DeleteHydrogens ();
      }
  } else {
    cout << "Warning: Cannot remove hydrogens. Resulting molecule would be empty!" << endl;
  }        

  conv.Write (&mol, &molstream2);

  outstring = molstream2.str ();

  // remove the trailling $$$$\n from the SDFile
  if (outstring.find ("$$$$\n", 0) != string::npos)
    {
      outstring = outstring.substr (0, outstring.length () - 5);
    }
  else if (outstring.find ("$$$$\r\n", 0) != string::npos)
    {
      outstring = outstring.substr (0, outstring.length () - 6);
    }

  tmpMolfile = strdup (outstring.c_str ());

  return (tmpMolfile);
}

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
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        cout << "no molecule loaded." << endl;
        continue;
      }
      cout << endl << "  Van der Waals energy = " << pFF->E_VDW() << " " << pFF->GetUnit() << endl << endl;
      continue;
    }
    
    if (EQn(commandline, "eeq", 3)) {
      if (mol.Empty()) {
        cout << "no molecule loaded." << endl;
        continue;
      }
      cout << endl << "  electrostatic energy = " << pFF->E_Electrostatic() << " " << pFF->GetUnit() << endl << endl;
      continue;
    }
    
    if (EQn(commandline, "addH", 4)) {
      int num1, num2;
      num1 = mol.NumAtoms();
      mol.AddHydrogens(false, true);
      num2 = mol.NumAtoms();
      cout << (num2 - num1) << " hydrogens added." << endl;
      
      if (!pFF->Setup(mol)) {
        cout << "error while initializing the force field for this molecule." <<endl;
        continue;
      }
      continue;
    }
    
    if (EQn(commandline, "delH", 4)) {
      int num1, num2;
      num1 = mol.NumAtoms();
      mol.DeleteHydrogens();
      num2 = mol.NumAtoms();
      cout << (num1 - num2) << " hydrogens deleted." << endl;
      
      if (!pFF->Setup(mol)) {
        cout << "error while initializing the force field for this molecule." <<endl;
        continue;
      }
      continue;
    }
    
    if (EQn(commandline, "gen", 3)) {
      //pFF->GenerateCoordinates();
      pFF->UpdateCoordinates(mol);
      continue;
    }
    
    if (EQn(commandline, "rs", 2)) {
      pFF->SystematicRotorSearch();
      pFF->UpdateCoordinates(mol);
      continue;
    }
    
    if (EQn(commandline, "nconf", 5)) {
      cout << endl << "  number of conformers = " << mol.NumConformers() << endl << endl;
      continue;
    }


    //
    // commands with parameters
    //
    tokenize(vs, commandline);

bool Database::readTreeSmi (string smi, Tid tid, Tid orig_tid, int line_nr) {

    OBMol mol;

    istringstream iss (smi, istringstream::in);
    OBConversion conv(&iss,&cout);

    // read the molecule
    conv.SetInAndOutFormats("SMI","SDF");

    if (!conv.ReadString(&mol,smi)) {
        cerr << "Error during conversion" << endl;
        return(0);
    }   

    if (!mol.DeleteHydrogens()) {
        cerr << "Unable to delete hydrogens" << endl;
        return(0);
    }

    // create and store new tree object
    DatabaseTreePtr tree = new DatabaseTree ( tid , orig_tid , line_nr );


    // SEGFAULT
    trees.push_back(tree);
    trees_map[orig_tid] = tree;

    int nodessize = 0, edgessize = 0;
    static vector<DatabaseTreeNode> nodes;
    static vector<vector<DatabaseTreeEdge> > edges;
    nodes.resize ( 0 );

//    cerr << "Atoms are (Type(ID)):" << endl;
    OBAtomIterator atom;
    mol.BeginAtom(atom);

	///////////
	// NODES //
	///////////

//    cerr << endl;
    do {

        //cerr << " " << (*atom)->GetType() << " (idx " << (*atom)->GetIdx() << ")" << endl;

        InputNodeLabel inputnodelabel=0;

        // set atom type as label
        // code for 'c' is set to -1 (aromatic carbon).
        if (fm::aromatic) {
            (*atom)->IsAromatic() ? inputnodelabel = (*atom)->GetAtomicNum()+150 : inputnodelabel = (*atom)->GetAtomicNum();
        }
        else inputnodelabel = (*atom)->GetAtomicNum();
        nodessize++;

        // Insert into map, using subsequent numbering for internal labels:
    	// node nr. 1, node nr. 2, ...
	    // Store direction inputlabel -> internal label
    	//cerr << "NodeLabelMap: insert " << inputnodelabel << " --> " << nodelabels.size() << endl;
        map_insert_pair ( nodelabelmap ) p = nodelabelmap.insert ( make_pair ( inputnodelabel, (nodelabels.size()) ) );

        // Store direction internal label -> node
        // if node has NOT been present, label it and set frequency to 1
        if ( p.second ) {
          vector_push_back ( DatabaseNodeLabel, nodelabels, nodelabel );
          nodelabel.inputlabel = inputnodelabel;
          nodelabel.occurrences.parent = NULL;
          nodelabel.occurrences.number = 1;
          nodelabel.lasttid = tid;
          //cerr << "Created node label " << nodelabel.inputlabel << " (freq " << nodelabel.frequency <<  ")" << endl;
        }
        // if node has been present, increase frequency
        else {
          DatabaseNodeLabel &nodelabel = nodelabels[p.first->second];
          if ( nodelabel.lasttid != tid ) nodelabel.frequency++;
          nodelabel.lasttid = tid;
          //cerr << "Updated node label " << nodelabel.inputlabel << " (freq " << nodelabel.frequency << ")" << endl;
        }

        // Tree node
        vector_push_back ( DatabaseTreeNode, nodes, node );
        node.nodelabel = p.first->second;                           // refer to nodelabel
    //    node.atom = (*atom);                                        // attach OB atom
        node.incycle = false;
        //cerr << "Created tree node for OB index " << node.atom->GetIdx() << " (nodelabel " << (int) node.nodelabel << ", nodes[] size " << nodessize << ")" << endl;

    } while (mol.NextAtom(atom));


    // copy nodes to tree and prepare edges storage size
    // edges[nodeid] gives the edges going out of node with id 'nodeid'
    tree->nodes.reserve ( nodessize );
    if ( edges.size () < (unsigned int) nodessize )
        edges.resize ( nodessize );					// edges stored per node
    for ( int i = 0; i < nodessize; i++ ) {
        edges[i].resize ( 0 );						// no edges yet
        tree->nodes.push_back ( nodes[i] );
    }
    //cerr << endl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  if (inFileArg == 0 || outFileArg == 0
      || (inFileArg < 0 && !gotInType)
      || (outFileArg < 0 && !gotOutType))
    usage();

  if (!gotInType)
    {
      if (extab.CanReadExtension(argv[inFileArg]))
	inFileType = extab.FilenameToType(argv[inFileArg]);
      else
	{
	  cerr << program_name << ": cannot read input format!" << endl;
	  usage();
	}
    }
  if (!gotOutType)
    {
      if (extab.CanWriteExtension(argv[outFileArg]))
	outFileType = extab.FilenameToType(argv[outFileArg]);
      else
	{
	  cerr << program_name << ": cannot write output format!" << endl;
	  usage();
	}
    }

  // Finally, we can do some work!
  OBMolVector moleculeList;
  ifstream inFileStream;
  bool usingStdin = false;
  bool canRead = true;
  int currentMol = 1;

  // read
  if (inFileArg > 0)
    {
      inFileStream.open(argv[inFileArg]);
      if (!inFileStream)
	{
	  cerr << program_name << ": cannot read input file!" << endl;
	  exit (-1);
	}
    }
  else
    usingStdin = true;

  while (canRead)
    {
      OBMol *mol = new OBMol(inFileType, outFileType);
      if (!usingStdin)
	fileFormat.ReadMolecule(inFileStream, *mol, argv[inFileArg]);
      else
	fileFormat.ReadMolecule(cin, *mol, "STDIN");

      if (mol->NumAtoms() != 0)
	{
	  // Perform any requested transformations
	  if (removeHydrogens)
	    mol->DeleteHydrogens();
	  if (addHydrogens)
	    mol->AddHydrogens(false, usePH);
	  if (centerCoords)
	    mol->Center();
      
	  if (currentMol >= firstMol && currentMol <= lastMol)
	    moleculeList.PushMol(mol);

	  if (!usingStdin && 
	      (inFileStream.peek() == EOF || !inFileStream.good()) )
	    canRead = false;
	  else if (usingStdin && (cin.peek() == EOF || !cin.good()) )
	    canRead = false;
	  else if (currentMol > lastMol)
	    canRead = false;
	}
      else // 0 atoms in this molecule! 
	{
	  //	  cerr << " error: has zero atoms! " << endl;
	  canRead = false;
	}

      currentMol++;
    }
 
  // write
  if (outFileArg > 0)
    {
      ofstream outFileStream(argv[outFileArg]);
      if (!outFileStream)
	{
	  cerr << program_name << ": cannot write to output file!" << endl;
	  exit (-1);
	}
      moleculeList.Write(outFileStream, formatOptions);
    }
  else
    moleculeList.Write(cout, formatOptions);

  return(0);
}