C++ Topology::SetParmBox方法代码示例

Exec::RetType Exec_ParmBox::Execute(CpptrajState& State, ArgList& argIn) {
  Box pbox;
  bool nobox = false;
  if ( argIn.hasKey("nobox") )
    nobox = true;
  else {
    pbox.SetX( argIn.getKeyDouble("x",0) );
    pbox.SetY( argIn.getKeyDouble("y",0) );
    pbox.SetZ( argIn.getKeyDouble("z",0) );
    pbox.SetAlpha( argIn.getKeyDouble("alpha",0) );
    pbox.SetBeta(  argIn.getKeyDouble("beta",0)  );
    pbox.SetGamma( argIn.getKeyDouble("gamma",0) );
  }
  Topology* parm = State.DSL().GetTopByIndex( argIn );
  if (parm == 0) return CpptrajState::ERR;
  if (nobox)
    mprintf("\tRemoving box information from parm %i:%s\n", parm->Pindex(), parm->c_str());
  else
    // Fill in missing parm box information from specified parm
    pbox.SetMissingInfo( parm->ParmBox() );
  if (argIn.hasKey("truncoct")) pbox.SetTruncOct();
  parm->SetParmBox( pbox );
  parm->ParmBox().PrintInfo();
  return CpptrajState::OK;
}

/** Read file as a Tripos Mol2 file. */
int Parm_Mol2::ReadParm(FileName const& fname, Topology &parmOut) {
  Mol2File infile;
  if (infile.OpenRead(fname)) return 1;
  mprintf("    Reading Mol2 file %s as topology file.\n",infile.Filename().base());
  // Get @<TRIPOS>MOLECULE information
  if (infile.ReadMolecule()) return 1;
  parmOut.SetParmName( infile.Mol2Title(), infile.Filename() );

  // Get @<TRIPOS>ATOM information
  if (infile.ScanTo( Mol2File::ATOM)) return 1;
  double XYZ[3];
  for (int atom=0; atom < infile.Mol2Natoms(); atom++) {
    if ( infile.Mol2XYZ(XYZ) ) return 1;
    parmOut.AddTopAtom( infile.Mol2Atom(), infile.Mol2Residue(), XYZ );
  }

  // Get @<TRIPOS>BOND information [optional]
  int at1 = 0;
  int at2 = 0;
  if (infile.ScanTo(Mol2File::BOND)==0) {
    for (int bond=0; bond < infile.Mol2Nbonds(); bond++) {
      if (infile.Mol2Bond(at1, at2)) return 1;
      // mol2 atom #s start from 1
      parmOut.AddBond(at1-1, at2-1);
    }
    needsBondSearch_ = false;
  } else {
    mprintf("      Mol2 file does not contain bond information.\n");
    needsBondSearch_ = true;
  }

  // No box
  parmOut.SetParmBox( Box() );

  mprintf("    Mol2 contains %i atoms, %i residues,\n", parmOut.Natom(),parmOut.Nres());
  //mprintf("    %i bonds to H, %i other bonds.\n", parmOut.NbondsWithH,parmOut.NbondsWithoutH);

  infile.CloseFile();

  return 0;
}

/** Read file as a Tinker file. */
int Parm_Tinker::ReadParm(FileName const& fname, Topology &parmOut) {
  TinkerFile infile;
  infile.SetTinkerName( fname );
  if (infile.OpenTinker()) return 1;
  mprintf("\tReading Tinker file %s as topology file.\n",infile.Filename().base());
  // Allocate memory for coordinates.
  double* Coords = new double[ infile.TinkerNatom() * 3 ];
  std::vector<int> Bonds;
  std::vector<Atom> Atoms = infile.ReadTinkerAtoms(Coords, Bonds);
  if (Atoms.empty()) return 1;
  // Use up to first 3 chars of title as residue name.
  std::string resname;
  for (std::string::const_iterator c = infile.TinkerTitle().begin();
                                   c != infile.TinkerTitle().end(); ++c)
    resname += *c;
  if (resname.size() > 3) resname.resize(3);
  Residue tinker_res( resname, 0, ' ', ' ' );
  // Put atoms into topology
  const double* XYZ = Coords;
  for (std::vector<Atom>::const_iterator atom = Atoms.begin();
                                         atom != Atoms.end();
                                       ++atom, XYZ += 3)
    parmOut.AddTopAtom( *atom, tinker_res, XYZ );
  delete[] Coords;
  // Add bond information
  for (std::vector<int>::const_iterator bond = Bonds.begin();
                                        bond != Bonds.end(); bond += 2)
    parmOut.AddBond( *bond, *(bond+1) );
  // Try to set up residue info based on bonds.
  if (parmOut.Setup_NoResInfo()) return 1;
  // Set topology box info.
  parmOut.SetParmBox( infile.TinkerBox() );
  parmOut.SetParmName( infile.TinkerTitle(), infile.Filename() );
  mprintf("\tTinker file contains %i atoms, %i residues,\n", parmOut.Natom(),parmOut.Nres());
  //mprintf("    %i bonds to H, %i other bonds.\n", parmOut.NbondsWithH,parmOut.NbondsWithoutH);

  infile.CloseFile();

  return 0;
}

