C++ ActionInit::TrajComm方法代码示例

// Action_Esander::Init()
Action::RetType Action_Esander::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
# ifdef MPI
  trajComm_ = init.TrajComm();
# endif
  SANDER_.SetDebug( debugIn );
  Init_ = init;
  // Get keywords
  outfile_ = init.DFL().AddDataFile( actionArgs.GetStringKey("out"), actionArgs );
  save_forces_ = actionArgs.hasKey("saveforces");
  ReferenceFrame REF = init.DSL().GetReferenceFrame( actionArgs );
  if (REF.error()) return Action::ERR;
  if (!REF.empty()) {
    refFrame_ = REF.Coord();
    currentParm_ = REF.ParmPtr();
  }
  if (SANDER_.SetInput( actionArgs )) return Action::ERR;
  // DataSet name and array
  setname_ = actionArgs.GetStringNext();
  if (setname_.empty())
    setname_ = init.DSL().GenerateDefaultName("ENE");
  Esets_.clear();
  Esets_.resize( (int)Energy_Sander::N_ENERGYTYPES, 0 );

  mprintf("    ESANDER: Calculating energy using Sander.\n");
  mprintf("\tTemporary topology file name is '%s'\n", SANDER_.TopFilename().full());
  if (save_forces_) mprintf("\tSaving force information to frame.\n");
  mprintf("\tReference for initialization");
  if (!REF.empty())
    mprintf(" is '%s'\n", REF.refName());
  else
    mprintf(" will be first frame.\n");
  return Action::OK;
}

Action::RetType Action_CreateCrd::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  // Keywords
  Topology* parm = init.DSL().GetTopology( actionArgs );
  if (parm == 0) {
    mprinterr("Error: createcrd: No parm files loaded.\n");
    return Action::ERR;
  }
  pindex_ = parm->Pindex();
  check_ = !actionArgs.hasKey("nocheck");
  // DataSet
  std::string setname = actionArgs.GetStringNext();
  if (setname == "_DEFAULTCRD_") {
    // Special case: Creation of COORDS DataSet has been requested by an
    //               analysis and should already be present.
    coords_ = (DataSet_Coords_CRD*)init.DSL().FindSetOfType(setname, DataSet::COORDS);
  } else 
    coords_ = (DataSet_Coords_CRD*)init.DSL().AddSet(DataSet::COORDS, setname, "CRD");
  if (coords_ == 0) return Action::ERR;
  // Do not set topology here since it may be modified later.

  mprintf("    CREATECRD: Saving coordinates from Top %s to \"%s\"\n",
          parm->c_str(), coords_->legend());
  if (!check_)
    mprintf("\tNot strictly enforcing that all frames have same # atoms.\n");
# ifdef MPI
  if (init.TrajComm().Size() > 1)
    mprintf("Warning: Synchronization of COORDS data sets over multiple threads is\n"
            "Warning:   experimental and may be slower than reading in via a single\n"
            "Warning:   thread. Users are encouraged to run benchmarks before\n"
            "Warning:   extensive usage.\n");
# endif
  return Action::OK;
}

// Action_Grid::Init()
Action::RetType Action_Grid::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  debug_ = debugIn;
  nframes_ = 0;
  // Get output filename
  std::string filename = actionArgs.GetStringKey("out");
  // Get grid options
  grid_ = GridInit( "GRID", actionArgs, init.DSL() );
  if (grid_ == 0) return Action::ERR;
# ifdef MPI
  if (ParallelGridInit(init.TrajComm(), grid_)) return Action::ERR;
# endif
  // Get extra options
  max_ = actionArgs.getKeyDouble("max", 0.80);
  madura_ = actionArgs.getKeyDouble("madura", 0);
  smooth_ = actionArgs.getKeyDouble("smoothdensity", 0);
  invert_ = actionArgs.hasKey("invert");
  pdbfile_ = init.DFL().AddCpptrajFile(actionArgs.GetStringKey("pdb"),"Grid PDB",DataFileList::PDB,true);
  density_ = actionArgs.getKeyDouble("density",0.033456);
  if (actionArgs.hasKey("normframe")) normalize_ = TO_FRAME;
  else if (actionArgs.hasKey("normdensity")) normalize_ = TO_DENSITY;
  else normalize_ = NONE;
  if (normalize_ != NONE && (smooth_ > 0.0 || madura_ > 0.0)) {
    mprinterr("Error: Normalize options are not compatible with smoothdensity/madura options.\n");
    init.DSL().RemoveSet( grid_ );
    return Action::ERR;
  }
  // Get mask
  std::string maskexpr = actionArgs.GetMaskNext();
  if (maskexpr.empty()) {
    mprinterr("Error: GRID: No mask specified.\n");
    init.DSL().RemoveSet( grid_ );
    return Action::ERR;
  }
  mask_.SetMaskString(maskexpr);

  // Setup output file
  // For backwards compat., if no 'out' assume next string is filename
  if (filename.empty() && actionArgs.Nargs() > 1 && !actionArgs.Marked(1))
    filename = actionArgs.GetStringNext();
  DataFile* outfile = init.DFL().AddDataFile(filename, actionArgs);
  if (outfile != 0) outfile->AddDataSet((DataSet*)grid_);

  // Info
  mprintf("    GRID:\n");
  GridInfo( *grid_ );
  if (outfile != 0) mprintf("\tGrid will be printed to file %s\n", outfile->DataFilename().full());
  mprintf("\tGrid data set: '%s'\n", grid_->legend());
  mprintf("\tMask expression: [%s]\n",mask_.MaskString());
  if (pdbfile_ != 0)
      mprintf("\tPseudo-PDB will be printed to %s\n", pdbfile_->Filename().full());
  if (normalize_ == TO_FRAME)
    mprintf("\tGrid will be normalized by number of frames.\n");
  else if (normalize_ == TO_DENSITY)
    mprintf("\tGrid will be normalized to a density of %g molecules/Ang^3.\n", density_);
  // TODO: print extra options

  return Action::OK;
}

