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

// Action_CheckStructure::SeparateSetup()
int Action_CheckStructure::SeparateSetup(Topology const& top, Box::BoxType btype, bool checkBonds)
{
  image_.SetupImaging( btype );
  bondList_.clear();
  // Set up masks
  if ( top.SetupIntegerMask( Mask1_ ) ) return 1;
  Mask1_.MaskInfo();
  if (Mask1_.None()) {
    mprinterr("Error: Mask '%s' has no atoms.\n", Mask1_.MaskString());
    return 1;
  }
  if (checkBonds) SetupBondList(Mask1_, top);
  if ( Mask2_.MaskStringSet() ) {
    if (top.SetupIntegerMask( Mask2_ ) ) return 1;
    Mask2_.MaskInfo();
    if (Mask2_.None()) {
      mprinterr("Error: Mask '%s' has no atoms.\n", Mask2_.MaskString());
      return 1;
    }
    int common = Mask1_.NumAtomsInCommon( Mask2_ );
    if (common > 0)
      mprintf("Warning: '%s' has %i atoms in common with '%s'. Some problems may be reported\n"
              "Warning:   more than once.\n", Mask1_.MaskString(), common, Mask2_.MaskString());
    // Outer mask should be the one with the most atoms.
    if ( Mask2_.Nselected() > Mask1_.Nselected() ) {
      OuterMask_ = Mask2_;
      InnerMask_ = Mask1_;
    } else {
      OuterMask_ = Mask1_;
      InnerMask_ = Mask2_;
    }
    if (checkBonds) SetupBondList(Mask2_, top);
  }
  return 0;
}

// ReferenceAction::SetRefMask()
int ReferenceAction::SetRefMask(Topology const& topIn, const char* call) {
  if (topIn.SetupIntegerMask( refMask_ )) return 1;
  mprintf("\tReference mask:");
  refMask_.BriefMaskInfo();
  mprintf("\n");
  if (refMask_.None()) {
    mprinterr("Error: %s: No reference atoms selected for parm %s, [%s]\n",
              call, topIn.c_str(), refMask_.MaskString());
    return 1;
  }
  selectedRef_.SetupFrameFromMask( refMask_, topIn.Atoms() );
  return 0;
}

/** Based on selected atoms in Mask1 (and optionally Mask2), set up
  * potential contact lists. Also set up atom/residue indices corresponding
  * to each potential contact.
  */
int Action_NativeContacts::SetupContactLists(Topology const& parmIn, Frame const& frameIn)
{
  // Setup first contact list
  if ( parmIn.SetupIntegerMask( Mask1_, frameIn ) ) return 1;
  if (!includeSolvent_) removeSelectedSolvent( parmIn, Mask1_ );
  Mask1_.MaskInfo();
  if (Mask1_.None())
  {
    mprinterr("Warning: Nothing selected for '%s'\n", Mask1_.MaskString());
    return 1;
  }
  if (debug_ > 0) DebugContactList( Mask1_, parmIn );
  contactIdx1_ = SetupContactIndices( Mask1_, parmIn );
  // Setup second contact list if necessary
  if ( Mask2_.MaskStringSet() ) {
    if (parmIn.SetupIntegerMask( Mask2_, frameIn ) ) return 1;
    if (!includeSolvent_) removeSelectedSolvent( parmIn, Mask2_ );
    Mask2_.MaskInfo();
    if (Mask2_.None())
    {
      mprinterr("Warning: Nothing selected for '%s'\n", Mask2_.MaskString());
      return 1;
    }
    // Warn if masks overlap
    int nOverlap = Mask1_.NumAtomsInCommon( Mask2_ );
    if (nOverlap > 0) {
      mprintf("Warning: Masks '%s' and '%s' overlap by %i atoms.\n"
              "Warning: Some contacts may be double-counted.\n", 
              Mask1_.MaskString(), Mask2_.MaskString(), nOverlap);
      if (mindist_ != 0)
        mprintf("Warning: Minimum distance will always be 0.0\n");
    }
    if (debug_ > 0) DebugContactList( Mask2_, parmIn );
    contactIdx2_ = SetupContactIndices( Mask2_, parmIn );
  }
  return 0;
}

