C++ OBMol类代码示例

  bool PDBFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
  {
    OBMol* pmol = dynamic_cast<OBMol*>(pOb);
    if(pmol==NULL)
      return false;

    //Define some references so we can use the old parameter names
    ostream &ofs = *pConv->GetOutStream();
    OBMol &mol = *pmol;

    unsigned int i;
    char buffer[BUFF_SIZE];
    char type_name[10], padded_name[10];
    char the_res[10];
    char the_chain = ' ';
    const char *element_name;
    int res_num;
    bool het=true;
    int model_num = 0;
    if (!pConv->IsLast() || pConv->GetOutputIndex() > 1)
      { // More than one molecule record
        model_num = pConv->GetOutputIndex(); // MODEL 1-based index
        snprintf(buffer, BUFF_SIZE, "MODEL %8d", model_num);
        ofs << buffer << endl;
      }

    // write back all fields (REMARKS, HELIX, SHEET, SITE, ...)
    bool compndWritten = false;
    bool authorWritten = false;
    std::vector<OBGenericData*> pairData = mol.GetAllData(OBGenericDataType::PairData);
    for (std::vector<OBGenericData*>::iterator data = pairData.begin(); data != pairData.end(); ++data) {
      OBPairData *pd = static_cast<OBPairData*>(*data);
      string attr = pd->GetAttribute();

      // filter to make sure we are writing pdb fields only
      if (attr != "HEADER" && attr != "OBSLTE" && attr != "TITLE" && attr != "SPLIT" &&
          attr != "CAVEAT" && attr != "COMPND" && attr != "SOURCE" && attr != "KEYWDS" &&
          attr != "EXPDTA" && attr != "NUMMDL" && attr != "MDLTYP" && attr != "AUTHOR" &&
          attr != "REVDAT" && attr != "SPRSDE" && attr != "JRNL" && attr != "REMARK" &&
          attr != "DBREF" && attr != "DBREF1" && attr != "DBREF2" && attr != "SEQADV" &&
          attr != "SEQRES" && attr != "MODRES" && attr != "HET" && attr != "HETNAM" &&
          attr != "HETSYN" && attr != "FORMUL" && attr != "HELIX" && attr != "SHEET" &&
          attr != "SSBOND" && attr != "LINK" && attr != "CISPEP" && attr != "SITE" &&
          attr != "ORIGX1" && attr != "ORIGX2" && attr != "ORIGX3" && attr != "SCALE1" &&
          attr != "SCALE2" && attr != "SCALE3" && attr != "MATRIX1" && attr != "MATRIX2" &&
          attr != "MATRIX3" && attr != "MODEL")
        continue;

      if (attr == "COMPND")
        compndWritten = true;
      if (attr == "AUTHOR")
        authorWritten = true;

      // compute spacing needed. HELIX, SITE, HET, ... are trimmed when reading
      int nSpacing = 6 - attr.size();
      for (int i = 0; i < nSpacing; ++i)
        attr += " ";


      std::string lines = pd->GetValue();
      string::size_type last = 0;
      string::size_type pos = lines.find('\n');
      while (last != string::npos) {
        string line = lines.substr(last, pos - last);
        if (pos == string::npos)
          last = string::npos;
        else
          last = pos + 1;
        pos = lines.find('\n', last);

        ofs << attr << line << endl;
      }
    }

    if (!compndWritten) {
      if (strlen(mol.GetTitle()) > 0)
        snprintf(buffer, BUFF_SIZE, "COMPND    %s ",mol.GetTitle());
      else
        snprintf(buffer, BUFF_SIZE, "COMPND    UNNAMED");
      ofs << buffer << endl;
    }

    if (!authorWritten) {
      snprintf(buffer, BUFF_SIZE, "AUTHOR    GENERATED BY OPEN BABEL %s",BABEL_VERSION);
      ofs << buffer << endl;
    }

    // Write CRYST1 record, containing unit cell parameters, space group
    // and Z value (supposed to be 1)
    if (pmol->HasData(OBGenericDataType::UnitCell))
      {
        OBUnitCell *pUC = (OBUnitCell*)pmol->GetData(OBGenericDataType::UnitCell);
        if(pUC->GetSpaceGroup()){
          string tmpHM=pUC->GetSpaceGroup()->GetHMName();
          // Do we have an extended HM symbol, with origin choice as ":1" or ":2" ? If so, remove it.
          size_t n=tmpHM.find(":");
          if(n!=string::npos) tmpHM=tmpHM.substr(0,n);
          snprintf(buffer, BUFF_SIZE,
                   "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s 1",
                   pUC->GetA(), pUC->GetB(), pUC->GetC(),
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	output.open(filepath, ios::out | ios::app ); // The file is open in append mode

	//              1     2       3        4        5          6          7        8         9          10        11    12
	output << "#METHOD: " << ff << " DATASET: " << filename.c_str() << std::endl;
	output << "#THREADS ENERGY MOL_MASS NUM_ATOMS NUM_ROTORS NUM_CONF TOT_TIME TIME_READ TIME_SETUP TIME_COMPUTE STEPS #MOL_NAME " << std::endl;

	// A second file is created to store information on the timings of specific parts
	// from the MMFF94 calculation routines. breakdown of time spent in calculations
	char filepath2[1100];
	std::ofstream output2;
	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();

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

  map<string,string>::const_iterator itr;

  // Is there an -s option?
  if(pConv->IsFirstInput())
  {
    _pOpIsoM = NULL; //assume no -s option
    itr = pmap->find("s");
    if(itr!=pmap->end())
    {
      //There is an -s option; check it is ok
      _pOpIsoM = static_cast<OpNewS*>(OBOp::FindType("s"));
      _stext = itr->second; //get its parameter(s)
      if(!_pOpIsoM || _stext.empty())
      {
        obErrorLog.ThrowError(__FUNCTION__,
        "No parameter on -s option, or its OBOp version is not loaded", obError);
        pConv->SetOneObjectOnly(); //to finish
        return false;
      }
    }
  }

  // If the output format is a 2D depiction format, then we should align
  // on the 2D coordinates and not the 3D coordinates (if present). This
  //means we need to generate the 2D coordinates at this point.
  if(pmol->GetDimension()==3 && (pConv->GetOutFormat()->Flags() & DEPICTION2D))
  {
    OBOp* pgen = OBOp::FindType("gen2D");
    if(pgen)
      pgen->Do(pmol);
  }

  // All molecules must have coordinates, so add them if 0D
  // They may be added again later when gen2D or gen3D is called, but they will be the same.
  // It would be better if this op was called after them, which would happen
  // if its name was alphabetically after "gen" (and before "s").
  if(pmol->GetDimension()==0)
  {
    //Will the coordinates be 2D or 3D?
    itr = pmap->find("gen3D");
    OBOp* pgen = (itr==pmap->end()) ? OBOp::FindType("gen2D") : OBOp::FindType("gen3D");
    if(pgen)
      pgen->Do(pmol);
  }

  //Do the alignment in 2D if the output format is svg, png etc. and there is no -xn option
  if(pmol->GetDimension()==3 && pConv && !pConv->IsOption("n"))
  {
    OBFormat* pOutFormat = pConv->GetOutFormat();
    if(pOutFormat->Flags() & DEPICTION2D)
    {
      OBOp* pgen = OBOp::FindType("gen2D");
      if(pgen)
        pgen->Do(pmol);
    }
  }

  if(pConv->IsFirstInput() || _refMol.NumAtoms()==0)
  {
    _refvec.clear();
    // Reference molecule is basically the first molecule
    _refMol = *pmol;
    if(!_pOpIsoM)
     //no -s option. Use a molecule reference.
     _align.SetRefMol(_refMol);
    else
    {
      //If there is a -s option, reference molecule has only those atoms that are matched
      //Call the -s option from here
      bool ret = _pOpIsoM->Do(pmol, _stext.c_str(), pmap, pConv);
      // Get the atoms that were matched
      vector<int> ats = _pOpIsoM->GetMatchAtoms();
      if(!ats.empty())
      {
        // Make a vector of the matching atom coordinates...
        for(vector<int>::iterator iter=ats.begin(); iter!=ats.end(); ++iter)
          _refvec.push_back((pmol->GetAtom(*iter))->GetVector());
        // ...and use a vector reference
        _align.SetRef(_refvec);
      }
      // Stop -s option being called normally, although it will still be called once
      //  in the DoOps loop already started for the current (first) molecule.
      pConv->RemoveOption("s",OBConversion::GENOPTIONS);
      if(!ret)
      {
        // the first molecule did not match the -s option so a reference molecule
        // could not be made. Keep trying.
        _refMol.Clear();
        //obErrorLog.ThrowError(__FUNCTION__, "The first molecule did not match the -s option\n"
        //  "so the reference structure was not derived from it", obWarning, onceOnly);
        return false; //not matched
      }
    }
  }
//.........这里部分代码省略.........

  void VibrationWidget::setMolecule(Molecule *molecule)
  {
    // update table
    ui.vibrationTable->clearContents();
      if (molecule == 0){
        ui.vibrationTable->setRowCount(0);
        ui.vibrationTable->horizontalHeader()->hide();
        return;
      }
    m_molecule = molecule;

    OBMol obmol = molecule->OBMol();
    m_vibrations = static_cast<OBVibrationData*>(obmol.GetData(OBGenericDataType::VibrationData));
    if (!m_vibrations) {
      ui.vibrationTable->setRowCount(0);
      ui.vibrationTable->horizontalHeader()->hide();
      //ui.exportButton->setEnabled(false);
      return;
    }

    ui.vibrationTable->horizontalHeader()->show();
    ui.vibrationTable->horizontalHeader()->setResizeMode(QHeaderView::Stretch);

    // OK, we have valid vibrations, so add them to the table
    vector<double> frequencies = m_vibrations->GetFrequencies();
    vector<double> intensities = m_vibrations->GetIntensities();
    m_frequencies = frequencies;
    m_intensities = intensities;
    vector<double> raman_activities = m_vibrations->GetRamanActivities();
    if (raman_activities.size() == 0) {
      //resize(274, height());
      ui.vibrationTable->setColumnCount(2);
      if(parentWidget())
        parentWidget()->setMinimumWidth(274);
    }
    else {
      //resize(310, height());
      ui.vibrationTable->setColumnCount(3);
      ui.vibrationTable->setHorizontalHeaderItem(2,
                                                 new QTableWidgetItem("Activity"));
      if(parentWidget())
        parentWidget()->setMinimumWidth(310);
    }

    // Generate an index vector to map sorted indicies to the old indices
    m_indexMap->clear();
    for (uint i = 0; i < frequencies.size(); i++)
      m_indexMap->push_back(i);

    // Setup progress dialog, just in case it takes longer than 2 seconds
    QProgressDialog prog(tr("Sorting %1 vibrations by frequency...")
                         .arg(frequencies.size()), "", 0, frequencies.size());
    prog.setWindowModality(Qt::WindowModal);
    prog.setMinimumDuration(2000);
    prog.setCancelButton(0);

    // Simple selection sort
    double tmp;
    int tmp_int;
    for (uint i = 0; i < frequencies.size(); i++) {
      for (uint j = i; j < frequencies.size(); j++) {
        if (i == j)
          continue; // Save a bit of time...
        if (frequencies.at(j) < frequencies.at(i)) {
          tmp = frequencies.at(j);
          frequencies.at(j) = frequencies.at(i);
          frequencies.at(i) = tmp;
          tmp = intensities.at(j);
          intensities.at(j) = intensities.at(i);
          intensities.at(i) = tmp;
          tmp_int = m_indexMap->at(j);
          m_indexMap->at(j) = m_indexMap->at(i);
          m_indexMap->at(i) = tmp_int;
        }
      }
      // Update progress bar
      prog.setValue(i);
    }

    ui.vibrationTable->setRowCount(frequencies.size());
    QString format("%1");

    for (unsigned int row = 0; row < frequencies.size(); ++row) {
        QTableWidgetItem *newFreq = new QTableWidgetItem(format.arg(frequencies[row], 0, 'f', 2));
        newFreq->setTextAlignment(Qt::AlignRight|Qt::AlignVCenter);
        // Some codes don't provide intensity data. Display "-" in place of intensities.
        QTableWidgetItem *newInten;
        if (row >= intensities.size()) {
            newInten = new QTableWidgetItem("-");
        }
        else {
            newInten = new QTableWidgetItem(format.arg(intensities[row], 0, 'f', 3));
        }
        newInten->setTextAlignment(Qt::AlignRight|Qt::AlignVCenter);
        ui.vibrationTable->setItem(row, 0, newFreq);
        ui.vibrationTable->setItem(row, 1, newInten);

        if (raman_activities.size() != 0) {
            QTableWidgetItem *newRaman;
            if (row >= raman_activities.size()) {
//.........这里部分代码省略.........

  bool COFFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
  {
    OBMol* pmol = dynamic_cast<OBMol*>(pOb);
    if(pmol==NULL)
      return false;
    stringstream errorMsg;

    ostream& ofs = *(pConv->GetOutStream());
    OBMol &mol = *pmol;

    ofs << "@CulgiVersion: 10.0.0\n"
      << "# Culgi Object File\n"
      << "# Generated by Open Babel\n\n"
      << "moleculekeys\n"
      << "\tname\n"
      << "\tid\n"
      << "end_moleculekeys\n\n"
      << "atomkeys\n"
      << "\tindex\n"
      << "\tname\n"
      << "\telement_type\n"
      << "\tforce_field_type\n"
      << "\tcharge\n"
      << "\tfixed\n"
      << "\tx\n"
      << "\ty\n"
      << "\tz\n"
      << "\tvelocity_x\n"
      << "\tvelocity_y\n"
      << "\tvelocity_z\n"
      << "\tvelocity_const_x\n"
      << "\tvelocity_const_y\n"
      << "\tvelocity_const_z\n"
      << "end_atomkeys\n\n"
      << "bondkeys\n"
      << "\tatom1_index\n"
      << "\tatom2_index\n"
      << "\tbond_order\n"
      << "end_bondkeys\n\n";

    FOR_ATOMS_OF_MOL(atom, mol)
    {
      if(atom->GetFormalCharge()!=0)
      {
        ofs << "formalqkeys\n"
          << "\tatom_index\n"
          << "\tformal_charge\n"
          << "end_formalqkeys\n\n";
        break;
      }
    }

    // Get molecule file name.  Some formats add file path
    // etc to the name, so we try to clean that up here
    const char *tit = pmol->GetTitle();
    std::string molname(tit);
    if(molname.empty())
      molname = pConv->GetTitle();
    size_t nPos = molname.find_last_of(".");
    if(nPos != std::string::npos)
      molname = molname.substr(0, nPos);
    nPos = molname.find_last_of("\\");
    if(nPos != std::string::npos)
      molname = molname.substr(nPos+1, molname.size());
    nPos = molname.find_last_of("/");
    if(nPos != std::string::npos)
      molname = molname.substr(nPos+1, molname.size());
    if(molname.empty())
      molname = "mol";
    ofs  << "molecule\t" << molname << "\t0" << endl;

    int i = 0;
    vector<string> names;
    vector<string> elems;
    vector<int> nelem;
    ostringstream sstream;
    FOR_ATOMS_OF_MOL(atom, mol)
    {
      sstream.str("");
      i++;
      string elem = OBElements::GetSymbol(atom->GetAtomicNum());
      // Culgi does not recognize atom type 'Xx' but does know 'X'.
      if(elem == "Xx")
        elem = "X";
      bool found = false;
      string ename;
      for( int j=0; j<elems.size(); j++)
      {
        if(elem == elems[j])
        {
          found = true;
          nelem[j]++;
          sstream.str("");
          sstream << elem << nelem[j];
          ename = sstream.str();
          names.push_back(ename);
          sstream.str("");
          break;
        }
      }
//.........这里部分代码省略.........

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

  bool OBResidueData::AssignBonds(OBMol &mol,OBBitVec &bv)
  {
    if (!_init)
      Init();

    OBAtom *a1,*a2;
    OBResidue *r1,*r2;
    vector<OBAtom*>::iterator i,j;
    vector3 v;

    int bo;
    string skipres = ""; // Residue Number to skip
    string rname = "";
    //assign residue bonds
    for (a1 = mol.BeginAtom(i);a1;a1 = mol.NextAtom(i))
      {
        r1 = a1->GetResidue();
        if (r1 == NULL) // atoms may not have residues
          continue;

        if (skipres.length() && r1->GetNumString() == skipres)
          continue;

        if (r1->GetName() != rname)
          {
            skipres = SetResName(r1->GetName()) ? "" : r1->GetNumString();
            rname = r1->GetName();
          }
        //assign bonds for each atom
        for (j=i,a2 = mol.NextAtom(j);a2;a2 = mol.NextAtom(j))
          {
            r2 = a2->GetResidue();
            if (r2 == NULL) // atoms may not have residues
              continue;

            if (r1->GetNumString() != r2->GetNumString())
              break;
            if (r1->GetName() != r2->GetName())
              break;
            if (r1->GetChain() != r2->GetChain())
              break; // Fixes PR#2889763 - Fabian

            if ((bo = LookupBO(r1->GetAtomID(a1),r2->GetAtomID(a2))))
              {
                // Suggested by Liu Zhiguo 2007-08-13
                // for predefined residues, don't perceive connection
                // by distance
                //                v = a1->GetVector() - a2->GetVector();
                //                if (v.length_2() < 3.5) //check by distance
                  mol.AddBond(a1->GetIdx(),a2->GetIdx(),bo);
              }
          }
      }

    int hyb;
    string type;

    //types and hybridization
    rname = ""; // name of current residue
    skipres = ""; // don't skip any residues right now
    for (a1 = mol.BeginAtom(i);a1;a1 = mol.NextAtom(i))
      {
        if (a1->GetAtomicNum() == OBElements::Oxygen && !a1->GetValence())
          {
            a1->SetType("O3");
            continue;
          }
        if (a1->GetAtomicNum() == OBElements::Hydrogen)
          {
            a1->SetType("H");
            continue;
          }

        //***valence rule for O-
        if (a1->GetAtomicNum() == OBElements::Oxygen && a1->GetValence() == 1)
          {
            OBBond *bond;
            bond = (OBBond*)*(a1->BeginBonds());
            if (bond->GetBO() == 2)
              {
                a1->SetType("O2");
                a1->SetHyb(2);
              }
            else if (bond->GetBO() == 1)
              {
                // Leave the protonation/deprotonation to phmodel.txt
                a1->SetType("O3");
                a1->SetHyb(3);
                // PR#3203039 -- Fix from Magnus Lundborg
                //                a1->SetFormalCharge(0);
              }
            continue;
          }

        r1 = a1->GetResidue();
        if (r1 == NULL) continue; // atoms may not have residues
        if (skipres.length() && r1->GetNumString() == skipres)
          continue;

        if (r1->GetName() != rname)
//.........这里部分代码省略.........

///////////////////////////////////////////////////////////////////////////////
//! \brief Set a tortional bond to a given angle
int main(int argc,char **argv)
{
  const char *Pattern=NULL;
  unsigned int i, t, errflg = 0;
  int c;
  char flags[255];
  string err;
  bool graphOutput=false;

  // parse the command line -- optional -a flag to change all matching torsions
  if (argc < 3 || argc > 4) {
    errflg++;
  } else {
    FileIn = argv[1];
    Pattern = "[!$(*#*)&!D1][email protected][!$(*#*)&!D1]";
    // Read the atom position
    c = sscanf(argv[2], "%d", &angleSum);
	angle = 360./angleSum;
    if (argc == 4)
	{
    		c = sscanf(argv[3], "%s", flags);
		int flagid=1;
    		while (flags[flagid]!=0)
			switch (flags[flagid++])
			{
			case 'g':
				graphOutput=true;
			case 'e':
				forceField=OBForceField::FindForceField("MMFF94");
				isEnergyCalcing=true;
				break;
    			}
 	}
  }
  if (errflg) {
    cerr << "Usage: rkrotate <filename> <angle> [options]" << endl;
    exit(-1);
  }

  // create pattern
  OBSmartsPattern sp;
  sp.Init(Pattern);

  OBFormat* format = conv.FormatFromExt(FileIn);
  if(!(format && conv.SetInAndOutFormats(format, format))) { //in and out formats same
    cerr << "obrotate: cannot read and/or write this file format!" << endl;
    exit (-1);
  } //...NF

  //Open the molecule file
  ifstream ifs;

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

  OBMol mol;
  vector< vector <int> > maplist;      // list of matched atoms
//  vector< vector <int> >::iterator m;  // and its iterators
  //   int tindex;
  
  // Set the angles
  for (;;) {
    mol.Clear();
    //NF      ifs >> mol;                   // Read molecule
    conv.Read(&mol,&ifs); //NF
    if (mol.Empty())
      break;

    if (sp.Match(mol)) {          
      // if match perform rotation
      maplist = sp.GetUMapList(); // get unique matches
      
      if (maplist.size() > 1)
        cerr << "obrotate: Found " << maplist.size() << " matches." << endl;

	energySheet=new MultiVector<double>(degrees=maplist.size(),angleSum);
	indexSheet=new int[maplist.size()];

      for (int EXO=0;EXO<maplist.size();++EXO)
	totalSum*=angleSum+EXO;
      // look at all the mapping atom but save only the last one.
	turnMol(mol,maplist,maplist.size()-1);
	
      if (graphOutput)
      {
	ofstream ofs("energyGraph.mlog");
	int ind[degrees];
	for (int i=0;i<degrees;++i)
		ind[i]=0;
	do
	{
		for (int i=0;i<degrees;++i)
			ofs<<ind[i]<<'\t';
		ofs<<energySheet->getVectorValue(ind)<<endl;
//.........这里部分代码省略.........

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

  char buffer[BUFF_SIZE];
  vector<OBSmartsPattern*> vsp;
  for (;ifs.getline(buffer,BUFF_SIZE);)
    {
      if (buffer[0] == '#') // skip comment line
        continue;

      OBSmartsPattern *sp = new OBSmartsPattern;

      if (sp->Init(buffer))
        vsp.push_back(sp);
      else
        delete sp;
    }

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

  ofs << vsp.size() << " patterns" << endl;
  std::ifstream mifs;
  if (!SafeOpen(mifs, smilestypes_file.c_str()))
    return;

  vector<int> vm;
  vector<vector<int> > mlist;
  vector<vector<int> >::iterator j;
  vector<OBSmartsPattern*>::iterator i;
  OBMol mol;
  OBConversion conv(&mifs, &cout);

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

  for (;mifs;)
    {
      mol.Clear();
      conv.Read(&mol);

      if (mol.Empty())
        continue;
      for (i = vsp.begin();i != vsp.end();i++)
        {
          (*i)->Match(mol);
          mlist = (*i)->GetMapList();
          for (j = mlist.begin();j != mlist.end();j++)
            {
              sprintf(buffer,"%3d",*(j->begin()));
              ofs << buffer;
            }
          ofs << endl;
        }
    }


  cerr << " SMARTS test results written successfully" << endl;
  return;
}

  bool EEMCharges::ComputeCharges(OBMol &mol)
  {
    mol.SetPartialChargesPerceived();

    if(_parameters.empty())
      _loadParameters();

    // Copied from spectrophore.cpp
    // CHI and ETA
    unsigned int _nAtoms = mol.NumAtoms();
    unsigned int dim(_nAtoms + 1);
    std::vector<double> CHI(dim);
    double** ETA = new double*[dim];
    for (unsigned int i = 0; i < dim; ++i)
      {
        ETA[i] = new double[dim];
      }
    double totalCharge(0.0);
    unsigned int i(0);
    double hardness;
    double electronegativity;
    for (OpenBabel::OBMolAtomIter atom(mol); atom; atom++, i++) {

      int n = atom->GetAtomicNum();
      int b = atom->HighestBondOrder();

      // Search for parameters for a particular atom type
      bool found = false;
      for(unsigned int j = 0; j < _parameters.size(); j++) {
        if((_parameters[j].Z == n && _parameters[j].bond_order == b) ||
            (_parameters[j].Z == n && _parameters[j].bond_order == - 1) ||
            (_parameters[j].Z == -1 && _parameters[j].bond_order == -1)) {

          electronegativity = _parameters[j].A;
          hardness = _parameters[j].B;
          found = true;
          break;
        }
      }

      if(!found) {
        std::stringstream ss;
        ss << "No parameters found for: " << etab.GetSymbol(n) << " " << b
           << ". EEM charges were not calculated for the molecule." << std::endl;
        obErrorLog.ThrowError(__FUNCTION__, ss.str(), obError);
        return false;
      }

      CHI[i] = -electronegativity;
      ETA[i][i] = hardness;

      // Adjust the total molecular charge
      totalCharge += atom->GetFormalCharge();
    }

    // Complete CHI
    CHI[_nAtoms] = totalCharge;

    // Complete ETA
    OBAtom *rAtom, *cAtom;
    for (unsigned int r = 0; r < _nAtoms; ++r)
      {
        rAtom = mol.GetAtom(r+1); // Atom index
        for (unsigned int c = r + 1; c < _nAtoms; ++c)
          {
            cAtom = mol.GetAtom(c+1); // Atom index
            ETA[r][c] = _kappa / cAtom->GetDistance(rAtom);
            ETA[c][r] = ETA[r][c];
          }
      }
    for (unsigned int i = 0; i < dim; ++i)
      {
        ETA[i][_nAtoms] = -1.0;
        ETA[_nAtoms][i] = +1.0;
      }
    ETA[_nAtoms][_nAtoms] = 0.0;

    // Solve the matrix equation
    _solveMatrix(ETA, &(CHI[0]), dim);    // CHI will contain the values

    OBAtom *atom;
    for (unsigned int i = 0; i < _nAtoms; ++i)
      {
        atom = mol.GetAtom(i+1); // atom index issue
        atom->SetPartialCharge(CHI[i]);
      }

    OBChargeModel::FillChargeVectors(mol);

    // Cleanup
    for(unsigned int i = 0; i < dim; i++)
      delete [] ETA[i];

    delete [] ETA;

    return true;
  }

  //! \return whether partial charges were successfully assigned to this molecule
  bool EQEqCharges::ComputeCharges(OBMol &mol)
  {
    int i, j, a, c, N = mol.NumAtoms();
    double cellVolume;
    VectorXf chi(N), J(N), b(N), x(N);
    MatrixXf J_ij(N, N), A(N, N);
    OBUnitCell *obuc;
    matrix3x3 unitcell, fourier;
    vector3 dx;
    int numNeighbors[3];
    OBAtom *atom;

    // If parameters have not yet been loaded, do that
    if (!_paramFileLoaded)
    {
      if (ParseParamFile())
      {
        _paramFileLoaded = true;
      } else
      {
        return false;
      }
    }

    // Calculate atomic properties based around their ionic charge
    for (i = 0; i < N; i++)
    {
      atom = mol.GetAtom(i + 1);
      a = atom->GetAtomicNum();
      c = _chargeCenter[a];

      // Fail if ionization data is missing for any atom in the molecule
      if (_ionizations[a][c + 1] == -1 || _ionizations[a][c] == -1 || a > TABLE_OF_ELEMENTS_SIZE)
      {
        obErrorLog.ThrowError(__FUNCTION__, "Insufficient ionization data for atoms in the given molecule. Update `data/eqeqIonizations.txt` with missing information and re-run this function.", obError);
        return false;
      }

      J(i) = _ionizations[a][c + 1] - _ionizations[a][c];
      chi(i) = 0.5 * (_ionizations[a][c + 1] + _ionizations[a][c]) - (a == 1? 0 : c * J(i));
    }

    // If a unit cell is defined, use the periodic Ewald calculation
    if (mol.HasData(OBGenericDataType::UnitCell))
    {
      // Get unit cell and calculate its Fourier transform + volume
      obuc = (OBUnitCell *) mol.GetData(OBGenericDataType::UnitCell);
      unitcell = obuc->GetCellMatrix();
      fourier = (2 * PI * unitcell.inverse()).transpose();
      cellVolume = obuc->GetCellVolume();

      // Get the number of radial unit cells to use in x, y, and z
      numNeighbors[0] = int(ceil(minCellLength / (2.0 * (obuc->GetA())))) - 1;
      numNeighbors[1] = int(ceil(minCellLength / (2.0 * (obuc->GetB())))) - 1;
      numNeighbors[2] = int(ceil(minCellLength / (2.0 * (obuc->GetC())))) - 1;

      for (i = 0; i < N; i++)
      {
        atom = mol.GetAtom(i + 1);
        for (j = 0; j < N; j++)
        {
          dx = atom->GetVector() - (mol.GetAtom(j + 1))->GetVector();
          J_ij(i, j) = GetPeriodicEwaldJij(J(i), J(j), dx, (i == j), unitcell, fourier, cellVolume, numNeighbors);
        }
      }
    // If no unit cell, use the simplified nonperiodic calculation
    } else
    {
      for (i = 0; i < N; i++)
      {
        atom = mol.GetAtom(i + 1);
        for (j = 0; j < N; j++)
        {
          J_ij(i, j) = GetNonperiodicJij(J(i), J(j), atom->GetDistance(j + 1), (i == j));
        }
        return false;
      }
    }

    // Formulate problem as A x = b, where x is the calculated partial charges
    // First equation is a simple overall balance: sum(Q) = 0
    A.row(0) = VectorXf::Ones(N);
    b(0) = 0;

    // Remaining equations are based off of the fact that, at equilibrium, the
    // energy of the system changes equally for a change in any charge:
    //     dE/dQ_1 = dE/dQ_2 = ... = dE/dQ_N
    A.block(1, 0, N - 1, N) = J_ij.block(0, 0, N - 1, N) - J_ij.block(1, 0, N - 1, N);
    b.tail(N - 1) = chi.tail(N - 1) - chi.head(N - 1);

    // The solution is a list of charges in the system
    x = A.colPivHouseholderQr().solve(b);

    // Now we are done calculating, pass all this back to OpenBabel molecule
    mol.SetPartialChargesPerceived();
    OBPairData *dp = new OBPairData;
    dp->SetAttribute("PartialCharges");
    dp->SetValue("EQEq");
    dp->SetOrigin(perceived);
//.........这里部分代码省略.........

  bool FastSearchFormat::ObtainTarget(OBConversion* pConv, vector<OBMol>& patternMols, const string& indexname)
  {
    //Obtains an OBMol from:
    // the filename in the -s option or
    // the SMARTS string in the -s option or
    // by converting the file in the -S or -aS options (deprecated).
    // If there is no -s -S or -aS option, information on the index file is displayed.

    OBMol patternMol;
    patternMol.SetIsPatternStructure();

    const char* p = pConv->IsOption("s",OBConversion::GENOPTIONS);

    bool OldSOption=false;
    //If no -s option, make OBMol from file in -S option or -aS option (both deprecated)
    if(!p)
    {
      p = pConv->IsOption("S",OBConversion::GENOPTIONS);
      if(!p)
        p = pConv->IsOption("S",OBConversion::INOPTIONS);//for GUI mainly
      OldSOption = true;
    }
    if(p)
    {
      vector<string> vec;
      tokenize(vec, p);

      //ignore leading ~ (not relevant to fastsearch)
      if(vec[0][0]=='~')
        vec[0].erase(0,1);

      if(vec.size()>1 && vec[1]=="exact")
        pConv->AddOption("e", OBConversion::INOPTIONS);

      OBConversion patternConv;
      OBFormat* pFormat;
      //Interpret as a filename if possible
      string& txt =vec [0];
      if( txt.empty() ||
          txt.find('.')==string::npos ||
          !(pFormat = patternConv.FormatFromExt(txt.c_str())) ||
          !patternConv.SetInFormat(pFormat) ||
          !patternConv.ReadFile(&patternMol, txt) ||
          patternMol.NumAtoms()==0)
        //if false, have a valid patternMol from a file
      {
        //is SMARTS/SMILES
        //Replace e.g. [#6] in SMARTS by C so that it can be converted as SMILES
        //for the fingerprint phase, but allow more generality in the SMARTS phase.
        for(;;)
        {
          string::size_type pos1, pos2;
          pos1 = txt.find("[#");
          if(pos1==string::npos)
            break;
          pos2 = txt.find(']');
          int atno;
          if(pos2!=string::npos &&  (atno = atoi(txt.substr(pos1+2, pos2-pos1-2).c_str())) && atno>0)
            txt.replace(pos1, pos2-pos1+1, etab.GetSymbol(atno));
          else
          {
            obErrorLog.ThrowError(__FUNCTION__,"Ill-formed [#n] atom in SMARTS", obError);
            return false;
          }
        }

        bool hasTildeBond;
        if( (hasTildeBond = (txt.find('~')!=string::npos)) ) // extra parens to indicate truth value
        {
          //Find ~ bonds and make versions of query molecule with a single and aromatic bonds
          //To avoid having to parse the SMILES here, replace ~ by $ (quadruple bond)
          //and then replace this in patternMol. Check first that there are no $ already
          //Sadly, isocynanides may have $ bonds.
          if(txt.find('$')!=string::npos)
          {
            obErrorLog.ThrowError(__FUNCTION__,
              "Cannot use ~ bonds in patterns with $ (quadruple) bonds.)", obError);
            return false;
          }
          replace(txt.begin(),txt.end(), '~' , '$');
        }

        //read as standard SMILES
        patternConv.SetInFormat("smi");
        if(!patternConv.ReadString(&patternMol, vec[0]))
        {
          obErrorLog.ThrowError(__FUNCTION__,"Cannot read the SMILES string",obError);
          return false;
        }
        if(hasTildeBond)
        {
          AddPattern(patternMols, patternMol, 0); //recursively add all combinations of tilde bond values
          return true;
        }
      }
    }

    if(OldSOption) //only when using deprecated -S and -aS options
    {
      //make -s option for later SMARTS test
//.........这里部分代码省略.........

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

  // The SMARTS and any other parameters are extracted on the first molecule
  // and stored in the static variables vec, inv. The parameter is cleared so that:
  // (a) the original -s option in transform.cpp is inactive, and
  // (b) the parsing does not have to be done again for multi-molecule files

  string txt(pmap->find(GetID())->second); // ID can be "s" or "v"
  static vector<string> vec;
  static bool inv;
  static int nPatternAtoms; //non-zero for exact matches
  static OBQuery* query;
  static vector<OBQuery*> queries;
  vector<OBQuery*>::iterator qiter;
  if(!txt.empty())
  {
    //Set up on first call
    tokenize(vec, txt);
    inv = GetID()[0]=='v';
    if(vec[0][0]=='~')
    {
      inv = true;
      vec[0].erase(0,1);
    }

    //Interpret as a filename if possible
    MakeQueriesFromMolInFile(queries, vec[0], &nPatternAtoms);

    if(vec.size()>1 && vec[1]=="exact")
    {
      if(queries.empty())
      {
        //Convert SMARTS to SMILES to count number of atoms
        OBConversion conv;
        OBMol patmol;
        if(!conv.SetInFormat("smi") || !conv.ReadString(&patmol, vec[0]))
        {
          obErrorLog.ThrowError(__FUNCTION__, "Cannot read the parameter of -s option, "
          "which has to be valid SMILES when the exact option is used.", obError, onceOnly);
          delete pmol;
          pConv->SetOneObjectOnly(); //stop conversion
          return false;
        }
        nPatternAtoms = patmol.NumHvyAtoms();
      }
    }
    else
      nPatternAtoms = 0;
    
    //disable old versions
    pConv->AddOption(GetID(), OBConversion::GENOPTIONS, "");
  }

  bool match;
  //These are a vector of each mapping, each containing atom indxs.
  vector<vector<int> > vecatomvec;
  vector<vector<int> >* pMappedAtoms = NULL;
  OBSmartsPattern sp;

  if(nPatternAtoms)
    if(pmol->NumHvyAtoms() != nPatternAtoms)
      return false;

  int imol=0; //index of mol in pattern file
  if(!queries.empty()) //filename supplied
  {
    //match is set true if any of the structures match - OR behaviour
    for(qiter=queries.begin();qiter!=queries.end();++qiter, ++imol) 
    {
      OBIsomorphismMapper* mapper = OBIsomorphismMapper::GetInstance(*qiter);
      OBIsomorphismMapper::Mappings mappings;
      mapper->MapUnique(pmol, mappings);
      if( (match = !mappings.empty()) ) // extra parens to indicate truth value
      {
        OBIsomorphismMapper::Mappings::iterator ita;
        OBIsomorphismMapper::Mapping::iterator itb;
        for(ita=mappings.begin(); ita!=mappings.end();++ita)//each mapping
        {
          vector<int> atomvec;
          for(itb=ita->begin(); itb!=ita->end();++itb)//each atom index
            atomvec.push_back(itb->second+1);
          vecatomvec.push_back(atomvec);
          atomvec.clear();
        }
        pMappedAtoms = &vecatomvec;
        break;
      }
    }
  }
  else //SMARTS supplied
  {
    if(!sp.Init(vec[0]))
    {
      string msg = vec[0] + " cannot be interpreted as either valid SMARTS "
        "or the name of a file with an extension known to OpenBabel "
        "that contains one or more pattern molecules.";
//.........这里部分代码省略.........

  static bool parseAtomRecord(char *buffer, OBMol &mol,int /*chainNum*/)
  /* ATOMFORMAT "(i5,1x,a4,a1,a3,1x,a1,i4,a1,3x,3f8.3,2f6.2,a2,a2)" */
  {
    string sbuf = &buffer[6];
    if (sbuf.size() < 48)
      return(false);

    bool hetatm = (EQn(buffer,"HETATM",6)) ? true : false;
    bool elementFound = false; // true if correct element found in col 77-78

    /* serial number */
    string serno = sbuf.substr(0,5);

    /* atom name */
    string atmid = sbuf.substr(6,4);

    /* chain */
    char chain = sbuf.substr(15,1)[0];

    /* element */
    string element = "  ";
    if (sbuf.size() > 71)
      {
        element = sbuf.substr(70,2);
        if (isalpha(element[1]))
          {
            if (element[0] == ' ')
              {
                element.erase(0, 1);
                elementFound = true;
              }
            else if (isalpha(element[0]))
              {
                elementFound = true;
              }
          }
      }

    if (!elementFound)
      {
        stringstream errorMsg;
        errorMsg << "WARNING: Problems reading a PDB file\n"
                 << "  Problems reading a HETATM or ATOM record.\n"
                 << "  According to the PDB specification,\n"
                 << "  columns 77-78 should contain the element symbol of an atom.\n"
                 << "  but OpenBabel found '" << element << "' (atom " << mol.NumAtoms()+1 << ")";
        obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obWarning);
      }

    // charge - optional
    string scharge;
    if (sbuf.size() > 73)
      {
        scharge = sbuf.substr(72,2);
      }

    //trim spaces on the right and left sides
    while (!atmid.empty() && atmid[0] == ' ')
      atmid = atmid.erase(0, 1);

    while (!atmid.empty() && atmid[atmid.size()-1] == ' ')
      atmid = atmid.substr(0,atmid.size()-1);

    /* residue name */
    string resname = sbuf.substr(11,3);
    if (resname == "   ")
      resname = "UNK";
    else
      {
        while (!resname.empty() && resname[0] == ' ')
          resname = resname.substr(1,resname.size()-1);

        while (!resname.empty() && resname[resname.size()-1] == ' ')
          resname = resname.substr(0,resname.size()-1);
      }

    string type;
    if (!elementFound) {
      // OK, we have to fall back to determining the element from the atom type
      // This is unreliable, but there's no other choice
      if (EQn(buffer,"ATOM",4)) {
        type = atmid.substr(0,2);
        if (isdigit(type[0])) {
          // sometimes non-standard files have, e.g 11HH
          if (!isdigit(type[1])) type = atmid.substr(1,1);
          else type = atmid.substr(2,1);
        } else if ((sbuf[6] == ' ' &&
                   strncasecmp(type.c_str(), "Zn", 2) != 0 &&
                   strncasecmp(type.c_str(), "Fe", 2) != 0) ||
                   isdigit(type[1]))	//type[1] is digit in Platon
          type = atmid.substr(0,1);     // one-character element


        if (resname.substr(0,2) == "AS" || resname[0] == 'N') {
          if (atmid == "AD1")
            type = "O";
          if (atmid == "AD2")
            type = "N";
        }
        if (resname.substr(0,3) == "HIS" || resname[0] == 'H') {
//.........这里部分代码省略.........

  bool ACRFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
  {
    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if(pmol==NULL)
      return false;

    istream& ifs = *pConv->GetInStream();

    pmol->BeginModify();

    /** Parse the input stream and use the OpenBabel API to populate the OBMol **/
    int id;
    char buf[BUFF_SIZE];
    int atoms, bonds, tmp;
    float scale, dtmp;
    bool atom_input = false, bond_input = false;
    string type;
    //int from, to;
    double X,Y,Z;
    vector<string> vs;

    // read in one at a time
    /* WARNING: Atom id starts from zero in Carine; not so in openbabel.
     * Solution: increment atom id's */

    while (true) {
      ifs.getline(buf, BUFF_SIZE);
      if (ifs.eof()) {
        break;
      }

      if (sscanf(buf, "General Scale=%f\n", &dtmp)) {
        scale = dtmp;
        continue;
      } else if (sscanf(buf, "Number of Atoms in Crystal=%d\n", &tmp)) {
        atoms = tmp;
        atom_input = true;

        // read table column names
        ifs.getline(buf, BUFF_SIZE);
        continue;
      } else if (sscanf(buf, "Number of Links in Crystal=%d\n", &tmp)) {
        atom_input = false;
        bond_input = true;
        bonds = tmp;

        // read table column names
        ifs.getline(buf, BUFF_SIZE);
        continue;
      } else if ( '#' == buf[0] || '\r' == buf[0] || '\n' == buf[0] ) {
        // between sections, in both windows and unix.
        continue;
      }
      tokenize(vs, buf, " \t\r\n");

      if (atom_input) {
	if (vs.size() < 9) return false; // timvdm 18/06/2008
        id = atoi((char*)vs[0].c_str()) + 1; // see warning above
        type = vs[1];
        X = atof((char*)vs[6].c_str())/scale;
        Y = atof((char*)vs[7].c_str())/scale;
        Z = atof((char*)vs[8].c_str())/scale;

        OBAtom* a = pmol->NewAtom();
        if (*(type.c_str()) != '*')
          a->SetAtomicNum(etab.GetAtomicNum(type.c_str()));
        a->SetVector(X,Y,Z);

      } else if (bond_input) {
	if (vs.size() < 2) return false; // timvdm 18/06/2008
        // add to pmol
        if (!pmol->AddBond(atoi((char*)vs[0].c_str()) + 1, atoi((char*)vs[1].c_str()) + 1,
                           1 /* bond order not specified in Carine, use PerceiveBondOrder later */))
          {
            obErrorLog.ThrowError(__FUNCTION__, "addition of bond between " + vs[0] + " and " + vs[1] + " failed", obError);
            return false;
          }
      }
    }

    /* got sanity? */
    if ( pmol->NumBonds() != bonds ) {
      // then we read a different number of bonds than those promised.
      obErrorLog.ThrowError(__FUNCTION__, "Number of bonds read does not match the number promised", obError);
      return false;
    } else if ( pmol->NumAtoms() != atoms ) {
      obErrorLog.ThrowError(__FUNCTION__, "Number of atoms read does not match the number promised", obError);
      return false;
    }

    pmol->PerceiveBondOrders();

    pmol->EndModify();

    return true;
  }