// Action_VelocityAutoCorr::Init()
Action::RetType Action_VelocityAutoCorr::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  if (actionArgs.hasKey("usevelocity")) {
    mprinterr("Error: The 'usevelocity' keyword is deprecated. Velocity information\n"
              "Error:   is now used by default if present. To force cpptraj to use\n"
              "Error:   coordinates to estimate velocities (not recommended) use the\n"
              "Error:   'usecoords' keyword.\n");
    return Action::ERR;
  }
  useVelInfo_ = !actionArgs.hasKey("usecoords");
  if (mask_.SetMaskString( actionArgs.GetMaskNext() )) return Action::ERR;
  DataFile* outfile =  init.DFL().AddDataFile( actionArgs.GetStringKey("out"), actionArgs );
  diffout_ = init.DFL().AddCpptrajFile( actionArgs.GetStringKey("diffout"),
                                        "VAC diffusion constants", DataFileList::TEXT, true );
  maxLag_ = actionArgs.getKeyInt("maxlag", -1);
  tstep_ = actionArgs.getKeyDouble("tstep", 1.0);
  useFFT_ = !actionArgs.hasKey("direct");
  normalize_ = actionArgs.hasKey("norm");
  // Set up output data set
  VAC_ = init.DSL().AddSet(DataSet::DOUBLE, actionArgs.GetStringNext(), "VAC");
  if (VAC_ == 0) return Action::ERR;
  // TODO: This should just be a scalar
  diffConst_ = init.DSL().AddSet(DataSet::DOUBLE,
                                 MetaData(VAC_->Meta().Name(), "D", MetaData::NOT_TS));
  if (diffConst_ == 0) return Action::ERR;
  if (outfile != 0) outfile->AddDataSet( VAC_ );
# ifdef MPI
  trajComm_ = init.TrajComm(); 
  if (trajComm_.Size() > 1 && !useVelInfo_)
    mprintf("\nWarning: When calculating velocities between consecutive frames,\n"
            "\nWarning:   'velocityautocorr' in parallel will not work correctly if\n"
            "\nWarning:   coordinates have been modified by previous actions (e.g. 'rms').\n\n");
  diffConst_->SetNeedsSync( false );
# endif
  mprintf("    VELOCITYAUTOCORR:\n"
          "\tCalculate velocity auto-correlation function for atoms in mask '%s'\n",
          mask_.MaskString());
  if (useVelInfo_)
    mprintf("\tUsing velocity information present in frames.\n");
  else
    mprintf("\tCalculating velocities between consecutive frames from coordinates.\n");
  if (outfile != 0)
    mprintf("\tOutput velocity autocorrelation function '%s' to '%s'\n", VAC_->legend(), 
            outfile->DataFilename().full());
  mprintf("\tWriting diffusion constants to '%s'\n", diffout_->Filename().full());
  if (maxLag_ < 1)
    mprintf("\tMaximum lag will be half total # of frames");
  else
    mprintf("\tMaximum lag is %i frames", maxLag_);
  mprintf(", time step between frames is %f ps\n", tstep_);
  if (useFFT_)
    mprintf("\tUsing FFT to calculate autocorrelation function.\n");
  else
    mprintf("\tUsing direct method to calculate autocorrelation function.\n");
  if (normalize_)
    mprintf("\tNormalizing autocorrelation function to 1.0\n");
  return Action::OK;
}

// Action_AtomicCorr::Init()
Action::RetType Action_AtomicCorr::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  debug_ = debugIn;
  outfile_ = init.DFL().AddDataFile( actionArgs.GetStringKey("out"), actionArgs );
  cut_ = actionArgs.getKeyDouble("cut", 0.0);
  if (cut_ < 0.0 || cut_ > 1.0) {
    mprinterr("Error: cut value must be between 0 and 1.\n");
    return Action::ERR;
  }
  min_ = actionArgs.getKeyInt("min",0);
  if (actionArgs.hasKey("byatom"))
    acorr_mode_ = ATOM;
  else if (actionArgs.hasKey("byres"))
    acorr_mode_ = RES;
  mask_.SetMaskString( actionArgs.GetMaskNext() );
  
  // Set up DataSet
  dset_ = init.DSL().AddSet( DataSet::MATRIX_FLT, actionArgs.GetStringNext(), "ACorr" );
  if (dset_ == 0) {
    mprinterr("Error: Could not allocate output data set.\n");
    return Action::ERR;
  }
  // Add DataSet to output file
  if (outfile_ != 0) outfile_->AddDataSet( dset_ );

  mprintf("    ATOMICCORR: Correlation of %s motions will be calculated for\n",
          ModeString[acorr_mode_]);
  mprintf("\tatoms in mask [%s]", mask_.MaskString());
  if (outfile_ != 0) mprintf(", output to file %s", outfile_->DataFilename().full());
  mprintf("\n\tData saved in set '%s'\n", dset_->legend());
  if (cut_ != 0)
    mprintf("\tOnly correlations greater than %.2f or less than -%.2f will be printed.\n",
            cut_,cut_);
  if (min_!=0)
    mprintf("\tOnly correlations for %ss > %i apart will be calculated.\n",
            ModeString[acorr_mode_],min_);
# ifdef MPI
  trajComm_ = init.TrajComm();
  if (trajComm_.Size() > 1)
    mprintf("\nWarning: 'atomiccorr' in parallel will not work correctly if coordinates have\n"
              "Warning:   been modified by previous actions (e.g. 'rms').\n\n");
  // Since matrix calc only happens in Print(), no sync needed.
  dset_->SetNeedsSync( false );
# endif
  return Action::OK;
}

// Action_VelocityAutoCorr::Init()
Action::RetType Action_VelocityAutoCorr::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  useVelInfo_ = actionArgs.hasKey("usevelocity");
  mask_.SetMaskString( actionArgs.GetMaskNext() );
  DataFile* outfile =  init.DFL().AddDataFile( actionArgs.GetStringKey("out"), actionArgs );
  maxLag_ = actionArgs.getKeyInt("maxlag", -1);
  tstep_ = actionArgs.getKeyDouble("tstep", 1.0);
  useFFT_ = !actionArgs.hasKey("direct");
  normalize_ = actionArgs.hasKey("norm");
  // Set up output data set
  VAC_ = init.DSL().AddSet(DataSet::DOUBLE, actionArgs.GetStringNext(), "VAC");
  if (VAC_ == 0) return Action::ERR;
  if (outfile != 0) outfile->AddDataSet( VAC_ );
