C++ ActionSetup::Nframes方法代码示例

// Action_Outtraj::Setup()
Action::RetType Action_Outtraj::Setup(ActionSetup& setup) {
  if (associatedParm_->Pindex() != setup.Top().Pindex())
    return Action::SKIP;
  if (!isSetup_) { // TODO: Trajout IsOpen?
    if (outtraj_.SetupTrajWrite(setup.TopAddress(), setup.CoordInfo(), setup.Nframes()))
      return Action::ERR;
    outtraj_.PrintInfo(0);
    isSetup_ = true;
  }
  return Action::OK;
}

/** Determine what atoms each mask pertains to for the current parm file.
  */
Action::RetType Action_ReplicateCell::Setup(ActionSetup& setup) {
  if (setup.Top().SetupIntegerMask( Mask1_ )) return Action::ERR;
  mprintf("\t%s (%i atoms)\n",Mask1_.MaskString(), Mask1_.Nselected());
  if (Mask1_.None()) {
    mprintf("Warning: One or both masks have no atoms.\n");
    return Action::SKIP;
  }
  // Set up imaging info for this parm
  image_.SetupImaging( setup.CoordInfo().TrajBox().Type() );
  if (!image_.ImagingEnabled()) {
    mprintf("Warning: Imaging cannot be performed for topology %s\n", setup.Top().c_str());
    return Action::SKIP;
  }
  // Create combined topology.
  if (combinedTop_.Natom() > 0) {
    // Topology already set up. Check that # atoms matches.
    if (Mask1_.Nselected() * ncopies_ != combinedTop_.Natom()) {
      mprintf("Warning: Unit cell can currently only be replicated for"
              " topologies with same # atoms.\n");
      return Action::SKIP;
    }
    // Otherwise assume top does not change.
  } else {
    // Set up topology and frame.
    Topology* stripParm = setup.Top().modifyStateByMask( Mask1_ );
    if (stripParm == 0) return Action::ERR;
    for (int cell = 0; cell != ncopies_; cell++)
      combinedTop_.AppendTop( *stripParm );
    combinedTop_.Brief("Combined parm:");
    delete stripParm;
    if (!parmfilename_.empty()) {
      ParmFile pfile;
      if (pfile.WriteTopology(combinedTop_, parmfilename_, ParmFile::UNKNOWN_PARM, 0)) {
        mprinterr("Error: Topology file %s not written.\n", parmfilename_.c_str());
        return Action::ERR;
      }
    }
    // Only coordinates for now. FIXME
    combinedFrame_.SetupFrameM(combinedTop_.Atoms());
    // Set up COORDS / output traj if necessary.
    if (coords_ != 0)
      coords_->CoordsSetup( combinedTop_, CoordinateInfo() );
    if (!trajfilename_.empty()) {
      if ( outtraj_.PrepareEnsembleTrajWrite(trajfilename_, trajArgs_,
                                             &combinedTop_, CoordinateInfo(),
                                             setup.Nframes(), TrajectoryFile::UNKNOWN_TRAJ,
                                             ensembleNum_) )
        return Action::ERR;
    }
  }

  return Action::OK;
}

Action::RetType Action_CreateCrd::Setup(ActionSetup& setup) {
  // Set COORDS topology now if not already set.
  if (setup.Top().Pindex() == pindex_ && coords_->Top().Natom() == 0) {
    coords_->CoordsSetup( setup.Top(), setup.CoordInfo() );
    // Estimate memory usage
    mprintf("\tEstimated memory usage (%i frames): %s\n",
            setup.Nframes(),
            ByteString(coords_->SizeInBytes(setup.Nframes()), BYTE_DECIMAL).c_str());
  }
  // If # atoms in currentParm does not match coords, warn user.
  if (setup.Top().Natom() != coords_->Top().Natom()) {
    if (check_) {
      mprinterr("Error: # atoms in current topology (%i) != # atoms in coords set \"%s\" (%i)\n",
                setup.Top().Natom(), coords_->legend(), coords_->Top().Natom());
      return Action::ERR;
    } else {
      mprintf("Warning: # atoms in current topology (%i) != # atoms in coords set \"%s\" (%i)\n"
              "Warning:   The resulting COORDS data set may have problems.\n",
              setup.Top().Natom(), coords_->legend(), coords_->Top().Natom());
    }
  }
  return Action::OK;
}

