C++ OBPairData::SetOrigin方法代码示例

  bool QEqCharges::ComputeCharges(OBMol &mol)
  {

     ///////////////////////////////////////////////////////////////////////////////
    //Some OpenBabel bookkeeping that I copied from the Gasteiger scheme
    mol.SetPartialChargesPerceived();

    OBPairData *dp = new OBPairData;
    dp->SetAttribute("PartialCharges");
    dp->SetValue("QEq");
    dp->SetOrigin(perceived);
    mol.SetData(dp);


    ///////////////////////////////////////////////////////////////////////////////
    //Read in atomic information from OpenBabel molecule and parameterize

    //Read in total number of atoms
    int i, N = mol.NumAtoms();

    Hardness = MatrixXd::Zero(N+1, N+1);
    Voltage = VectorXd::Zero(N+1);
    Electronegativity = VectorXd::Zero(N);
    VectorXd BasisSet = VectorXd::Zero(N);

    Vector3d Parameters;

    FOR_ATOMS_OF_MOL(atom, mol)
    {
       	Parameters = GetParameters(atom->GetAtomicNum(), atom->GetFormalCharge());
	i = atom->GetIdx() - 1;

	if (Parameters[0] == 0.)
        {
		stringstream msg;
		msg << "Some QEq Parameters not found!" << endl
	            << "Parameters not found for atom no. " << i+1 << endl
		    << "Atom will be ignored in the charge computation.";
		obErrorLog.ThrowError(__FUNCTION__, msg.str(), obError);
	}

	Electronegativity[i] = Parameters[0];
	Hardness(i,i) = Parameters[1];
	BasisSet[i] = Parameters[2];
    }

//.........这里部分代码省略.........
      } else if (strstr(buffer,"MOPAC") != NULL) {
        // MOPAC output
        formatName = "mopout";
        break;
      } else if (strstr(buffer,"Program PWSCF") != NULL) {
        // PWSCF
        formatName = "pwscf";
        break;
      } else if (strstr(buffer,"Welcome to Q-Chem") != NULL) {
        // Q-Chem output
        formatName = "qcout";
        break;
      } else if (strstr(buffer,"Amsterdam Density Functional") != NULL) {
        // ADF output
        // Determine the kind of ADF output
        while (ifs.getline(buffer, BUFF_SIZE)) {
          if (strstr(buffer, "|     A D F     |") != NULL) {
            formatName = "adfout";
            break;
          } else if (strstr(buffer, "|     B A N D     |") != NULL) {
            formatName = "adfband";
            break;
          } else if (strstr(buffer, "|     D F T B     |") != NULL) {
            formatName = "adfdftb";
            break;
          } else if (strstr(buffer, "DFTB Engine") != NULL) {
            // "|     D F T B     |" is no longer printed in ADF 2018
            // Hopefully, "DFTB Engine" will work fine...
            formatName = "adfdftb";
            break;
          }
        }
        break;
      } else if (strstr(buffer,"Northwest Computational Chemistry") != NULL) {
        // NWChem output
        formatName = "nwo";
        break;
      } else if (strstr(buffer,"MPQC: Massively Parallel Quantum Chemistry") != NULL) {
        // MPQC output
        formatName = "mpqc";
        break;
      } else if (strstr(buffer,"PROGRAM SYSTEM MOLPRO") != NULL) {
        // MOLPRO output
        formatName = "mpo";
        break;
      } else if ((strstr(buffer,"Schrodinger, Inc.") != NULL) &&
                 (strstr(buffer,"Jaguar") != NULL)) {
        // Jaguar
        formatName = "jout";
        break;
      } else if (strstr(buffer, "ABINIT") != NULL) {
        // Abinit
        formatName = "abinit";
        break;
      } else if (strstr(buffer, "ACES2") != NULL) {
        // ACESII
        formatName = "acesout";
        break;
      } else if (strstr(buffer, "CRYSTAL06") != NULL ||
                 strstr(buffer, "CRYSTAL09") != NULL) {
        // CRYSTAL09
        formatName = "c09out";
        break;
      } else if (strstr(buffer, "* O   R   C   A *") != NULL) {
        // ORCA
        formatName = "orca";
        break;
      } else if (strstr(buffer, "WELCOME TO SIESTA") != NULL) {
        // SIESTA
        formatName = "siesta";
        break;
      }
    }

    // if we assigned something above, let's try to find it
    if (formatName.length())
      pFormat = pConv->FindFormat(formatName);

    if (pFormat) {
      ifs.seekg (0, ios::beg); // reset the stream to the beginning
      ifs.clear();
      bool success = pFormat->ReadMolecule(pOb, pConv);

      // Tag the molecule with the format (e.g., if a program wants to know the kind of "out" or "log" file)
      // We have to do this *after* ReadMolecule returns, or the data might be cleared
      if (pOb) {
        OBPairData *dp = new OBPairData;
        dp->SetAttribute("File Format");
        dp->SetValue(formatName);
        dp->SetOrigin(fileformatInput);
        pOb->SetData(dp);
      }

      return success;
    }

    obErrorLog.ThrowError(__FUNCTION__,
                          "Problems reading an output file: Could not determine the format of this file. Please report it to the openbabel-discuss @ lists.sourceforge.net mailing list.", obError);
    return(false); // we couldn't figure out the format
  }

