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

/* FUNCTION: newMolecule */
OBMol* NXMoleculeSet::newMolecule() {
	OBMol* molecule = new OBMol();
	NXMoleculeData* moleculeData = new NXMoleculeData();
	moleculeData->SetIdx(NextMoleculeIndex);
	NextMoleculeIndex++;
	molecule->SetData(moleculeData);
	molecules.push_back(molecule);
	return molecule;
}

// Delete hydrogens should not remove charged or isotopic hydrogens or [H][H] or [Cu][H][Cu]
// or hydrogens with assigned atom classes
void test_Issue178_DeleteHydrogens()
{
  OBConversion conv;
  conv.SetInFormat("smi");
  OBMol mol;
  // Test DeleteHydrogens() and DeleteNonPolarHydrogens()
  static const char *smi[] = { "C[H]", "[H][H]", "C[1H]", "C[H]C", "C[H+]" };
  int numHs[] = { 0, 2, 1, 1, 1 };
  for (int i = 0; i < 5; ++i) {
    for (int j = 0; j < 2; ++j) {
      conv.ReadString(&mol, smi[i]);
      if (j == 0)
        mol.DeleteHydrogens();
      else
        mol.DeleteNonPolarHydrogens();
      int myNumHs = 0;
      FOR_ATOMS_OF_MOL(atom, mol)
        if (atom->IsHydrogen())
          myNumHs++;
      OB_COMPARE(myNumHs, numHs[i]);
    }
  }
  // Test DeletePolarHydrogens()
  static const char *smiB[] = { "N[H]", "[H][H]", "N[1H]", "N[H]C", "N[H+]" };
  int numHsB[] = { 0, 2, 1, 1, 1 };
  for (int i = 0; i < 5; ++i) {
    conv.ReadString(&mol, smiB[i]);
    mol.DeletePolarHydrogens();
    int myNumHs = 0;
    FOR_ATOMS_OF_MOL(atom, mol)
      if (atom->IsHydrogen())
        myNumHs++;
    OB_COMPARE(myNumHs, numHsB[i]);
  }
  // Test atom class
  // Currently, the SMILES parser does not retain atom classes for hydrogens on reading so...
  conv.ReadString(&mol, "C[H]");
  OBAtomClassData *ac = new OBAtomClassData;
  ac->Add(2, 99); // Assign the hydrogen (atom 2) a class of 99
  mol.SetData(ac);
  mol.DeleteHydrogens();
  int myNumHs = 0;
  FOR_ATOMS_OF_MOL(atom, mol)
    if (atom->IsHydrogen())
      myNumHs++;
  OB_COMPARE(myNumHs, 1);
}

  bool QEqCharges::ComputeCharges(OBMol &mol)
  {

     ///////////////////////////////////////////////////////////////////////////////
    //Some OpenBabel bookkeeping that I copied from the Gasteiger scheme
    mol.SetPartialChargesPerceived();

    OBPairData *dp = new OBPairData;
    dp->SetAttribute("PartialCharges");
    dp->SetValue("QEq");
    dp->SetOrigin(perceived);
    mol.SetData(dp);


    ///////////////////////////////////////////////////////////////////////////////
    //Read in atomic information from OpenBabel molecule and parameterize

    //Read in total number of atoms
    int i, N = mol.NumAtoms();

    Hardness = MatrixXd::Zero(N+1, N+1);
    Voltage = VectorXd::Zero(N+1);
    Electronegativity = VectorXd::Zero(N);
    VectorXd BasisSet = VectorXd::Zero(N);

    Vector3d Parameters;

    FOR_ATOMS_OF_MOL(atom, mol)
    {
       	Parameters = GetParameters(atom->GetAtomicNum(), atom->GetFormalCharge());
	i = atom->GetIdx() - 1;

	if (Parameters[0] == 0.)
        {
		stringstream msg;
		msg << "Some QEq Parameters not found!" << endl
	            << "Parameters not found for atom no. " << i+1 << endl
		    << "Atom will be ignored in the charge computation.";
		obErrorLog.ThrowError(__FUNCTION__, msg.str(), obError);
	}

	Electronegativity[i] = Parameters[0];
	Hardness(i,i) = Parameters[1];
	BasisSet[i] = Parameters[2];
    }

bool ThermoFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
{
  OBMol* pmol = pOb->CastAndClear<OBMol>();
  if(!pmol)
    return false;
  bool stopOnEnd = pConv->IsOption("e",OBConversion::INOPTIONS)!=NULL;
  pmol->SetDimension(0);
  OBNasaThermoData* pND = new OBNasaThermoData; //to store rate constant data
  pND->SetOrigin(fileformatInput);
  pmol->SetData(pND);

  istream &ifs = *pConv->GetInStream();

  double DefaultMidT = 1500;
  char ln[BUFF_SIZE];
  unsigned int i;

  //find line with 1 in col 80
  do
  {
    if(!ifs.getline(ln,BUFF_SIZE) || stopOnEnd && !strncasecmp(ln,"END",3))
      return false;
  }while(ln[79]!='1');

  char phase, nam[25], dum[7], elname[3];
  elname[2]=0;
  int elnum;
  double Coeff[14];

  sscanf(ln,"%18s%6s",nam,dum);
  pmol->SetTitle(nam);
  char* p=ln+24;
  if(ln[80]=='&')
  {
    //Reaction Design extension
    p+=20;
    string line;
    if(!getline(ifs,line))return false;
    vector<string> toks;
    tokenize(toks,line," \t\n\r");
    for(i=0;i<toks.size();i+=2)
    {
      OBAtom atom;
      atom.SetAtomicNum(etab.GetAtomicNum(toks[i].c_str()));
      elnum = atoi(toks[i+1].c_str());
      atom.ForceNoH();
      for(;elnum>0;--elnum)
        pmol->AddAtom(atom);
    }
  }
  else
  {
    for(i=0;i<4;i++,p+=5)
    {
      char snum[4]={0,0,0,0};//Was problem with F   10   0 reading as ten
      sscanf(p,"%c%c%c%c%c",elname,elname+1,snum,snum+1,snum+2);
      elnum=atoi(snum);
      if(elname[0]!=' ' && elname[0]!='0')
      {
        if(elname[1]==' ')
          elname[1]=0;
        OBAtom atom;
        atom.SetAtomicNum(etab.GetAtomicNum(elname));
        atom.ForceNoH();
        for(;elnum>0;--elnum)
          pmol->AddAtom(atom);
      }
    }
  }
  double LoT, HiT, MidT=0;
  /* int nc = */sscanf(p,"%c%10lf%10lf10%lf",&phase, &LoT, &HiT, &MidT);
  pND->SetPhase(phase);
  pND->SetLoT(LoT);
  pND->SetHiT(HiT);
  if(MidT>HiT || MidT<LoT)
    MidT=DefaultMidT;
  pND->SetMidT(MidT);
  if (!ifs.getline(ln, BUFF_SIZE)) return false;
  p=ln;
  for(i=0;i<5;i++,p+=15)
    sscanf(p,"%15lf",&Coeff[i]);
  if (!ifs.getline(ln, BUFF_SIZE)) return false;
  p=ln;
  for(i=5;i<10;i++,p+=15)
    sscanf(p,"%15lf",&Coeff[i]);
  if (!ifs.getline(ln, BUFF_SIZE)) return false;
  p=ln;
  for(i=10;i<14;i++,p+=15)
    sscanf(p,"%15lf",&Coeff[i]);

  for(i=0;i<14;++i)
    pND->SetCoeff(i, Coeff[i]);

  pmol->AssignSpinMultiplicity();
  return true;
}

  /*! Makes a new OBMol on the heap by combining two molecules according to the rule below. 
    If both have OBGenericData of the same type, or OBPairData with the
    same attribute,  the version from pFirst is used.
    Returns a pointer to a new OBMol which will need deleting by the calling program
    (probably by being sent to an output format). 
    If the molecules cannot be regarded as being the same structure a NULL
    pointer is returned and an error message logged.
    
    pFirst and pSecond and the objects they point to are not changed. (const
    modifiers difficult because class OBMol not designed appropriately)
    
    Combining molecules: rules for each of the three parts
    Title:
    Use the title of pFirst unless it has none, when use that of pSecond.
    Warning if neither molecule has a title.
    
    Structure
    - a structure with atoms replaces one with no atoms
    - a structure with bonds replaces one with no bonds,
    provided the formula is the same, else an error.
    - structures with atoms and bonds are compared by InChI; error if not the same. 
    - a structure with 3D coordinates replaces one with 2D coordinates
    - a structure with 2D coordinates replace one with 0D coordinates
    
    OBGenericData
    OBPairData
  */
  OBMol* OBMoleculeFormat::MakeCombinedMolecule(OBMol* pFirst, OBMol* pSecond)
  {
    //Decide on which OBMol provides the new title
    string title("No title");
    if(*pFirst->GetTitle()!=0)
      title = pFirst->GetTitle();
    else
      {
        if(*pSecond->GetTitle()!=0)
          title = pSecond->GetTitle();
        else
          obErrorLog.ThrowError(__FUNCTION__,"Combined molecule has no title", obWarning);
      }

    //Decide on which OBMol provides the new structure
    bool swap=false;
    if(pFirst->NumAtoms()==0 && pSecond->NumAtoms()!=0)
      swap=true;
    else if(pSecond->NumAtoms()!=0)
      {
        if(pFirst->GetSpacedFormula()!=pSecond->GetSpacedFormula())
          {
            obErrorLog.ThrowError(__FUNCTION__, 
                                  "Molecules with name = " + title + " have different formula",obError);
            return NULL;
          }
        else
          {
            if(pSecond->NumBonds()!=0 && pFirst->NumBonds()==0)
              swap=true;
            else
              {
                //Compare by inchi; error if different NOT YET IMPLEMENTED
                //Use the one with the higher dimension
                if(pSecond->GetDimension() > pFirst->GetDimension())
                  swap=true;
              }
          }
      }

    OBMol* pNewMol = new OBMol;
    pNewMol->SetTitle(title);

    OBMol* pMain = swap ? pSecond : pFirst;
    OBMol* pOther = swap ? pFirst : pSecond;
    
    *pNewMol = *pMain; //Now copies all data 

    //Copy some OBGenericData from the OBMol which did not provide the structure
    vector<OBGenericData*>::iterator igd;
    for(igd=pOther->BeginData();igd!=pOther->EndData();++igd)
      {
        //copy only if not already data of the same type from molecule already copied to pNewMol
        unsigned datatype = (*igd)->GetDataType();
        OBGenericData* pData = pNewMol->GetData(datatype);
        if(datatype==OBGenericDataType::PairData)
          {
            if(pData->GetAttribute() == (*igd)->GetAttribute())
              continue;
          }
        else if(pNewMol->GetData(datatype)!=NULL)
          continue;

        OBGenericData* pCopiedData = (*igd)->Clone(pNewMol);
        pNewMol->SetData(pCopiedData);
      }
    return pNewMol;
  }

//.........这里部分代码省略.........
    vector3 origin(x, y, z);

    // now three lines with the x, y, and z axes
    vector<vector3> axes;
    for (unsigned int i = 0; i < 3; ++i) {
      if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "delta", 5))
        return false;
      else {
        tokenize(vs, buffer);
        if (vs.size() != 4)
          return false;
        x = atof(vs[1].c_str());
        y = atof(vs[2].c_str());
        z = atof(vs[3].c_str());
        axes.push_back(vector3(x, y, z));
      }
    }

    // Two remaining header lines before the data:
    /*
      object 2 class gridconnections counts nx ny nz
      object 3 class array type double rank 0 times n data follows
    */
    if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "object", 6))
      return false;
    if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "object", 6))
      return false;

    pmol->BeginModify();
    pmol->SetDimension(3);

    OBGridData *gd = new OBGridData;
    gd->SetAttribute("OpenDX");

    // get all values as one vector<double>
    char *endptr;
    vector<double> values;
    int n = voxels[0]*voxels[1]*voxels[2];
    int line = 0;
    values.reserve(n);
    while (ifs.getline(buffer, BUFF_SIZE))
    {
      ++line;
      if (EQn(buffer, "attribute", 9))
        break; // we're finished with reading data -- although we should probably have a voxel check in here too

      tokenize(vs, buffer);
      if (vs.size() == 0)
      {
        errorMsg << "Problem reading the OpenDX grid file: cannot"
                 << " read line " << line
                 << ", there does not appear to be any data in it.\n"
                 << buffer << "\n";
        obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obError);
        return false;
      }

      for (unsigned int l = 0; l < vs.size(); ++l)
      {
        values.push_back(strtod(static_cast<const char*>(vs[l].c_str()), &endptr));
      }
    }

    gd->SetNumberOfPoints(voxels[0], voxels[1], voxels[2]);
    gd->SetLimits(origin, axes[0], axes[1], axes[2]);
    gd->SetUnit(OBGridData::ANGSTROM);
    gd->SetOrigin(fileformatInput); // i.e., is this data from a file or determined by Open Babel
    gd->SetValues(values); // set the values
    pmol->SetData(gd); // store the grids in the OBMol
    pmol->EndModify();

    // Trailing control lines
    /*
     attribute "dep" string "positions"
     object "regular positions regular connections" class field
     component "positions" value 1
     component "connections" value 2
     component "data" value 3
    */
    if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "object", 6))
      return false;
    if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "component", 9))
      return false;
    if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "component", 9))
      return false;
    if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "component", 9))
      return false;

    // clean out any remaining blank lines
    std::streampos ipos;
    do
    {
      ipos = ifs.tellg();
      ifs.getline(buffer,BUFF_SIZE);
    }
    while(strlen(buffer) == 0 && !ifs.eof() );
    ifs.seekg(ipos);

    return true;
  }