// Action_Outtraj::Setup()
Action::RetType Action_Outtraj::Setup(ActionSetup& setup) {
  if (!isActive_ || associatedParm_->Pindex() != setup.Top().Pindex()) {
    mprintf("\tOutput trajectory not active for topology '%s'\n", setup.Top().c_str());
    return Action::SKIP;
  }
  if (!isSetup_) { // TODO: Trajout IsOpen?
    if (outtraj_.SetupTrajWrite(setup.TopAddress(), setup.CoordInfo(), setup.Nframes()))
      return Action::ERR;
    mprintf("      "); //TODO this is a kludge; PrintInfo should be a string.
    outtraj_.PrintInfo(0);
    mprintf("\tHas %s\n", outtraj_.Traj().CoordInfo().InfoString().c_str());
    isSetup_ = true;
  }
  return Action::OK;
}

//.........这里部分代码省略.........
                 *ba < setup.Top().Res(*res).LastAtom() )
            {
              if ( !selected[ *ba - res_first_atom ] &&
                   setup.Top()[ *ba ].Element() == Atom::HYDROGEN )
                heavyAtomGroup.push_back( *ba );
            }
          if (!heavyAtomGroup.empty())
            potentialSites.push_back( Site(*res, heavyAtomGroup) );
        }
      }
    }
    mprintf("\t%zu potential NOE sites:\n", potentialSites.size());
    for (SiteArray::const_iterator site = potentialSites.begin();
                                   site != potentialSites.end(); ++site)
    {
      mprintf("  %u\tRes %i:", site - potentialSites.begin(), site->ResNum()+1);
      for (unsigned int idx = 0; idx != site->Nindices(); ++idx)
        mprintf(" %s", setup.Top().TruncAtomNameNum( site->Idx(idx) ).c_str());
      mprintf("\n");
    }
    if (noeArray_.empty()) {
      size_t siteArraySize = 0;
      // Set up all potential NOE pairs. Keep track of size.
      for (SiteArray::const_iterator site1 = potentialSites.begin();
                                     site1 != potentialSites.end(); ++site1)
      {
        for (SiteArray::const_iterator site2 = site1 + 1;
                                       site2 != potentialSites.end(); ++site2)
        {
          if (site1->ResNum() != site2->ResNum()) {
            std::string legend = site1->SiteLegend(setup.Top()) + "--" +
                                 site2->SiteLegend(setup.Top());
            DataSet* ds = 0;
            if (series_) {
              ds = masterDSL_->AddSet(DataSet::FLOAT,
                                      MetaData(setname_, "foundNOE", noeArray_.size()));
              if (ds == 0) return Action::ERR;
              // Construct a data set name.
              ds->SetLegend(legend);
            }
            noeArray_.push_back( NOEtype(*site1, *site2, ds, legend) );
            siteArraySize += (2 * sizeof(int) * site1->Nindices()) +
                             (2 * sizeof(int) * site2->Nindices());
          }
        }
      }
      numNoePairs_ = noeArray_.size();
      size_t siteSize = sizeof(int) + (2 * sizeof(Iarray)) + sizeof(Site);
      size_t noeSize = (2 * siteSize) + sizeof(DataSet*) + sizeof(double) 
                       + sizeof(NOEtype);
      if (series_) noeSize += sizeof(std::vector<float>);
      size_t noeArraySize = (noeSize * numNoePairs_) + siteArraySize;
      if (series_)
        noeArraySize += (setup.Nframes() * numNoePairs_ * sizeof(float));
      mprintf("\t%zu potential NOE pairs. Estimated memory usage is %s\n",
              numNoePairs_, ByteString(noeArraySize, BYTE_DECIMAL).c_str());
    } else if (numNoePairs_ != potentialSites.size()) {
      mprinterr("Warning: Found NOE matrix has already been set up for %zu potential\n"
                "Warning:   NOEs, but %zu NOEs currently found.\n", numNoePairs_,
                potentialSites.size());
      return Action::SKIP;
    }
  }
  // ---------------------------------------------
  // Set up NOEs specified on the command line
  if (!Pairs_.empty()) {
    if (!specifiedNOEs_.empty()) {
      mprintf("Warning: Specifying NOEs currently only works with first topology used.\n");
      return Action::SKIP;
    }
    for (MaskPairArray::iterator mp = Pairs_.begin(); mp != Pairs_.end(); mp++) {
      if (setup.Top().SetupIntegerMask( mp->first )) return Action::ERR;
      int res1 = CheckSameResidue(setup.Top(), mp->first);
      if (res1 < 0) continue;
      if (setup.Top().SetupIntegerMask( mp->second )) return Action::ERR;
      int res2 = CheckSameResidue(setup.Top(), mp->second);
      if (res2 < 0) continue;
      Site site1( res1, mp->first.Selected() );
      Site site2( res2, mp->second.Selected() );
      std::string legend = site1.SiteLegend(setup.Top()) + "--" +
                           site2.SiteLegend(setup.Top());
      DataSet* ds = 0;
      if (series_) {
        ds = masterDSL_->AddSet(DataSet::FLOAT,
                                MetaData(setname_, "specNOE", specifiedNOEs_.size()));
        if (ds == 0) return Action::ERR;
        ds->SetLegend(legend);
      }
      specifiedNOEs_.push_back( NOEtype(site1, site2, ds, legend) );
    }
  } 
  // Set up imaging info for this parm
  Image_.SetupImaging( setup.CoordInfo().TrajBox().Type() );
  if (Image_.ImagingEnabled())
    mprintf("\tImaged.\n");
  else
    mprintf("\tImaging off.\n");

  return Action::OK;  
}