Exec::RetType Exec_ParmStrip::Execute(CpptrajState& State, ArgList& argIn) {
    Topology* parm = State.DSL().GetTopByIndex( argIn );
  if (parm == 0) return CpptrajState::ERR;
  // Check if this topology has already been used to set up an input
  // trajectory, as this will break the traj read.
  bool topology_in_use = false;
  const char* fname = 0;
  for (TrajinList::trajin_it tIn = State.InputTrajList().trajin_begin();
                             tIn != State.InputTrajList().trajin_end(); ++tIn)
    if ( (*tIn)->Traj().Parm() == parm ) {
      topology_in_use = true;
      fname = (*tIn)->Traj().Filename().full();
      break;
    }
  if (!topology_in_use) {
    for (TrajinList::ensemble_it eIn = State.InputTrajList().ensemble_begin();
                                 eIn != State.InputTrajList().ensemble_end(); ++eIn)
      if ( (*eIn)->Traj().Parm() == parm ) {
        topology_in_use = true;
        fname = (*eIn)->Traj().Filename().full();
        break;
      }
  }
  if (topology_in_use) {
    mprinterr("Error: Topology '%s' has already been used to set up trajectory '%s'.\n"
              "Error:   To strip this topology use the 'strip' action.\n",
              parm->c_str(), fname);
    return CpptrajState::ERR;
  }
  AtomMask tempMask( argIn.GetMaskNext() );
  // Since want to keep atoms outside mask, invert selection
  tempMask.InvertMaskExpression();
  if (parm->SetupIntegerMask( tempMask )) return CpptrajState::ERR;
  mprintf("\tStripping atoms in mask [%s] (%i) from %s\n",tempMask.MaskString(),
           parm->Natom() - tempMask.Nselected(), parm->c_str());
  Topology* tempParm = parm->modifyStateByMask(tempMask);
  if (tempParm==0) {
    mprinterr("Error: %s: Could not strip parm.\n", argIn.Command());
    return CpptrajState::ERR;
  } else {
    // Replace parm with stripped version
    *parm = *tempParm;
    parm->Brief("Stripped parm:");
    delete tempParm;
  }
  return CpptrajState::OK;
}

/** Process a mask from the command line. */
int Cpptraj::ProcessMask( Sarray const& topFiles, Sarray const& refFiles,
                          std::string const& maskexpr,
                          bool verbose, bool residue ) const
{
  SetWorldSilent(true);
  if (topFiles.empty()) {
    mprinterr("Error: No topology file specified.\n");
    return 1;
  }
  ParmFile pfile;
  Topology parm;
  if (pfile.ReadTopology(parm, topFiles[0], State_.Debug())) return 1;
  if (!refFiles.empty()) {
    DataSet_Coords_REF refCoords;
    if (refCoords.LoadRefFromFile( refFiles[0], parm, State_.Debug())) return 1;
    parm.SetDistMaskRef( refCoords.RefFrame() );
  }
  if (!verbose) {
    AtomMask tempMask( maskexpr );
    if (parm.SetupIntegerMask( tempMask )) return 1;
    loudPrintf("Selected=");
    if (residue) {
      int res = -1;
      for (AtomMask::const_iterator atom = tempMask.begin(); 
                                    atom != tempMask.end(); ++atom)
      {
        if (parm[*atom].ResNum() > res) {
          loudPrintf(" %i", parm[*atom].ResNum()+1);
          res = parm[*atom].ResNum();
        }
      }
    } else
      for (AtomMask::const_iterator atom = tempMask.begin();
                                    atom != tempMask.end(); ++atom)
        loudPrintf(" %i", *atom + 1);
    loudPrintf("\n");
  } else {
    if (residue)
      parm.PrintResidueInfo( maskexpr );
    else
      parm.PrintAtomInfo( maskexpr );
  }
  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;
}

/** Like the strip action, closest will modify the current parm keeping info
  * for atoms in mask plus the closestWaters solvent molecules. Set up the
  * vector of MolDist objects, one for every solvent molecule in the original
  * parm file. Atom masks for each solvent molecule will be set up.
  */