# ifdef MPI
  trajComm_ = init.TrajComm(); 
  if (trajComm_.Size() > 1 && !useVelInfo_)
    mprintf("\nWarning: When calculating velocities between consecutive frames,\n"
            "\nWarning:   'velocityautocorr' in parallel will not work correctly if\n"
            "\nWarning:   coordinates have been modified by previous actions (e.g. 'rms').\n\n");
# endif
  mprintf("    VELOCITYAUTOCORR:\n"
          "\tCalculate velocity auto-correlation function for atoms in mask '%s'\n",
          mask_.MaskString());
  if (useVelInfo_)
    mprintf("\tUsing velocity information present in frames.\n");
  else
    mprintf("\tCalculating velocities between consecutive frames.\n");
  if (outfile != 0)
    mprintf("\tOutput data set '%s' to '%s'\n", VAC_->legend(), 
            outfile->DataFilename().full());
  if (maxLag_ < 1)
    mprintf("\tMaximum lag will be half total # of frames");
  else
    mprintf("\tMaximum lag is %i frames", maxLag_);
  mprintf(", time step is %f ps\n", tstep_);
  if (useFFT_)
    mprintf("\tUsing FFT to calculate autocorrelation function.\n");
  else
    mprintf("\tUsing direct method to calculate autocorrelation function.\n");
  if (normalize_)
    mprintf("\tNormalizing autocorrelation function to 1.0\n");
  return Action::OK;
}

// Action_GridFreeEnergy::init()
Action::RetType Action_GridFreeEnergy::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  // Get output filename
  DataFile* outfile = init.DFL().AddDataFile(actionArgs.GetStringNext(), actionArgs);
  if (outfile == 0) {
    mprinterr("Error: GridFreeEnergy: no output filename specified.\n");
    return Action::ERR;
  }
  // Get grid options (<nx> <dx> <ny> <dy> <nz> <dz> [box|origin] [negative])
  grid_ = GridInit( "GridFreeEnergy", actionArgs, init.DSL() );
  if (grid_ == 0) return Action::ERR;
# ifdef MPI
  if (ParallelGridInit(init.TrajComm(), grid_)) return Action::ERR;
# endif
  //grid_.PrintXplor( filename_, "", "REMARKS Change in Free energy from bulk solvent with bin normalisation of " + integerToString(currentLargestVoxelOccupancyCount) );

  // Get mask
  std::string maskexpr = actionArgs.GetMaskNext();
  if (maskexpr.empty()) {
    mprinterr("Error: GridFreeEnergy: No mask specified.\n");
    init.DSL().RemoveSet( grid_ );
    return Action::ERR;
  }
  mask_.SetMaskString(maskexpr);

  // Get extra args
  tempInKevin_ = actionArgs.getKeyDouble("temp", 293.0);
  outfile->AddDataSet( grid_ );

  // Info
  mprintf("Warning: DNAIONTRACKER is experimental code!\n");
  mprintf("    GridFreeEnergy\n");
  GridInfo( *grid_ );
  mprintf("\tGrid will be printed to file %s\n",outfile->DataFilename().full());
  mprintf("\tMask expression: [%s]\n",mask_.MaskString());
  mprintf("\tTemp is : %f K\n",tempInKevin_);

  // Allocate grid
  //if (GridAllocate()) return 1;

  return Action::OK;
}

// Action_Density::Init()
Action::RetType Action_Density::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
# ifdef MPI
  trajComm_ = init.TrajComm();