// Action_LESsplit::Setup()
Action::RetType Action_LESsplit::Setup(ActionSetup& setup) {
  if ( !setup.Top().LES().HasLES() ) {
    mprintf("Warning: No LES parameters in '%s', skipping.\n", setup.Top().c_str());
    return Action::SKIP;
  }
  if (lesParm_ == 0) { // First time setup
    // Set up masks for all copies
    lesMasks_.clear();
    lesMasks_.resize( setup.Top().LES().Ncopies() );
    unsigned int atom = 0;
    for (LES_Array::const_iterator les = setup.Top().LES().Array().begin();
                                   les != setup.Top().LES().Array().end(); ++les, ++atom)
    {
      // Copy 0 is in all copies
      if ( les->Copy() == 0 ) {
        for (MaskArray::iterator mask = lesMasks_.begin(); mask != lesMasks_.end(); ++mask)
          mask->AddAtom( atom );
      } else
        lesMasks_[ les->Copy() - 1 ].AddAtom( atom );
    }
    for (unsigned int i = 0; i < lesMasks_.size(); i++) {
      mprintf("\t%i atoms in LES copy %u\n", lesMasks_[i].Nselected(), i+1);
      if ( lesMasks_[i].Nselected() != lesMasks_[0].Nselected() ) {
        mprinterr("Error: Currently all LES copies MUST have same # atoms.\n");
        return Action::ERR;
      }
    }
    // Create topology for first copy
    lesParm_ = setup.Top().modifyStateByMask( lesMasks_[0] );
    if (lesParm_ == 0) return Action::ERR;
    // Set up frames to hold individual copies
    lesFrames_.resize( lesMasks_.size() );
    lesFrames_.SetupFrames(lesParm_->Atoms(), setup.CoordInfo());
    lesPtrs_.resize( lesMasks_.size() );
    for (unsigned int i = 0; i != lesMasks_.size(); i++)
      lesPtrs_[i] = &lesFrames_[i];
    if (lesSplit_) {
      // Set up output ensemble FIXME check overwrites TODO combine init/setup?
      if (lesTraj_.InitEnsembleWrite(trajfilename_, trajArgs_, lesMasks_.size(),
                                     TrajectoryFile::UNKNOWN_TRAJ))
        return Action::ERR;
      if (lesTraj_.SetupEnsembleWrite(lesParm_, setup.CoordInfo(), setup.Nframes()))
         return Action::ERR;
      lesTraj_.PrintInfo(0);
    }
    if (lesAverage_) {
      // For average only care about coords.
      avgFrame_.SetupFrame( lesParm_->Natom() );
      if (avgTraj_.PrepareTrajWrite( avgfilename_, trajArgs_, lesParm_,
                                     CoordinateInfo(), setup.Nframes(),
                                     TrajectoryFile::UNKNOWN_TRAJ ))
        return Action::ERR;
      avgTraj_.PrintInfo(0);
    }
  } else {
    if (lesParm_->Pindex() != setup.Top().Pindex()) {
      mprintf("Warning: Already set up for LES parm '%s'. Skipping '%s'\n",
              lesParm_->c_str(), setup.Top().c_str());
      return Action::SKIP;
    }
  }

  return Action::OK;
}

