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

void testTetrahedralStereo1()
{
  cout << "testTetrahedralStereo1()" << endl;
  // read a smiles string
  OBMol mol;
  OBConversion conv;
  OB_REQUIRE( conv.SetInFormat("smi") );
  cout << "smiles: C[[email protected]](O)N" << endl;
  OB_REQUIRE( conv.ReadString(&mol, "C[[email protected]](O)N") );

  // get the stereo data
  OB_REQUIRE( mol.HasData(OBGenericDataType::StereoData) );
  std::vector<OBGenericData *> stereoData = mol.GetAllData(OBGenericDataType::StereoData);
  OB_REQUIRE( stereoData.size() == 1 );

  // convert to tetrahedral data
  OB_REQUIRE( ((OBStereoBase*)stereoData[0])->GetType() == OBStereo::Tetrahedral );
  OBTetrahedralStereo *ts = dynamic_cast<OBTetrahedralStereo*>(stereoData[0]);
  OB_REQUIRE( ts );

  // print the configuration
  cout << *ts << endl;

  // construct a valid configuration here
  //
  // C[[email protected]](O)N
  // 0 1 2  3 4  <- ids
  //
  OBTetrahedralStereo::Config cfg(1, 0, OBStereo::MakeRefs(4, 3, 2), OBStereo::Clockwise);

  // compare stereochemistry
  OB_REQUIRE( ts->GetConfig() == cfg );

  cout << endl;
}

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

bool WriteReport(ostream &ofs,OBMol &mol)
{
  char buffer[BUFF_SIZE];
  ofs << "FILENAME: " << mol.GetTitle() << endl;
  ofs << "MASS: ";
  sprintf(buffer, "%5.4f", mol.GetMolWt());
  ofs << buffer << endl;
  ofs << "EXACT MASS: ";
  sprintf(buffer, "%5.7f", mol.GetExactMass());
  ofs << buffer << endl;
  if (mol.GetTotalCharge() != 0)
    {
      ofs << "TOTAL CHARGE: ";
      sprintf(buffer, "%d", mol.GetTotalCharge());
      ofs << buffer << endl;
    }
  if (mol.GetTotalSpinMultiplicity() != 1)
    {
      ofs << "TOTAL SPIN: ";
      sprintf(buffer, "%d", mol.GetTotalSpinMultiplicity());
      ofs << buffer << endl;
    }
  ofs << "INTERATOMIC DISTANCES" << endl;
  WriteDistanceMatrix(ofs, mol);
  ofs << endl << endl << "ATOMIC CHARGES" << endl;
  WriteCharges(ofs, mol);
  ofs << endl << endl << "BOND ANGLES" << endl;
  WriteAngles(ofs, mol);
  ofs << endl << endl << "TORSION ANGLES" << endl;
  WriteTorsions(ofs, mol);
  if (mol.IsChiral())
    {
      ofs << endl << endl << "CHIRAL ATOMS" << endl;
      WriteChiral(ofs, mol);
    }
  if (mol.HasData(obCommentData)) {
    ofs << endl << endl << "COMMENTS" << endl;
    OBCommentData *cd = (OBCommentData*)mol.GetData(obCommentData);
    ofs << cd->GetData() << endl;
  }
  ofs << endl << endl;
  return(true);
}

