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

//.........这里部分代码省略.........
          }
        else if (strstr(buffer, "natom")) {
          tokenize(vs, buffer);
          if (vs.size() != 2)
            continue;
          natom = atoi(vs[1].c_str());
        }
        else if (strstr(buffer, "rprim")) {
          numTranslationVectors = 0;
          int column;
          for (int i = 0; i < 3; ++i) {
            tokenize(vs, buffer);
            if (vs.size() < 3)
              break;

            // first line, rprim takes up a token
            if (i == 0)
              column = 1;
            else
              column = 0;

            x = atof((char*)vs[column].c_str()) * BOHR_TO_ANGSTROM;
            y = atof((char*)vs[column+1].c_str()) * BOHR_TO_ANGSTROM;
            z = atof((char*)vs[column+2].c_str()) * BOHR_TO_ANGSTROM;
            translationVectors[numTranslationVectors++].Set(x, y,z);
            ifs.getline(buffer,BUFF_SIZE);
          }
        }
        else if (strstr(buffer, "Symmetries")) {
          tokenize(vs, buffer, "()");
          // Should be something like (#160)
          symmetryCode = atoi(vs[1].substr(1).c_str());
        }
        else if (strstr(buffer, "typat")) {
          tokenize(vs, buffer);
          atomTypes.clear();
          for (unsigned int i = 1; i < vs.size(); ++i) {
            atomTypes.push_back(atoi(vs[i].c_str()));
          }
        }
        else if (strstr(buffer, "znucl")) {
          tokenize(vs, buffer);
          // make sure znucl is first token
          if (vs[0] != "znucl")
            continue;
          // push back the remaining tokens into atomicNumbers
          atomicNumbers.clear();
          atomicNumbers.push_back(0); // abinit starts typat with 1
          for (unsigned int i = 1; i < vs.size(); ++i)
            atomicNumbers.push_back(int(atof(vs[i].c_str())));
        }
        // xangst
        // forces
      }

    for (int i = 0; i < natom; ++i) {
      atom = mol.NewAtom();
      //set atomic number
      int idx = atom->GetIdx();
      int type = atomTypes[idx - 1];
      atom->SetAtomicNum(atomicNumbers[type]);
      // we set the coordinates by conformers in another loop
    }

    mol.EndModify();

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

    // Attach unit cell translation vectors if found
    if (numTranslationVectors > 0) {
      OBUnitCell* uc = new OBUnitCell;
      uc->SetData(acell[0] * translationVectors[0], acell[1] * translationVectors[1], acell[2] * translationVectors[2]);
      uc->SetOrigin(fileformatInput);
      if (symmetryCode)
        uc->SetSpaceGroup(symmetryCode);
      mol.SetData(uc);
    }

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

