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

int Parm_CharmmPsf::WriteParm(FileName const& fname, Topology const& parm) {
  // TODO: CMAP etc info
  CpptrajFile outfile;
  if (outfile.OpenWrite(fname)) return 1;
  // Write PSF
  outfile.Printf("PSF\n\n");
  // Write title
  std::string titleOut = parm.ParmName();
  titleOut.resize(78);
  outfile.Printf("%8i !NTITLE\n* %-78s\n\n", 1, titleOut.c_str());
  // Write NATOM section
  outfile.Printf("%8i !NATOM\n", parm.Natom());
  unsigned int idx = 1;
  // Make fake segment ids for now.
  char segid[2];
  segid[0] = 'A';
  segid[1] = '\0';
  mprintf("Warning: Assigning single letter segment IDs.\n");
  int currentMol = 0;
  bool inSolvent = false;
  for (Topology::atom_iterator atom = parm.begin(); atom != parm.end(); ++atom, ++idx) {
    int resnum = atom->ResNum();
    if (atom->MolNum() != currentMol) {
      if (!inSolvent) {
        inSolvent = parm.Mol(atom->MolNum()).IsSolvent();
        currentMol = atom->MolNum();
        segid[0]++;
      } else
        inSolvent = parm.Mol(atom->MolNum()).IsSolvent();
    }
    // TODO: Print type name for xplor-like PSF
    int typeindex = atom->TypeIndex() + 1;
    // If type begins with digit, assume charmm numbers were read as
    // type. Currently Amber types all begin with letters.
    if (isdigit(atom->Type()[0]))
      typeindex = convertToInteger( *(atom->Type()) );
    // ATOM# SEGID RES# RES ATNAME ATTYPE CHRG MASS (REST OF COLUMNS ARE LIKELY FOR CMAP AND CHEQ)
    outfile.Printf("%8i %-4s %-4i %-4s %-4s %4i %14.6G %9g  %10i\n", idx, segid,
                   parm.Res(resnum).OriginalResNum(), parm.Res(resnum).c_str(),
                   atom->c_str(), typeindex, atom->Charge(),
                   atom->Mass(), 0);
  }
  outfile.Printf("\n");
  // Write NBOND section
  outfile.Printf("%8u !NBOND: bonds\n", parm.Bonds().size() + parm.BondsH().size());
  idx = 1;
  for (BondArray::const_iterator bond = parm.BondsH().begin();
                                 bond != parm.BondsH().end(); ++bond, ++idx)
  {
    outfile.Printf("%8i%8i", bond->A1()+1, bond->A2()+1);
    if ((idx % 4)==0) outfile.Printf("\n"); 
  }
  for (BondArray::const_iterator bond = parm.Bonds().begin();
                                 bond != parm.Bonds().end(); ++bond, ++idx)
  {
    outfile.Printf("%8i%8i", bond->A1()+1, bond->A2()+1);
    if ((idx % 4)==0) outfile.Printf("\n"); 
  }
  if ((idx % 4)!=0) outfile.Printf("\n");
  outfile.Printf("\n");
  // Write NTHETA section
  outfile.Printf("%8u !NTHETA: angles\n", parm.Angles().size() + parm.AnglesH().size());
  idx = 1;
  for (AngleArray::const_iterator ang = parm.AnglesH().begin();
                                  ang != parm.AnglesH().end(); ++ang, ++idx)
  {
    outfile.Printf("%8i%8i%8i", ang->A1()+1, ang->A2()+1, ang->A3()+1);
    if ((idx % 3)==0) outfile.Printf("\n");
  }
  for (AngleArray::const_iterator ang = parm.Angles().begin();
                                  ang != parm.Angles().end(); ++ang, ++idx)
  {
    outfile.Printf("%8i%8i%8i", ang->A1()+1, ang->A2()+1, ang->A3()+1);
    if ((idx % 3)==0) outfile.Printf("\n");
  }
  if ((idx % 3)==0) outfile.Printf("\n");
  outfile.Printf("\n");
  // Write out NPHI section
  outfile.Printf("%8u !NPHI: dihedrals\n", parm.Dihedrals().size() + parm.DihedralsH().size());
  idx = 1;
  for (DihedralArray::const_iterator dih = parm.DihedralsH().begin();
                                     dih != parm.DihedralsH().end(); ++dih, ++idx)
  {
    outfile.Printf("%8i%8i%8i%8i", dih->A1()+1, dih->A2()+1, dih->A3()+1, dih->A4()+1);
    if ((idx % 2)==0) outfile.Printf("\n");
  }
  for (DihedralArray::const_iterator dih = parm.Dihedrals().begin();
                                     dih != parm.Dihedrals().end(); ++dih, ++idx)
  {
    outfile.Printf("%8i%8i%8i%8i", dih->A1()+1, dih->A2()+1, dih->A3()+1, dih->A4()+1);
    if ((idx % 2)==0) outfile.Printf("\n");
  }
  if ((idx % 2)==0) outfile.Printf("\n");
  outfile.Printf("\n");

  outfile.CloseFile();
  return 0;
}