# endif
  DataFile* outfile = init.DFL().AddDataFile(actionArgs.GetStringKey("out"), actionArgs);

  std::string dsname = actionArgs.GetStringKey("name");

  if (actionArgs.hasKey("x") ) {
    axis_ = DX;
    area_coord_[0] = DY;
    area_coord_[1] = DZ;
  } else if (actionArgs.hasKey("y") ) {
    axis_ = DY;
    area_coord_[0] = DX;
    area_coord_[1] = DZ;
  } else if (actionArgs.hasKey("z") ) {
    axis_ = DZ;
    area_coord_[0] = DX;
    area_coord_[1] = DY;
  }

  property_ = NUMBER;
  if      (actionArgs.hasKey("number") )   property_ = NUMBER;
  else if (actionArgs.hasKey("mass") )     property_ = MASS;
  else if (actionArgs.hasKey("charge") )   property_ = CHARGE;
  else if (actionArgs.hasKey("electron") ) property_ = ELECTRON;

  binType_ = CENTER;
  if      (actionArgs.hasKey("bincenter")) binType_ = CENTER;
  else if (actionArgs.hasKey("binedge")  ) binType_ = EDGE;

  delta_ = actionArgs.getKeyDouble("delta", 0.01);
  if (delta_ <= 0) {
    mprinterr("Error: Delta must be > 0.0\n");
    return Action::ERR;
  }

  // for compatibility with ptraj, ignored because we rely on the atom code to
  // do the right thing, see Atom.{h,cpp}
  if (actionArgs.hasKey("efile"))
    mprintf("Warning: The 'efile' keyword is deprecated.\n");

  // read the rest of the command line as a series of masks
  std::string maskstr;

  unsigned int idx = 1;
  while ( (maskstr = actionArgs.GetMaskNext() ) != emptystring) {
    masks_.push_back( AtomMask(maskstr) );
    if (dsname.empty())
      dsname = init.DSL().GenerateDefaultName("DENSITY");
    MetaData MD(dsname, "avg", idx);
    MD.SetTimeSeries( MetaData::NOT_TS );
    // Hold average density
    DataSet* ads = init.DSL().AddSet( DataSet::DOUBLE, MD );
    if (ads == 0) return Action::ERR;
    ads->SetLegend( NoWhitespace(masks_.back().MaskExpression()) );
    AvSets_.push_back( ads );
    if (outfile != 0) outfile->AddDataSet( ads );
    // Hold SD density
    MD.SetAspect("sd");
    DataSet* sds = init.DSL().AddSet( DataSet::DOUBLE, MD );
    if (sds == 0) return Action::ERR;
    sds->SetLegend( NoWhitespace("sd(" + masks_.back().MaskExpression() + ")") );
    SdSets_.push_back( sds );
    if (outfile != 0) outfile->AddDataSet( sds );
#   ifdef MPI
    ads->SetNeedsSync( false ); // Populated in Print()
    sds->SetNeedsSync( false );
#   endif
    idx++;
  }
  if (masks_.empty()) {
    // If no masks assume we want total system density.
    if (dsname.empty())
      dsname = actionArgs.GetStringNext();
    density_ = init.DSL().AddSet(DataSet::DOUBLE, dsname, "DENSITY");
    if (density_ == 0) return Action::ERR;
    if (outfile != 0) outfile->AddDataSet( density_ );
    image_.InitImaging( true );
    // Hijack delta for storing sum of masses
    delta_ = 0.0;
  } else {
    // Density selected by mask(s) along an axis
    density_ = 0;
    histograms_.resize(masks_.size() );
  }

  mprintf("    DENSITY:");
  if (density_ == 0) {
    const char* binStr[] = {"center", "edge"};
    mprintf(" Determining %s density for %zu masks.\n", PropertyStr_[property_], masks_.size());
    mprintf("\troutine version: %s\n", ROUTINE_VERSION_STRING);
    mprintf("\tDelta is %f\n", delta_);
    mprintf("\tAxis is %s\n", AxisStr_[axis_]);
    mprintf("\tData set name is '%s'\n", dsname.c_str());
    mprintf("\tData set aspect [avg] holds the mean, aspect [sd] holds standard deviation.\n");
    mprintf("\tBin coordinates will be to bin %s.\n", binStr[binType_]);
//.........这里部分代码省略.........

// -----------------------------------------------------------------------------
// Action_NMRrst::Init()
Action::RetType Action_NMRrst::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
# ifdef MPI
  if (init.TrajComm().Size() > 1) {
    mprinterr("Error: 'nmrrst' action does not work with > 1 thread (%i threads currently).\n",
              init.TrajComm().Size());
    return Action::ERR;
  }
# endif
  debug_ = debugIn;
  // Get Keywords
  Image_.InitImaging( !(actionArgs.hasKey("noimage")) );
  useMass_ = !(actionArgs.hasKey("geom"));
  findNOEs_ = actionArgs.hasKey("findnoes");
  findOutput_ = init.DFL().AddCpptrajFile(actionArgs.GetStringKey("findout"), "Found NOEs",
                                    DataFileList::TEXT, true);
  specOutput_ = init.DFL().AddCpptrajFile(actionArgs.GetStringKey("specout"), "Specified NOEs",
                                    DataFileList::TEXT, true);
  if (findOutput_ == 0 || specOutput_ == 0) return Action::ERR;
  resOffset_ = actionArgs.getKeyInt("resoffset", 0);
  DataFile* outfile = init.DFL().AddDataFile( actionArgs.GetStringKey("out"), actionArgs );
  max_cut_ = actionArgs.getKeyDouble("cut", 6.0);
  strong_cut_ = actionArgs.getKeyDouble("strongcut", 2.9);
  medium_cut_ = actionArgs.getKeyDouble("mediumcut", 3.5);
  weak_cut_ = actionArgs.getKeyDouble("weakcut", 5.0);
  series_ = actionArgs.hasKey("series");
  std::string rstfilename = actionArgs.GetStringKey("file");
  viewrst_ = actionArgs.GetStringKey("viewrst");
  setname_ = actionArgs.GetStringKey("name");
  if (setname_.empty())
    setname_ = init.DSL().GenerateDefaultName("NMR");
  nframes_ = 0;

  // Atom Mask
  Mask_.SetMaskString( actionArgs.GetMaskNext() );

  // Pairs specified on command line.
  std::string pair1 = actionArgs.GetStringKey("pair");
  while (!pair1.empty()) {
    std::string pair2 = actionArgs.GetStringNext();
    if (pair2.empty()) {
      mprinterr("Error: Only one mask specified for pair (%s)\n", pair1.c_str());
      return Action::ERR;
    }
    Pairs_.push_back( MaskPairType(AtomMask(pair1), AtomMask(pair2)) );
    pair1 = actionArgs.GetStringKey("pair");
  }

  // Check that something will be done
  if (!findNOEs_ && rstfilename.empty() && Pairs_.empty()) {
    mprinterr("Error: Must specify restraint file, 'pair', and/or 'findnoes'.\n");
    return Action::ERR;
  }

  // Read in NMR restraints.
  if (!rstfilename.empty()) {
    if (ReadNmrRestraints( rstfilename )) return Action::ERR;
  }

  // Set up distances.
  int num_noe = 1;
  for (noeDataArray::iterator noe = NOEs_.begin(); noe != NOEs_.end(); ++noe, ++num_noe) {
    // Translate any ambiguous atom names
    TranslateAmbiguous( noe->aName1_ ); 
    TranslateAmbiguous( noe->aName2_ );
    // Create mask expressions from resnum/atom name
    noe->dMask1_.SetMaskString( MaskExpression( noe->resNum1_, noe->aName1_ ) ); 
    noe->dMask2_.SetMaskString( MaskExpression( noe->resNum2_, noe->aName2_ ) ); 
    // Dataset to store distances
    AssociatedData_NOE noeData(noe->bound_, noe->boundh_, noe->rexp_);
    MetaData md(setname_, "NOE", num_noe);
    md.SetLegend( noe->dMask1_.MaskExpression() + " and " + noe->dMask2_.MaskExpression());
    md.SetScalarMode( MetaData::M_DISTANCE );
    md.SetScalarType( MetaData::NOE );
    noe->dist_ = init.DSL().AddSet(DataSet::DOUBLE, md);
    if (noe->dist_==0) return Action::ERR;
    noe->dist_->AssociateData( &noeData );
    // Add dataset to data file
    if (outfile != 0) outfile->AddDataSet( noe->dist_ );
  }

  masterDSL_ = init.DslPtr();
 
  mprintf("Warning: *** THIS ACTION IS EXPERIMENTAL. ***\n");
  mprintf("    NMRRST: %zu NOEs from NMR restraint file.\n", NOEs_.size());
  mprintf("\tShifting residue numbers in restraint file by %i\n", resOffset_);
  // DEBUG - print NOEs
  for (noeDataArray::const_iterator noe = NOEs_.begin(); noe != NOEs_.end(); ++noe)
    mprintf("\t'%s'  %f < %f < %f\n", noe->dist_->legend(),
            noe->bound_, noe->rexp_, noe->boundh_);
  if (findNOEs_) {
    mprintf("\tSearching for potential NOEs. Max cutoff is %g Ang.\n", max_cut_);
    mprintf("\tNOE distance criteria (Ang.): S= %g, M= %g, W= %g\n",
            strong_cut_, medium_cut_, weak_cut_);
    if (series_)
      mprintf("\tDistance data for NOEs less than cutoff will be saved as '%s[foundNOE]'.\n",
              setname_.c_str());
    mprintf("\tFound NOEs will be written to '%s'\n", findOutput_->Filename().full());
//.........这里部分代码省略.........

// Action_ReplicateCell::Init()
Action::RetType Action_ReplicateCell::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  // Require imaging.
  image_.InitImaging( true );
  // Set up output traj
  std::string trajfilename = actionArgs.GetStringKey("out");
  parmfilename_ = actionArgs.GetStringKey("parmout");
  bool setAll = actionArgs.hasKey("all");
  std::string dsname = actionArgs.GetStringKey("name");
  if (!dsname.empty()) {
    coords_ = (DataSet_Coords*)init.DSL().AddSet(DataSet::COORDS, dsname, "RCELL");
    if (coords_ == 0) return Action::ERR;
  }
  if (trajfilename.empty() && coords_ == 0) {
    mprinterr("Error: Either 'out <traj filename> or 'name <dsname>' must be specified.\n");
    return Action::ERR;
  }
  // Get Mask
  Mask1_.SetMaskString( actionArgs.GetMaskNext() );

  // Determine which directions to set
  if (setAll) {
    for (int ix = -1; ix < 2; ix++)
      for (int iy = -1; iy < 2; iy++)
        for (int iz = -1; iz < 2; iz++) {
          directionArray_.push_back( ix );
          directionArray_.push_back( iy );
          directionArray_.push_back( iz );
        }
  } else {
    std::string dirstring = actionArgs.GetStringKey("dir");
    while (!dirstring.empty()) {
      std::vector<int> ixyz(3, -2);
      std::vector<int>::iterator iptr = ixyz.begin();
      for (std::string::const_iterator c = dirstring.begin();
                                       c != dirstring.end(); ++c)
      {
        if (iptr == ixyz.end()) {
          mprinterr("Error: 'dir' string has too many characters.\n");
          return Action::ERR;
        }
        int sign = 1;
        if      (*c == '+') ++c;
        else if (*c == '-') { sign = -1; ++c; }
        
        if (isdigit( *c ))
          *iptr = toDigit( *c ) * sign;
        else {
          mprinterr("Error: illegal character '%c' in 'dir' string '%s'; only numbers allowed.\n",
                    *c, dirstring.c_str());
          return Action::ERR;
        }
        ++iptr;
      }
      //mprintf("DEBUG: %s = %i %i %i\n", dirstring.c_str(), ixyz[0], ixyz[1], ixyz[2]);
      directionArray_.push_back( ixyz[0] );
      directionArray_.push_back( ixyz[1] );
      directionArray_.push_back( ixyz[2] );
      dirstring = actionArgs.GetStringKey("dir");
    }
  }
  ncopies_ = (int)(directionArray_.size() / 3);
  if (ncopies_ < 1) {
    mprinterr("Error: No directions (or 'all') specified.\n");
    return Action::ERR;
  }
  // Initialize output trajectory with remaining arguments
  if (!trajfilename.empty()) {
    outtraj_.SetDebug( debugIn );
    if ( outtraj_.InitEnsembleTrajWrite(trajfilename, actionArgs.RemainingArgs(),
                                        TrajectoryFile::UNKNOWN_TRAJ, init.DSL().EnsembleNum()) )
      return Action::ERR;
    writeTraj_ = true;
#   ifdef MPI
    outtraj_.SetTrajComm( init.TrajComm() );
#   endif
  } else
    writeTraj_ = false;

  mprintf("    REPLICATE CELL: Replicating cell in %i directions:\n", ncopies_);
  mprintf("\t\t X  Y  Z\n");
  for (unsigned int i = 0; i != directionArray_.size(); i += 3)
    mprintf("\t\t%2i %2i %2i\n", directionArray_[i], 
            directionArray_[i+1], directionArray_[i+2]);
  mprintf("\tUsing atoms in mask '%s'\n", Mask1_.MaskString());
  if (writeTraj_)
    mprintf("\tWriting to trajectory %s\n", outtraj_.Traj().Filename().full());
  if (!parmfilename_.empty())
    mprintf("\tWriting topology %s\n", parmfilename_.c_str());
  if (coords_ != 0)
    mprintf("\tSaving coords to data set %s\n", coords_->legend());

  return Action::OK;
}

// Action_CreateReservoir::Init()
Action::RetType Action_CreateReservoir::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
# ifndef BINTRAJ
  mprinterr("Error: NetCDF reservoir requires NetCDF support. Recompile with -DBINTRAJ.\n");
  return Action::ERR;
# endif
# ifdef MPI
  if (init.TrajComm().Size() > 1) {
    mprinterr("Error: 'createreservoir' action does not work with > 1 process (%i processes currently).\n", init.TrajComm().Size());
    return Action::ERR;
  }
# endif
  // Get keywords
  filename_.SetFileName( actionArgs.GetStringNext() );
  if (filename_.empty()) {
    mprinterr("Error: createreservoir: No filename specified.\n");
    return Action::ERR;
  }
  reservoirT_ = actionArgs.getKeyDouble("temp0", -1.0);
  if (reservoirT_ < 0.0) {
    mprinterr("Error: Reservoir temperature must be specified and cannot be < 0.0\n");
    return Action::ERR;
  }
  iseed_ = actionArgs.getKeyInt("iseed", 0);
  if (iseed_ < 1) {
    mprinterr("Error: Reservoir random seed must be specified and > 0\n");
    return Action::ERR;
  }
  useVelocity_ = !actionArgs.hasKey("novelocity");
  useForce_ = !actionArgs.hasKey("noforce");
  // Get parm for reservoir traj
  original_trajparm_ = init.DSL().GetTopology( actionArgs );
  if (original_trajparm_ == 0) {
    mprinterr("Error: createreservoir: no topology.\n");
    return Action::ERR;
  }
  // Get energy data set
  std::string eneDsname = actionArgs.GetStringKey("ene");
  DataSet* dstmp = init.DSL().GetDataSet( eneDsname );
  if (dstmp == 0) {
    mprinterr("Error: could not get energy data set %s\n", eneDsname.c_str());
    return Action::ERR;
  }
  if (dstmp->Type() != DataSet::FLOAT &&
      dstmp->Type() != DataSet::DOUBLE &&
      dstmp->Type() != DataSet::XYMESH)
  {
    mprinterr("Error: energy data set %s must be type FLOAT, DOUBLE, or XYMESH.\n",
              dstmp->legend());
    return Action::ERR;
  }
  if (dstmp->Ndim() != 1) {
    mprinterr("Error: energy data set is not 1D (%zu)\n", dstmp->Ndim());
    return Action::ERR;
  }
  ene_ = static_cast<DataSet_1D*>( dstmp );
  // Get bin data set
  std::string binDSname = actionArgs.GetStringKey("bin");
  if (!binDSname.empty()) {
    dstmp = init.DSL().GetDataSet( binDSname );
    if (dstmp == 0) {
      mprinterr("Error: could not get bin data set %s\n", binDSname.c_str());
      return Action::ERR;
    } else if (dstmp->Ndim() != 1) {
      mprinterr("Error: bin data set must be one dimensional.\n");
      return Action::ERR;
    }
    bin_ = static_cast<DataSet_1D*>( dstmp );
  }
  trajIsOpen_ = false;
  nframes_ = 0;
  // Setup output reservoir file
  reservoir_.SetDebug( debugIn );
  // Set title
  title_ = actionArgs.GetStringKey("title");
  if (title_.empty())
    title_.assign("Cpptraj Generated structure reservoir");

  mprintf("    CREATERESERVOIR: '%s', energy data '%s'", filename_.full(), ene_->legend());
  if (bin_ != 0)
    mprintf(", bin data '%s'", bin_->legend());
  mprintf("\n\tTitle: %s\n", title_.c_str());
  mprintf("\tReservoir temperature= %.2f, random seed= %i\n", reservoirT_, iseed_);
  if (useVelocity_)
    mprintf("\tVelocities will be written to reservoir if present.\n");
  else
    mprintf("\tVelocities will not be written to reservoir.\n");
  if (useForce_)
    mprintf("\tForces will be written to reservoir if present.\n");
  else
    mprintf("\tForces will not be written to reservoir.\n");
  mprintf("\tTopology: %s\n", original_trajparm_->c_str());
  return Action::OK;
}

