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

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

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

  conv.Read (&mol, &molstream1);

  mol.AddHydrogens (polaronly != 0, correct4PH != 0);

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

// Reading an InChI and then adding hydrogens messed up the structure
void test_Issue134_InChI_addH()
{
  OBConversion conv;
  conv.SetInFormat("inchi");
  OBMol mol;
  conv.ReadString(&mol, "InChI=1S/C2H7NO/c1-2(3)4/h2,4H,3H2,1H3/t2-/m0/s1");
  OB_ASSERT(!mol.HasData(OBGenericDataType::VirtualBondData));
  mol.AddHydrogens();
  conv.SetOutFormat("smi");
  std::string res = conv.WriteString(&mol, true);
  OB_COMPARE(res, "C[[email protected]@H](N)O");
}

  void HydrogensCommand::redo()
  {
    if (m_SelectedList.size() == 0) {
      switch(m_action) {
      case AddHydrogens:
        m_molecule->addHydrogens();
        break;
      case AddHydrogensPH:
        {
          OBMol obmol = m_molecule->OBMol();
          obmol.UnsetFlag(OB_PH_CORRECTED_MOL);
          FOR_ATOMS_OF_MOL (a, obmol)
            a->SetFormalCharge(0.0);
          obmol.SetAutomaticFormalCharge(true);
          obmol.AddHydrogens(false, true, m_pH);
          m_molecule->setOBMol(&obmol);
          break;
        }
      case RemoveHydrogens:
        m_molecule->removeHydrogens();
        break;
      }
    }
    else { // user selected some atoms, only operate on those

      foreach(unsigned long id, m_SelectedList.subList(Primitive::AtomType))
      {
        Atom *atom = m_molecule->atomById(id);
        if(atom)
        {
          switch(m_action) {
            case AddHydrogens:
              m_molecule->addHydrogens(atom);
              break;
            case RemoveHydrogens:
              m_molecule->removeHydrogens(atom);
              break;
            default:
              break;
          }
        }
      }
    } // end adding to selected atoms
    m_molecule->update();
  }

  void OpConfab::Run(OBConversion* pConv, OBMol* pmol)
  {
    OBMol mol = *pmol;
    
    N++;
    cout << "**Molecule " << N << endl << "..title = " << mol.GetTitle() << endl;
    cout << "..number of rotatable bonds = " << mol.NumRotors() << endl;
    mol.AddHydrogens();
    bool success = pff->Setup(mol);
    if (!success) {
      cout << "!!Cannot set up forcefield for this molecule\n"
           << "!!Skipping\n" << endl;
      return;
    }
    pff->DiverseConfGen(rmsd_cutoff, conf_cutoff, energy_cutoff, verbose);

    pff->GetConformers(mol);
    int nconfs = include_original ? mol.NumConformers() : mol.NumConformers() - 1;
    unsigned int c = include_original ? 0 : 1;

    // If mol.NumRotors is 0 and originals have not been included, then nconfs
    // may be 0. Here, if nconfs is 0, we include the original input conformer
    if (nconfs == 0) {
      nconfs = mol.NumConformers();
      c = 0;
    }

    cout << "..generated " << nconfs << " conformers" << endl;

    for (; c < mol.NumConformers(); ++c) {
      mol.SetConformer(c);
      if(!pConv->GetOutFormat()->WriteMolecule(&mol, pConv))
        break;
    }
    cout << endl;

  }

