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

extern "C" char *
ob_mol_to_canonical_smiles (char *molfile, int omit_iso_and_chiral_markings)
{
  OBMol mol;
  OBConversion conv;
  string tmpStr (molfile);
  string outstring;
  istringstream molstream (tmpStr);
  ostringstream smilesstream;
  char *tmpSmiles;

  conv.SetInAndOutFormats ("MDL", "CAN");
  conv.AddOption ("n", OBConversion::OUTOPTIONS);
  //conv.AddOption ("c", OBConversion::OUTOPTIONS);
  if (omit_iso_and_chiral_markings != 0)
    {
      conv.AddOption ("i", OBConversion::OUTOPTIONS);
    }

  conv.Read (&mol, &molstream);

  if (mol.Empty ())
    return NULL;

  conv.Write (&mol, &smilesstream);

  outstring = smilesstream.str ();

  outstring = outstring.substr (0, outstring.length () - 1);

  tmpSmiles = strdup (outstring.c_str ());

  return (tmpSmiles);
}

void GenerateFormulaReference()
{
  std::ifstream ifs;
  if (!SafeOpen(ifs, smilestypes_file.c_str()))
    return;

  std::ofstream ofs;
  if (!SafeOpen(ofs, results_file.c_str()))
    return;

  OBMol mol;
  OBConversion conv(&ifs, &cout);

  if(! conv.SetInAndOutFormats("SMI","SMI"))
    {
      cerr << "SMILES format is not loaded" << endl;
      return;
    }

  for (;ifs;)
    {
      mol.Clear();
      conv.Read(&mol);
      if (mol.Empty())
        continue;

      //write out formula, molecular weight and exact mass
      ofs << mol.GetFormula() << " " << mol.GetMolWt() << " " 
          << mol.GetExactMass() << endl;
    }

	cerr << " Molecular formula results written successfully" << endl;
  return;
}

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

  if (argc != 2) {
    cerr << " Usage: " << program_name << " <input file>\n";
    exit(-1);
  }
  else {
      FileIn  = argv[1];
  }

  // Find Input filetype
  OBConversion conv(&cin, &cout);
  OBFormat *inFormat = conv.FormatFromExt(FileIn);

  if (!inFormat || !conv.SetInFormat(inFormat)) {
    cerr << program_name << ": cannot read input format!" << endl;
    exit (-1);
  }
  // If we can't also use this for an output format, use XYZ
  if (!conv.SetOutFormat(inFormat))
    conv.SetOutFormat(conv.FindFormat("xyz"));

  ifstream ifs;

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

  OBMol mol;
  OBPointGroup pg;

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

      // not needed by OBPointGroup, but useful for external programs
      pg.Setup(&mol);
      cerr << "Point Group: " << pg.IdentifyPointGroup() << endl;
      pg.Symmetrize(&mol);

      conv.Write(&mol, &cout);

    } // end for loop

  return(1);
}