/** Find potential symmetric atoms. All residues up to the last selected
  * residue are considered.
  */
int SymmetricRmsdCalc::SetupSymmRMSD(Topology const& topIn, AtomMask const& tgtMask, bool remapIn)
{
  // Allocate space for remapping selected atoms in target frame. This will
  // also put the correct masses in based on the mask.
  tgtRemap_.SetupFrameFromMask(tgtMask, topIn.Atoms());
  // Create map of original atom numbers to selected indices
  Iarray SelectedIdx( topIn.Natom(), -1 );
  int tgtIdx = 0;
  for (int originalAtom = 0; originalAtom != topIn.Natom(); ++originalAtom)
    if ( originalAtom == tgtMask[tgtIdx] )
      SelectedIdx[originalAtom] = tgtIdx++;
  if (debug_ > 0) {
    mprintf("DEBUG: Original atom -> Selected Index mapping:\n");
    for (int originalAtom = 0; originalAtom != topIn.Natom(); ++originalAtom)
      mprintf("\t%8i -> %8i\n", originalAtom + 1, SelectedIdx[originalAtom] + 1);
  }
  // Create initial 1 to 1 atom map for all selected atoms; indices in 
  // SymmetricAtomIndices will correspond to positions in AMap.
  AMap_.resize( tgtRemap_.Natom() );
  // Determine last selected residue.
  int last_res = topIn[tgtMask.back()].ResNum() + 1;
  mprintf("\tResidues up to %s will be considered for symmetry correction.\n",
          topIn.TruncResNameNum(last_res-1).c_str());
  // In each residue, determine which selected atoms are symmetric.
  SymmetricAtomIndices_.clear();
  AtomMap resmap;
  if (debug_ > 1) resmap.SetDebug(1);
  for (int res = 0; res < last_res; ++res) {
    AtomMap::AtomIndexArray residue_SymmetricGroups;
    if (resmap.SymmetricAtoms(topIn, residue_SymmetricGroups, res)) {
      mprinterr("Error: Finding symmetric atoms in residue '%s'\n",
                topIn.TruncResNameNum(res).c_str());
      return 1;
    }
    if (!residue_SymmetricGroups.empty()) {
      // Which atoms in symmetric groups are selected?
      bool resHasSelectedSymmAtoms = false;
      for (AtomMap::AtomIndexArray::const_iterator symmGroup = residue_SymmetricGroups.begin();
                                                   symmGroup != residue_SymmetricGroups.end();
                                                 ++symmGroup)
      {
        Iarray selectedAtomIndices;
        for (Iarray::const_iterator atnum = symmGroup->begin();
                                    atnum != symmGroup->end(); ++atnum)
        {
          if ( SelectedIdx[*atnum] != -1 )
            selectedAtomIndices.push_back( SelectedIdx[*atnum] ); // Store tgtMask indices
        }
        if (!selectedAtomIndices.empty()) {
          SymmetricAtomIndices_.push_back( selectedAtomIndices );
          resHasSelectedSymmAtoms = true;
        }
      }
      // If remapping and not all atoms in a residue are selected, warn user.
      // TODO: Should they just be considered even if not selected?
      if (remapIn && resHasSelectedSymmAtoms) {
        for (int atom = topIn.Res(res).FirstAtom(); atom != topIn.Res(res).LastAtom(); ++atom)
          if (SelectedIdx[atom] == -1) {
            mprintf("Warning: Not all atoms selected in residue '%s'. Re-mapped\n"
                    "Warning:   structures may appear distorted.\n", 
                    topIn.TruncResNameNum(res).c_str());
            break;
          }
      }
    }
  }
  if (debug_ > 0) {
    mprintf("DEBUG: Potential Symmetric Atom Groups:\n");
    for (AtomIndexArray::const_iterator symmatoms = SymmetricAtomIndices_.begin();
                                        symmatoms != SymmetricAtomIndices_.end();
                                        ++symmatoms)
    {
      mprintf("\t%8u) ", symmatoms - SymmetricAtomIndices_.begin());
      for (Iarray::const_iterator atom = symmatoms->begin();
                                  atom != symmatoms->end(); ++atom)
        mprintf(" %s(%i)", topIn.AtomMaskName(tgtMask[*atom]).c_str(), tgtMask[*atom] + 1);
      mprintf("\n");
    } 
  }
  return 0;
}