void genericSmilesCanonicalTest(const std::string &smiles)
{
  cout << "Testing generic smiles <-> canonical smiles" << endl;
  // read a smiles string
  OBMol mol;
  OBConversion conv;
  OB_REQUIRE( conv.SetInFormat("smi") );
  OB_REQUIRE( conv.SetOutFormat("can") );
  cout << "smiles: " << smiles << endl;
  // read a smiles string
  OB_REQUIRE( conv.ReadString(&mol, smiles) );

  // store the stereo data for the smiles string using unique symmetry ids
  std::vector<OBTetrahedralStereo::Config> tetrahedral1;
  std::vector<OBCisTransStereo::Config> cistrans1;
  std::vector<OBSquarePlanarStereo::Config> squareplanar1;

  // get the stereo data
  OB_ASSERT( mol.HasData(OBGenericDataType::StereoData) );
  std::vector<OBGenericData *> stereoData = mol.GetAllData(OBGenericDataType::StereoData);

  std::vector<unsigned int> canlbls;
  std::vector<unsigned int> symclasses;
  OBGraphSym gs1(&mol);
  gs1.GetSymmetry(symclasses);
  CanonicalLabels(&mol, symclasses, canlbls);
  cout << "mol.NumAtoms = " << mol.NumAtoms() << endl;
  for (std::vector<OBGenericData*>::iterator data = stereoData.begin(); data != stereoData.end(); ++data) {
    if (((OBStereoBase*)*data)->GetType() == OBStereo::Tetrahedral) {
      // convert to tetrahedral data
      OBTetrahedralStereo *ts = dynamic_cast<OBTetrahedralStereo*>(*data);
      OB_REQUIRE( ts );
      OB_ASSERT( ts->IsValid() );
      if (!ts->IsValid())
        continue;

      OBTetrahedralStereo::Config config = ts->GetConfig();
      // convert atom ids to symmetry ids
     if (mol.GetAtomById(config.center))
        config.center = canlbls.at( mol.GetAtomById(config.center)->GetIdx() - 1 );
      if (mol.GetAtomById(config.from))
        config.from = canlbls.at( mol.GetAtomById(config.from)->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[0]))
        config.refs[0] = canlbls.at( mol.GetAtomById(config.refs[0])->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[1]))
        config.refs[1] = canlbls.at( mol.GetAtomById(config.refs[1])->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[2]))
        config.refs[2] = canlbls.at( mol.GetAtomById(config.refs[2])->GetIdx() - 1 );
      cout << "Config with symmetry ids: " << config << endl;
      tetrahedral1.push_back(config);
    } else
    if (((OBStereoBase*)*data)->GetType() == OBStereo::CisTrans) {
      // convert to tetrahedral data
      OBCisTransStereo *ct = dynamic_cast<OBCisTransStereo*>(*data);
      OB_REQUIRE( ct );
      OB_ASSERT( ct->IsValid() );

      OBCisTransStereo::Config config = ct->GetConfig();
      // convert atom ids to symmetry ids
      config.begin = canlbls.at( mol.GetAtomById(config.begin)->GetIdx() - 1 );
      config.end = canlbls.at( mol.GetAtomById(config.end)->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[0]))
        config.refs[0] = canlbls.at( mol.GetAtomById(config.refs[0])->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[1]))
        config.refs[1] = canlbls.at( mol.GetAtomById(config.refs[1])->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[2]))
        config.refs[2] = canlbls.at( mol.GetAtomById(config.refs[2])->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[3]))
        config.refs[3] = canlbls.at( mol.GetAtomById(config.refs[3])->GetIdx() - 1 );
      cout << "Config with symmetry ids: " << config << endl;
      cistrans1.push_back(config);
    } else
    if (((OBStereoBase*)*data)->GetType() == OBStereo::SquarePlanar) {
      // convert to tetrahedral data
      OBSquarePlanarStereo *sp = dynamic_cast<OBSquarePlanarStereo*>(*data);
      OB_REQUIRE( sp );
      OB_ASSERT( sp->IsValid() );
      if (!sp->IsValid())
        continue;

      OBSquarePlanarStereo::Config config = sp->GetConfig();
      // convert atom ids to symmetry ids
     if (mol.GetAtomById(config.center))
        config.center = canlbls.at( mol.GetAtomById(config.center)->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[0]))
        config.refs[0] = canlbls.at( mol.GetAtomById(config.refs[0])->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[1]))
        config.refs[1] = canlbls.at( mol.GetAtomById(config.refs[1])->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[2]))
        config.refs[2] = canlbls.at( mol.GetAtomById(config.refs[2])->GetIdx() - 1 );
      if (mol.GetAtomById(config.refs[3]))
        config.refs[3] = canlbls.at( mol.GetAtomById(config.refs[3])->GetIdx() - 1 );
      cout << "Config with symmetry ids: " << config << endl;
      squareplanar1.push_back(config);
    }


  }

  // write to can smiles