//.........这里部分代码省略.........
    if (!getline(ifs,ln))
      return false; //creator
    if (!getline(ifs, ln))
      return false; //comment
    // Originally the comment was added to the reaction via:
    //     pmol->SetComment(Trim(ln));

    if (!getline(ifs, ln))
      return false; // num reactants, products, and optionally agents

    unsigned int nReactants = 0, nProducts = 0, nAgents = 0;
    bool ok = ParseComponent(ln.c_str() + 0, &nReactants);
    if (!ok)
      return false;
    ok = ParseComponent(ln.c_str() + 3, &nProducts);
    if (!ok)
      return false;
    if (ln[6] != '\0') { // optional agents
      ok = ParseComponent(ln.c_str() + 6, &nAgents);
      if (!ok)
        return false;
    }

    if(nReactants + nProducts + nAgents)
    {
      //Read the first $MOL. The others are read at the end of the previous MOL
      if(!getline(ifs, ln))
        return false;
      if(Trim(ln).find("$MOL")==string::npos)
        return false;
    }

    OBReactionFacade rxnfacade(pmol);

    // Note: If we supported it, we could read each of the rxn components directly
    // into the returned OBMol instead of having to do a copy. Unfortunately,
    // this isn't possible at the moment (MOL format will need some work first).
    // Here is some example code to do it:
    //
    //unsigned int old_numatoms = 0;
    //unsigned int compid = 1;
    //for (int i = 0; i<nReactants; i++)
    //{
    //  //Read a MOL file	using the same OBConversion object but with a different format
    //  if (!pMolFormat->ReadMolecule(pmol, pConv))
    //    obErrorLog.ThrowError(__FUNCTION__, "Failed to read a reactant", obWarning);
    //  unsigned int numatoms = pmol->NumAtoms();
    //  for (unsigned int idx = old_numatoms + 1; idx <= numatoms; ++idx) {
    //    OBAtom* atom = pmol->GetAtom(idx);
    //    rxnfacade.SetRole(atom, REACTANT);
    //    rxnfacade.SetComponentId(atom, compid);
    //  }
    //  old_numatoms = numatoms;
    //  compid++;
    //}

    const char* type[3] = {"a reactant", "a product", "an agent"};
    OBReactionRole role;
    unsigned int num_components;
    for(unsigned int N=0; N<3; N++) {
      switch(N) {
      case 0:
        role = REACTANT;
        num_components = nReactants;
        break;
      case 1:
        role = PRODUCT;
        num_components = nProducts;
        break;
      case 2:
        role = AGENT;
        num_components = nAgents;
        break;
      }
      for (int i=0; i<num_components; i++)
      {
        //Read a MOL file	using the same OBConversion object but with a different format
        OBMol mol;
        if (!pMolFormat->ReadMolecule(&mol, pConv)) {
          std::string error = "Failed to read ";
          error += type[N];
          obErrorLog.ThrowError(__FUNCTION__, error, obWarning);
          continue;
        }
        if (mol.NumAtoms() == 0) {
          OBAtom* dummy = mol.NewAtom(); // Treat the empty OBMol as having a single dummy atom
          OBPairData *pd = new OBPairData();
          pd->SetAttribute("rxndummy");
          pd->SetValue("");
          pd->SetOrigin(fileformatInput);
          dummy->SetData(pd);
        }

        rxnfacade.AddComponent(&mol, role);
      }
    }

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

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

  //All molecules
  if(pmol->GetDimension()!= _refMol.GetDimension())
  {
    stringstream ss;
    ss << "The molecule" << pmol->GetTitle()
       << " does not have the same dimensions as the reference molecule "
       << _refMol.GetTitle() << " and is ignored.";
       obErrorLog.ThrowError(__FUNCTION__, ss.str().c_str(), obError);
    return false;
  }

  if(_pOpIsoM) //Using -s option
  {
    //Ignore mol if it does not pass -s option
    if(!_pOpIsoM->Do(pmol, "", pmap, pConv)) // "" means will use existing parameters
      return false;

    // Get the atoms equivalent to those in ref molecule
    vector<int> ats = _pOpIsoM->GetMatchAtoms();

    // Make a vector of their coordinates and get the centroid
    vector<vector3> vec;
    vector3 centroid;
    for(vector<int>::iterator iter=ats.begin(); iter!=ats.end(); ++iter) {
      vector3 v = pmol->GetAtom(*iter)->GetVector();
      centroid += v;
      vec.push_back(v);
    }
    centroid /= vec.size();

    // Do the alignment
    _align.SetTarget(vec);
    if(!_align.Align())
      return false;

    // Get the centroid of the reference atoms
    vector3 ref_centroid;
    for(vector<vector3>::iterator iter=_refvec.begin(); iter!=_refvec.end(); ++iter)
      ref_centroid += *iter;
    ref_centroid /= _refvec.size();

    //subtract the centroid, rotate the target molecule, then add the centroid
    matrix3x3 rotmatrix = _align.GetRotMatrix();
    for (unsigned int i = 1; i <= pmol->NumAtoms(); ++i)
    {
      vector3 tmpvec = pmol->GetAtom(i)->GetVector();
      tmpvec -= centroid;
      tmpvec *= rotmatrix; //apply the rotation
      tmpvec += ref_centroid;
      pmol->GetAtom(i)->SetVector(tmpvec);
    }
  }
  else //Not using -s option)
  {
    _align.SetTargetMol(*pmol);
    if(!_align.Align())
      return false;
    _align.UpdateCoords(pmol);
  }

  //Save rmsd as a property
  OBPairData* dp = new OBPairData;
  dp->SetAttribute("rmsd");
  double val = _align.GetRMSD();
  if(val<1e-12)
    val = 0.0;
  dp->SetValue(toString(val));
  dp->SetOrigin(local);
  pmol->SetData(dp);

  return true;
}