// This function will call the Babel library to add 
// hydrogens to the residues
void PDB::addHydrogensToPair(AminoAcid& a, AminoAcid& b, int cd1, int cd2)
{
  OBMol mol;
  string addedH;
  istringstream tempss;
  bool ligand;
  if(b.atom[0]->line.find("HETATM") != string::npos)
    {
      ligand = true;
    }
  else
    {
      ligand = false;
    }

  // This section is just to suppress all of the 
  // warning message that aren't important to us
  {
    OBConversion apiConv;
    OBFormat* pAPI = OBConversion::FindFormat("obapi");
    if(pAPI)
      {
        apiConv.SetOutFormat(pAPI);
        apiConv.AddOption("errorlevel", OBConversion::GENOPTIONS, "0");
        apiConv.Write(NULL, &std::cout);
      }
  }

  // Now, let's pack up the information into a string
  string packedFile="";
  for(unsigned int i=0; i < a.altlocs[cd1].size(); i++)
    {
      if( !a.altlocs[cd1][i]->skip )
        {
          packedFile += a.altlocs[cd1][i]->line + "\n";
        }
    }
  
  int cd2_al = cd2;
  if(b.residue == "ASP" || b.residue == "GLU")
    {
      cd2_al = cd2%(b.altlocs.size());
    }

  for(unsigned int i=0; i < b.altlocs[cd2_al].size(); i++)
    {
      if( !b.altlocs[cd2_al][i]->skip )
        {
            packedFile += b.altlocs[cd2_al][i]->line + "\n";
        }
    }
  packedFile += a.makeConect(cd1);
  packedFile += b.makeConect(cd2_al);


  // Now, let's set up some Babel information
  // First, we get the PDB format to tell
  // Babel how to read the information and 
  // how to output it
  OBFormat* pdbformat = this->conv.FindFormat("pdb");
  this->conv.SetInFormat(pdbformat);
  this->conv.SetOutFormat(pdbformat);

  // Here is where Babel reads everything
  // and adds hydrogens to the pair
  // TO ADD: option to set pH
  this->conv.ReadString(&mol,packedFile);
  mol.AddHydrogens(false,true,PH_LEVEL);

  // Let's write the newly written hydrogens to 
  // a string and parse it
  addedH = this->conv.WriteString(&mol);
  tempss.str(addedH);

  // This ensures that the ligand hydrogens are labeled as
  // HETATM instead of ATOM just for the sake of STAAR. 
  // This may be wrong, but it should be fine since we are
  // stripping out that information later when we write 
  // the GAMESS inp files
  if( ligand )
    {
      string line;
      string f = "";
      while( getline(tempss,line) )
        {
          if( line.find(b.residue) != string::npos )
            {
              line.replace(0,6,"HETATM");
            }
          f += line + "\n";
        }
      tempss.seekg(ios_base::beg);
      tempss.clear();
      tempss.str(f);
    }

  this->failure = false;
  this->parsePDB(tempss,99999.99);
//.........这里部分代码省略.........

int main(int argc,char **argv)
{
  char *program_name= argv[0];
  int c;
  int verbose = 0;
  bool hydrogens = false;
  string basename, filename = "", option, option2, ff = "";

  if (argc < 2) {
    cout << "Usage: obenergy [options] <filename>" << endl;
    cout << endl;
    cout << "options:      description:" << endl;
    cout << endl;
    cout << "  -v          verbose: print out indivual energy interactions" << endl;
    cout << endl;
    cout << "  -h          add hydrogens before calculating energy" << endl;
    cout << endl;
    cout << "  -ff ffid    select a forcefield" << endl;
    cout << endl;
    cout << "              available forcefields:" << endl;
    cout << endl;
    OBPlugin::List("forcefields", "verbose");
    exit(-1);
  } else {
    int ifile = 1;
    for (int i = 1; i < argc; i++) {
      option = argv[i];
      
      if (option == "-v") {
        verbose = 1;
        ifile++;
        break;
      }

      if (option == "-h") {
        hydrogens = true;
        ifile++;
      }

      if ((option == "-ff") && (argc > (i+1))) {
        ff = argv[i+1];
        ifile += 2;
      }
    }
    
    basename = filename = argv[ifile];
    size_t extPos = filename.rfind('.');

    if (extPos!= string::npos) {
      basename = filename.substr(0, extPos);
    }


  }

  // Find Input filetype
  OBConversion conv;
  OBFormat *format_in = conv.FormatFromExt(filename.c_str());
    
  if (!format_in || !conv.SetInFormat(format_in)) {
    cerr << program_name << ": cannot read input format!" << endl;
    exit (-1);
  }

  ifstream ifs;
  ofstream ofs;

  // Read the file
  ifs.open(filename.c_str());
  if (!ifs) {
    cerr << program_name << ": cannot read input file!" << endl;
    exit (-1);
  }

  OBForceField* pFF = OBForceField::FindForceField(ff);
  if (!pFF) {
    cerr << program_name << ": could not find forcefield '" << ff << "'." <<endl;
    exit (-1);
  }
  pFF->SetLogFile(&cout);
  if (verbose)
    pFF->SetLogLevel(OBFF_LOGLVL_HIGH);
  else
    pFF->SetLogLevel(OBFF_LOGLVL_MEDIUM);

  OBMol mol;
  double energy;
  for (c=1;;c++) {
    mol.Clear();
    if (!conv.Read(&mol, &ifs))
      break;
    if (mol.Empty())
      break;

    if (hydrogens)
      mol.AddHydrogens();
       
    if (!pFF->Setup(mol)) {
      cerr << program_name << ": could not setup force field." << endl;
      exit (-1);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  //           1      2       3        4         5     6     7      8       9         10         11      12
  output2 << "#METHOD: " << ff << " THREADS: " << nthreads << " DATASET: " << filename.c_str() << std::endl;
  output2 << "#E_BOND E_ANGLE E_STRBND E_TORSION E_OOP E_VDW E_ELEC N_ATOMS PAIRS_VDW PAIRS_ELEC MEM_VDW MEM_ELEC " << std::endl;

	// A third file is created to store information on the memory allocation of data structures
	// from the MMFF94 calculation routines. breakdown of memory allocated for each calculation type
	char filepath3[1100];
	std::ofstream output3;
	sprintf(filepath3, "%s/%s_%s_t%d_malloc.mat", cwd, statsfile, ff.c_str(), nthreads);
	std::cout << "Writing memory allocation breakdown detail file to: " << filepath3 << std::endl;
	output3.open(filepath3, ios::out | ios::app ); // The file is open in append mode

	//           1     2      3       4        5         6     7     8      9      10      11       12        13    14    15
	output3 << "#METHOD: " << ff << " THREADS: " << nthreads << " DATASET: " << filename.c_str() << std::endl;
	output3 << "#ATOMS M_BOND M_ANGLE M_STRBND M_TORSION M_OOP M_VDW M_ELEC C_BOND C_ANGLE C_STRBND C_TORSION C_OOP C_VDW C_ELEC" << std::endl;

  double bondCalcTime, angleCalcTime, strbndCalcTime, torsionCalcTime, oopCalcTime, vdwCalcTime, electrostaticCalcTime;
  int numPairsVDW, numPairsElectrostatic;

  OBMol mol;
  double energy;
  for (c=1;;c++) {
    mol.Clear();

    totalTimer.start();
    readTimer.start();

    if (!conv.Read(&mol, &ifs))
      break;
    if (mol.Empty())
      break;

    if (hydrogens)
      mol.AddHydrogens();
       
    readTime = readTimer.get();
    setupTimer.start();


    if (!pFF->Setup(mol)) {
      cerr << program_name << ": could not setup force field." << endl;
      exit (-1);
    }

    setupTime = setupTimer.get();
    computeTimer.start();
    
    energy = pFF->Energy(false);

    computeTime = computeTimer.get();
    totalTime = totalTimer.get();

    // THREADS  ENERGY  MOL_MASS  NUM_ATOMS  NUM_ROTORS  NUM_CONF  TOT_TIME  TIME_READ  TIME_SETUP  TIME_COMPUTE STEPS  #MOL_NAME
    output << nthreads << " " << energy << " " << mol.GetExactMass() << " " << mol.NumAtoms()
    		<< " " << mol.NumRotors() << "  " <<  mol.NumConformers() << " "
    		<< totalTime <<  " " << readTime << " " << " " << setupTime << " " << computeTime << " "
    		<< totalSteps << " #" << mol.GetTitle()  // comment added to avoid errors when reading matrix in Octave
    		<< std::endl;

    map<string, double> timings = pFF->getTimings();
    map<string, size_t> memalloc = pFF->getAllocatedMemory();
    MapKeys mk;
    // 1      2       3        4         5     6     7      8       9         10         11      12
    // E_BOND E_ANGLE E_STRBND E_TORSION E_OOP E_VDW E_ELEC N_ATOMS PAIRS_VDW PAIRS_ELEC MEM_VDW MEM_ELEC
    output2 << timings[mk.TIME_BOND_CALCULATIONS] << " "  // 1
    		<< timings[mk.TIME_ANGLE_CALCULATIONS] << " " // 2

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

  if (argc != 1)
    {
      cout << "Usage: conversion" << endl;
      cout << " Unit tests for OBConversion " << endl;
      return(-1);
    }

  cout << "# Unit tests for OBConversion \n";

  // the number of tests for "prove"
  cout << "1..9\n";

  cout << "ok 1\n"; // for loading tests

  OBMol obMol;
  OBConversion obConversion;
  obConversion.SetInAndOutFormats("smi", "mdl");
  cout << "ok 2\n";

  obConversion.ReadString(&obMol, "C1=CC=CS1");
  cout << "ok 3\n";

  if (obMol.NumAtoms() == 5) {
    cout << "ok 4\n";
  } else {
    cout << "not ok 4\n";
  }

  obMol.AddHydrogens();
  if (obMol.NumAtoms() == 9) {
    cout << "ok 5\n";
  } else {
    cout << "not ok 5\n";
  }

  if ( (obConversion.WriteString(&obMol)).length() > 0)
    cout << "ok 6\n";
  else
    cout << "not ok 6\n";

  // PR#1474265
  obConversion.WriteFile(&obMol, "test.mdl");
  ifstream ifs("test.mdl");
  if (ifs.good())
    cout << "ok 7\n";
  else
    cout << "not ok 7\n";

  // PR#143577
  obConversion.SetInFormat("mdl");
  obConversion.ReadFile(&obMol, "test.mdl");
  if ( remove("test.mdl") != -1)
    cout << "ok 8\n";
  else
    cout << "not ok 8\n";
  
  // gzip input
  // gzip output

  // multi-molecule reading
  // PR#1465586
  // aromatics.smi
  // attype.00.smi

  //ReadFile()
  //Read()
  //WriteString()
  // GetOutputIndex()
  // IsLast

  //ReadString()
  //IsFirstInput
  //Read()

  // splitting
  
  // splitting using gzip-input
  // PR#1357705
  
  // size 0 input
  // PR#1250900
  
  // RegisterFormat
  // FindFormat
  // FormatFromExt
  // FormatFromMIME
  // GetNextFormat
  // GetDefaultFormat

  // BatchFileName
  // IncrementedFileName

  // option handling
  // AddOption
  // IsOption
//.........这里部分代码省略.........

int main(int argc,char **argv)
{
  char *program_name= argv[0];
  int c;
  char *FileIn = NULL;

  if (argc != 2)
    {
      string err = "Usage: ";
      err += program_name;
      err += " <filename>\n"
      "Output format:\n"
        "name NAME\n"
        "formula  FORMULA\n"
        "mol_weight MOLECULAR_WEIGHT\n"
        "exact_mass ISOTOPIC MASS\n"
        "canonical_SMILES STRING\n"
        "InChI  STRING\n"
        "num_atoms  NUM\n"
        "num_bonds  NUM\n"
        "num_residues  NUM\n"
	"num_rotors NUM\n"
        "sequence RESIDUE_SEQUENCE\n"
        "num_rings NUMBER_OF_RING_(SSSR)\n"
        "logP   NUM\n"
        "PSA    POLAR_SURFACE_AREA\n"
        "MR     MOLAR REFRACTIVITY";
      err += "$$$$";
//      ThrowError(err); wasn't being output because error level too low
      cerr << err; //Why not do directly
      exit(-1);
    }
  else
    {
      FileIn  = argv[1];
    }

  // Find Input filetype
  OBConversion conv;
  OBFormat *format = conv.FormatFromExt(FileIn);
    
  if (!format || !conv.SetInFormat(format))
    {
      cerr << program_name << ": cannot read input format!" << endl;
      exit (-1);
    }

  ifstream ifs;

  // Read the file
  ifs.open(FileIn);
  if (!ifs)
    {
      cerr << program_name << ": cannot read input file!" << endl;
      exit (-1);
    }
  
  OBMol mol;
  OBFormat *canSMIFormat = conv.FindFormat("can");
  OBFormat *inchiFormat = conv.FindFormat("inchi");


  ////////////////////////////////////////////////////////////////////////////
  // List of properties
  // Name
  // Molecular weight (Standard molar mass given by IUPAC atomic masses)
  // Number of rings : the size of the smallest set of smallest rings (SSSR)
  
  //.....ADD YOURS HERE.....
  
  for (c = 1;; ++c)
    {
      mol.Clear();
      conv.Read(&mol, &ifs);
      if (mol.Empty())
        break;
      
      if (!mol.HasHydrogensAdded())
        mol.AddHydrogens();
      // Print the properties
      if (strlen(mol.GetTitle()) != 0)
        cout << "name             " << mol.GetTitle() << endl;
      else 
        cout << "name             " << FileIn << " " << c << endl;

      cout << "formula          " << mol.GetFormula() << endl;
      cout << "mol_weight       " << mol.GetMolWt() << endl;
      cout << "exact_mass       " << mol.GetExactMass() << endl;

      string smilesString = "-";
      if (canSMIFormat) {
        conv.SetOutFormat(canSMIFormat);
        smilesString = conv.WriteString(&mol);
        if ( smilesString.length() == 0 )
        {
          smilesString = "-";
        }
      }
      cout << "canonical_SMILES " << smilesString << endl;

//.........这里部分代码省略.........

///////////////////////////////////////////////////////////////////////////////
//! \brief  Generate rough 3D coordinates for SMILES (or other 0D files).
//
int main(int argc,char **argv)
{
  char *program_name= argv[0];
  int c;
  string basename, filename = "", option, option2, ff = "MMFF94";

  list<string> argl(argv+1, argv+argc);

  list<string>::iterator optff = find(argl.begin(), argl.end(), "-ff");
  if (optff != argl.end()) {
    list<string>::iterator optffarg = optff;
    ++optffarg;

    if (optffarg != argl.end()) {
      ff = *optffarg;

      argl.erase(optff,++optffarg);
    } else {
      argl.erase(optff);
    }
  }

  if (argl.empty()) {
    cout << "Usage: obgen <filename> [options]" << endl;
    cout << endl;
    cout << "options:      description:" << endl;
    cout << endl;
    cout << "  -ff         select a forcefield" << endl;
    cout << endl;
    OBPlugin::List("forcefields", "verbose");
    exit(-1);
  }

  basename = filename = *argl.begin();
  size_t extPos = filename.rfind('.');

  if (extPos!= string::npos) {
    basename = filename.substr(0, extPos);
  }

  // Find Input filetype
  OBConversion conv;
  OBFormat *format_in = conv.FormatFromExt(filename.c_str());
  OBFormat *format_out = conv.FindFormat("sdf");

  if (!format_in || !format_out || !conv.SetInAndOutFormats(format_in, format_out)) {
    cerr << program_name << ": cannot read input/output format!" << endl;
    exit (-1);
  }

  ifstream ifs;
  ofstream ofs;

  // Read the file
  ifs.open(filename.c_str());
  if (!ifs) {
    cerr << program_name << ": cannot read input file!" << endl;
    exit (-1);
  }

  OBMol mol;

  for (c=1;;c++) {
      mol.Clear();
      if (!conv.Read(&mol, &ifs))
        break;
      if (mol.Empty())
        break;

      OBForceField* pFF = OBForceField::FindForceField(ff);
      if (!pFF) {
        cerr << program_name << ": could not find forcefield '" << ff << "'." <<endl;
        exit (-1);
      }

      //mol.AddHydrogens(false, true); // hydrogens must be added before Setup(mol) is called

      pFF->SetLogFile(&cerr);
      pFF->SetLogLevel(OBFF_LOGLVL_LOW);

      //pFF->GenerateCoordinates();
      OBBuilder builder;
      builder.Build(mol);

      mol.AddHydrogens(false, true); // hydrogens must be added before Setup(mol) is called
      if (!pFF->Setup(mol)) {
        cerr << program_name << ": could not setup force field." << endl;
        exit (-1);
      }

      pFF->SteepestDescent(500, 1.0e-4);
      pFF->WeightedRotorSearch(250, 50);
      pFF->SteepestDescent(500, 1.0e-6);

      pFF->UpdateCoordinates(mol);
      //pFF->ValidateGradients();
      //pFF->SetLogLevel(OBFF_LOGLVL_HIGH);
//.........这里部分代码省略.........

  bool OpConformer::Do(OBBase* pOb, const char* OptionText, OpMap* pmap, OBConversion*)
  {
    OBMol* pmol = dynamic_cast<OBMol*>(pOb);
    if(!pmol)
      return false;
    pmol->AddHydrogens(false, false);

    OpMap::const_iterator iter;
    bool log = false;
    bool systematic = false;
    bool random = false;
    bool weighted = false;
    int numConformers = 30;

    iter = pmap->find("log");
    if(iter!=pmap->end())
      log=true;

    iter = pmap->find("nconf");
    if(iter!=pmap->end())
      getInteger(iter->second, numConformers);

    iter = pmap->find("systematic");
    if(iter!=pmap->end())
      systematic = true;

    iter = pmap->find("random");
    if(iter!=pmap->end())
      random = true;

    iter = pmap->find("weighted");
    if(iter!=pmap->end())
      weighted = true;

    if (systematic || random || weighted) {
      std::string ff = "MMFF94";
      iter = pmap->find("ff");
      if(iter!=pmap->end())
        ff = iter->second;
      OBForceField* pFF = OBForceField::FindForceField(ff);

      // set some force field variables
      pFF->SetLogFile(&clog);
      pFF->SetLogLevel(log ? OBFF_LOGLVL_MEDIUM : OBFF_LOGLVL_NONE);

      if (!pFF->Setup(*pmol)) {
        cerr  << "Could not setup force field." << endl;
        return false;
      }
    } else {
      int numChildren = 5;
      int mutability = 5;
      int convergence = 25;
      std::string score = "rmsd";

      iter = pmap->find("children");
      if(iter!=pmap->end())
        getInteger(iter->second, numChildren);

      iter = pmap->find("mutability");
      if(iter!=pmap->end())
        getInteger(iter->second, mutability);

      iter = pmap->find("convergence");
      if(iter!=pmap->end())
        getInteger(iter->second, convergence);

      iter = pmap->find("score");
      if(iter!=pmap->end())
        score = iter->second;

      OBConformerSearch cs;
      if (score == "energy")
        cs.SetScore(new OBEnergyConformerScore);

      if (cs.Setup(*pmol, numConformers, numChildren, mutability, convergence)) {
        cs.Search();
        cs.GetConformers(*pmol);
      }
    }

    return true;
  }

  bool OpConformer::Do(OBBase* pOb, const char* OptionText, OpMap* pmap, OBConversion*)
  {
    OBMol* pmol = dynamic_cast<OBMol*>(pOb);
    if(!pmol)
      return false;
    pmol->AddHydrogens(false, false);

    OpMap::const_iterator iter;
    bool log = false;
    bool systematic = false;
    bool random = false;
    bool weighted = false;
    bool fast = false;
    int numConformers = 30;

    iter = pmap->find("log");
    if(iter!=pmap->end())
      log=true;

    iter = pmap->find("nconf");
    if(iter!=pmap->end())
      getInteger(iter->second, numConformers);

    iter = pmap->find("systematic");
    if(iter!=pmap->end())
      systematic = true;

    iter = pmap->find("random");
    if(iter!=pmap->end())
      random = true;

    iter = pmap->find("fast");
    if(iter!=pmap->end())
      fast = true;

    iter = pmap->find("weighted");
    if(iter!=pmap->end())
      weighted = true;

    if (systematic || random || fast || weighted) {
      std::string ff = "MMFF94";
      iter = pmap->find("ff");
      if(iter!=pmap->end())
        ff = iter->second;
      OBForceField* pFF = OBForceField::FindForceField(ff);

      // set some force field variables
      pFF->SetLogFile(&clog);
      pFF->SetLogLevel(log ? OBFF_LOGLVL_MEDIUM : OBFF_LOGLVL_NONE);

      // Add cut-offs for faster conformer searching
      // Generally people will perform further optimization on a final conformer
      pFF->EnableCutOff(true);
      pFF->SetVDWCutOff(10.0);
      pFF->SetElectrostaticCutOff(20.0);
      pFF->SetUpdateFrequency(10); // delay updates of non-bonded distances

      if (!pFF->Setup(*pmol)) {
        cerr  << "Could not setup force field." << endl;
        return false;
      }

      // Perform search
      if (systematic) {
        pFF->SystematicRotorSearch(10); // 10 steepest-descent forcfield steps per conformer
      } else if (fast) {
        pFF->FastRotorSearch(true); // permute rotors
      } else if (random) {
        pFF->RandomRotorSearch(numConformers, 10);
      } else if (weighted) {
        pFF->WeightedRotorSearch(numConformers, 10);
      }
      pFF->GetConformers(*pmol);

    } else { // GA-based searching
      int numChildren = 5;
      int mutability = 5;
      int convergence = 5;
      std::string score = "rmsd";

      iter = pmap->find("children");
      if(iter!=pmap->end())
        getInteger(iter->second, numChildren);

      iter = pmap->find("mutability");
      if(iter!=pmap->end())
        getInteger(iter->second, mutability);

      iter = pmap->find("convergence");
      if(iter!=pmap->end())
        getInteger(iter->second, convergence);

      iter = pmap->find("score");
      if(iter!=pmap->end())
        score = iter->second;

      OBConformerSearch cs;
      if (score == "energy")
        cs.SetScore(new OBEnergyConformerScore);
      else if (score == "mine" || score == "minenergy")
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    {
      cout << "Bail out! SMILES format is not loaded" << endl;
      return -1;
    }

  for (;mifs;)
    {
      mol.Clear();
      conv.Read(&mol);
      if (mol.Empty())
        continue;
      if (!rifs.getline(buffer,BUFF_SIZE))
        {
          cout << "Bail out! error reading reference data" << endl;
          return -1; // test failed
        }

      tokenize(vs,buffer);
      if (vs.size() != 3)
        {
          cout << "Bail out! Reference data has incorrect format" << endl;
          return -1; // test failed
        }

      if (vs[0] != mol.GetFormula())
        {
          cout << "not ok " << ++currentTest << " # molecular formula incorrect"
               << " for molecule " << mol.GetTitle() << "\n";
        }
      else
        cout << "ok " << ++currentTest << " # molecular formula\n";

      if ( fabs(atof(vs[1].c_str()) - mol.GetMolWt() ) > 1.0e-3)
        {
          cout << "not ok " << ++currentTest << " # molecular weight incorrect"
               << " for molecule " << mol.GetTitle() << "\n";
          cout << "# Expected " << atof(vs[1].c_str()) << " found " <<
            mol.GetMolWt() << "\n";
        }
      else
        cout << "ok " << ++currentTest << " # molecular weight\n";

      if ( fabs(atof(vs[2].c_str()) - mol.GetExactMass() ) > 1.0e-3)
        {
          cout << "not ok " << ++currentTest << " # exact mass incorrect"
               << " for molecule " << mol.GetTitle() << "\n";
          cout << "# Expected " << atof(vs[2].c_str()) << " found " <<
            mol.GetExactMass() << "\n";
        }
      else
        cout << "ok " << ++currentTest << " # molecular exact mass\n";


      // now after adding explict hydrogens -- should be identical
      //  since we'll add hydrogens that were implicit before

      // PR#1485580
      mol.AddHydrogens();

      if (vs[0] != mol.GetFormula())
        {
          cout << "not ok " << ++currentTest << " # molecular formula incorrect"
               << " for hydrogen-added molecule " << mol.GetTitle() << "\n";
        }
      else
        cout << "ok " << ++currentTest << " # molecular hydrogen-added formula\n";

      if ( fabs(atof(vs[1].c_str()) - mol.GetMolWt() ) > 1.0e-3)
        {
          cout << "not ok " << ++currentTest << " # molecular weight incorrect"
               << " for hydrogen-added molecule " << mol.GetTitle() << "\n";
          cout << "# Expected " << atof(vs[1].c_str()) << " found " <<
            mol.GetMolWt() << "\n";
          cout << "# Difference " << fabs(atof(vs[1].c_str()) - mol.GetMolWt())
               << "\n";
        }
      else
        cout << "ok " << ++currentTest << " # molecule + hydrogens weight\n";

      if ( fabs(atof(vs[2].c_str()) - mol.GetExactMass() ) > 1.0e-3)
        {
          cout << "not ok " << ++currentTest << " # exact mass incorrect"
               << " for hydrogen-added molecule " << mol.GetTitle() << "\n";
          cout << "# Expected " << atof(vs[2].c_str()) << " found " <<
            mol.GetExactMass() << "\n";
          cout << "# Difference " << fabs(atof(vs[2].c_str()) - mol.GetExactMass())
               << "\n";	
        }
      else
        cout << "ok " << ++currentTest << " # molecular exact mass"
             << " after hydrogen addition\n";

    }

  // return number of tests run
  cout << "1.." << currentTest << endl;

  // Passed tests
  return 0;
}

void formula_test()
{
#ifdef TESTDATADIR
  string testdatadir = TESTDATADIR;
  string results_file = testdatadir + "formularesults.txt";
  string smilestypes_file = testdatadir + "attype.00.smi";
#else
  string results_file = "files/formularesults.txt";
  string smilestypes_file = "files/attype.00.smi";
#endif

  cout << "# Testing molecular formulas..." << endl;

  std::ifstream mifs;
  BOOST_REQUIRE_MESSAGE( SafeOpen(mifs, smilestypes_file.c_str()), "Bail out! Cannot read file " );

  std::ifstream rifs;
  BOOST_REQUIRE_MESSAGE( SafeOpen(rifs, results_file.c_str()), "Bail out! Cannot read file " );

  char buffer[BUFF_SIZE];
  char message[BUFF_SIZE];
  vector<string> vs;
  OBMol mol;
  OBConversion conv(&mifs, &cout);
  unsigned int currentTest = 0;
  // double mass;

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

      tokenize(vs,buffer);
      BOOST_REQUIRE_MESSAGE( vs.size() == 3, "Bail out! Reference data has incorrect format" );

      currentTest++;
      snprintf(message, BUFF_SIZE, "not ok %d # molecular formula incorrect for molecule %s"
		                   " # Expected %s, found %s", currentTest, mol.GetTitle(),
				   vs[0].c_str(), mol.GetFormula().c_str());
      BOOST_CHECK_MESSAGE(vs[0] == mol.GetFormula(), message );

      currentTest++;
      snprintf(message, BUFF_SIZE, "not ok %d # molecular weight incorrect for molecule %s"
		                   " # Expected %f, found %f ", currentTest, mol.GetTitle(),
				   atof(vs[1].c_str()), mol.GetMolWt());
      BOOST_CHECK_MESSAGE( fabs(atof(vs[1].c_str()) - mol.GetMolWt() ) < 1.0e-3, message );

      currentTest++;
      snprintf(message, BUFF_SIZE, "not ok %d # exact mass incorrect for molecule %s"
		                   " # Expected %f, found %f ", currentTest, mol.GetTitle(),
				   atof(vs[2].c_str()), mol.GetExactMass());
      BOOST_CHECK_MESSAGE( fabs(atof(vs[2].c_str()) - mol.GetExactMass() ) < 1.0e-3, message );

      // now after adding explict hydrogens -- should be identical
      //  since we'll add hydrogens that were implicit before

      // PR#1485580
      BOOST_CHECK( mol.AddHydrogens() );

      currentTest++;
      snprintf(message, BUFF_SIZE, "not ok %d # molecular formula incorrect for "
		                   "hydrogen-added molecule %s", currentTest, mol.GetTitle());
      BOOST_CHECK_MESSAGE(vs[0] == mol.GetFormula(), message);

      currentTest++;
      snprintf(message, BUFF_SIZE, "not ok %d # molecular weight incorrect for hydrogen-added "
		                   "molecule %s # Expected %f, found %f ", currentTest, mol.GetTitle(),
				   atof(vs[1].c_str()), mol.GetMolWt());
      BOOST_CHECK_MESSAGE( fabs(atof(vs[1].c_str()) - mol.GetMolWt() ) < 1.0e-3, message);

      currentTest++;
      snprintf(message, BUFF_SIZE, "not ok %d # exact mass  incorrect for hydrogen-added "
		                   "molecule %s # Expected %f, found %f ", currentTest, mol.GetTitle(),
				   atof(vs[2].c_str()), mol.GetExactMass());
      BOOST_CHECK_MESSAGE( fabs(atof(vs[2].c_str()) - mol.GetExactMass() ) < 1.0e-3, message);

    }

}

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