/** Determine from selected mask atoms which dihedrals will be rotated. */
Action::RetType Action_DihedralScan::Setup(ActionSetup& setup) {
  DihedralScanType dst;
  // If range is empty (i.e. no resrange arg given) look through all 
  // solute residues.
  Range actualRange;
  if (resRange_.Empty())
    actualRange = setup.Top().SoluteResidues();
  else 
    actualRange = resRange_;
  // Search for dihedrals
  if (dihSearch_.FindDihedrals(setup.Top(), actualRange))
    return Action::ERR;
  // For each found dihedral, set up mask of atoms that will move upon 
  // rotation. Also set up mask of atoms in this residue that will not
  // move, including atom2.
  if (debug_>0)
    mprintf("DEBUG: Dihedrals:\n");
  for (DihedralSearch::mask_it dih = dihSearch_.begin();
                               dih != dihSearch_.end(); ++dih)
  {
    dst.checkAtoms.clear();
    // Set mask of atoms that will move during dihedral rotation.
    dst.Rmask = DihedralSearch::MovingAtoms(setup.Top(), dih->A1(), dih->A2());
    // If randomly rotating angles, check for atoms that are in the same
    // residue as A1 but will not move. They need to be checked for clashes
    // since further rotations will not help them.
    if (mode_ == RANDOM && check_for_clashes_) {
      CharMask cMask( dst.Rmask.ConvertToCharMask(), dst.Rmask.Nselected() );
      int a1res = setup.Top()[dih->A1()].ResNum();
      for (int maskatom = setup.Top().Res(a1res).FirstAtom();
               maskatom < setup.Top().Res(a1res).LastAtom(); ++maskatom)
        if (!cMask.AtomInCharMask(maskatom))
          dst.checkAtoms.push_back( maskatom );
      dst.checkAtoms.push_back(dih->A1()); // TODO: Does this need to be added first?
      // Since only the second atom and atoms it is bonded to move during 
      // rotation, base the check on the residue of the second atom.
      dst.resnum = a1res;
    }
    dst.atom0 = dih->A0(); // FIXME: This duplicates info
    dst.atom1 = dih->A1();
    dst.atom2 = dih->A2();
    dst.atom3 = dih->A3();
    BB_dihedrals_.push_back(dst);
    // DEBUG: List dihedral info.
    if (debug_ > 0) {
      mprintf("\t%s-%s-%s-%s\n", 
              setup.Top().TruncResAtomName(dih->A0()).c_str(),
              setup.Top().TruncResAtomName(dih->A1()).c_str(),
              setup.Top().TruncResAtomName(dih->A2()).c_str(),
              setup.Top().TruncResAtomName(dih->A3()).c_str() );
      if (debug_ > 1 && mode_ == RANDOM && check_for_clashes_) {
        mprintf("\t\tCheckAtoms=");
        for (std::vector<int>::const_iterator ca = dst.checkAtoms.begin();
                                              ca != dst.checkAtoms.end(); ++ca)
          mprintf(" %i", *ca + 1);
        mprintf("\n");
      }
      if (debug_ > 2) {
        mprintf("\t\t");
        dst.Rmask.PrintMaskAtoms("Rmask:");
      }
    }
  }

  // Set up CheckStructure for this parm (false = nobondcheck)
  if (checkStructure_.SeparateSetup(setup.Top(),
                                    setup.CoordInfo().TrajBox().Type(), false) != Action::OK)
    return Action::ERR;

  // Set the overall max number of rotations to try
  max_rotations_ = (int) BB_dihedrals_.size();
  max_rotations_ *= max_factor_;

  // Set up simple structure check. First step is coarse; check distances 
  // between a certain atom in each residue (first, COM, CA, some other atom?)
  // to see if residues are in each others neighborhood. Second step is to 
  // check the atoms in each close residue.
  if (check_for_clashes_) {
    ResidueCheckType rct;
    int res = 0;
    for (Topology::res_iterator residue = setup.Top().ResStart();
                                residue != setup.Top().ResEnd(); ++residue)
    {
      rct.resnum = res++;
      rct.start = residue->FirstAtom();
      rct.stop = residue->LastAtom();
      rct.checkatom = rct.start;
      ResCheck_.push_back(rct);
    }
  }

  if (!outfilename_.empty() && CurrentParm_ == 0) // FIXME: Correct frames for # of rotations
    outtraj_.SetupTrajWrite(setup.TopAddress(), setup.CoordInfo(), setup.Nframes());

  CurrentParm_ = setup.TopAddress();
  return Action::OK;  
}