void GenerateCharges()
{
  std::ifstream ifs;
  if (!SafeOpen(ifs, molecules_file.c_str()))
    return;

  std::ofstream rofs;
  if (!SafeOpen(rofs, results_file.c_str()))
    return;

  std::ofstream dofs;
  if (!SafeOpen(dofs, dipole_file.c_str()))
    return;

  OBMol mol;
  OBConversion conv(&ifs, &cout);
  char buffer[BUFF_SIZE];
  
  if(! conv.SetInAndOutFormats("SDF","SDF"))
    {
      cerr << "SDF format is not loaded" << endl;
      return;
    }

  OBChargeModel *pCM = OBChargeModel::FindType("mmff94");

  if (pCM == NULL) {
    cerr << "Cannot load charge model!" << endl;
    return;
  }

  std::vector<double> partialCharges;
  vector3 dipoleMoment;
  for (;ifs;)
    {
      mol.Clear();
      conv.Read(&mol);
      if (mol.Empty())
        continue;

      if (pCM->ComputeCharges(mol)) {
        partialCharges = pCM->GetPartialCharges();
      }

      // write out the dipole moment
      dipoleMoment = pCM->GetDipoleMoment(mol);
      sprintf(buffer, "%15.5f%15.5f%15.5f\n", dipoleMoment.x(), dipoleMoment.y(), dipoleMoment.z());
      dofs << buffer;
      
      // and write all the partial charges
      FOR_ATOMS_OF_MOL(atom, mol) {
        sprintf(buffer, "%15.5f\n", atom->GetPartialCharge());
        rofs << buffer;
      }
    }

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

  if (argc != 2) {
    cerr << " Usage: " << program_name << " <input file>\n";
    exit(-1);
  }
  else {
      FileIn  = argv[1];
  }

  // Find Input filetype
  OBConversion conv;
  OBFormat *inFormat = conv.FormatFromExt(FileIn);
    
  if (!inFormat || !conv.SetInFormat(inFormat)) {
    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;
  OBPointGroup pg;

  for (c = 1;; ++c)
    {
      mol.Clear();
      conv.Read(&mol, &ifs);
      if (mol.Empty())
        break;
      
      // not needed by OBPointGroup, but useful for external programs
      mol.Center();
      mol.ToInertialFrame();

      pg.Setup(&mol);
      cout << "Point Group: " << pg.IdentifyPointGroup() << endl;

    } // end for loop
  
  return(1);
}

void GenerateEnergies()
{
  std::ifstream ifs;
  if (!SafeOpen(ifs, molecules_file.c_str()))
    return;

  std::ofstream ofs;
  if (!SafeOpen(ofs, results_file.c_str()))
    return;

  OBMol mol;
  OBConversion conv(&ifs, &cout);
  char buffer[BUFF_SIZE];
  
  if(! conv.SetInAndOutFormats("SDF","SDF"))
    {
      cerr << "SDF format is not loaded" << endl;
      return;
    }

  OBForceField* pFF = OBForceField::FindForceField("Ghemical");

  if (pFF == NULL) {
    cerr << "Cannot load force field!" << endl;
    return;
  }

  pFF->SetLogFile(&cout);
  pFF->SetLogLevel(OBFF_LOGLVL_NONE);

  for (;ifs;)
    {
      mol.Clear();
      conv.Read(&mol);
      if (mol.Empty())
        continue;

      if (!pFF->Setup(mol)) {
        cerr << "Could not setup force field on molecule: " << mol.GetTitle() << endl;
        return;
      }
      
      // Don't compute gradients
      sprintf(buffer, "%15.5f\n", pFF->Energy(false));
      ofs << buffer;
    }

	cerr << " Ghemical force field energies written successfully" << endl;
  return;
}

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)
    {
      if (strncmp(argv[1], "-g", 2))
        {
          cout << "Usage: formula" << endl;
          cout << "   Tests Open Babel molecular formula, weight, and exact mass." << endl;
          return 0;
        }
      else
        {
          GenerateFormulaReference();
          return 0;
        }
    }

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

  std::ifstream mifs;
  if (!SafeOpen(mifs, smilestypes_file.c_str()))
    {
      cout << "Bail out! Cannot read file " << smilestypes_file << endl;
      return -1; // test failed
    }

  std::ifstream rifs;
  if (!SafeOpen(rifs, results_file.c_str()))
    {
      cout << "Bail out! Cannot read file " << results_file << endl;
      return -1; // test failed
    }

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

  if(! conv.SetInAndOutFormats("SMI","SMI"))
    {
      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

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

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

//.........这里部分代码省略.........
	sprintf(filepath2, "%s/%s_%s_t%d_compute.mat", cwd, statsfile, ff.c_str(), nthreads);
	std::cout << "Writing method breakdown detail file to: " << filepath2 << std::endl;
	output2.open(filepath2, ios::out | ios::app ); // The file is open in append mode

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

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

      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
  OBFormat *format_in = conv.FormatFromExt(filename.c_str());

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

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

  // set some force field variables
  pFF->SetLogFile(&cerr);
  pFF->SetLogLevel(OBFF_LOGLVL_LOW);
  pFF->SetVDWCutOff(rvdw);
  pFF->SetElectrostaticCutOff(rele);
  pFF->SetUpdateFrequency(freq);
  pFF->EnableCutOff(cut);
  if (newton)
    pFF->SetLineSearchType(LineSearchType::Newton2Num);

  OBMol mol;

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

    bool done = true;
    OBStopwatch timer;
    timer.Start();
    if (sd) {
      pFF->SteepestDescentInitialize(steps, crit);
    } else {
      pFF->ConjugateGradientsInitialize(steps, crit);
    }

    unsigned int totalSteps = 1;
    while (done) {
      if (sd)
        done = pFF->SteepestDescentTakeNSteps(1);
      else
        done = pFF->ConjugateGradientsTakeNSteps(1);
      totalSteps++;

      if (pFF->DetectExplosion()) {
        cerr << "explosion has occured!" << endl;
        conv.Write(&mol, &cout);
        return(1);
      } else
        pFF->GetCoordinates(mol);
    }
    double timeElapsed = timer.Elapsed();

    pFF->GetCoordinates(mol);

    conv.Write(&mol, &cout);
    cerr << "Time: " << timeElapsed << "seconds. Iterations per second: " <<  double(totalSteps) / timeElapsed << endl;
  } // end for loop

  return(0);
}

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

  // Define location of file formats for testing