// Action_Density::Init()
Action::RetType Action_Density::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
# ifdef MPI
  if (init.TrajComm().Size() > 1) {
    mprinterr("Error: 'density' action does not work with > 1 thread (%i threads currently).\n",
              init.TrajComm().Size());
    return Action::ERR;
  }
# endif

  DataFile* outfile = init.DFL().AddDataFile(actionArgs.GetStringKey("out"), actionArgs);

  std::string dsname = actionArgs.GetStringKey("name");
  if (dsname.empty())
    dsname = init.DSL().GenerateDefaultName("DENSITY");

  if (actionArgs.hasKey("x") ) {
    axis_ = DX;
    area_coord_[0] = DY;
    area_coord_[1] = DZ;
  } else if (actionArgs.hasKey("y") ) {
    axis_ = DY;
    area_coord_[0] = DX;
    area_coord_[1] = DZ;
  } else if (actionArgs.hasKey("z") ) {
    axis_ = DZ;
    area_coord_[0] = DX;
    area_coord_[1] = DY;
  }

  property_ = NUMBER;
  if (actionArgs.hasKey("number") )   property_ = NUMBER;
  if (actionArgs.hasKey("mass") )     property_ = MASS;
  if (actionArgs.hasKey("charge") )   property_ = CHARGE;
  if (actionArgs.hasKey("electron") ) property_ = ELECTRON;

  delta_ = actionArgs.getKeyDouble("delta", 0.01);

  // for compatibility with ptraj, ignored because we rely on the atom code to
  // do the right thing, see Atom.{h,cpp}
  if (actionArgs.hasKey("efile"))
    mprintf("Warning: The 'efile' keyword is deprecated.\n");

  // read the rest of the command line as a series of masks
  std::string maskstr;

  unsigned int idx = 1;
  while ( (maskstr = actionArgs.GetMaskNext() ) != emptystring) {
    masks_.push_back( AtomMask(maskstr) );
    MetaData MD(dsname, "avg", idx);
    MD.SetTimeSeries( MetaData::NOT_TS );
    // Hold average density
    DataSet* ads = init.DSL().AddSet( DataSet::DOUBLE, MD );
    if (ads == 0) return Action::ERR;
    ads->SetLegend( NoWhitespace(masks_.back().MaskExpression()) );
    AvSets_.push_back( ads );
    if (outfile != 0) outfile->AddDataSet( ads );
    // Hold SD density
    MD.SetAspect("sd");
    DataSet* sds = init.DSL().AddSet( DataSet::DOUBLE, MD );
    if (sds == 0) return Action::ERR;
    sds->SetLegend( NoWhitespace("sd(" + masks_.back().MaskExpression() + ")") );
    SdSets_.push_back( sds );
    if (outfile != 0) outfile->AddDataSet( sds );
#   ifdef MPI
    ads->SetNeedsSync( false ); // Populated in Print()
    sds->SetNeedsSync( false );
#   endif
    idx++;
  }
  if (masks_.empty()) {
    mprinterr("Error: No masks specified.\n");
    return Action::ERR;
  }

  minus_histograms_.resize(masks_.size() );
  plus_histograms_.resize(masks_.size() );

  mprintf("    DENSITY: Determining %s density for %zu masks.\n", PropertyStr_[property_],
          masks_.size());
  mprintf("\troutine version: %s\n", ROUTINE_VERSION_STRING);
  mprintf("\tDelta is %f\n", delta_);
  mprintf("\tAxis is %s\n", AxisStr_[axis_]);
  mprintf("\tData set name is '%s'\n", dsname.c_str());
  if (outfile != 0)
    mprintf("\tOutput to '%s'\n", outfile->DataFilename().full());

  return Action::OK;
}