//.........这里部分代码省略.........
  }
  // Does this CIF contain multiple models?
  int Nmodels = 0;
  int model_col = block.ColumnIndex("pdbx_PDB_model_num");
  if (model_col != -1) {
    line = block.end();
    --line;
    Nmodels = convertToInteger( (*line)[model_col] );
    if (Nmodels > 1)
      mprintf("Warning: CIF '%s' contains %i models. Using first model for topology.\n", 
              fname.full(), Nmodels);
  }
  // Get essential columns
  int COL[NENTRY];
  for (int i = 0; i < (int)NENTRY; i++) {
    COL[i] = block.ColumnIndex(Entries[i]);
    if (COL[i] == -1) {
      mprinterr("Error: In CIF file '%s' could not find entry '%s' in block '%s'\n",
                fname.full(), Entries[i], block.Header().c_str());
      return 1;
    }
    if (debug_>0) mprintf("DEBUG: '%s' column = %i\n", Entries[i], COL[i]);
  }
  // Get optional columns
  int occ_col = block.ColumnIndex("occupancy");
  int bfac_col = block.ColumnIndex("B_iso_or_equiv");
  int icode_col = block.ColumnIndex("pdbx_PDB_ins_code");
  int altloc_col = block.ColumnIndex("label_alt_id");
  std::vector<AtomExtra> extra;

  // Loop over all atom sites
  int current_res = 0;
  double XYZ[3];
  double occupancy = 1.0;
  double bfactor = 0.0;
  char altloc = ' ';
  char icode;
  icode = ' ';
  Frame Coords;
  for (line = block.begin(); line != block.end(); ++line) {
    // If more than 1 model check if we are done.
    if (Nmodels > 1) {
      if ( convertToInteger( (*line)[model_col] ) > 1 )
        break;
    }
    if (occ_col != -1) occupancy = convertToDouble( (*line)[ occ_col ] );
    if (bfac_col != -1) bfactor = convertToDouble( (*line)[ bfac_col ] );
    if (altloc_col != -1) altloc = (*line)[ altloc_col ][0];
    // '.' altloc means blank?
    if (altloc == '.') altloc = ' ';
    extra.push_back( AtomExtra(occupancy, bfactor, altloc) );
    if (icode_col != -1) {
      icode = (*line)[ icode_col ][0];
      // '?' icode means blank
      if (icode == '?') icode = ' ';
    }
    XYZ[0] = convertToDouble( (*line)[ COL[X] ] );
    XYZ[1] = convertToDouble( (*line)[ COL[Y] ] );
    XYZ[2] = convertToDouble( (*line)[ COL[Z] ] );
    NameType currentResName( (*line)[ COL[RNAME] ] );
    // It seems that in some CIF files, there doesnt have to be a residue
    // number. Check if residue name has changed.
    if ( (*line)[ COL[RNUM] ][0] == '.' ) {
      Topology::res_iterator lastResidue = TopIn.ResEnd();
      --lastResidue;
      if ( currentResName != (*lastResidue).Name() )
        current_res = TopIn.Nres() + 1;
    } else
      current_res = convertToInteger( (*line)[ COL[RNUM] ] );
    TopIn.AddTopAtom( Atom((*line)[ COL[ANAME] ], "  "),
                      Residue(currentResName, current_res, icode,
                              (*line)[ COL[CHAINID] ][0]) );
    Coords.AddXYZ( XYZ );
  }
  if (TopIn.SetExtraAtomInfo( 0, extra )) return 1;
  // Search for bonds // FIXME nobondsearch?
  BondSearch( TopIn, Coords, Offset_, debug_ );
  // Get title. 
  CIFfile::DataBlock const& entryblock = infile.GetDataBlock("_entry");
  std::string ciftitle;
  if (!entryblock.empty())
    ciftitle = entryblock.Data("id");
  TopIn.SetParmName( ciftitle, infile.CIFname() );
  // Get unit cell parameters if present.
  CIFfile::DataBlock const& cellblock = infile.GetDataBlock("_cell");
  if (!cellblock.empty()) {
    double cif_box[6];
    cif_box[0] = convertToDouble( cellblock.Data("length_a") );
    cif_box[1] = convertToDouble( cellblock.Data("length_b") );
    cif_box[2] = convertToDouble( cellblock.Data("length_c") );
    cif_box[3] = convertToDouble( cellblock.Data("angle_alpha") );
    cif_box[4] = convertToDouble( cellblock.Data("angle_beta" ) );
    cif_box[5] = convertToDouble( cellblock.Data("angle_gamma") );
    mprintf("\tRead cell info from CIF: a=%g b=%g c=%g alpha=%g beta=%g gamma=%g\n",
              cif_box[0], cif_box[1], cif_box[2], cif_box[3], cif_box[4], cif_box[5]);
    TopIn.SetParmBox( Box(cif_box) ); 
  }
  
  return 0;
}