//.........这里部分代码省略.........
      //If there is a -s option, reference molecule has only those atoms that are matched
      //Call the -s option from here
      bool ret = _pOpIsoM->Do(pmol, _stext.c_str(), pmap, pConv);
      // Get the atoms that were matched
      vector<int> ats = _pOpIsoM->GetMatchAtoms();
      if(!ats.empty())
      {
        // Make a vector of the matching atom coordinates...
        for(vector<int>::iterator iter=ats.begin(); iter!=ats.end(); ++iter)
          _refvec.push_back((pmol->GetAtom(*iter))->GetVector());
        // ...and use a vector reference
        _align.SetRef(_refvec);
      }
      // Stop -s option being called normally, although it will still be called once
      //  in the DoOps loop already started for the current (first) molecule.
      pConv->RemoveOption("s",OBConversion::GENOPTIONS);
      if(!ret)
      {
        // the first molecule did not match the -s option so a reference molecule
        // could not be made. Keep trying.
        _refMol.Clear();
        //obErrorLog.ThrowError(__FUNCTION__, "The first molecule did not match the -s option\n"
        //  "so the reference structure was not derived from it", obWarning, onceOnly);
        return false; //not matched
      }
    }
  }

  //All molecules
  if(pmol->GetDimension()!= _refMol.GetDimension())
  {
    stringstream ss;
    ss << "The molecule" << pmol->GetTitle()
       << " does not have the same dimensions as the reference molecule "
       << _refMol.GetTitle() << " and is ignored.";
       obErrorLog.ThrowError(__FUNCTION__, ss.str().c_str(), obError);
    return false;
  }

  if(_pOpIsoM) //Using -s option
  {
    //Ignore mol if it does not pass -s option
    if(!_pOpIsoM->Do(pmol, "", pmap, pConv)) // "" means will use existing parameters
      return false;

    // Get the atoms equivalent to those in ref molecule
    vector<int> ats = _pOpIsoM->GetMatchAtoms();

    // Make a vector of their coordinates and get the centroid
    vector<vector3> vec;
    vector3 centroid;
    for(vector<int>::iterator iter=ats.begin(); iter!=ats.end(); ++iter) {
      vector3 v = pmol->GetAtom(*iter)->GetVector();
      centroid += v;
      vec.push_back(v);
    }
    centroid /= vec.size();

    // Do the alignment
    _align.SetTarget(vec);
    if(!_align.Align())
      return false;

    // Get the centroid of the reference atoms
    vector3 ref_centroid;
    for(vector<vector3>::iterator iter=_refvec.begin(); iter!=_refvec.end(); ++iter)
      ref_centroid += *iter;
    ref_centroid /= _refvec.size();

    //subtract the centroid, rotate the target molecule, then add the centroid
    matrix3x3 rotmatrix = _align.GetRotMatrix();
    for (unsigned int i = 1; i <= pmol->NumAtoms(); ++i)
    {
      vector3 tmpvec = pmol->GetAtom(i)->GetVector();
      tmpvec -= centroid;
      tmpvec *= rotmatrix; //apply the rotation
      tmpvec += ref_centroid;
      pmol->GetAtom(i)->SetVector(tmpvec);
    }
  }
  else //Not using -s option)
  {
    _align.SetTargetMol(*pmol);
    if(!_align.Align())
      return false;
    _align.UpdateCoords(pmol);
  }

  //Save rmsd as a property
  OBPairData* dp = new OBPairData;
  dp->SetAttribute("rmsd");
  double val = _align.GetRMSD();
  if(val<1e-12)
    val = 0.0;
  dp->SetValue(toString(val));
  dp->SetOrigin(local);
  pmol->SetData(dp);

  return true;
}