// Action_FilterByData::Init()
Action::RetType Action_FilterByData::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  maxmin_ = init.DSL().AddSet( DataSet::INTEGER, actionArgs.GetStringKey("name"), "Filter" );
  if (maxmin_ == 0) return Action::ERR;
  DataFile* maxminfile = init.DFL().AddDataFile( actionArgs.GetStringKey("out"), actionArgs );
  if (maxminfile != 0)
    maxminfile->AddDataSet( maxmin_ );
  // Get min and max args.
  while (actionArgs.Contains("min"))
    Min_.push_back( actionArgs.getKeyDouble("min", 0.0) );
  while (actionArgs.Contains("max"))
    Max_.push_back( actionArgs.getKeyDouble("max", 0.0) );
  if (Min_.empty()) {
    mprinterr("Error: At least one 'min' arg must be specified.\n");
    return Action::ERR;
  }
  if (Max_.empty()) {
    mprinterr("Error: At least one 'max' arg must be specified.\n");
    return Action::ERR;
  }
  if (Min_.size() != Max_.size()) {
    mprinterr("Error: # of 'min' args (%zu) != # of 'max' args (%zu)\n",
              Min_.size(), Max_.size());
    return Action::ERR;
  }
  // Get DataSets from remaining arguments
  Dsets_.AddSetsFromArgs( actionArgs.RemainingArgs(), init.DSL() );

  if (Dsets_.empty()) {
    mprinterr("Error: No data sets specified.\n");
    return Action::ERR;
  }

  if ( Dsets_.size() < Min_.size() ) {
    mprinterr("Error: More 'min'/'max' args (%zu) than data sets (%zu).\n",
              Min_.size(), Dsets_.size());
    return Action::ERR;
  }
  if ( Dsets_.size() > Min_.size() ) {
    unsigned int Nremaining = Dsets_.size() - Min_.size();
    double useMin = Min_.back();
    double useMax = Max_.back();
    mprintf("Warning: More data sets than 'min'/'max' args.\n"
            "Warning:  Using min=%f and max=%f for last %zu data sets.\n",
            useMin, useMax, Nremaining);
    for (unsigned int ds = 0; ds < Nremaining; ++ds) {
      Min_.push_back( useMin );
      Max_.push_back( useMax );
    }
  }

  mprintf("    FILTER: Filtering out frames using %zu data sets.\n", Dsets_.size());
  for (unsigned int ds = 0; ds < Dsets_.size(); ds++)
    mprintf("\t%.4f < '%s' < %.4f\n", Min_[ds], Dsets_[ds]->legend(), Max_[ds]);
  if (maxminfile != 0)
    mprintf("\tFilter frame info will be written to %s\n", maxminfile->DataFilename().full());
