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

  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

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

  QUndoCommand* UnitCellExtension::performAction(QAction *, GLWidget *widget)
  {
    // FIXME: this is bad mmmkay
    m_widget = widget;

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

    OBUnitCell *uc = m_molecule->OBUnitCell();

    if (uc == NULL) {
      // show warning and ask if the user wants to create a unit cell
      // (otherwise this extension isn't very useful)

      QMessageBox::StandardButton ret;
      ret = QMessageBox::warning(qobject_cast<QWidget*>(parent()),
                                 tr("Avogadro"),
                                 tr("This document is currently an isolated molecule.\n\n"
                                    "Do you want to create a crystal unit cell?"),
                                 QMessageBox::Yes
                                 | QMessageBox::No);
      if (ret == QMessageBox::Yes) {
        // Set some initial data (e.g., a box about the size of the molecule)
        // and one unit cell in each direction
        uc = new OBUnitCell;
        double estimatedSize = widget->radius() + 2.0;
        uc->SetData(estimatedSize, estimatedSize, estimatedSize,
                    90.0, 90.0, 90.0);
        m_molecule->setOBUnitCell(uc);

        widget->setUnitCells(1, 1, 1);
        widget->setUnitCellColor(m_color);
      } else { // do nothing -- user picked "Cancel"
        return NULL;
      }

    } // end if (existing unit cell or create a new one)

    // Don't emit signals while we update these
    disconnect(m_dialog, SIGNAL(unitCellDisplayChanged(int, int, int)),
            this, SLOT(unitCellDisplayChanged(int, int, int)));
    disconnect(m_dialog, SIGNAL(unitCellParametersChanged(double, double, double, double, double, double)),
            this, SLOT(unitCellParametersChanged(double, double, double, double, double, double)));

    m_dialog->aCells(widget->aCells());
    m_dialog->bCells(widget->bCells());
    m_dialog->cCells(widget->cCells());

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

    // OK, now we can handle signal/slots
    connect(m_dialog, SIGNAL(unitCellDisplayChanged(int, int, int)),
            this, SLOT(unitCellDisplayChanged(int, int, int)));
    connect(m_dialog, SIGNAL(unitCellParametersChanged(double, double, double, double, double, double)),
            this, SLOT(unitCellParametersChanged(double, double, double, double, double, double)));

    m_dialog->show();

    return NULL;
  }

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

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

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

