C++ OBUnitCell::GetSpaceGroup方法代码示例

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

  void UnitCellExtension::fillUnitCell()
  {
    /* 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 == NULL)
      return;

    const SpaceGroup *sg = uc->GetSpaceGroup(); // the actual space group and transformations for this unit cell

    if (!sg) {
      QMessageBox::warning(qobject_cast<QWidget*>(parent()),
                           tr("Avogadro"),
                           tr("This unit cell does not have an associated spacegroup."));
      return;
    }

    // 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
      #ifdef OPENBABEL_IS_NEWER_THAN_2_2_99
      uniqueV = uc->CartesianToFractional(uniqueV);
      #else
      uniqueV *= uc->GetFractionalMatrix();
      #endif
      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;

        addAtom = mol.NewAtom();
        addAtom->Duplicate(atom);
        #ifdef OPENBABEL_IS_NEWER_THAN_2_2_99
        addAtom->SetVector(uc->FractionalToCartesian(updatedCoordinate));
        #else
        addAtom->SetVector(uc->GetOrthoMatrix() * updatedCoordinate);
        #endif
      } // end loop of transformed atoms

      // Put the original atom into the proper space in the unit cell too
      #ifdef OPENBABEL_IS_NEWER_THAN_2_2_99
      atom->SetVector(uc->FractionalToCartesian(uniqueV));
      #else
      atom->SetVector(uc->GetOrthoMatrix() * uniqueV);
      #endif
    } // end loop of atoms

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

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

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

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