//.........这里部分代码省略.........
      {
        string title(GetTitle());
        title += itr->second;
        SetTitle(title.c_str());
      }

/*    itr = pOptions->find("addformula"); //Appends tab + formula to title
    if(itr!=pOptions->end())
      {
        string title(GetTitle());
        title += '\t' + GetSpacedFormula(1,"");//actually unspaced
        SetTitle(title.c_str());
      }
*/
    //Add an extra property to the molecule.
    //Parameter has atrribute and value separated by a space
    itr = pOptions->find("property");
    if(itr!=pOptions->end())
      {
        string txt(itr->second);
        string::size_type pos = txt.find(' ');
        if(pos==string::npos)
          {
            obErrorLog.ThrowError(__FUNCTION__, "Missing property value", obError);
            ret=false;
          }
        else
          {
            string attr(txt.substr(0,pos)), val(txt.substr(pos+1));
            //Update value if it already exists
            OBPairData* dp = dynamic_cast<OBPairData*>(GetData(attr));
            if(dp) {
              dp->SetValue(val);
              dp->SetOrigin(userInput);
            }
            else {
              // Pair did not exist; make new one
              dp = new OBPairData;
              dp->SetAttribute(attr);
              dp->SetValue(val);
              dp->SetOrigin(userInput);
              SetData(dp);
            }
          }
      }

    itr = pOptions->find("add");  //adds new properties from descriptors in list
    if(itr!=pOptions->end())
      OBDescriptor::AddProperties(this, itr->second);

    itr = pOptions->find("delete"); //deletes the specified properties
    if(itr!=pOptions->end())
      OBDescriptor::DeleteProperties(this, itr->second);

    itr = pOptions->find("append"); //Appends values of descriptors or properties to title
    if(itr!=pOptions->end())
      {
        string title(GetTitle());
        title += OBDescriptor::GetValues(this, itr->second);
        if(ispunct(title[0]))
          title[0]=' ';//a leading punct char is used only as a separator, not at start
        SetTitle(Trim(title).c_str());
      }