// ----- M A I N ---------------------------------------------------------------
int main(int argc, char** argv) {
  SetWorldSilent(true); // No STDOUT output from cpptraj routines.
  std::string topname, crdname, title, bres, pqr, sybyltype, writeconect;
  std::string aatm(" pdbatom"), ter_opt(" terbyres"), box(" sg \"P 1\"");
  TrajectoryFile::TrajFormatType fmt = TrajectoryFile::PDBFILE;
  bool ctr_origin = false;
  bool useExtendedInfo = false;
  int res_offset = 0;
  int debug = 0;
  int numSoloArgs = 0;
  for (int i = 1; i < argc; ++i) {
    std::string arg( argv[i] );
    if (arg == "-p" && i+1 != argc && topname.empty()) // Topology
      topname = std::string( argv[++i] );
    else if (arg == "-c" && i+1 != argc && crdname.empty()) // Coords
      crdname = std::string( argv[++i] );
    else if (arg == "-tit" && i+1 != argc && title.empty()) // Title
      title = " title " + std::string( argv[++i] );
    else if (arg == "-offset" && i+1 != argc) // Residue # offset
      res_offset = convertToInteger( argv[++i] );
    else if ((arg == "-d" || arg == "--debug") && i+1 != argc) // Debug level
      debug = convertToInteger( argv[++i] );
    else if (arg == "-h" || arg == "--help") { // Help
      Help(argv[0], true);
      return 0;
    } else if (arg == "-v" || arg == "--version") { // Version info
      WriteVersion();
      return 0;
    } else if (arg == "-aatm") // Amber atom names, include extra pts
      aatm.assign(" include_ep");
    else if (arg == "-sybyl") // Amber atom types to SYBYL
      sybyltype.assign(" sybyltype");
    else if (arg == "-conect") // Write CONECT records from bond info
      writeconect.assign(" conect");
    else if (arg == "-ep") // PDB atom names, include extra pts
      aatm.append(" include_ep");
    else if (arg == "-bres") // PDB residue names
      bres.assign(" pdbres");
    else if (arg == "-ext") // Use extended PDB info from Topology
      useExtendedInfo = true;
    else if (arg == "-ctr")  // Center on origin
      ctr_origin = true;
    else if (arg == "-noter") // No TER cards
      ter_opt.assign(" noter");
    else if (arg == "-nobox") // No CRYST1 record
      box.assign(" nobox");
    else if (arg == "-pqr") { // Charge/Radii in occ/bfactor cols
      pqr.assign(" dumpq");
      ++numSoloArgs;
    } else if (arg == "-mol2") { // output as mol2
      fmt = TrajectoryFile::MOL2FILE;
      ++numSoloArgs;
    } else if (Unsupported(arg)) {
      mprinterr("Error: Option '%s' is not yet supported.\n\n", arg.c_str());
      return 1;
    } else {
      mprinterr("Error: Unrecognized option '%s'\n", arg.c_str());
      Help(argv[0], false);
      return 1;
    }
  }
  if (debug > 0) WriteVersion();
  // Check command line for errors.
  if (topname.empty()) topname.assign("prmtop");
  if (debug > 0 && crdname.empty())
    mprinterr("| Reading Amber restart from STDIN\n");
  if (numSoloArgs > 1) {
    mprinterr("Error: Only one alternate output format option may be specified (found %i)\n",
              numSoloArgs);
    Help(argv[0], true);
    return 1;
  }
  if (!sybyltype.empty() && fmt != TrajectoryFile::MOL2FILE) {
    mprinterr("Warning: -sybyl is only valid for MOL2 file output.\n");
    sybyltype.clear();
  }
  if (debug > 0) {
    mprinterr("Warning: debug is %i; debug info will be written to STDOUT.\n", debug);
    SetWorldSilent(false);
  }
  // Topology
  ParmFile pfile;
  Topology parm;
  if (pfile.ReadTopology(parm, topname, debug)) {
    if (topname == "prmtop") Help(argv[0], false);
    return 1;
  }
  if (!useExtendedInfo)
    parm.ResetPDBinfo();
  if (res_offset != 0)
    for (int r = 0; r < parm.Nres(); r++)
      parm.SetRes(r).SetOriginalNum( parm.Res(r).OriginalResNum() + res_offset );
  ArgList trajArgs;
  // Input coords
  Frame TrajFrame;
  CoordinateInfo cInfo;
  if (!crdname.empty()) {
    Trajin_Single trajin;
    if (trajin.SetupTrajRead(crdname, trajArgs, &parm)) return 1;
//.........这里部分代码省略.........