C++ OBUnitCell类代码示例

  bool MDFFFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
  {
    //No surprises in this routine, cartesian coordinates are written out
    //and if at least a single atom has information about constraints,
    //then selective dynamics is used and the info is written out.
    //The atoms are ordered according to their atomic number so that the
    //output looks nice, this can be reversed by using command line flag "-xw".
    //
    OBMol* pmol = dynamic_cast<OBMol*>(pOb);
    if (pmol == NULL) {
      return false;
    }

    ostream& ofs = *pConv->GetOutStream();
    OBMol mol(*pmol);
    
    if(mol.HasData(OBGenericDataType::UnitCell))
    {
      OBUnitCell *uc = static_cast<OBUnitCell*>(mol.GetData(OBGenericDataType::UnitCell));
      uc->FillUnitCell(&mol);
    }            

    char buffer[BUFF_SIZE];
    OBUnitCell *uc = NULL;
    vector<vector3> cell;

    const char * sortAtoms     = pConv->IsOption("w", OBConversion::OUTOPTIONS);
    const char * sortAtomsList = pConv->IsOption("u", OBConversion::OUTOPTIONS);
    const char * writeIONS     = pConv->IsOption("i", OBConversion::OUTOPTIONS);    

    // Create a list of ids. These may be sorted by atomic number depending
    // on the value of keepOrder.
    
    map<int, int> indl; 
    
    if (sortAtoms != NULL) 
    {
      indl.clear();
      for(int i = 0; i < 200; i++)
        indl[i] = i;
    }
    
    if (sortAtomsList != NULL) 
    {
      indl.clear();
      vector<string> vs;
      tokenize(vs, sortAtomsList);
      for(int i = 0; i < vs.size(); i++)
        indl[etab.GetAtomicNum(vs[i].c_str())] = i;
    }
    
    map<aindx, OBAtom *> amap;
    
    FOR_ATOMS_OF_MOL(atom, mol) 
    {
      aindx ndx;
      ndx.index_param = indl[atom->GetAtomicNum()];
      ndx.atom_index = atom->GetIndex();
      amap[ndx] = &(*atom);
    }

void testSpaceGroupClean()
{
  // See https://github.com/openbabel/openbabel/pull/254
  OBConversion conv;
  OBMol mol;
  conv.SetInFormat("cif");
  conv.SetOutFormat("pdb");
  conv.ReadFile(&mol, GetFilename("test02.cif"));
  OBUnitCell* pUC = (OBUnitCell*)mol.GetData(OBGenericDataType::UnitCell);
  const SpaceGroup* pSG = pUC->GetSpaceGroup();
  SpaceGroup* sg = new SpaceGroup(*pSG);
  pSG = SpaceGroup::Find(sg);
  OB_ASSERT( pSG != NULL );

  // Check also for errors and warnings
  string summary = obErrorLog.GetMessageSummary();
  OB_ASSERT( summary.find("error") == string::npos);
  OB_ASSERT( summary.find("warning") == string::npos);

  OB_ASSERT( pSG->GetId() == 166 );

  string pdb = conv.WriteString(&mol);
  pdb = conv.WriteString(&mol);

  OB_ASSERT(pdb.find("H -3 m") != string::npos);
}

  void UnitCellExtension::unitCellParametersChanged(double a, double b, double c,
                                                    double alpha, double beta, double gamma)
  {
    if (m_molecule) {
      OBUnitCell *uc = m_molecule->OBUnitCell();
      if (uc == NULL) // huh? strange, we lost our unit cell, just return
        return;

      uc->SetData(a, b, c, alpha, beta, gamma);
      m_molecule->setOBUnitCell(uc);
      m_molecule->update();

      if (m_widget)
        m_widget->update();

    } // end if molecule
  } // end parameters changed

void testAlternativeOrigin()
{
  // See https://github.com/openbabel/openbabel/pull/1558
  OBConversion conv;
  OBMol mol;
  conv.SetInFormat("cif");
  conv.ReadFile(&mol, GetFilename("test04.cif"));
  OBUnitCell* pUC = (OBUnitCell*)mol.GetData(OBGenericDataType::UnitCell);
  const SpaceGroup* pSG = pUC->GetSpaceGroup();
  SpaceGroup* sg = new SpaceGroup(*pSG);
  pSG = SpaceGroup::Find(sg);

  string summary = obErrorLog.GetMessageSummary();
  OB_ASSERT( summary.find("warning") == string::npos);
  OB_ASSERT( pSG != NULL );
  OB_ASSERT( pSG->GetOriginAlternative() == 1);
}

  void UnitCellExtension::setMolecule(Molecule *molecule)
  {
    m_molecule = molecule;

    if (m_molecule == NULL || m_dialog == NULL)
      return; // nothing we can do

    OBUnitCell *uc = m_molecule->OBUnitCell();
    if (!uc)
      return; // no unit cell

    // We don't want to send signals while we update
    disconnect(m_dialog, SIGNAL(unitCellParametersChanged(double, double, double, double, double, double)),
               this, SLOT(unitCellParametersChanged(double, double, double, double, double, double)));

    m_dialog->aLength(uc->GetA());
    m_dialog->bLength(uc->GetB());
    m_dialog->cLength(uc->GetC());

    m_dialog->alpha(uc->GetAlpha());
    m_dialog->beta(uc->GetBeta());
    m_dialog->gamma(uc->GetGamma());

    // reconnect
    connect(m_dialog, SIGNAL(unitCellParametersChanged(double, double, double, double, double, double)),
            this, SLOT(unitCellParametersChanged(double, double, double, double, double, double)));
  }

void testDecayToP1()
{
  // See https://github.com/openbabel/openbabel/pull/261
  OBConversion conv;
  OBMol mol;
  conv.SetInFormat("cif");
  conv.ReadFile(&mol, GetFilename("test03.cif"));
  OBUnitCell* pUC = (OBUnitCell*)mol.GetData(OBGenericDataType::UnitCell);
  const SpaceGroup* pSG = pUC->GetSpaceGroup();
  SpaceGroup* sg = new SpaceGroup(*pSG);
  pSG = SpaceGroup::Find(sg);
  OB_ASSERT( pSG != NULL );

  // Check also for errors and warnings
  string summary = obErrorLog.GetMessageSummary();
  OB_ASSERT( summary.find("2 warnings") != string::npos);

  OB_ASSERT( pSG->GetId() == 1 );
}

  bool FreeFormFractionalFormat::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;

    char buffer[BUFF_SIZE];
    OBUnitCell *uc = NULL;

    ofs << mol.GetTitle() << endl;

    if (!mol.HasData(OBGenericDataType::UnitCell))
      ofs << "   1.00000   1.00000   1.00000  90.00000  90.00000  90.00000\n";
    else
      {
        uc = (OBUnitCell*)mol.GetData(OBGenericDataType::UnitCell);
        snprintf(buffer, BUFF_SIZE,
                 "%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f",
                 uc->GetA(), uc->GetB(), uc->GetC(),
                 uc->GetAlpha() , uc->GetBeta(), uc->GetGamma());
        ofs << buffer << "\n";
      }

    vector3 v;
    FOR_ATOMS_OF_MOL(atom, mol)
      {
        v = atom->GetVector();
        if (uc != NULL)
          v = uc->CartesianToFractional(v);

        snprintf(buffer, BUFF_SIZE, "%s %10.5f%10.5f%10.5f",
                 etab.GetSymbol(atom->GetAtomicNum()),
                 v.x(),
                 v.y(),
                 v.z());
        ofs << buffer << endl;
      }

  bool CacaoFormat::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;

    OBAtom *atom;
    char buffer[BUFF_SIZE];
    vector<OBAtom*>::iterator i;

    snprintf(buffer, BUFF_SIZE, "%s\n",mol.GetTitle());
    ofs << buffer;
    snprintf(buffer, BUFF_SIZE, "%3d   DIST  0  0  0\n",mol.NumAtoms());
    ofs << buffer;

    if (!mol.HasData(OBGenericDataType::UnitCell))
      ofs << "CELL 1.,1.,1.,90.,90.,90.\n";
    else
      {
        OBUnitCell *uc = (OBUnitCell*)mol.GetData(OBGenericDataType::UnitCell);
        snprintf(buffer, BUFF_SIZE, "CELL %f,%f,%f,%f,%f,%f\n",
                 uc->GetA(), uc->GetB(), uc->GetC(),
                 uc->GetAlpha(), uc->GetBeta(), uc->GetGamma());
        ofs << buffer;
      }

    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      {
        snprintf(buffer,BUFF_SIZE,"%2s %7.4f, %7.4f, %7.4f\n",
                 etab.GetSymbol(atom->GetAtomicNum()),
                 atom->x(),
                 atom->y(),
                 atom->z());
        ofs << buffer;
      }

    return(true);
  }

  bool BGFFormat::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;
    mol.SetTitle( pConv->GetTitle()); //default title is the filename
    mol.BeginModify();

    char buffer[BUFF_SIZE];
    char tmp[16],tmptyp[16];
    vector<string> vs;

    while (ifs.getline(buffer,BUFF_SIZE)) {
      if (EQn(buffer,"CRYSTX",6)) {
        // Parse unit cell
        tokenize(vs,buffer," \n\t,");
        if (vs.size() != 7)
          continue; // something strange

        double A, B, C, Alpha, Beta, Gamma;
        A = atof(vs[1].c_str());
        B = atof(vs[2].c_str());
        C = atof(vs[3].c_str());
        Alpha = atof(vs[4].c_str());
        Beta  = atof(vs[5].c_str());
        Gamma = atof(vs[6].c_str());
        OBUnitCell *uc = new OBUnitCell;
        uc->SetOrigin(fileformatInput);
        uc->SetData(A, B, C, Alpha, Beta, Gamma);
        mol.SetData(uc);
      } else if (EQn(buffer,"FORMAT",6))
        break;
    }

    ttab.SetFromType("DRE");
    ttab.SetToType("INT");
    OBAtom *atom;
    double x,y,z,chrg;
    for (;;)
      {
        if (!ifs.getline(buffer,BUFF_SIZE))
          break;
        if (EQn(buffer,"FORMAT",6))
          break;

        sscanf(buffer,"%*s %*s %*s %*s %*s %*s %lf %lf %lf %15s %*s %*s %lf",
               &x,&y,&z,
               tmptyp,
               &chrg);
        atom = mol.NewAtom();

        ttab.Translate(tmp,tmptyp);
        atom->SetType(tmp);

        CleanAtomType(tmptyp);
        atom->SetAtomicNum(etab.GetAtomicNum(tmptyp));

        atom->SetVector(x,y,z);
      }
    unsigned int i;
    vector<int> vtmp;
    vector<vector<int> > vcon;
    vector<vector<int> > vord;

    for (i = 0; i < mol.NumAtoms();i++)
      {
        vcon.push_back(vtmp);
        vord.push_back(vtmp);
      }

    unsigned int bgn;
    for (;;)
      {
        if (!ifs.getline(buffer,BUFF_SIZE) || EQn(buffer,"END",3))
          break;

        tokenize(vs,buffer);
        if (vs.empty() || vs.size() < 3 || vs.size() > 10)
          continue;

        if (EQn(buffer,"CONECT",6))
          {
            bgn = atoi((char*)vs[1].c_str()) - 1;
            if (bgn < 1 || bgn > mol.NumAtoms())
              continue;
            for (i = 2;i < vs.size();i++)
              {
                vcon[bgn].push_back(atoi((char*)vs[i].c_str()));
                vord[bgn].push_back(1);
              }
          }
        else
          if (EQn(buffer,"ORDER",5))
            {
              bgn = atoi((char*)vs[1].c_str()) - 1;
//.........这里部分代码省略.........

  void SuperCellExtension::fillCell()
  {
    /* Change coords back to inverse space, apply the space group transforms
     *  then change coords back to real space
     */
    if (!m_molecule)
      return;

    OBUnitCell *uc = m_molecule->OBUnitCell();
    if (!uc) {
      qDebug() << "No unit cell found - fillCell() returning...";
      return;
    }

    const SpaceGroup *sg = uc->GetSpaceGroup(); // the actual space group and transformations for this unit cell
    if (sg) {
      qDebug() << "Space group:" << sg->GetId();// << sg->GetHMName();
      // We operate on a copy of the Avogadro molecule
      // For each atom, we loop through:
      // * convert the coords back to inverse space
      // * apply the transformations
      // * create new (duplicate) atoms
      OBMol mol = m_molecule->OBMol();
      vector3 uniqueV, newV;
      list<vector3> transformedVectors; // list of symmetry-defined copies of the atom
      list<vector3>::iterator transformIterator, duplicateIterator;
      vector3 updatedCoordinate;
      bool foundDuplicate;

      OBAtom *addAtom;
      QList<OBAtom*> atoms; // keep the current list of unique atoms -- don't double-create
      list<vector3> coordinates; // all coordinates to prevent duplicates
      FOR_ATOMS_OF_MOL(atom, mol)
        atoms.push_back(&(*atom));

      foreach(OBAtom *atom, atoms) {
        uniqueV = atom->GetVector();
        // Assert: won't crash because we already ensure uc != NULL
        uniqueV = uc->CartesianToFractional(uniqueV);
        uniqueV = transformedFractionalCoordinate(uniqueV);
        coordinates.push_back(uniqueV);

        transformedVectors = sg->Transform(uniqueV);
        for (transformIterator = transformedVectors.begin();
             transformIterator != transformedVectors.end(); ++transformIterator) {
          // coordinates are in reciprocal space -- check if it's in the unit cell
          // if not, transform it in place
          updatedCoordinate = transformedFractionalCoordinate(*transformIterator);
          foundDuplicate = false;

          // Check if the transformed coordinate is a duplicate of an atom
          for (duplicateIterator = coordinates.begin();
               duplicateIterator != coordinates.end(); ++duplicateIterator) {
            if (duplicateIterator->distSq(updatedCoordinate) < 1.0e-4) {
              foundDuplicate = true;
              break;
            }
          }
          if (foundDuplicate)
            continue;

          coordinates.push_back(updatedCoordinate); // make sure to check the new atom for dupes
          addAtom = mol.NewAtom();
          addAtom->Duplicate(atom);
          addAtom->SetVector(uc->FractionalToCartesian(updatedCoordinate));
        } // end loop of transformed atoms

        // Put the original atom into the proper space in the unit cell too
        atom->SetVector(uc->FractionalToCartesian(uniqueV));
      } // end loop of atoms

      m_molecule->setOBMol(&mol);
      qDebug() << "Spacegroups done...";

      // Need a fresh pointer to the new unit cell - setOBMol is invalidating
      // the old one. This should be cleaned up to use a more permanent data
      // structure.
      uc = m_molecule->OBUnitCell();
      uc->SetSpaceGroup(1);
    }

  bool CSSRFormat::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;

    char buffer[BUFF_SIZE];

    if (!mol.HasData(OBGenericDataType::UnitCell))
      {
        snprintf(buffer, BUFF_SIZE,
                 " REFERENCE STRUCTURE = 00000   A,B,C =%8.3f%8.3f%8.3f",
                 1.0,1.0,1.0);
        ofs << buffer << endl;
        snprintf(buffer, BUFF_SIZE,
                 "   ALPHA,BETA,GAMMA =%8.3f%8.3f%8.3f    SPGR =    P1"
                 , 90.0f, 90.0f, 90.0f);
        ofs << buffer << endl;
      }
    else
      {
        OBUnitCell *uc = (OBUnitCell*)mol.GetData(OBGenericDataType::UnitCell);
        snprintf(buffer, BUFF_SIZE,
                 " REFERENCE STRUCTURE = 00000   A,B,C =%8.3f%8.3f%8.3f",
                 uc->GetA(), uc->GetB(), uc->GetC());
        ofs << buffer << endl;
        snprintf(buffer, BUFF_SIZE,
                 "   ALPHA,BETA,GAMMA =%8.3f%8.3f%8.3f    SPGR =    P1",
                 uc->GetAlpha() , uc->GetBeta(), uc->GetGamma());
        ofs << buffer << endl;
      }

    snprintf(buffer, BUFF_SIZE, "%4d   1 %s\n",mol.NumAtoms(), mol.GetTitle());
    ofs << buffer << endl << endl;

    OBAtom *atom,*nbr;
    vector<OBAtom*>::iterator i;
    vector<OBBond*>::iterator j;
    vector<int> vtmp(106,0);
    int bonds;

    for(atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      {
        //assign_pdb_number(pdb_types,atom->GetIdx());
        vtmp[atom->GetAtomicNum()]++;
        snprintf(buffer, BUFF_SIZE, "%4d%2s%-3d  %9.5f %9.5f %9.5f ",
                 atom->GetIdx(),
                 etab.GetSymbol(atom->GetAtomicNum()),
                 vtmp[atom->GetAtomicNum()],
                 atom->x(),
                 atom->y(),
                 atom->z());
        ofs << buffer;
        bonds = 0;
        for (nbr = atom->BeginNbrAtom(j); nbr; nbr = atom->NextNbrAtom(j))
          {
            if (bonds > 8) break;
            snprintf(buffer, BUFF_SIZE, "%4d",nbr->GetIdx());
            ofs << buffer;
            bonds++;
          }
        for (; bonds < 8; bonds ++)
          {
            snprintf(buffer, BUFF_SIZE, "%4d",0);
            ofs << buffer;
          }
        snprintf(buffer, BUFF_SIZE, " %7.3f%4d", atom->GetPartialCharge(), 1);
        ofs << buffer << endl;
      }

    return(true);
  }

//.........这里部分代码省略.........
          while (vs.size() >= 4) {
            if (!createdAtoms) {
              atom = mol.NewAtom();
              //set atomic number
              atomicNum = OBElements::GetAtomicNum(vs[0].c_str());
              if (atomicNum == 0) {
                atomicNum = atoi(vs[0].c_str());
              }
              atom->SetAtomicNum(atomicNum);
            }
            x = atof((char*)vs[1].c_str());
            y = atof((char*)vs[2].c_str());
            z = atof((char*)vs[3].c_str());
            atomPositions.push_back(vector3(x, y, z)); // we may have a movie or animation

            ifs.getline(buffer, BUFF_SIZE);
            tokenize(vs, buffer);
          }
          createdAtoms = true; // don't run NewAtom() anymore
        }
        else if ( strstr(buffer, "PRIMVEC")
                 || strstr(buffer, "CONVVEC") ) {
          // translation vectors
          numTranslationVectors = 0; // if we have an animation
          while (numTranslationVectors < 3 && ifs.getline(buffer,BUFF_SIZE)) {
            tokenize(vs,buffer); // we really need to check that it's 3 entries only
            if (vs.size() < 3) return false; // timvdm 18/06/2008
            x = atof((char*)vs[0].c_str());
            y = atof((char*)vs[1].c_str());
            z = atof((char*)vs[2].c_str());
            translationVectors[numTranslationVectors++].Set(x, y, z);
          }
        }
        else if (strstr(buffer, "PRIMCOORD") != NULL) {
          // read the coordinates
          ifs.getline(buffer, BUFF_SIZE);
          tokenize(vs, buffer);
          if (vs.size() < 2) return false;
          int numAtoms = atoi(vs[0].c_str());
          for (int a = 0; a < numAtoms; ++a) {
            if (!ifs.getline(buffer,BUFF_SIZE))
              break;
            tokenize(vs,buffer);
            if (vs.size() < 4)
              break;

            if (!createdAtoms) {
              atom = mol.NewAtom();
              //set atomic number
              atomicNum = OBElements::GetAtomicNum(vs[0].c_str());
              if (atomicNum == 0) {
                atomicNum = atoi(vs[0].c_str());
              }
              atom->SetAtomicNum(atomicNum);
            }
            x = atof((char*)vs[1].c_str());
            y = atof((char*)vs[2].c_str());
            z = atof((char*)vs[3].c_str());
            atomPositions.push_back(vector3(x, y, z));
          }
        }
      }

    mol.EndModify();

    int natom = mol.NumAtoms();
    int numConformers = atomPositions.size() / natom;
    for (int i = 0; i < numConformers; ++i) {
      double *coordinates = new double[natom * 3];
      for (int j = 0; j < natom; ++j) {
        vector3 currentPosition = atomPositions[i*natom + j];
        coordinates[j*3] = currentPosition.x();
        coordinates[j*3 + 1] = currentPosition.y();
        coordinates[j*3 + 2] = currentPosition.z();
      }
      mol.AddConformer(coordinates);
    }
    // Delete first conformer, created by EndModify, bunch of 0s
    mol.DeleteConformer(0);
    // Set geometry to last one
    mol.SetConformer(mol.NumConformers() - 1);

    if (!pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.ConnectTheDots();
    if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.PerceiveBondOrders();

    // Add final properties
    mol.SetTitle(title);
    if (numTranslationVectors == 3) {
      OBUnitCell *uc = new OBUnitCell;
      uc->SetOrigin(fileformatInput);
      uc->SetData(translationVectors[0],
                  translationVectors[1],
                  translationVectors[2]);
      mol.SetData(uc);
    }

    return(true);
  }

//.........这里部分代码省略.........
    // but remove last geometry -- it's a duplicate
    if (vconf.size() > 1)
      vconf.pop_back();
    mol.SetConformers(vconf);
    mol.SetConformer(mol.NumConformers() - 1);
    // Copy the conformer data too
    confData->SetDimension(confDimensions);
    confData->SetEnergies(confEnergies);
    confData->SetForces(confForces);
    mol.SetData(confData);

    // Attach orbital data, if there is any
    if (orbitals.size() > 0)
      {
        OBOrbitalData *od = new OBOrbitalData;
        if (aHOMO == bHOMO) {
          od->LoadClosedShellOrbitals(orbitals, symmetries, aHOMO);
        } else {
          // we have to separate the alpha and beta vectors
          std::vector<double>      betaOrbitals;
          std::vector<std::string> betaSymmetries;
          unsigned int initialSize = orbitals.size();
          for (unsigned int i = betaStart; i < initialSize; ++i) {
            betaOrbitals.push_back(orbitals[i]);
            if (symmetries.size() > 0)
              betaSymmetries.push_back(symmetries[i]);
          }
          // ok, now erase the end elements of orbitals and symmetries
          for (unsigned int i = betaStart; i < initialSize; ++i) {
            orbitals.pop_back();
            if (symmetries.size() > 0)
              symmetries.pop_back();
          }
          // and load the alphas and betas
          od->LoadAlphaOrbitals(orbitals, symmetries, aHOMO);
          od->LoadBetaOrbitals(betaOrbitals, betaSymmetries, bHOMO);
        }
        od->SetOrigin(fileformatInput);
        mol.SetData(od);
      }

    //Attach vibrational data, if there is any, to molecule
    if(Frequencies.size()>0)
    {
      OBVibrationData* vd = new OBVibrationData;
      vd->SetData(Lx, Frequencies, Intensities);
      vd->SetOrigin(fileformatInput);
      mol.SetData(vd);
    }
    //Attach rotational data, if there is any, to molecule
    if(RotConsts[0]!=0.0)
    {
      OBRotationData* rd = new OBRotationData;
      rd->SetData(RotorType, RotConsts, RotSymNum);
      rd->SetOrigin(fileformatInput);
      mol.SetData(rd);
    }
    // Attach unit cell translation vectors if found
    if (numTranslationVectors > 0) {
      OBUnitCell* uc = new OBUnitCell;
      uc->SetData(translationVectors[0], translationVectors[1], translationVectors[2]);
      uc->SetOrigin(fileformatInput);
      mol.SetData(uc);
    }
    //Attach electronic transition data, if there is any, to molecule
    if(Forces.size() > 0 && Forces.size() == Wavelengths.size())
    {
      OBElectronicTransitionData* etd = new OBElectronicTransitionData;
      etd->SetData(Wavelengths, Forces);
      if (EDipole.size() == Forces.size())
        etd->SetEDipole(EDipole);
      if (RotatoryStrengthsLength.size() == Forces.size())
        etd->SetRotatoryStrengthsLength(RotatoryStrengthsLength);
      if (RotatoryStrengthsVelocity.size() == Forces.size())
        etd->SetRotatoryStrengthsVelocity(RotatoryStrengthsVelocity);
      etd->SetOrigin(fileformatInput);
      mol.SetData(etd);
    }

    if (!pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.ConnectTheDots();
    if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.PerceiveBondOrders();

    if (hasPartialCharges) {
      mol.SetPartialChargesPerceived();

      // Annotate that partial charges come from Mulliken
      OBPairData *dp = new OBPairData;
      dp->SetAttribute("PartialCharges");
      dp->SetValue(chargeModel); // Mulliken, ESP, etc.
      dp->SetOrigin(fileformatInput);
      mol.SetData(dp);
    }
    mol.SetTotalCharge(charge);
    mol.SetTotalSpinMultiplicity(spin);

    mol.SetTitle(title);
    return(true);
  }

//.........这里部分代码省略.........
    int spaceGroup = -1;
    if (location != string::npos) {
      // e.g., "#230"
      string spaceGroupNumber = readTitle.substr(location + 1, 4); // +1 to skip #
      string::size_type nonNumber = spaceGroupNumber.find_first_not_of("0123456789");
      if (nonNumber != string::npos)
        spaceGroupNumber.erase(nonNumber);
      // Finally get the space group from the file
      spaceGroup = atoi(spaceGroupNumber.c_str());
    }

    location = readTitle.find_first_not_of(" \t\n\r");
    // Is there non-whitespace
    if (location != string::npos)
      mol.SetTitle(readTitle);
    else
      mol.SetTitle(defaultTitle);

    vector3 v1, v2, v3;
    double x,y,z;
    vector<string> vs;
    OBAtom *atom;
    int atomicNum;
    bool setCellVectors = false;
    // go through remaining lines, particularly looking for cell vectors
    mol.BeginModify();

    while(ifs.peek() != EOF && ifs.good()) {
      ifs.getline(buffer,BUFF_SIZE);
      if (strstr(buffer, "Primitive vectors")) {
        // three lines: a(#) = .. .. ..
        ifs.getline(buffer, BUFF_SIZE);
        tokenize(vs,buffer);
        if (vs.size() != 5)
          continue;
        v1.SetX(atof(vs[2].c_str()));
        v1.SetY(atof(vs[3].c_str()));
        v1.SetZ(atof(vs[4].c_str()));

        ifs.getline(buffer, BUFF_SIZE);
        tokenize(vs,buffer);
        if (vs.size() != 5)
          continue;
        v2.SetX(atof(vs[2].c_str()));
        v2.SetY(atof(vs[3].c_str()));
        v2.SetZ(atof(vs[4].c_str()));

        ifs.getline(buffer, BUFF_SIZE);
        tokenize(vs,buffer);
        if (vs.size() != 5)
          continue;
        v3.SetX(atof(vs[2].c_str()));
        v3.SetY(atof(vs[3].c_str()));
        v3.SetZ(atof(vs[4].c_str()));

        setCellVectors = true;
      }
      if (strstr(buffer, "Basis Vectors:")) {
        ifs.getline(buffer,BUFF_SIZE);  // column titles
        ifs.getline(buffer,BUFF_SIZE); // blank
        // real atomic data
        while (ifs.getline(buffer,BUFF_SIZE) && strlen(buffer)) {
          tokenize(vs,buffer);
          if (vs.size() != 7)
            break;

          atom = mol.NewAtom();
          // check to see if first column is number or element symbol
          // (PCModel has files of the form X Y Z symbol)
          atomicNum = OBElements::GetAtomicNum(vs[0].c_str());
          x = atof(vs[4].c_str());
          y = atof(vs[5].c_str());
          z = atof(vs[6].c_str());
          atom->SetVector(x,y,z);
          atom->SetAtomicNum(atomicNum);
        }
      }
    }

    if (setCellVectors) {
      OBUnitCell* uc = new OBUnitCell;
      uc->SetData(v1, v2, v3);
      uc->SetOrigin(fileformatInput);
      mol.SetData(uc);

      // hopefully, we have a space group too
      if (spaceGroup != -1) {
        uc->SetSpaceGroup(spaceGroup);
      }
    }

    if (!pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.ConnectTheDots();
    if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.PerceiveBondOrders();

    mol.EndModify();

    return(true);
  }

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