# ifdef MPI
  if (init.TrajComm().Size() > 1)
    mprintf("Warning: Trajectories written after 'filter' may have issues if\n"
            "Warning:   the number of threads writing is > 1 (currently %i threads)\n",
            init.TrajComm().Size());
# endif
  return Action::OK;
}

// Action_Outtraj::Init()
Action::RetType Action_Outtraj::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  // Set up output traj
  std::string trajfilename = actionArgs.GetStringNext();
  if (trajfilename.empty()) {
    mprinterr("Error: No filename given.\nError: Usage: ");
    Help();
    return Action::ERR;
  }
  associatedParm_ = init.DSL().GetTopology(actionArgs);
  if (associatedParm_ == 0) {
    mprinterr("Error: Could not get associated topology for %s\n",trajfilename.c_str());
    return Action::ERR;
  }
  std::string rangeArg = actionArgs.GetStringKey("onlymembers");
  if (rangeArg.empty())
    isActive_ = true;
  else {
    Range members;
    if (members.SetRange( rangeArg )) return Action::ERR;
    isActive_ = members.InRange( init.DSL().EnsembleNum() );
  }
  // If maxmin, get the name of the dataset as well as the max and min values.
  double lastmin = 0.0;
  double lastmax = 0.0;
  while ( actionArgs.Contains("maxmin") ) {
    std::string datasetName = actionArgs.GetStringKey("maxmin");
    if (!datasetName.empty()) {
      DataSet* dset = init.DSL().GetDataSet(datasetName);
      if (dset==0) {
        mprintf("Error: maxmin: Could not get dataset %s\n",datasetName.c_str());
        return Action::ERR;
      } else {
        // Currently only allow int, float, or double datasets
        if (dset->Type() != DataSet::INTEGER &&
            dset->Type() != DataSet::FLOAT &&
            dset->Type() != DataSet::DOUBLE) 
        {
          mprinterr("Error: maxmin: Only int, float, or double dataset (%s) supported.\n",
                  datasetName.c_str());
          return Action::ERR;
        }
        Dsets_.push_back( (DataSet_1D*)dset );
        Max_.push_back( actionArgs.getKeyDouble("max",lastmax) );
        Min_.push_back( actionArgs.getKeyDouble("min",lastmin) );
        lastmax = Max_.back();
        lastmin = Min_.back();
      }
    } else {
      mprinterr("Error: maxmin Usage: maxmin <setname> max <max> min <min>\n");
      return Action::ERR;
    }
  }