//.........这里部分代码省略.........
        _paramFileLoaded = true;
      } else
      {
        return false;
      }
    }

    // Calculate atomic properties based around their ionic charge
    for (i = 0; i < N; i++)
    {
      atom = mol.GetAtom(i + 1);
      a = atom->GetAtomicNum();
      c = _chargeCenter[a];

      // Fail if ionization data is missing for any atom in the molecule
      if (_ionizations[a][c + 1] == -1 || _ionizations[a][c] == -1 || a > TABLE_OF_ELEMENTS_SIZE)
      {
        obErrorLog.ThrowError(__FUNCTION__, "Insufficient ionization data for atoms in the given molecule. Update `data/eqeqIonizations.txt` with missing information and re-run this function.", obError);
        return false;
      }

      J(i) = _ionizations[a][c + 1] - _ionizations[a][c];
      chi(i) = 0.5 * (_ionizations[a][c + 1] + _ionizations[a][c]) - (a == 1? 0 : c * J(i));
    }

    // If a unit cell is defined, use the periodic Ewald calculation
    if (mol.HasData(OBGenericDataType::UnitCell))
    {
      // Get unit cell and calculate its Fourier transform + volume
      obuc = (OBUnitCell *) mol.GetData(OBGenericDataType::UnitCell);
      unitcell = obuc->GetCellMatrix();
      fourier = (2 * PI * unitcell.inverse()).transpose();
      cellVolume = obuc->GetCellVolume();

      // Get the number of radial unit cells to use in x, y, and z
      numNeighbors[0] = int(ceil(minCellLength / (2.0 * (obuc->GetA())))) - 1;
      numNeighbors[1] = int(ceil(minCellLength / (2.0 * (obuc->GetB())))) - 1;
      numNeighbors[2] = int(ceil(minCellLength / (2.0 * (obuc->GetC())))) - 1;

      for (i = 0; i < N; i++)
      {
        atom = mol.GetAtom(i + 1);
        for (j = 0; j < N; j++)
        {
          dx = atom->GetVector() - (mol.GetAtom(j + 1))->GetVector();
          J_ij(i, j) = GetPeriodicEwaldJij(J(i), J(j), dx, (i == j), unitcell, fourier, cellVolume, numNeighbors);
        }
      }
    // If no unit cell, use the simplified nonperiodic calculation
    } else
    {
      for (i = 0; i < N; i++)
      {
        atom = mol.GetAtom(i + 1);
        for (j = 0; j < N; j++)
        {
          J_ij(i, j) = GetNonperiodicJij(J(i), J(j), atom->GetDistance(j + 1), (i == j));
        }
        return false;
      }
    }

    // Formulate problem as A x = b, where x is the calculated partial charges
    // First equation is a simple overall balance: sum(Q) = 0
    A.row(0) = VectorXf::Ones(N);
    b(0) = 0;

    // Remaining equations are based off of the fact that, at equilibrium, the
    // energy of the system changes equally for a change in any charge:
    //     dE/dQ_1 = dE/dQ_2 = ... = dE/dQ_N
    A.block(1, 0, N - 1, N) = J_ij.block(0, 0, N - 1, N) - J_ij.block(1, 0, N - 1, N);
    b.tail(N - 1) = chi.tail(N - 1) - chi.head(N - 1);

    // The solution is a list of charges in the system
    x = A.colPivHouseholderQr().solve(b);

    // Now we are done calculating, pass all this back to OpenBabel molecule
    mol.SetPartialChargesPerceived();
    OBPairData *dp = new OBPairData;
    dp->SetAttribute("PartialCharges");
    dp->SetValue("EQEq");
    dp->SetOrigin(perceived);
    mol.SetData(dp);

    m_partialCharges.clear();
    m_partialCharges.reserve(N);
    m_formalCharges.clear();
    m_formalCharges.reserve(N);

    for (i = 0; i < N; i ++)
    {
      atom = mol.GetAtom(i + 1);
      atom->SetPartialCharge(x(i));
      m_partialCharges.push_back(x(i));
      m_formalCharges.push_back(atom->GetFormalCharge());
    }

    obErrorLog.ThrowError(__FUNCTION__, "EQEq charges successfully assigned.", obInfo);
    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;

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