Action_Closest::RetType Action_Closest::Setup(Topology const& topIn, CoordinateInfo const& cInfoIn)
{
  // If there are no solvent molecules this action is not valid.
  if (topIn.Nsolvent()==0) {
    mprintf("Warning: Parm %s does not contain solvent.\n",topIn.c_str());
    return Action_Closest::SKIP;
  }
  // If # solvent to keep >= solvent in this parm the action is not valid.
  if (closestWaters_ >= topIn.Nsolvent()) {
    mprintf("Warning: # solvent to keep (%i) >= # solvent molecules in '%s' (%i)\n",
            closestWaters_, topIn.c_str(), topIn.Nsolvent());
    return Action_Closest::SKIP;
  }
  image_.SetupImaging( cInfoIn.TrajBox().Type() );
  if (image_.ImagingEnabled())
    mprintf("\tDistances will be imaged.\n");
  else
    mprintf("\tImaging off.\n"); 
  // LOOP OVER MOLECULES
  // 1: Check that all solvent molecules contain same # atoms. Solvent 
  //    molecules must be identical for the command to work properly; 
  //    the prmtop strip occurs only once so the solvent params become fixed.
  // 2: Set up a mask for all solvent molecules.
  SolventMols_.clear();
  // NOTE: May not be necessary to init 'solvent'
  MolDist solvent;
  solvent.D = 0.0;
  solvent.mol = 0;
  SolventMols_.resize(topIn.Nsolvent(), solvent);
  std::vector<MolDist>::iterator mdist = SolventMols_.begin();
  // 3: Set up the soluteMask for all non-solvent molecules.
  int molnum = 1;
  int nclosest = 0;
  int NsolventAtoms = -1;
  for (Topology::mol_iterator Mol = topIn.MolStart();
                              Mol != topIn.MolEnd(); ++Mol)
  {
    if ( Mol->IsSolvent() ) {
      // Solvent, check for same # of atoms.
      if (NsolventAtoms == -1)
        NsolventAtoms = Mol->NumAtoms();
      else if ( NsolventAtoms != Mol->NumAtoms() ) {
        mprinterr("Error: Solvent molecules in '%s' are not of uniform size.\n"
                  "Error:   First solvent mol = %i atoms, solvent mol %i = %i atoms.\n",
                  topIn.c_str(), NsolventAtoms, molnum, (*Mol).NumAtoms());
        return Action_Closest::ERR;
      }
      // mol here is the output molecule number which is why it starts from 1.
      mdist->mol = molnum;
      // Solvent molecule mask
      mdist->mask.AddAtomRange( Mol->BeginAtom(), Mol->EndAtom() );
      // Atoms in the solvent molecule to actually calculate distances to.
      if (firstAtom_) {
        mdist->solventAtoms.assign(1, Mol->BeginAtom() );
      } else {
        mdist->solventAtoms.clear();
        mdist->solventAtoms.reserve( Mol->NumAtoms() );
        for (int svatom = Mol->BeginAtom(); svatom < Mol->EndAtom(); svatom++)
          mdist->solventAtoms.push_back( svatom );
      }
      if (debug_ > 0) {
        mprintf("DEBUG:\tSet up mol %i:", mdist->mol); // DEBUG
        mdist->mask.PrintMaskAtoms("solvent"); // DEBUG
        mprintf("\n"); // DEBUG
      }
      ++mdist;
    }
    ++molnum;
  }

  // Setup distance atom mask
  // NOTE: Should ensure that no solvent atoms are selected!
  if ( topIn.SetupIntegerMask(distanceMask_) ) return Action_Closest::ERR;
  if (distanceMask_.None()) {
    mprintf("Warning: Distance mask '%s' contains no atoms.\n",
            distanceMask_.MaskString());
    return Action_Closest::SKIP;
  }
  distanceMask_.MaskInfo();

  // Check the total number of solvent atoms to be kept.
  NsolventAtoms *= closestWaters_;
  mprintf("\tKeeping %i solvent atoms.\n",NsolventAtoms);
  if (NsolventAtoms < 1) {
    mprintf("Warning: # of solvent atoms to be kept is < 1.\n");
    return Action_Closest::SKIP;
  }
  NsolventMolecules_ = (int)SolventMols_.size();

  return Action_Closest::OK;
}