// Reading Gaussian output has been tested for G98 and G03 to some degree
// If you have problems (or examples of older output), please contact
// the [email protected] mailing list and/or post a bug
bool GaussianOutputFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
{
    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if(pmol==NULL)
        return false;

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();
    OBMol &mol = *pmol;
    const char* title = pConv->GetTitle();

    char buffer[BUFF_SIZE];
    string str,str1,str2,thermo_method;
    double x,y,z;
    OBAtom *atom;
    vector<string> vs,vs2;
    int total_charge = 0;
    unsigned int spin_multiplicity = 1;
    bool hasPartialCharges = false;
    string chargeModel; // descriptor for charges (e.g. "Mulliken")

    // Variable for G2/G3/G4 etc. calculations
    double ezpe,Hcorr,Gcorr,E0,CV;
    bool ezpe_set=false,Hcorr_set=false,Gcorr_set=false,E0_set=false,CV_set=false;
    double temperature = 0; /* Kelvin */
    std::vector<double> Scomponents;
    // Electrostatic potential
    OBFreeGrid *esp = NULL;

    // coordinates of all steps
    // Set conformers to all coordinates we adopted
    std::vector<double*> vconf; // index of all frames/conformers
    std::vector<double> coordinates; // coordinates in each frame
    int natoms = 0; // number of atoms -- ensure we don't go to a new job with a different molecule

    // OBConformerData stores information about multiple steps
    // we can change attribute later if needed (e.g., IRC)
    OBConformerData *confData = new OBConformerData();
    confData->SetOrigin(fileformatInput);
    std::vector<unsigned short> confDimensions = confData->GetDimension(); // to be fair, set these all to 3D
    std::vector<double>         confEnergies   = confData->GetEnergies();
    std::vector< std::vector< vector3 > > confForces = confData->GetForces();

    //Vibrational data
    std::vector< std::vector< vector3 > > Lx;
    std::vector<double> Frequencies, Intensities;
    //Rotational data
    std::vector<double> RotConsts(3);
    int RotSymNum=1;
    OBRotationData::RType RotorType = OBRotationData::UNKNOWN;

    // Translation vectors (if present)
    vector3 translationVectors[3];
    int numTranslationVectors = 0;

    //Electronic Excitation data
    std::vector<double> Forces, Wavelengths, EDipole,
        RotatoryStrengthsVelocity, RotatoryStrengthsLength;

    // Orbital data
    std::vector<double> orbitals;
    std::vector<std::string> symmetries;
    int aHOMO, bHOMO, betaStart;
    aHOMO = bHOMO = betaStart = -1;

    int i=0;
    bool no_symmetry=false;
    char coords_type[25];

    //Prescan file to find second instance of "orientation:"
    //This will be the kind of coords used in the chk/fchk file
    //Unless the "nosym" keyword has been requested
    while (ifs.getline(buffer,BUFF_SIZE))
    {
        if (strstr(buffer,"Symmetry turned off by external request.") != NULL)
        {
            // The "nosym" keyword has been requested
            no_symmetry = true;
        }
        if (strstr(buffer,"orientation:") !=NULL)
        {
            i++;
            tokenize (vs, buffer);
            // gotta check what types of orientation are present
            strncpy (coords_type, vs[0].c_str(), 24);
            strcat (coords_type, " orientation:");
        }
        if ((no_symmetry && i==1) || i==2)
            break;
    }
    // Reset end-of-file pointers etc.
    ifs.clear();
    ifs.seekg(0);  //rewind

    mol.BeginModify();
    while (ifs.getline(buffer,BUFF_SIZE))
    {
//.........这里部分代码省略.........

  // Reading Gaussian output has been tested for G98 and G03 to some degree
  // If you have problems (or examples of older output), please contact
  // the [email protected] mailing list and/or post a bug
  bool GaussianOutputFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
  {
    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if(pmol==NULL)
      return false;

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();
    OBMol &mol = *pmol;
    const char* title = pConv->GetTitle();

    char buffer[BUFF_SIZE];
    string str,str1;
    double x,y,z;
    OBAtom *atom;
    vector<string> vs;
    int charge = 0;
    unsigned int spin = 1;
    bool hasPartialCharges = false;
    string chargeModel; // descriptor for charges (e.g. "Mulliken")

    // coordinates of all steps
    // Set conformers to all coordinates we adopted
    std::vector<double*> vconf; // index of all frames/conformers
    std::vector<double> coordinates; // coordinates in each frame
    int natoms = 0; // number of atoms -- ensure we don't go to a new job with a different molecule

    // OBConformerData stores information about multiple steps
    // we can change attribute later if needed (e.g., IRC)
    OBConformerData *confData = new OBConformerData();
    confData->SetOrigin(fileformatInput);
    std::vector<unsigned short> confDimensions = confData->GetDimension(); // to be fair, set these all to 3D
    std::vector<double>         confEnergies   = confData->GetEnergies();
    std::vector< std::vector< vector3 > > confForces = confData->GetForces();

    //Vibrational data
    std::vector< std::vector< vector3 > > Lx;
    std::vector<double> Frequencies, Intensities;
    //Rotational data
    std::vector<double> RotConsts(3);
    int RotSymNum=1;
    OBRotationData::RType RotorType;

    // Translation vectors (if present)
    vector3 translationVectors[3];
    int numTranslationVectors = 0;

    //Electronic Excitation data
    std::vector<double> Forces, Wavelengths, EDipole,
      RotatoryStrengthsVelocity, RotatoryStrengthsLength;

    // Orbital data
    std::vector<double> orbitals;
    std::vector<std::string> symmetries;
    int aHOMO, bHOMO, betaStart;

    //Put some metadata into OBCommentData
    string comment("Gaussian ");
    ifs.getline(buffer,BUFF_SIZE);
    if(*buffer)
    {
      comment += strchr(buffer,'=')+2;
      comment += "";
      for(unsigned i=0; i<115, ifs; ++i)
      {
        ifs.getline(buffer,BUFF_SIZE);
        if(buffer[1]=='#')
        {
          //the line describing the method
          comment += buffer;
          OBCommentData *cd = new OBCommentData;
          cd->SetData(comment);
          cd->SetOrigin(fileformatInput);
          mol.SetData(cd);
          break;
        }
      }
    }

    int i=0;
    bool no_symmetry=false;
    char coords_type[25];

    //Prescan file to find second instance of "orientation:"
    //This will be the kind of coords used in the chk/fchk file
    //Unless the "nosym" keyword has been requested
    while (ifs.getline(buffer,BUFF_SIZE))
      {
        if (strstr(buffer,"Symmetry turned off by external request.") != NULL)
          {
            // The "nosym" keyword has been requested
            no_symmetry = true;
          }
        if (strstr(buffer,"orientation:") !=NULL)
          {
            i++;
            tokenize (vs, buffer);
//.........这里部分代码省略.........

  bool MOL2Format::ReadMolecule(OBBase* pOb, OBConversion* pConv)
  {

    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if(pmol==NULL)
      return false;

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

    //Old code follows...
    bool foundAtomLine = false;
    char buffer[BUFF_SIZE];
    char *comment = NULL;
    string str,str1;
    vector<string> vstr;
    int len;

    mol.BeginModify();

    for (;;)
      {
        if (!ifs.getline(buffer,BUFF_SIZE))
          return(false);
        if (EQn(buffer,"@<TRIPOS>MOLECULE",17))
          break;
      }

    // OK, just read MOLECULE line
    int lcount;
    int natoms,nbonds;
    bool hasPartialCharges = true;
    for (lcount=0;;lcount++)
      {
        if (!ifs.getline(buffer,BUFF_SIZE))
          return(false);
        if (EQn(buffer,"@<TRIPOS>ATOM",13))
          {
            foundAtomLine = true;
            break;
          }

        if (lcount == 0)
          {
            tokenize(vstr,buffer);
            if (!vstr.empty())
              mol.SetTitle(buffer);
          }
        else if (lcount == 1)
          sscanf(buffer,"%d%d",&natoms,&nbonds);
        else if (lcount == 3) // charge descriptions
          {
            // Annotate origin of partial charges
            OBPairData *dp = new OBPairData;
            dp->SetAttribute("PartialCharges");
            dp->SetValue(buffer);
            dp->SetOrigin(fileformatInput);
            mol.SetData(dp);

            if (strncasecmp(buffer, "NO_CHARGES", 10) == 0)
              hasPartialCharges = false;
          }
        else if (lcount == 4) //energy (?)
          {
            tokenize(vstr,buffer);
            if (!vstr.empty() && vstr.size() == 3)
              if (vstr[0] == "Energy")
                mol.SetEnergy(atof(vstr[2].c_str()));
          }
        else if (lcount == 5) //comment
          {
            if ( buffer[0] )
              {
                len = (int) strlen(buffer)+1;
                //! @todo allow better multi-line comments
                // which don't allow ill-formed data to consume memory
                // Thanks to Andrew Dalke for the pointer
                if (comment != NULL)
                  delete [] comment;
                comment = new char [len];
                memcpy(comment,buffer,len);
              }
          }
      }

    if (!foundAtomLine)
      {
        mol.EndModify();
        mol.Clear();
        obErrorLog.ThrowError(__FUNCTION__, "Unable to read Mol2 format file. No atoms found.", obWarning);
        return(false);
      }

    mol.ReserveAtoms(natoms);

    int i;
    vector3 v;
    OBAtom atom;
    double x,y,z,pcharge;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            mol.SetTitle(defaultTitle);
          }
      }
    else
      return(false);

    mol.BeginModify();
    mol.ReserveAtoms(natoms);
    connections.resize(natoms+1);

    /***********************************************************************/

    // Get Type Bonds, BondOrder, X, Y, Z

    double x,y,z;
    vector3 v;
    char temp_type[10];
    int i,j;
    double charge;
    OBAtom atom;

    ttab.SetFromType("MMD");
    for (i = 1; i <= natoms; i++)
      {
        if (!ifs.getline(buffer,BUFF_SIZE))
          break;

        int end[6], order[6];

        sscanf(buffer,"%9s%d%d%d%d%d%d%d%d%d%d%d%d%lf%lf%lf",
               temp_type,&end[0],&order[0],&end[1],&order[1],&end[2],&order[2],
               &end[3], &order[3], &end[4], &order[4], &end[5], &order[5],
               &x, &y, &z);

        pair<int,int> tmp;
        for ( j = 0 ; j <=5 ; j++ )
          {
            if ( end[j] > 0  && end[j] > i)
              {
                tmp.first = end[j];
                tmp.second = order[j];
                connections[i].push_back(tmp);
              }
          }

        v.SetX(x);
        v.SetY(y);
        v.SetZ(z);
        atom.SetVector(v);

        string str = temp_type,str1;
        ttab.SetToType("ATN");
        ttab.Translate(str1,str);
        atom.SetAtomicNum(atoi(str1.c_str()));
        ttab.SetToType("INT");
        ttab.Translate(str1,str);
        atom.SetType(str1);

        // stuff for optional fields

        buffer[109]='\0';
        sscanf(&buffer[101],"%lf", &charge);
        atom.SetPartialCharge(charge);
        mol.AddAtom(atom);
      }

    for (i = 1; i <= natoms; i++)
      for (j = 0; j < (signed)connections[i].size(); j++)
        mol.AddBond(i, connections[i][j].first, connections[i][j].second);

    mol.EndModify();
    mol.SetPartialChargesPerceived();
    // Annotate origin of partial charges
    OBPairData *dp = new OBPairData;
    dp->SetAttribute("PartialCharges");
    dp->SetValue("MACROMODEL");
    dp->SetOrigin(fileformatInput);
    mol.SetData(dp);

    OBBond *bond;
    vector<OBBond*>::iterator bi;
    for (bond = mol.BeginBond(bi);bond;bond = mol.NextBond(bi))
      if (bond->GetBondOrder() == 5 && !bond->IsInRing())
        bond->SetBondOrder(1);

    if ( natoms != (signed)mol.NumAtoms() )
      return(false);

    // clean out remaining blank lines
    std::streampos ipos;
    do
    {
      ipos = ifs.tellg();
      ifs.getline(buffer,BUFF_SIZE);
    }
    while(strlen(buffer) == 0 && !ifs.eof() );
    ifs.seekg(ipos);

    return(true);
  }