#ifdef FORMATDIR
    char env[BUFF_SIZE];
    snprintf(env, BUFF_SIZE, "BABEL_LIBDIR=%s", FORMATDIR);
    putenv(env);
#endif

  if (argc != 1)
    {
      if (strncmp(argv[1], "-g", 2))
        {
          cout << "Usage: charge-mmff94" << endl;
          return 0;
        }
      else
        {
          GenerateCharges();
          return 0;
        }
    }

  cout << "# Testing MMFF94 Charge Model..." << endl;

  std::ifstream mifs;
  if (!SafeOpen(mifs, molecules_file.c_str()))
    {
      cout << "Bail out! Cannot read file " << molecules_file << endl;
      return -1; // test failed
    }

  std::ifstream rifs;
  if (!SafeOpen(rifs, results_file.c_str()))
    {
      cout << "Bail out! Cannot read file " << results_file << endl;
      return -1; // test failed
    }

  std::ifstream difs;
  if (!SafeOpen(difs, dipole_file.c_str()))
    {
      cout << "Bail out! Cannot read file " << dipole_file << endl;
      return -1; // test failed
    }

  char buffer[BUFF_SIZE];
  vector<string> vs;
  OBMol mol;
  OBConversion conv(&mifs, &cout);
  unsigned int currentTest = 0;
  vector3 dipoleMoment, result;
  
  std::vector<double> partialCharges;

  if(! conv.SetInAndOutFormats("SDF","SDF"))
    {
      cout << "Bail out! SDF format is not loaded" << endl;
      return -1; // test failed
    }
    
  OBChargeModel *pCM = OBChargeModel::FindType("mmff94");

  if (pCM == NULL) {
    cerr << "Bail out! Cannot load charge model!" << endl;
    return -1; // test failed
  }

  while(mifs)
    {
      mol.Clear();
      conv.Read(&mol);
      if (mol.Empty())
        continue;
      if (!difs.getline(buffer,BUFF_SIZE))
        {
          cout << "Bail out! error reading reference data" << endl;
          return -1; // test failed
        }
        
      if (!pCM->ComputeCharges(mol)) {
        cout << "Bail out! could not compute charges on " << mol.GetTitle() << endl;
        return -1; // test failed
      }
      partialCharges = pCM->GetPartialCharges();

      // compare the calculated energy to our reference data
      tokenize(vs, buffer);
      if (vs.size() < 3)
        return -1;

      dipoleMoment.SetX(atof(vs[0].c_str()));
      dipoleMoment.SetY(atof(vs[1].c_str()));
      dipoleMoment.SetZ(atof(vs[2].c_str()));
      result = pCM->GetDipoleMoment(mol) - dipoleMoment;
                        
      if ( fabs(result.length_2()) > 1.0e-4)
//.........这里部分代码省略.........

int main(int argc,char **argv)
{
  char c;
  int transX, transY, transZ, centerIdx;
  transX = transY = transZ = centerIdx = -1;

  char *FileIn = NULL;
  char *program_name = argv[0];
  //   char *iext;

  OBConversion conv(&cin,&cout);
  OBFormat *pFormat = conv.FindFormat("smi"); // default format is SMILES

  // Still need to add:
  // rotate X, Y, or Z by set angle
  // translate X, Y, or Z by set amount
    
  // Parse options
  while ((c = getopt(argc, argv, "x:y:z:c:")) != -1)
    {
      switch (c)
        {
        case 'c': /// atom to be centered

          c = sscanf(optarg, "%d", &centerIdx);
          if (c != 1 )
            {
              cerr << program_name << ": unable to parse -c option" << endl;
              exit (-1);
            }
          break;

        case 'x': /// atom to be centered

          c = sscanf(optarg, "%d", &transX);
          if (c != 1 )
            {
              cerr << program_name << ": unable to parse -x option" << endl;
              exit (-1);
            }
          break;

        case 'y': /// atom to y-axis

          c = sscanf(optarg, "%d", &transY);
          if (c != 1 )
            {
              cerr << program_name << ": unable to parse -y option" << endl;
              exit (-1);
            }
          break;

        case 'z': /// atom to z-axis

          c = sscanf(optarg, "%d", &transZ);
          if (c != 1 )
            {
              cerr << program_name << ": unable to parse -z option" << endl;
              exit (-1);
            }
          break;

        }
    }
    
  ifstream ifs;
  FileIn  = argv[optind];
  if (FileIn != NULL)
    {
      // Read the file
      ifs.open(FileIn);
      if (!ifs)
        {
          cerr << program_name << ": cannot read input file!" << endl;
          exit (-1);
        }
      conv.SetInStream(&ifs);
	
      // Find Input filetype
      pFormat = conv.FormatFromExt(FileIn);
      if (pFormat == NULL)
        {
          cerr << program_name << ": cannot read input format!" << endl;
          return (-1);
        }
    }

  if (! conv.SetInAndOutFormats(pFormat, pFormat))
    {
      cerr << program_name << ": cannot read or write to this file format" << endl;
      return (-1);
    }

  OBMol mol;
  mol.Clear();
  conv.Read(&mol);
  if (mol.Empty())
    return(1);

  vector3 v;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	
	
      // Find Input filetype
      if (pFormat == NULL) {
          pFormat = conv.FormatFromExt(FileIn);
          if (pFormat == NULL)
            {
              cerr << program_name << ": cannot read input format!" << endl;
              return (-1);
            }
      }
    }

  if (! conv.SetInAndOutFormats(pFormat, pFormat))
    {
      cerr << program_name << ": cannot read or write to this file format" << endl;
      return (-1);
    }

  // Match the SMART
  OBSmartsPattern sp;
  vector< vector <int> > maplist;      // list of matched atoms
  sp.Init(Pattern);

  OBMol mol;

  bool impossible_match;

  // Search for pattern
  for (c=0;;)
    {
      mol.Clear();
      conv.Read(&mol);
      if (mol.Empty())
        break;


      ////////////////////////////////////////////////////////////////
      // Do not loose time trying to match the pattern if the matching
      // is impossible.
      // It is impossible to make a full match if the number of atoms is
      // different
      if (full )
        impossible_match = (sp.NumAtoms() == mol.NumHvyAtoms()) ? false : true;
      else
        impossible_match = false;

      if (impossible_match)
        { // -> avoid useless SMART matching attempt
          if (invert)
            {
              if (!count)
                {
                  if ( name_only )
                    cout << mol.GetTitle() << endl;
                  else
                    conv.Write(&mol, &cout);
                }
              numMatching++;
            }
          continue;
        }


      ////////////////////////////////////////////////////////////////
      // perform SMART matching

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

void mmff94_validate()
{
  OBForceField* pFF = OBForceField::FindForceField("MMFF94");

  OBConversion conv;
  OBFormat *format_in = conv.FindFormat("mol2");
  vector<string> vs;
  vector<int> types;
  vector<double> fcharges, pcharges;
  vector<double> bond_lengths;
  char buffer[BUFF_SIZE], _logbuf[BUFF_SIZE];
  bool molfound, atomfound, bondfound, fchgfound, pchgfound;
  double etot, ebond, eangle, eoop, estbn, etor, evdw, eeq;
  double termcount; //1=bond, 2=angle, 3=strbnd, 4=torsion, 5=oop
  int n = 0;

  BOOST_REQUIRE_MESSAGE( format_in && conv.SetInFormat(format_in), "Could not set mol2 input format" );

  ifstream ifs, ifs2;
  ofstream ofs;

  ifs.open("MMFF94_dative.mol2");
  BOOST_REQUIRE_MESSAGE( ifs, "Could not open ./MMFF94_dative.mol2" );
 
  ifs2.open("MMFF94_opti.log");
  BOOST_REQUIRE_MESSAGE( ifs2, "Could not open ./MMFF94_opti.log" );
    
  ofs.open("MMFF94_openbabel.log");
  BOOST_REQUIRE_MESSAGE( ofs, "Could not open ./MMFF94_openbabel.log" );
    
  pFF->SetLogFile(&ofs);
  pFF->SetLogLevel(OBFF_LOGLVL_HIGH);
   
  OBMol mol;
  for (unsigned int c=1;; c++) {
    mol.Clear();
    types.clear();
    fcharges.clear();
    pcharges.clear();
    bond_lengths.clear();

    if (!conv.Read(&mol, &ifs))
      break;
    if (mol.Empty())
      break;
     
    BOOST_CHECK_MESSAGE( pFF->Setup(mol), "Could not setup calculations (missing parameters...)" );

    pFF->GetAtomTypes(mol);
    //pFF->GetFormalCharges(mol);
    pFF->GetPartialCharges(mol);

    termcount = 0;
    molfound = false;
    atomfound = false;
    bondfound = false;
    fchgfound = false;
    pchgfound = false;

    // Parse log file for types, charges, energies, ..
    while (ifs2.getline(buffer, 150)) {
      tokenize(vs, buffer);
      if (vs.size() == 0) {
        bondfound = false;
        continue;
      }
	
      string str(buffer);
      if (string::npos != str.find(mol.GetTitle(),0))
        molfound = true;

      // read atom types
      if (atomfound) {
        if (n) {
          types.push_back(atoi(vs[2].c_str()));
          types.push_back(atoi(vs[5].c_str()));
          types.push_back(atoi(vs[8].c_str()));
          types.push_back(atoi(vs[11].c_str()));
        } else {
          if (vs.size() > 2)
            types.push_back(atoi(vs[2].c_str()));
          if (vs.size() > 5)
            types.push_back(atoi(vs[5].c_str()));
          if (vs.size() > 8)
            types.push_back(atoi(vs[8].c_str()));
   
          atomfound = false;
        }
        n--;
      }
        
      // read formal charges
      if (fchgfound) {
        if (n) {
          fcharges.push_back(atof(vs[2].c_str()));
          fcharges.push_back(atof(vs[5].c_str()));
          fcharges.push_back(atof(vs[8].c_str()));
          fcharges.push_back(atof(vs[11].c_str()));
        } else {
          if (vs.size() > 2)
//.........这里部分代码省略.........