//.........这里部分代码省略.........
            int atomicNum = atoi(vs[3].c_str());
            if (atomicNum == 0)
              atomicNum = 1; // hydrogen ?
            if (atomicNum <= 0 || atomicNum > etab.GetNumberOfElements())
              continue;

            // valid element, so create the atom
            atom = mol.NewAtom();
            atom->SetAtomicNum(atomicNum);
            continue;
          }
          else if (strstr(buffer, "XYZ") != NULL) {
            tokenize(vs, buffer);
            // size should be 6 -- need a test here
            if (vs.size() != 6) return false; // timvdm 18/06/2008
            vs[3].erase(0,1); // remove ( character
            vs[5].erase(vs[5].length()-2, 2); // remove trailing )) characters
            atom->SetVector(atof(vs[3].c_str()),
                            atof(vs[4].c_str()),
                            atof(vs[5].c_str()));
            continue;
          }
        } // end of atom records

        // bond information
        if (bondRecord) {
          if (strstr(buffer, "Atom1") != NULL) {
            tokenize(vs, buffer);
            if (vs.size() < 4) return false; // timvdm 18/06/2008
            vs[3].erase(vs[3].length()-1,1);
            startBondAtom = atoi(vs[3].c_str());
            continue;
          }
          else if (strstr(buffer, "Atom2") != NULL) {
            tokenize(vs, buffer);
            if (vs.size() < 4) return false; // timvdm 18/06/2008
            vs[3].erase(vs[3].length()-1,1);
            endBondAtom = atoi(vs[3].c_str());
            continue;
          }
          else if (strstr(buffer, "Type") != NULL) {
            tokenize(vs, buffer);
            if (vs.size() < 4) return false; // timvdm 18/06/2008
            vs[3].erase(vs[3].length()-1,1);
            bondOrder = atoi(vs[3].c_str());
            if (bondOrder == 4) // triple bond?
              bondOrder = 3;
            else if (bondOrder == 8) // aromatic?
              bondOrder = 5;
            else if (bondOrder != 2) // 1 OK, 2 OK, others unknown
              bondOrder = 1;
            continue;
          }
        }

        // ending a "tag" -- a lone ")" on a line
        if (strstr(buffer,")") != NULL && strstr(buffer, "(") == NULL) {
          openParens--;
          if (atomRecord) {
            atomRecord = false;
          }
          if (bondRecord) {
            // Bond records appear to be questionable
            mol.AddBond(startBondAtom - 1, endBondAtom - 1, bondOrder);
            bondRecord = false;
          }

          if (openParens == 0) {
            ifs.getline(buffer, BUFF_SIZE);
            break; // closed this molecule
          }
        }
      }

    mol.EndModify();

    // clean out any remaining blank lines
    while(ifs.peek() != EOF && ifs.good() &&
          (ifs.peek() == '\n' || ifs.peek() == '\r'))
      ifs.getline(buffer,BUFF_SIZE);

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

    if (numTranslationVectors > 0) {
      OBUnitCell* uc = new OBUnitCell;
      uc->SetData(translationVectors[0], translationVectors[1], translationVectors[2]);
      uc->SetOrigin(fileformatInput);
      if (setSpaceGroup) {
        uc->SetSpaceGroup(sg);
      }
      mol.SetData(uc);
    }

    return(true);
  }

  bool CARFormat::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();

    bool hasPartialCharges = false;
    char buffer[BUFF_SIZE];
    string str;
    double x,y,z;
    OBAtom *atom;
    vector<string> vs;

    mol.BeginModify();

    while (ifs.getline(buffer,BUFF_SIZE))
      {
        if(strstr(buffer,"end") != NULL)
          {
            if (mol.NumAtoms() > 0) // we've already read in a molecule, so exit
              break;
            // else, we hit the end of the previous molecular system
            // (in a multimolecule file)
            ifs.getline(buffer,BUFF_SIZE); // title
            ifs.getline(buffer,BUFF_SIZE); // DATE
          }

        if (strncmp(buffer, "!BIOSYM", 7) == 0)
          {
            continue;
          }

        if(strstr(buffer,"PBC") != NULL)
          {
            if(strstr(buffer,"ON") != NULL)
              {
                ifs.getline(buffer,BUFF_SIZE); // title
                ifs.getline(buffer,BUFF_SIZE); // DATE
                ifs.getline(buffer,BUFF_SIZE); // PBC a b c alpha beta gamma SG

                // parse cell parameters
                tokenize(vs,buffer);
                if (vs.size() == 8)
                  {
                    //parse cell values
                    double A,B,C,Alpha,Beta,Gamma;
                    A = atof((char*)vs[1].c_str());
                    B = atof((char*)vs[2].c_str());
                    C = atof((char*)vs[3].c_str());
                    Alpha = atof((char*)vs[4].c_str());
                    Beta  = atof((char*)vs[5].c_str());
                    Gamma = atof((char*)vs[6].c_str());
                    OBUnitCell *uc = new OBUnitCell;
                    uc->SetOrigin(fileformatInput);
                    uc->SetData(A, B, C, Alpha, Beta, Gamma);
                    uc->SetSpaceGroup(vs[7]);
                    mol.SetData(uc);
                  }
              }
            else // PBC=OFF
              {
                ifs.getline(buffer,BUFF_SIZE); // title
                ifs.getline(buffer,BUFF_SIZE); // !DATE
              }
            continue;
          } // PBC

        // reading real data!
        tokenize(vs,buffer);
        if (vs.size() < 8) {
          break;
        }

        atom = mol.NewAtom();
      
        atom->SetAtomicNum(etab.GetAtomicNum(vs[7].c_str()));
        x = atof((char*)vs[1].c_str());
        y = atof((char*)vs[2].c_str());
        z = atof((char*)vs[3].c_str());
        atom->SetVector(x,y,z);
      
        // vs[0] contains atom label
        // vs[4] contains "type of residue containing atom"
        // vs[5] contains "residue sequence name"
        // vs[6] contains "potential type of atom"
      
        if (vs.size() == 9)
          {
            atom->SetPartialCharge(atof((char*)vs[8].c_str()));
            hasPartialCharges = true;
          }
      }
    
    if (!pConv->IsOption("b",OBConversion::INOPTIONS))
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      }
    }
    // 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]));
        }
      }
    }
  }
  
  // Set spacegroup to P1, since we generated all symmetrics
  pUC->SetSpaceGroup("P1");
/*  
  list<vector3> transformedVectors; // list of symmetry-defined copies of the atom
  vector3 uniqueV, newV, updatedCoordinate;
    list<vector3> coordinates; // all coordinates to prevent duplicates

    vector3 uniqueV, newV, updatedCoordinate;
    list<vector3> transformedVectors; // list of symmetry-defined copies of the atom
    list<vector3>::iterator transformIterator, duplicateIterator;
    OBAtom *newAtom;
    list<OBAtom*> atoms; // keep the current list of unique atoms -- don't double-create
    list<vector3> coordinates; // all coordinates to prevent duplicates
    bool foundDuplicate;
    FOR_ATOMS_OF_MOL(atom, *mol)
      atoms.push_back(&(*atom));

    list<OBAtom*>::iterator i;
    for (i = atoms.begin(); i != atoms.end(); ++i) {
      uniqueV = (*i)->GetVector();
      uniqueV = 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
        newAtom = mol->NewAtom();
        newAtom->Duplicate(*i);
        newAtom->SetVector(FractionalToCartesian(updatedCoordinate));
      } // end loop of transformed atoms
      (*i)->SetVector(FractionalToCartesian(uniqueV));
    } // end loop of atoms
*/
  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);
  }

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