//.........这里部分代码省略.........

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

  //! \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 OBMoleculeFormat::ReadChemObjectImpl(OBConversion* pConv, OBFormat* pFormat)
  {
    std::istream *ifs = pConv->GetInStream();
    if (!ifs || !ifs->good())
      return false;

    OBMol* pmol = new OBMol;

    std::string auditMsg = "OpenBabel::Read molecule ";
    std::string description(pFormat->Description());
    auditMsg += description.substr(0,description.find('\n'));
    obErrorLog.ThrowError(__FUNCTION__,
                          auditMsg,
                          obAuditMsg);

    if(pConv->IsOption("C",OBConversion::GENOPTIONS))
      return DeferMolOutput(pmol, pConv, pFormat);

    bool ret=true;
   if(pConv->IsOption("separate",OBConversion::GENOPTIONS))
   {
     //On first call, separate molecule and put fragments in MolArray.
     //On subsequent calls, remove a fragment from MolArray and send it for writing
     //Done this way so that each fragment can be written to its own file (with -m option)
     if(!StoredMolsReady)
     {
       while(ret) //do all the molecules in the file
       {
         ret = pFormat->ReadMolecule(pmol,pConv);

         if(ret && (pmol->NumAtoms() > 0 || (pFormat->Flags()&ZEROATOMSOK)))
         {
           vector<OBMol> SepArray = pmol->Separate(); //use un-transformed molecule
           //Add an appropriate title to each fragment
           if(SepArray.size()>1)
             for (unsigned int i=0; i<SepArray.size(); ++i)
             {
               stringstream ss;
               ss << pmol->GetTitle() << '#' << i+1;
               string title = ss.str();
               SepArray[i].SetTitle(title);
             }
           else
              SepArray[0].SetTitle(pmol->GetTitle());

           copy(SepArray.begin(),SepArray.end(),back_inserter(MolArray));
         }
       }
       reverse(MolArray.begin(),MolArray.end());
       StoredMolsReady = true;
       //Clear the flags of the input stream(which may have found eof) to ensure will
       //try to read anothe molecule and allow the stored ones to be sent for output.
       pConv->GetInStream()->clear();
     }

     if(MolArray.empty()) //normal end of fragments
       ret =false;
     else
     {
       // Copying is needed because the OBMol passed to AddChemObject will be deleted.
       // The OBMol in the vector is deleted here.
       OBMol* pMolCopy = new OBMol( MolArray.back());
       MolArray.pop_back();
       ret = pConv->AddChemObject(
           pMolCopy->DoTransformations(pConv->GetOptions(OBConversion::GENOPTIONS), pConv))!=0;
     }
     if(!ret)
       StoredMolsReady = false;

     delete pmol;
     return ret;
   }

    ret=pFormat->ReadMolecule(pmol,pConv);

    OBMol* ptmol = NULL;
    //Molecule is valid if it has some atoms
    //or the format allows zero-atom molecules and it has a title or properties
    if(ret && (pmol->NumAtoms() > 0 
      || (pFormat->Flags()&ZEROATOMSOK && (*pmol->GetTitle() || pmol->HasData(1)))))
    {
      ptmol = static_cast<OBMol*>(pmol->DoTransformations(pConv->GetOptions(OBConversion::GENOPTIONS),pConv));
      if(ptmol && (pConv->IsOption("j",OBConversion::GENOPTIONS)
                || pConv->IsOption("join",OBConversion::GENOPTIONS)))
      {
        //With j option, accumulate all mols in one stored in this class
        if(pConv->IsFirstInput())
          _jmol = new OBMol;
        pConv->AddChemObject(_jmol);
        //will be discarded in WriteChemObjectImpl until the last input mol. This complication
        //is needed to allow joined molecules to be from different files. pOb1 in AddChem Object
        //is zeroed at the end of a file and _jmol is in danger of not being output.
        *_jmol += *ptmol;
        delete ptmol;
        return true;
      }
    }
    else
      delete pmol;

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

bool OpFillUC::Do(OBBase* pOb, const char* OptionText, OpMap* pOptions, OBConversion* pConv)
{
  OBMol* pmol = dynamic_cast<OBMol*>(pOb);
  if(!pmol)
    return false;
  
  if (!(pmol->HasData(OBGenericDataType::UnitCell)))
  {
    obErrorLog.ThrowError(__FUNCTION__, "Cannot fill unit cell without a unit cell !" , obWarning);
    return false;
  }
  OBUnitCell *pUC = (OBUnitCell*)pmol->GetData(OBGenericDataType::UnitCell);
  const SpaceGroup* pSG = pUC->GetSpaceGroup();
  if (pSG == NULL)
  {
    obErrorLog.ThrowError(__FUNCTION__, "Cannot fill unit cell without spacegroup information !" , obWarning);
    return false;
  }
  // Now loop over all symmetry operations, and generate symmetric atoms one at a time
  // Avoid creating overlapping atoms (duplicate), and bring back atoms within the unit cell
  // using two options:
  // "--fillUC strict": keep only atoms that are strictly inside the unit cell 
  //                    (fractionnal coordinates 0<= <1)
  // "--fillUC keepconnect": generate symmetrics of the molecule, and translate 
  //                         it back in the unit cell if necessary

  std::map<OBAtom*,std::vector<vector3> > vatoms;// key: original atoms, value=all generated symmetrics
  FOR_ATOMS_OF_MOL(atom, *pmol)
      vatoms[&(*atom)]=std::vector<vector3>();
  
  for(std::map<OBAtom*,std::vector<vector3> >:: iterator atom=vatoms.begin();
      atom!=vatoms.end();++atom){
    vector3 orig = atom->first->GetVector();
    orig = pUC->CartesianToFractional(orig);// To fractionnal coordinates
    
    // Loop over symmetry operators
    transform3dIterator ti;
    const transform3d *t = pSG->BeginTransform(ti);
    while(t){
      atom->second.push_back ( (transform3d)(*t) * orig);
      t = pSG->NextTransform(ti);
    }
  }
  if(0==strncasecmp(OptionText, "keepconnect", 11)){
    // First, bring back all symmetrical molecules back in the UC
    for(unsigned int i=0;i<vatoms.begin()->second.size();++i){
      vector3 ccoord(0,0,0);//geometrical center
      for(std::map<OBAtom*,std::vector<vector3> >:: iterator atom=vatoms.begin();
        atom!=vatoms.end();++atom){
        ccoord+=atom->second[i];
      }
      ccoord/=vatoms.size();
      ccoord=transformedFractionalCoordinate2(ccoord)-ccoord;
      for(std::map<OBAtom*,std::vector<vector3> >:: iterator atom=vatoms.begin();
        atom!=vatoms.end();++atom){
        atom->second[i]+=ccoord;
      }
    }
    // Now add atoms that are not duplicates
    for(std::map<OBAtom*,std::vector<vector3> >:: iterator atom=vatoms.begin();
        atom!=vatoms.end();++atom){
      for(unsigned int i=1;i<atom->second.size();++i){
        bool foundDuplicate = false;
        for(unsigned int j=0;j<i;++j){
          if(atom->second[i].distSq(atom->second[j])<1e-4){
            foundDuplicate=true;
            break;
          }
        }
        if(!foundDuplicate){
          OBAtom *newAtom = pmol->NewAtom();
          newAtom->Duplicate(atom->first);
          newAtom->SetVector( pUC->FractionalToCartesian(atom->second[i]));
        }
      }
    }
  }
  else{
    if(0!=strncasecmp(OptionText, "strict", 6))
      obErrorLog.ThrowError(__FUNCTION__, "fillUC: lacking \"strict\n or \"keepconnect\" option, using strict" , obWarning);
    for(std::map<OBAtom*,std::vector<vector3> >:: iterator atom=vatoms.begin();
        atom!=vatoms.end();++atom){
      // Bring back within unit cell
      for(unsigned int i=0;i<atom->second.size();++i){
        atom->second[i]=transformedFractionalCoordinate2(atom->second[i]);
      }
      for(unsigned int i=1;i<atom->second.size();++i){
        bool foundDuplicate = false;
        for(unsigned int j=0;j<i;++j){
          if(atom->second[i].distSq(atom->second[j])<1e-4){
            foundDuplicate=true;
            break;
          }
        }
        if(!foundDuplicate){
          OBAtom *newAtom = pmol->NewAtom();
          newAtom->Duplicate(atom->first);
          newAtom->SetVector( pUC->FractionalToCartesian(atom->second[i]));
        }
      }
//.........这里部分代码省略.........

extern "C" int readgrid_(float *vdata, int *nx, int *ny, int *nz,
  float *x0, float *y0, float *z0, float *xx, float *yy, float *zz,
  char *file, int fsize) {

  float dx,dy,dz;
  int gsize;
  float v, vmin, vmax, vavg;
  int navg = 0;

  OBMol mol;
  OBConversion conv;
  conv.SetInFormat("cube");
 
  std::string fname = file;
  int blank = fname.find(" ");
  //printf("%d,'%s'\n", blank, fname.substr(0,blank).c_str());
  conv.ReadFile(&mol, fname.substr(0,blank).c_str());
  //cout << mol.NumAtoms() << " atoms." << endl;
  if (mol.HasData(OBGenericDataType::GridData)) {
    std::vector<OBGenericData*> grids = mol.GetAllData(OBGenericDataType::GridData);
    OBGridData *grid = dynamic_cast<OBGridData *> (grids[0]);
    gsize = grid->GetNumberOfPoints();
    grid->GetNumberOfPoints(*nx, *ny, *nz);
    vector3 origin = grid->GetOriginVector();
    *x0 = origin[0]; *y0 = origin[1]; *z0 = origin[2];
    vector3 maxv = grid->GetMaxVector();
    *xx = maxv[0]; *yy = maxv[1]; *zz = maxv[2];
    dx=(*xx-*x0)/(*nx-1);
    dy=(*yy-*y0)/(*ny-1);
    dz=(*zz-*z0)/(*nz-1);
    printf("%s %d=%d*%d*%d\n", grid->GetAttribute().c_str(), gsize, *nx, *ny, *nz);
    printf("%f %f\n", grid->GetMinValue(), grid->GetMaxValue());
    printf("%f,%f,%f\n", *x0,*y0,*z0);
    printf("%f,%f,%f\n", *xx,*yy,*zz);
    printf("%f,%f,%f\n", dx,dy,dz);
 
    //vdata = (float *)calloc(gsize, sizeof(float));
/* this is for fortran, so reverse sense of slowest/fastest moving dimensions */
    //for (int i=0; i<*nx; ++i) {
    vmin = grid->GetValue(0,0,0);
    vmax = vmin;
    for (int i=0; i<*nz; ++i) {
      for (int j=0; j<*ny; ++j) {
        //for (int k=0; k<*nz; ++k) {
        for (int k=0; k<*nx; ++k) {
          //*vdata++ = grid->GetValue(i,j,k);
          v = grid->GetValue(k,j,i);
          if (v < 1e30 && v > -1e30) {
            if (v < vmin) vmin = v;
            if (v > vmax) vmax = v;
            ++navg;
            vavg += v;
          } else {
            v = vmax * 1000; // just a guess
          }
          *vdata++ = v;
        }
      }
    }
  }
  vavg = vavg/navg;
  printf("min/avg/max = %f/%f/%f\n", vmin, vavg, vmax);
  return 0;
}

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

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

    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);
	
        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(),
                 pUC->GetAlpha(), pUC->GetBeta(), pUC->GetGamma(),
                 pUC->GetSpaceGroup() ?
                 pUC->GetSpaceGroup()->GetHMName().c_str() : "P1");
        ofs << buffer << endl;
      }

    // before we write any records, we should check to see if any coord < -1000
    // which will cause errors in the formatting

    double minX, minY, minZ;
    minX = minY = minZ = -999.0f;
    FOR_ATOMS_OF_MOL(a, mol)
      {
        if (a->GetX() < minX)
          minX = a->GetX();
        if (a->GetY() < minY)
          minY = a->GetY();
        if (a->GetZ() < minZ)
          minZ = a->GetZ();
      }
    vector3 transV = VZero;
    if (minX < -999.0)
      transV.SetX(-1.0*minX - 900.0);
    if (minY < -999.0)
      transV.SetY(-1.0*minY - 900.0);
    if (minZ < -999.0)
      transV.SetZ(-1.0*minZ - 900.0);

    // if minX, minY, or minZ was never changed, shift will be 0.0f
    // otherwise, move enough so that smallest coord is > -999.0f
    mol.Translate(transV);

    OBAtom *atom;
    OBResidue *res;
    for (i = 1; i <= mol.NumAtoms(); i++)
      {
        atom = mol.GetAtom(i);
        strncpy(type_name, etab.GetSymbol(atom->GetAtomicNum()), sizeof(type_name));
        type_name[sizeof(type_name) - 1] = '\0';

        //two char. elements are on position 13 and 14 one char. start at 14
        if (strlen(type_name) > 1)
          type_name[1] = toupper(type_name[1]);
        else
          {
            char tmp[10];
            strncpy(tmp, type_name, 10);
            snprintf(type_name, sizeof(type_name), " %-3s", tmp);
          }

        if ( (res = atom->GetResidue()) != 0 )
          {
            het = res->IsHetAtom(atom);
            snprintf(the_res,4,"%s",(char*)res->GetName().c_str());
            snprintf(type_name,5,"%s",(char*)res->GetAtomID(atom).c_str());
            the_chain = res->GetChain();
//.........这里部分代码省略.........