# ifdef MPI
  trajComm_ = init.TrajComm();
  if (trajComm_.Size() > 1 && !Dsets_.empty()) {
    mprinterr("Error: outtraj 'maxmin' currently does not work when using > 1 thread\n"
              "Error:   to write trajectory (currently %i threads)\n", trajComm_.Size());
    return Action::ERR;
  }
  outtraj_.SetTrajComm( trajComm_ );
# endif
  // Initialize output trajectory with remaining arguments
  if (isActive_) {
    outtraj_.SetDebug(debugIn);
    if ( outtraj_.InitEnsembleTrajWrite(trajfilename, actionArgs.RemainingArgs(),
                                        TrajectoryFile::UNKNOWN_TRAJ, init.DSL().EnsembleNum()) )
      return Action::ERR;
  }
  isSetup_ = false;

  mprintf("    OUTTRAJ: Writing frames associated with topology '%s'\n", associatedParm_->c_str());
  if (!rangeArg.empty())
    mprintf("\tonlymembers: Only writing members %s\n", rangeArg.c_str());
  for (unsigned int ds = 0; ds < Dsets_.size(); ++ds)
    mprintf("\tmaxmin: Printing trajectory frames based on %g <= %s <= %g\n",
            Min_[ds], Dsets_[ds]->legend(), Max_[ds]);

  return Action::OK;
} 

// Action_Matrix::Init()
Action::RetType Action_Matrix::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
# ifdef MPI
  trajComm_ = init.TrajComm();
# endif
  debug_ = debugIn;
  // Get Keywords
  std::string outfilename = actionArgs.GetStringKey("out");
  // Get start/stop/offset
  if (InitFrameCounter(actionArgs)) return Action::ERR;
  // Determine matrix type
  MetaData::scalarType mtype = MetaData::DIST;
  if (actionArgs.hasKey("distcovar"))
    mtype = MetaData::DISTCOVAR;
  else if (actionArgs.hasKey("mwcovar"))
    mtype = MetaData::MWCOVAR;
  else if (actionArgs.hasKey("dist"))
    mtype = MetaData::DIST;
  else if (actionArgs.hasKey("covar"))
    mtype = MetaData::COVAR;
  else if (actionArgs.hasKey("correl"))
    mtype = MetaData::CORREL;
  else if (actionArgs.hasKey("idea"))
    mtype = MetaData::IDEA;
  else if (actionArgs.hasKey("ired"))
    mtype = MetaData::IREDMAT;
  else if (actionArgs.hasKey("dihcovar"))
    mtype = MetaData::DIHCOVAR;
  // Output type
  if (actionArgs.hasKey("byres"))
    outtype_ = BYRESIDUE;
  else if (actionArgs.hasKey("bymask"))
    outtype_ = BYMASK;
  else if (actionArgs.hasKey("byatom"))
    outtype_ = BYATOM;
  else
    outtype_ = BYATOM;
  // Check if output type is valid for matrix type
  if ( outtype_ != BYATOM && (mtype == MetaData::COVAR || 
                              mtype == MetaData::MWCOVAR || 
                              mtype == MetaData::IREDMAT ) )
  {
    mprinterr("Error: matrix: for COVAR, MWCOVAR, or IRED matrix only byatom output possible\n");
    return Action::ERR;
  }
  // Get matrix name
  std::string name = actionArgs.GetStringKey("name");
  // UseMass
  useMass_ = actionArgs.hasKey("mass");
  // NOTE: Determine matrix kind here so subsequent Actions/Analyses know about it.
  DataSet_2D::MatrixKindType mkind = DataSet_2D::HALF;
  if (mtype == MetaData::IREDMAT) { // IRED matrix
    // Setup IRED vectors and determine Legendre order
    order_ = actionArgs.getKeyInt("order",1);
    if (order_ <= 0) {
      mprinterr("Error: matrix: order parameter <= 0, ignoring command\n");
      return Action::ERR;
    }
    for ( DataSetList::const_iterator DS = init.DSL().begin(); DS != init.DSL().end(); ++DS) {
      if ( (*DS)->Type() == DataSet::VECTOR && (*DS)->Meta().ScalarType() == MetaData::IREDVEC )
        IredVectors_.push_back( (DataSet_Vector*)*DS );
    }
    if (IredVectors_.empty()) {
      mprinterr("Error: matrix: no vectors defined for IRED\n");
      return Action::ERR;
    }
  } else if (mtype == MetaData::DIHCOVAR) { // Dihedral Covariance
    // Get data set mask for dihedral covariance
    DihedralSets_.clear();
    DihedralSets_.AddTorsionSets( init.DSL().GetMultipleSets( actionArgs.GetStringKey("dihedrals") ) );
    if ( DihedralSets_.empty() ) {
      mprinterr("Error: No valid data sets found.\n");
      return Action::ERR;
    }
  } else {
    // Get masks if not IRED/DIHCOVAR
    mask1_.SetMaskString( actionArgs.GetMaskNext() );
    std::string maskexpr = actionArgs.GetMaskNext();
    if (!maskexpr.empty()) useMask2_ = true;
    if ( useMask2_ && (mtype == MetaData::IDEA || mtype == MetaData::DISTCOVAR) )
    {
      mprinterr("Error: Mask 2 [%s] specified but not used for %s matrix\n",
                maskexpr.c_str(), MetaData::TypeString(mtype));
      useMask2_ = false;
      return Action::ERR;
    }
    if (useMask2_) {
      mask2_.SetMaskString( maskexpr );
      mkind = DataSet_2D::FULL;
    }
  }

  // Create matrix BYATOM DataSet
  Mat_ = (DataSet_MatrixDbl*)init.DSL().AddSet(DataSet::MATRIX_DBL,
                                         MetaData(name, MetaData::M_MATRIX, mtype), "Mat");
  if (Mat_ == 0) return Action::ERR;
  // NOTE: Type/Kind is set here so subsequent analyses/actions know about it.
  Mat_->SetMatrixKind( mkind );
  // Set default precision for backwards compat.
//.........这里部分代码省略.........