//.........这里部分代码省略.........
    }
    // Reset end-of-file pointers etc.
    ifs.clear();
    ifs.seekg(0);  //rewind

    mol.BeginModify();
    while (ifs.getline(buffer,BUFF_SIZE))
    {
        if(strstr(buffer, "Entering Gaussian") != NULL)
        {
            //Put some metadata into OBCommentData
            string comment("Gaussian ");

            if(NULL != strchr(buffer,'='))
            {
                comment += strchr(buffer,'=')+2;
                comment += "";
                for(unsigned i=0; i<115 && ifs; ++i)
                {
                    ifs.getline(buffer,BUFF_SIZE);
                    if(strstr(buffer,"Revision") != NULL)
                    {
                        if (buffer[strlen(buffer)-1] == ',')
                        {
                            buffer[strlen(buffer)-1] = '\0';
                        }
                        add_unique_pairdata_to_mol(&mol,"program",buffer,0);
                    }
                    else if(buffer[1]=='#')
                    {
                        //the line describing the method
                        comment += buffer;
                        OBCommentData *cd = new OBCommentData;
                        cd->SetData(comment);
                        cd->SetOrigin(fileformatInput);
                        mol.SetData(cd);

                        tokenize(vs,buffer);
                        if (vs.size() > 1)
                        {
                            char *str = strdup(vs[1].c_str());
                            char *ptr = strchr(str,'/');

                            if (NULL != ptr)
                            {
                                *ptr = ' ';
                                add_unique_pairdata_to_mol(&mol,"basis",ptr,0);
                                *ptr = '\0';
                                add_unique_pairdata_to_mol(&mol,"method",str,0);
                            }
                        }

                        break;
                    }
                }
            }
        }

        else if (strstr(buffer,"Multiplicity") != NULL)
        {
            tokenize(vs, buffer, " \t\n");
            if (vs.size() == 6)
            {
                total_charge = atoi(vs[2].c_str());
                spin_multiplicity = atoi(vs[5].c_str());
            }

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

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

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

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

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

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

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

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

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

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

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

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

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

    char buffer[BUFF_SIZE], tag[BUFF_SIZE];
    double x,y,z;
    double alat = 1.0;
    vector<string> vs;
    matrix3x3 ortho;
    int atomicNum;
    OBUnitCell *cell = new OBUnitCell();
    bool hasEnthalpy=false;
    double enthalpy, pv;

    pmol->BeginModify();

    while (ifs.getline(buffer,BUFF_SIZE)) {

      // Older version of pwscf may use this for alat
      if (strstr(buffer, "lattice parameter (a_0)")) {
        tokenize(vs, buffer);
        alat = atof(vs.at(4).c_str());
      }

      // Newer versions will use this for alat instead
      if (strstr(buffer, "lattice parameter (alat)")) {
        tokenize(vs, buffer);
        alat = atof(vs.at(4).c_str());
      }

      // Unit cell info
      // Newer versions will also say "CELL_PARAMETERS" to complain that no
      // units were specified
      if (strstr(buffer, "CELL_PARAMETERS") &&
          !strstr(buffer, "no units specified in CELL_PARAMETERS card")) {
        // Discover units
        double conv = 1.0;
        tokenize(vs, buffer);

        if (strstr(vs[1].c_str(), "alat")) {
          conv = alat * BOHR_TO_ANGSTROM;
        }
        else if (strstr(vs[1].c_str(), "bohr")) {
          conv = BOHR_TO_ANGSTROM;
        }
        // Add others if needed

        double v11, v12, v13,
          v21, v22, v23,
          v31, v32, v33;

        ifs.getline(buffer,BUFF_SIZE); // v1
        tokenize(vs, buffer);
        v11 = atof(vs.at(0).c_str()) * conv;
        v12 = atof(vs.at(1).c_str()) * conv;
        v13 = atof(vs.at(2).c_str()) * conv;

        ifs.getline(buffer,BUFF_SIZE); // v2
        tokenize(vs, buffer);
        v21 = atof(vs.at(0).c_str()) * conv;
        v22 = atof(vs.at(1).c_str()) * conv;
        v23 = atof(vs.at(2).c_str()) * conv;

        ifs.getline(buffer,BUFF_SIZE); // v3
        tokenize(vs, buffer);
        v31 = atof(vs.at(0).c_str()) * conv;
        v32 = atof(vs.at(1).c_str()) * conv;
        v33 = atof(vs.at(2).c_str()) * conv;

        // Build unit cell
        cell->SetData(vector3(v11,v12,v13),
                      vector3(v21,v22,v23),
                      vector3(v31,v32,v33));
      }

      // Unit cell info (for non-variable cell calcs)
      if (strstr(buffer, "crystal axes: (cart. coord. in units of a_0)") ||
          strstr(buffer, "crystal axes: (cart. coord. in units of alat)")) {
        double conv = alat * BOHR_TO_ANGSTROM;
        double v11, v12, v13,
          v21, v22, v23,
          v31, v32, v33;

        ifs.getline(buffer,BUFF_SIZE); // v1
        tokenize(vs, buffer);
        v11 = atof(vs.at(3).c_str()) * conv;
        v12 = atof(vs.at(4).c_str()) * conv;
        v13 = atof(vs.at(5).c_str()) * conv;

        ifs.getline(buffer,BUFF_SIZE); // v2
        tokenize(vs, buffer);
        v21 = atof(vs.at(3).c_str()) * conv;
        v22 = atof(vs.at(4).c_str()) * conv;
        v23 = atof(vs.at(5).c_str()) * conv;

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