// Analysis_Modes::Setup()
Analysis::RetType Analysis_Modes::Setup(ArgList& analyzeArgs, AnalysisSetup& setup, int debugIn)
{
  debug_ = debugIn;
  // Analysis type
  if (analyzeArgs.hasKey("fluct"))
    type_ = FLUCT;
  else if (analyzeArgs.hasKey("displ"))
    type_ = DISPLACE;
  else if (analyzeArgs.hasKey("corr"))
    type_ = CORR;
  else if (analyzeArgs.Contains("trajout"))
    type_ = TRAJ;
  else if (analyzeArgs.hasKey("eigenval"))
    type_ = EIGENVAL;
  else if (analyzeArgs.hasKey("rmsip"))
    type_ = RMSIP;
  else {
    mprinterr("Error: No analysis type specified.\n");
    return Analysis::ERR;
  }

  // Get modes name
  std::string modesfile = analyzeArgs.GetStringKey("name");
  if (modesfile.empty()) {
    // Check for deprecated args
    CheckDeprecated(analyzeArgs, modesfile, "file");
    CheckDeprecated(analyzeArgs, modesfile, "stack");
    if (modesfile.empty()) {
      mprinterr("Error: No 'name <modes data set name>' argument given.\n");
      return Analysis::ERR;
    }
  }
  // Get second modes name for RMSIP
  std::string modesfile2 = analyzeArgs.GetStringKey("name2");
  if (type_ == RMSIP) {
    if (modesfile2.empty()) {
      mprinterr("Error: 'rmsip' requires second modes data 'name2 <modes>'\n");
      return Analysis::ERR;
    }
  } else
    modesfile2.clear();
  // Get topology for TRAJ/CORR
  Topology* analyzeParm = setup.DSL().GetTopology( analyzeArgs ); 

  if (type_ == TRAJ ) {
    // Get trajectory format args for projected traj
    beg_ = analyzeArgs.getKeyInt("beg",1) - 1; // Args start at 1
    std::string tOutName = analyzeArgs.GetStringKey("trajout");
    if (tOutName.empty()) {
      mprinterr("Error: Require output trajectory filename, 'trajout <name>'\n");
      return Analysis::ERR;
    }
    TrajectoryFile::TrajFormatType tOutFmt = TrajectoryFile::UNKNOWN_TRAJ;
    if ( analyzeArgs.Contains("trajoutfmt") )
      tOutFmt = TrajectoryFile::GetFormatFromString( analyzeArgs.GetStringKey("trajoutfmt") );
    if (analyzeParm == 0) {
      mprinterr("Error: Could not get topology for output trajectory.\n");
      return Analysis::ERR;
    }
    AtomMask tOutMask( analyzeArgs.GetStringKey("trajoutmask") );
    if ( analyzeParm->SetupIntegerMask( tOutMask ) || tOutMask.None() ) {
      mprinterr("Error: Could not setup output trajectory mask.\n");
      return Analysis::ERR;
    }
    tOutMask.MaskInfo();
    // Strip topology to match mask.
    if (tOutParm_ != 0) delete tOutParm_;
    tOutParm_ = analyzeParm->modifyStateByMask( tOutMask );
    if (tOutParm_ == 0) {
      mprinterr("Error: Could not create topology to match mask.\n");
      return Analysis::ERR;
    }
    // Setup output traj
    if (trajout_.InitTrajWrite( tOutName, ArgList(), tOutFmt ) != 0) {
      mprinterr("Error: Could not init output trajectory.\n");
      return Analysis::ERR;
    }
    // Get min and max for PC
    pcmin_ = analyzeArgs.getKeyDouble("pcmin", -10.0);
    pcmax_ = analyzeArgs.getKeyDouble("pcmax",  10.0);
    if (pcmax_ < pcmin_ || pcmax_ - pcmin_ < Constants::SMALL) {
      mprinterr("Error: pcmin must be less than pcmax\n");
      return Analysis::ERR;
    }
    tMode_ = analyzeArgs.getKeyInt("tmode", 1);
  } else {
    // Args for everything else
    beg_ = analyzeArgs.getKeyInt("beg",7) - 1; // Args start at 1
    bose_ = analyzeArgs.hasKey("bose");
    calcAll_ = analyzeArgs.hasKey("calcall");
  }
  end_ = analyzeArgs.getKeyInt("end", 50);
  factor_ = analyzeArgs.getKeyDouble("factor",1.0);
  std::string setname = analyzeArgs.GetStringKey("setname");

  // Check if modes name exists on the stack
  modinfo_ = (DataSet_Modes*)setup.DSL().FindSetOfType( modesfile, DataSet::MODES );
  if (modinfo_ == 0) {
    mprinterr("Error: '%s' not found: %s\n", modesfile.c_str(), DataSet_Modes::DeprecateFileMsg);
//.........这里部分代码省略.........

/** Set up variable with value. In this case allow any amount of whitespace,
  * so re-tokenize the original argument line (minus the command).
  */
CpptrajState::RetType
  Control_Set::SetupControl(CpptrajState& State, ArgList& argIn, Varray& CurrentVars)
{
  ArgList remaining = argIn.RemainingArgs();
  size_t pos0 = remaining.ArgLineStr().find_first_of("=");
  if (pos0 == std::string::npos) {
    mprinterr("Error: Expected <var>=<value>\n");
    return CpptrajState::ERR;
  }
  size_t pos1 = pos0;
  bool append = false;
  if (pos0 > 0 && remaining.ArgLineStr()[pos0-1] == '+') {
    pos0--;
    append = true;
  }
  std::string variable = NoWhitespace( remaining.ArgLineStr().substr(0, pos0) );
  if (variable.empty()) {
    mprinterr("Error: No variable name.\n");
    return CpptrajState::ERR;
  }
  ArgList equals( NoLeadingWhitespace(remaining.ArgLineStr().substr(pos1+1)) );
  std::string value;
  if (equals.Contains("inmask")) {
    AtomMask mask( equals.GetStringKey("inmask") );
    Topology* top = State.DSL().GetTopByIndex( equals );
    if (top == 0) return CpptrajState::ERR;
    if (top->SetupIntegerMask( mask )) return CpptrajState::ERR;
    if (equals.hasKey("atoms"))
      value = integerToString( mask.Nselected() );
    else if (equals.hasKey("residues")) {
      int curRes = -1;
      int nres = 0;
      for (AtomMask::const_iterator at = mask.begin(); at != mask.end(); ++at) {
        int res = (*top)[*at].ResNum();
        if (res != curRes) {
          nres++;
          curRes = res;
        }
      }
      value = integerToString( nres );
    } else if (equals.hasKey("molecules")) {
      int curMol = -1;
      int nmol = 0;
      for (AtomMask::const_iterator at = mask.begin(); at != mask.end(); ++at) {
        int mol = (*top)[*at].MolNum();
        if (mol != curMol) {
          nmol++;
          curMol = mol;
        }
      }
      value = integerToString( nmol );
    } else {
      mprinterr("Error: Expected 'atoms', 'residues', or 'molecules'.\n");
      return CpptrajState::ERR;
    }
  } else if (equals.hasKey("trajinframes")) {
    value = integerToString(State.InputTrajList().MaxFrames());
  } else
    value = equals.ArgLineStr();
  if (append)
    CurrentVars.AppendVariable( "$" + variable, value );
  else
    CurrentVars.UpdateVariable( "$" + variable, value );
  mprintf("\tVariable '%s' set to '%s'\n", variable.c_str(), value.c_str());
  for (int iarg = 0; iarg < argIn.Nargs(); iarg++)
    argIn.MarkArg( iarg );
  return CpptrajState::OK;
}

/** Set up each mask/integer loop. */
int ControlBlock_For::SetupBlock(CpptrajState& State, ArgList& argIn) {
  mprintf("    Setting up 'for' loop.\n");
  Vars_.clear();
  Topology* currentTop = 0;
  static const char* TypeStr[] = { "ATOMS ", "RESIDUES ", "MOLECULES ",
                                   "MOL_FIRST_RES ", "MOL_LAST_RES " };
  static const char* OpStr[] = {"+=", "-=", "<", ">"};
  description_.assign("for (");
  int MaxIterations = -1;
  int iarg = 0;
  while (iarg < argIn.Nargs())
  {
    // Advance to next unmarked argument.
    while (iarg < argIn.Nargs() && argIn.Marked(iarg)) iarg++;
    if (iarg == argIn.Nargs()) break;
    // Determine 'for' type
    ForType ftype = UNKNOWN;
    bool isMaskFor = true;
    int argToMark = iarg;
    if      ( argIn[iarg] == "atoms"       ) ftype = ATOMS;
    else if ( argIn[iarg] == "residues"    ) ftype = RESIDUES;
    else if ( argIn[iarg] == "molecules"   ) ftype = MOLECULES;
    else if ( argIn[iarg] == "molfirstres" ) ftype = MOLFIRSTRES;
    else if ( argIn[iarg] == "mollastres"  ) ftype = MOLLASTRES;
    else if ( argIn[iarg].find(";") != std::string::npos ) {
      isMaskFor = false;
      ftype = INTEGER;
    }
    // If type is still unknown, check for list.
    if (ftype == UNKNOWN) {
      if (iarg+1 < argIn.Nargs() && argIn[iarg+1] == "in") {
        ftype = LIST;
        isMaskFor = false;
        argToMark = iarg+1;
      }
    }
    // Exit if type could not be determined.
    if (ftype == UNKNOWN) {
      mprinterr("Error: for loop type not specfied.\n");
      return 1;
    }
    argIn.MarkArg(argToMark);
    Vars_.push_back( LoopVar() );
    LoopVar& MH = Vars_.back();
    int Niterations = -1;
    // Set up for specific type
    if (description_ != "for (") description_.append(", ");
    // -------------------------------------------
    if (isMaskFor)
    {
      // {atoms|residues|molecules} <var> inmask <mask> [TOP KEYWORDS]
      if (argIn[iarg+2] != "inmask") {
        mprinterr("Error: Expected 'inmask', got %s\n", argIn[iarg+2].c_str());
        return 1;
      }
      AtomMask currentMask;
      if (currentMask.SetMaskString( argIn.GetStringKey("inmask") )) return 1;
      MH.varType_ = ftype;
      Topology* top = State.DSL().GetTopByIndex( argIn );
      if (top != 0) currentTop = top;
      if (currentTop == 0) return 1;
      MH.varname_ = argIn.GetStringNext();
      if (MH.varname_.empty()) {
        mprinterr("Error: 'for inmask': missing variable name.\n");
        return 1;
      }
      MH.varname_ = "$" + MH.varname_;
      // Set up mask
      if (currentTop->SetupIntegerMask( currentMask )) return 1;
      currentMask.MaskInfo();
      if (currentMask.None()) return 1;
      // Set up indices
      if (MH.varType_ == ATOMS)
        MH.Idxs_ = currentMask.Selected();
      else if (MH.varType_ == RESIDUES) {
        int curRes = -1;
        for (AtomMask::const_iterator at = currentMask.begin(); at != currentMask.end(); ++at) {
          int res = (*currentTop)[*at].ResNum();
          if (res != curRes) {
            MH.Idxs_.push_back( res );
            curRes = res;
          }
        }
      } else if (MH.varType_ == MOLECULES ||
                 MH.varType_ == MOLFIRSTRES ||
                 MH.varType_ == MOLLASTRES)
      {
        int curMol = -1;
        for (AtomMask::const_iterator at = currentMask.begin(); at != currentMask.end(); ++at) {
          int mol = (*currentTop)[*at].MolNum();
          if (mol != curMol) {
            if (MH.varType_ == MOLECULES)
              MH.Idxs_.push_back( mol );
            else {
              int res;
              if (MH.varType_ == MOLFIRSTRES)
                res = (*currentTop)[ currentTop->Mol( mol ).BeginAtom() ].ResNum();
              else // MOLLASTRES
                res = (*currentTop)[ currentTop->Mol( mol ).EndAtom()-1 ].ResNum();
//.........这里部分代码省略.........

// Analysis_Matrix::Setup()
Analysis::RetType Analysis_Matrix::Setup(ArgList& analyzeArgs, AnalysisSetup& setup, int debugIn)
{
#ifdef NO_MATHLIB
  mprinterr("Error: Compiled without LAPACK routines.\n");
  return Analysis::ERR;
#else
  // Get matrix name
  std::string mname = analyzeArgs.GetStringNext();
  if (mname.empty()) {
    mprinterr("Error: Missing matrix name (first argument).\n");
    return Analysis::ERR;
  }
  // Find matrix in DataSetList.
  matrix_ = (DataSet_2D*)setup.DSL().FindSetOfType( mname, DataSet::MATRIX_DBL );
  if (matrix_ == 0)
    matrix_ = (DataSet_2D*)setup.DSL().FindSetOfType( mname, DataSet::MATRIX_FLT );
  if (matrix_ == 0) {
    mprinterr("Error: Could not find matrix named %s\n",mname.c_str());
    return Analysis::ERR;
  }
  // Check that matrix is symmetric (half-matrix incl. diagonal).
  if (matrix_->MatrixKind() != DataSet_2D::HALF) {
    mprinterr("Error: Only works for symmetric matrices (i.e. no mask2)\n");
    return Analysis::ERR;
  }
  
  // nmwiz flag
  nmwizopt_ = analyzeArgs.hasKey("nmwiz");
  if (nmwizopt_) { 
    nmwizvecs_ = analyzeArgs.getKeyInt("nmwizvecs", 20);
    if (nmwizvecs_ < 1) {
      mprinterr("Error: nmwizvecs must be >= 1\n");
      return Analysis::ERR;
    }
    nmwizfile_ = setup.DFL().AddCpptrajFile(analyzeArgs.GetStringKey("nmwizfile"), "NMwiz output",
                                       DataFileList::TEXT, true);
    Topology* parmIn = setup.DSL().GetTopology( analyzeArgs); // TODO: Include with matrix
    if (parmIn == 0) {
      mprinterr("Error: nmwiz: No topology specified.\n");
      return Analysis::ERR;
    }
    AtomMask nmwizMask( analyzeArgs.GetStringKey("nmwizmask") );
    if (parmIn->SetupIntegerMask( nmwizMask )) return Analysis::ERR;
    nmwizMask.MaskInfo();
    Topology* nparm = parmIn->partialModifyStateByMask( nmwizMask );
    if (nparm == 0) return Analysis::ERR;
    nmwizParm_ = *nparm;
    delete nparm;
    nmwizParm_.Brief("nmwiz topology");
  }
  
  // Filenames
  DataFile* outfile = setup.DFL().AddDataFile( analyzeArgs.GetStringKey("out"), analyzeArgs);
  // Thermo flag
  thermopt_ = analyzeArgs.hasKey("thermo");
  if (thermopt_) {
    outthermo_ = setup.DFL().AddCpptrajFile(analyzeArgs.GetStringKey("outthermo"), "'thermo' output",
                                       DataFileList::TEXT, true);
    if (outthermo_ == 0) return Analysis::ERR;
  }
  thermo_temp_ = analyzeArgs.getKeyDouble("temp", 298.15);
  if (thermopt_ && matrix_->Meta().ScalarType() != MetaData::MWCOVAR) {
    mprinterr("Error: Parameter 'thermo' only works for mass-weighted covariance matrix ('mwcovar').\n");
    return Analysis::ERR;
  }
  // Number of eigenvectors; allow "0" only in case of 'thermo'. -1 means 'All'
  nevec_ = analyzeArgs.getKeyInt("vecs",-1);
  if (nevec_ == 0 && !thermopt_) {
    mprintf("Warning: # of eigenvectors specified is 0 and 'thermo' not specified.\n");
    mprintf("Warning: Specify # eigenvectors with 'vecs <#>'. Setting to All.\n");
    nevec_ = -1;
  }
  // Reduce flag
  reduce_ = analyzeArgs.hasKey("reduce");
  // Set up DataSet_Modes. Set Modes DataSet type to be same as input matrix. 
  MetaData md( analyzeArgs.GetStringKey("name") );
  md.SetScalarType( matrix_->Meta().ScalarType() );
  modes_ = (DataSet_Modes*)setup.DSL().AddSet( DataSet::MODES, md, "Modes" );
  if (modes_==0) return Analysis::ERR;
  if (outfile != 0) outfile->AddDataSet( modes_ );

  // Print Status
  mprintf("    DIAGMATRIX: Diagonalizing matrix %s",matrix_->legend());
  if (outfile != 0)
    mprintf(" and writing modes to %s", outfile->DataFilename().full());
  if (nevec_ > 0)
    mprintf("\n\tCalculating %i eigenvectors.\n", nevec_);
  else if (nevec_ == 0)
    mprintf("\n\tNot calculating eigenvectors.\n");
  else
    mprintf("\n\tCalculating all eigenvectors.\n");
  if (thermopt_)
    mprintf("\tCalculating thermodynamic data at %.2f K, output to %s\n",
            thermo_temp_, outthermo_->Filename().full());
  if (nmwizopt_)
    mprintf("\tWriting %i modes to NMWiz file %s", nmwizvecs_, nmwizfile_->Filename().full());
  if (nevec_>0 && reduce_)
    mprintf("\tEigenvectors will be reduced\n");
  mprintf("\tStoring modes with name: %s\n", modes_->Meta().Name().c_str());
//.........这里部分代码省略.........

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