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

/** 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;
}

/** An atom pair list consists of 2 values for each entry, a beginning
  * index and ending index. For molecules and residues this is the first
  * and just beyond the last atom; for atoms it is just the atom itself
  * twice.
  */
Image::PairType Image::CreatePairList(Topology const& Parm, Mode modeIn,
                                       std::string const& maskExpression)
{
  PairType atomPairs;
  // Set up mask based on desired imaging mode.
  if ( modeIn == BYMOL || modeIn == BYRES ) {
    CharMask cmask( maskExpression );
    if ( Parm.SetupCharMask( cmask ) ) return atomPairs;
    cmask.MaskInfo();
    if (cmask.None()) return atomPairs;
    // Set up atom range for each entity to be imaged.
    if (modeIn == BYMOL) {
      atomPairs.reserve( Parm.Nmol()*2 );
      for (Topology::mol_iterator mol = Parm.MolStart();
                                  mol != Parm.MolEnd(); ++mol)
        CheckRange( atomPairs, cmask, mol->BeginAtom(), mol->EndAtom());
    } else { // BYRES
      atomPairs.reserve( Parm.Nres()*2 );
      for (Topology::res_iterator residue = Parm.ResStart();
                                  residue != Parm.ResEnd(); ++residue)
        CheckRange( atomPairs, cmask, residue->FirstAtom(), residue->LastAtom() );
    }
  } else { // BYATOM
    AtomMask imask( maskExpression );
    if ( Parm.SetupIntegerMask( imask ) ) return atomPairs;
    imask.MaskInfo();
    if (imask.None()) return atomPairs;
    atomPairs.reserve( Parm.Natom()*2 );
    for (AtomMask::const_iterator atom = imask.begin(); atom != imask.end(); ++atom) {
      atomPairs.push_back(  *atom    );
      atomPairs.push_back( (*atom)+1 );
    }
  }
//  mprintf("\tNumber of %ss to be imaged is %zu based on mask '%s'\n",
//           ModeString[modeIn], atomPairs.size()/2, maskIn.MaskString());
  return atomPairs;
}

// Parm_CIF::ReadParm()
int Parm_CIF::ReadParm(FileName const& fname, Topology &TopIn) {
  CIFfile infile;
  CIFfile::DataBlock::data_it line;

  if (infile.Read( fname, debug_ )) return 1;
  CIFfile::DataBlock const& block = infile.GetDataBlock("_atom_site");
  if (block.empty()) {
    mprinterr("Error: CIF data block '_atom_site' not found.\n");
    return 1;
  }
  // 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") );
//.........这里部分代码省略.........

// Exec_PermuteDihedrals::RandomizeAngles()
void Exec_PermuteDihedrals::RandomizeAngles(Frame& currentFrame, Topology const& topIn) {
  Matrix_3x3 rotationMatrix;
# ifdef DEBUG_PERMUTEDIHEDRALS
  // DEBUG
  int debugframenum=0;
  Trajout_Single DebugTraj;
  DebugTraj.PrepareTrajWrite("debugtraj.nc",ArgList(),(Topology*)&topIn,
                             CoordinateInfo(), BB_dihedrals_.size()*max_factor_,
                             TrajectoryFile::AMBERNETCDF);
  DebugTraj.WriteSingle(debugframenum++,currentFrame);
# endif
  int next_resnum;
  int bestLoop = 0;
  int number_of_rotations = 0;
  // Set max number of rotations to try.
  int max_rotations = (int)BB_dihedrals_.size();
  max_rotations *= max_factor_;

  // Loop over all dihedrals
  std::vector<PermuteDihedralsType>::const_iterator next_dih = BB_dihedrals_.begin();
  next_dih++;
  for (std::vector<PermuteDihedralsType>::const_iterator dih = BB_dihedrals_.begin();
                                                     dih != BB_dihedrals_.end(); 
                                                     ++dih, ++next_dih)
  {
    ++number_of_rotations;
    // Get the residue atom of the next dihedral. Residues up to and
    // including this residue will be checked for bad clashes 
    if (next_dih != BB_dihedrals_.end()) 
      next_resnum = next_dih->resnum;
    else
      next_resnum = dih->resnum - 1;
    // Set axis of rotation
    Vec3 axisOfRotation = currentFrame.SetAxisOfRotation(dih->atom1, dih->atom2);
    // Generate random value to rotate by in radians
    // Guaranteed to rotate by at least 1 degree.
    // NOTE: could potentially rotate 360 - prevent?
    // FIXME: Just use 2PI and rn_gen, get everything in radians
    double theta_in_degrees = ((int)(RN_.rn_gen()*100000) % 360) + 1;
    double theta_in_radians = theta_in_degrees * Constants::DEGRAD;
    // Calculate rotation matrix for random theta
    rotationMatrix.CalcRotationMatrix(axisOfRotation, theta_in_radians);
    int loop_count = 0;
    double clash = 0;
    double bestClash = 0;
    if (debug_>0) mprintf("DEBUG: Rotating dihedral %zu res %8i:\n", dih - BB_dihedrals_.begin(),
                          dih->resnum+1);
    bool rotate_dihedral = true;
    while (rotate_dihedral) {
      if (debug_>0) {
        mprintf("\t%8i %12s %12s, +%.2lf degrees (%i).\n",dih->resnum+1,
                topIn.AtomMaskName(dih->atom1).c_str(),
                topIn.AtomMaskName(dih->atom2).c_str(),
                theta_in_degrees,loop_count);
      }
      // Rotate around axis
      currentFrame.Rotate(rotationMatrix, dih->Rmask);
#     ifdef DEBUG_PERMUTEDIHEDRALS
      // DEBUG
      DebugTraj.WriteSingle(debugframenum++,currentFrame);
#     endif
      // If we dont care about sterics exit here
      if (!check_for_clashes_) break;
      // Check resulting structure for issues
      int checkresidue;
      if (!checkAllResidues_)
        checkresidue = CheckResidue(currentFrame, topIn, *dih, next_resnum, clash);
      else
        checkresidue = CheckResidue(currentFrame, topIn, *dih, topIn.Nres(), clash);
      if (checkresidue==0)
        rotate_dihedral = false;
      else if (checkresidue==-1) {
        if (dih - BB_dihedrals_.begin() < 2) {
          mprinterr("Error: Cannot backtrack; initial structure already has clashes.\n");
          number_of_rotations = max_rotations + 1;
        } else {
          dih--; //  0
          dih--; // -1
          next_dih = dih;
          next_dih++;
          if (debug_>0)
            mprintf("\tCannot resolve clash with further rotations, trying previous again.\n");
        }
        break;
      }
      if (clash > bestClash) {bestClash = clash; bestLoop = loop_count;}
      //n_problems = CheckResidues( currentFrame, second_atom );
      //if (n_problems > -1) {
      //  mprintf("%i\tCheckResidues: %i problems.\n",frameNum,n_problems);
      //  rotate_dihedral = false;
      //} else if (loop_count==0) {
      if (loop_count==0 && rotate_dihedral) {
        if (debug_>0)
          mprintf("\tTrying dihedral increments of +%i\n",increment_);
        // Instead of a new random dihedral, try increments
        theta_in_degrees = (double)increment_;
        theta_in_radians = theta_in_degrees * Constants::DEGRAD;
        // Calculate rotation matrix for new theta
        rotationMatrix.CalcRotationMatrix(axisOfRotation, theta_in_radians);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // Advance to <natom> !NATOM
  int natom = FindTag(tag, "!NATOM", 6, infile);
  if (debug_>0) mprintf("\tPSF: !NATOM tag found, natom=%i\n", natom);
  // If no atoms, probably issue with PSF file
  if (natom < 1) {
    mprinterr("Error: No atoms in PSF file.\n");
    return 1;
  }
  // Read the next natom lines
  int psfresnum = 0;
  char psfresname[6];
  char psfname[6];
  char psftype[6];
  double psfcharge;
  double psfmass;
  for (int atom=0; atom < natom; atom++) {
    if ( (buffer=infile.NextLine()) == 0 ) {
      mprinterr("Error: ReadParmPSF(): Reading atom %i\n",atom+1);
      return 1;
    }
    // Read line
    // ATOM# SEGID RES# RES ATNAME ATTYPE CHRG MASS (REST OF COLUMNS ARE LIKELY FOR CMAP AND CHEQ)
    sscanf(buffer,"%*i %*s %i %s %s %s %lf %lf",&psfresnum, psfresname, 
           psfname, psftype, &psfcharge, &psfmass);
    parmOut.AddTopAtom( Atom( psfname, psfcharge, psfmass, psftype), 
                        Residue( psfresname, psfresnum, ' ', ' '), 0 );
  } // END loop over atoms 
  // Advance to <nbond> !NBOND
  int bondatoms[9];
  int nbond = FindTag(tag, "!NBOND", 6, infile);
  if (nbond > 0) {
    if (debug_>0) mprintf("\tPSF: !NBOND tag found, nbond=%i\n", nbond);
    int nlines = nbond / 4;
    if ( (nbond % 4) != 0) nlines++;
    for (int bondline=0; bondline < nlines; bondline++) {
      if ( (buffer=infile.NextLine()) == 0 ) {
        mprinterr("Error: ReadParmPSF(): Reading bond line %i\n",bondline+1);
        return 1;
      }
      // Each line has 4 pairs of atom numbers
      int nbondsread = sscanf(buffer,"%i %i %i %i %i %i %i %i",bondatoms,bondatoms+1,
                              bondatoms+2,bondatoms+3, bondatoms+4,bondatoms+5,
                              bondatoms+6,bondatoms+7);
      // NOTE: Charmm atom nums start from 1
      for (int bondidx=0; bondidx < nbondsread; bondidx+=2)
        parmOut.AddBond(bondatoms[bondidx]-1, bondatoms[bondidx+1]-1);
    }
  } else
    mprintf("Warning: PSF has no bonds.\n");
  // Advance to <nangles> !NTHETA
  int nangle = FindTag(tag, "!NTHETA", 7, infile);
  if (nangle > 0) {
    if (debug_>0) mprintf("\tPSF: !NTHETA tag found, nangle=%i\n", nangle);
    int nlines = nangle / 3;
    if ( (nangle % 3) != 0) nlines++;
    for (int angleline=0; angleline < nlines; angleline++) {
      if ( (buffer=infile.NextLine()) == 0) {
        mprinterr("Error: Reading angle line %i\n", angleline+1);
        return 1;
      }
      // Each line has 3 groups of 3 atom numbers
      int nanglesread = sscanf(buffer,"%i %i %i %i %i %i %i %i %i",bondatoms,bondatoms+1,
                              bondatoms+2,bondatoms+3, bondatoms+4,bondatoms+5,
                              bondatoms+6,bondatoms+7, bondatoms+8);
      for (int angleidx=0; angleidx < nanglesread; angleidx += 3)
        parmOut.AddAngle( bondatoms[angleidx  ]-1,
                          bondatoms[angleidx+1]-1,
                          bondatoms[angleidx+2]-1 );
    }
  } else
    mprintf("Warning: PSF has no angles.\n");
  // Advance to <ndihedrals> !NPHI
  int ndihedral = FindTag(tag, "!NPHI", 5, infile);
  if (ndihedral > 0) {
    if (debug_>0) mprintf("\tPSF: !NPHI tag found, ndihedral=%i\n", ndihedral);
    int nlines = ndihedral / 2;
    if ( (ndihedral % 2) != 0) nlines++;
    for (int dihline = 0; dihline < nlines; dihline++) {
      if ( (buffer=infile.NextLine()) == 0) {
        mprinterr("Error: Reading dihedral line %i\n", dihline+1);
        return 1;
      }
      // Each line has 2 groups of 4 atom numbers
      int ndihread = sscanf(buffer,"%i %i %i %i %i %i %i %i",bondatoms,bondatoms+1,
                              bondatoms+2,bondatoms+3, bondatoms+4,bondatoms+5,
                              bondatoms+6,bondatoms+7);
      for (int dihidx=0; dihidx < ndihread; dihidx += 4)
        parmOut.AddDihedral( bondatoms[dihidx  ]-1,
                             bondatoms[dihidx+1]-1,
                             bondatoms[dihidx+2]-1,
                             bondatoms[dihidx+3]-1 );
    }
  } else
    mprintf("Warning: PSF has no dihedrals.\n");
  mprintf("\tPSF contains %i atoms, %i residues.\n", parmOut.Natom(), parmOut.Nres());

  infile.CloseFile();

  return 0;
}

/** Perform setup required for per residue rmsd calculation.
  * Need to set up a target mask, reference mask, and dataset for each
  * residue specified in ResRange.
  * NOTE: Residues in the range arguments from user start at 1, internal
  *       res nums start from 0.
  */
int Action_Rmsd::perResSetup(Topology const& currentParm, Topology const& refParm) {
  Range tgt_range; // Selected target residues
  Range ref_range; // Selected reference residues

  // If no target range previously specified do all solute residues
  if (TgtRange_.Empty()) { 
    tgt_range = currentParm.SoluteResidues();
    tgt_range.ShiftBy(1); // To match user range arg which would start from 1
  } else
    tgt_range = TgtRange_;
  // If the reference range is empty, set it to match the target range
  if (RefRange_.Empty()) 
    ref_range = tgt_range;
  else
    ref_range = RefRange_;
  // Check that the number of reference residues matches number of target residues
  if (tgt_range.Size() != ref_range.Size()) {
    mprintf("Warning: Number of target residues %i does not match\n"
            "Warning:   number of reference residues %i.\n",
            tgt_range.Size(), ref_range.Size());
    return 1;
  }

  // Setup a dataset, target mask, and reference mask for each residue.
  int maxNatom = 0;
  Range::const_iterator ref_it = ref_range.begin();
  MetaData md(rmsd_->Meta().Name(), "res");
  md.SetScalarMode( MetaData::M_RMS );
  for (Range::const_iterator tgt_it = tgt_range.begin();
                             tgt_it != tgt_range.end(); ++tgt_it, ++ref_it)
  {
    int tgtRes = *tgt_it;
    int refRes = *ref_it;
    // Check if either the residue num or the reference residue num out of range.
    if ( tgtRes < 1 || tgtRes > currentParm.Nres()) {
      mprintf("Warning: Specified residue # %i is out of range.\n", tgtRes);
      continue;
    }
    if ( refRes < 1 || refRes > refParm.Nres() ) {
      mprintf("Warning: Specified reference residue # %i is out of range.\n", refRes);
      continue;
    }
    // Check if a perResType has been set for this residue # yet.
    perResArray::iterator PerRes;
    for (PerRes = ResidueRMS_.begin(); PerRes != ResidueRMS_.end(); ++PerRes)
      if ( PerRes->data_->Meta().Idx() == tgtRes ) break;
    // If necessary, create perResType for residue
    if (PerRes == ResidueRMS_.end()) {
      perResType p;
      md.SetIdx( tgtRes );
      md.SetLegend( currentParm.TruncResNameNum(tgtRes-1) );
      p.data_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::DOUBLE, md);
      if (p.data_ == 0) {
        mprinterr("Internal Error: Could not set up per residue data set.\n");
        return 2;
      }
      if (perresout_ != 0) perresout_->AddDataSet( p.data_ );
      // Setup mask strings. Note that masks are based off user residue nums
      p.tgtResMask_.SetMaskString(":" + integerToString(tgtRes) + perresmask_);
      p.refResMask_.SetMaskString(":" + integerToString(refRes) + perresmask_);
      ResidueRMS_.push_back( p );
      PerRes = ResidueRMS_.end() - 1;
    }
    PerRes->isActive_ = false;
    // Setup the reference mask
    if (refParm.SetupIntegerMask(PerRes->refResMask_)) {
      mprintf("Warning: Could not setup reference mask for residue %i\n",refRes);
      continue;
    }
    if (PerRes->refResMask_.None()) {
      mprintf("Warning: No atoms selected for reference residue %i\n",refRes);
      continue;
    }
    // Setup the target mask
    if (currentParm.SetupIntegerMask(PerRes->tgtResMask_)) {
      mprintf("Warning: Could not setup target mask for residue %i\n",tgtRes);
      continue;
    }
    if (PerRes->tgtResMask_.None()) {
      mprintf("Warning: No atoms selected for target residue %i\n",tgtRes);
      continue;
    }
    // Check that # atoms in target and reference masks match
    if (PerRes->tgtResMask_.Nselected() != PerRes->refResMask_.Nselected()) {
      mprintf("Warning: Res %i: # atoms in Tgt (%i) != # atoms in Ref (%i)\n",
              tgtRes, PerRes->tgtResMask_.Nselected(), PerRes->refResMask_.Nselected());
      if (debug_ > 0) {
        mprintf("    Target Atoms:\n");
        for (AtomMask::const_iterator t = PerRes->tgtResMask_.begin();
                                      t != PerRes->tgtResMask_.end(); ++t)
          mprintf("\t%s\n", currentParm.AtomMaskName(*t).c_str());
        mprintf("    Ref Atoms:\n");
        for (AtomMask::const_iterator r = PerRes->refResMask_.begin();
                                      r != PerRes->refResMask_.end(); ++r)
//.........这里部分代码省略.........

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