C++ AtomMask::Nselected方法代码示例

/** Set up the exclusion list based on the given mask and parm.
  * \return the total number of interactions, -1 on error.
  */
int Action_Pairwise::SetupNonbondParm(AtomMask const& maskIn, Topology const& ParmIn) {
  // Check if LJ parameters present - need at least 2 atoms for it to matter.
  if (ParmIn.Natom() > 1 && !ParmIn.Nonbond().HasNonbond()) {
    mprinterr("Error: Topology does not have LJ information.\n");
    return -1;
  }

  // Determine the actual number of pairwise interactions that will be calcd.
  // This is ((N^2 - N) / 2) - SUM[ #excluded atoms]
  int N_interactions = ((maskIn.Nselected() * maskIn.Nselected()) - maskIn.Nselected()) / 2;
  for (AtomMask::const_iterator at = maskIn.begin(); at != maskIn.end(); ++at)
    N_interactions -= ParmIn[ *at ].Nexcluded();

  // DEBUG - Print total number of interactions for this parm
  mprintf("\t%i interactions for this parm.\n",N_interactions);

  // DEBUG - Print exclusion list for each atom
  /*for (unsigned int atom = 0; atom < exclusionList.size(); atom++) {
    mprintf("\t%8u:",atom + 1);
    for (std::vector<int>::iterator eat = exclusionList[atom].begin();
                                    eat != exclusionList[atom].end();
                                    eat++)
    {
      mprintf(" %i",*eat + 1);
    }
    mprintf("\n");
  }*/
  return N_interactions;
}

/** Write file containing only cut atoms and energies as charges. */
int Action_Pairwise::WriteCutFrame(int frameNum, Topology const& Parm, AtomMask const& CutMask, 
                                   Darray const& CutCharges,
                                   Frame const& frame, std::string const& outfilename) 
{
  if (CutMask.Nselected() != (int)CutCharges.size()) {
    mprinterr("Error: WriteCutFrame: # of charges (%u) != # mask atoms (%i)\n",
              CutCharges.size(), CutMask.Nselected());
    return 1;
  }
  Frame CutFrame(frame, CutMask);
  Topology* CutParm = Parm.modifyStateByMask( CutMask );
  if (CutParm == 0) return 1;
  // Set new charges
  for (int i = 0; i != CutParm->Natom(); i++)
    CutParm->SetAtom(i).SetCharge( CutCharges[i] );
  int err = 0;
  Trajout_Single tout;
  if (tout.PrepareTrajWrite(outfilename, "multi", CutParm, CoordinateInfo(), 1,
                            TrajectoryFile::MOL2FILE))
  {
    mprinterr("Error: Could not set up cut mol2 file %s\n", outfilename.c_str());
    err = 1;
  } else {
    tout.WriteSingle(frameNum, CutFrame);
    tout.EndTraj();
  }
  delete CutParm;
  return err;
}

// Action_Matrix::MaskToMatResArray()
Action_Matrix::MatResArray Action_Matrix::MaskToMatResArray(Topology const& currentParm,
                                                            AtomMask const& mask) const
{
  MatResArray residues;
  int currentResNum = -1;
  matrix_res blank_res;
  for (int idx = 0; idx != mask.Nselected(); idx++)
  {
    int atom1 = mask[idx];
    int resNum = currentParm[atom1].ResNum();
    if (resNum != currentResNum) {
      residues.push_back( blank_res );
      residues.back().resnum_ = resNum;
      currentResNum = resNum;
    }
    residues.back().maskIdxs_.push_back( idx );
  }
  if (debug_ > 0) {
    mprintf("DEBUG: BYRES: MASK '%s'\n", mask.MaskString());
    for (MatResArray::const_iterator res = residues.begin(); res != residues.end(); ++res) {
      mprintf("\tRes %i:", res->resnum_+1);
      for (Iarray::const_iterator it = res->maskIdxs_.begin(); it != res->maskIdxs_.end(); ++it)
        mprintf(" %i (%i)", mask[*it] + 1, *it);
      mprintf("\n");
    }
  }
  return residues;
}

// -----------------------------------------------------------------------------
void Image::WrapToCell0(std::vector<double>& CoordsIn, Frame const& frmIn,
                        AtomMask const& maskIn,
                        Matrix_3x3 const& ucell, Matrix_3x3 const& recip)
{
  double* uFrac = &CoordsIn[0];
  int nUatoms = maskIn.Nselected();
  int idx;
  double* result;
  const double* XYZ;
# ifdef _OPENMP
# pragma omp parallel private(idx, result, XYZ)
  {
# pragma omp for
# endif
  for (idx = 0; idx < nUatoms; idx++)
  {
    result = uFrac + idx*3;
    XYZ = frmIn.XYZ( maskIn[idx] );
    // Convert to frac coords
    recip.TimesVec( result, XYZ );
    // Wrap to primary unit cell
    result[0] = result[0] - floor(result[0]);
    result[1] = result[1] - floor(result[1]);
    result[2] = result[2] - floor(result[2]);
    // Convert back to Cartesian
    ucell.TransposeMult( result, result );
  }
# ifdef _OPENMP
  } // END pragma omp parallel
# endif
}

// Action_Matrix::FillMassArray()
Action_Matrix::Darray Action_Matrix::FillMassArray(Topology const& currentParm, 
                                                   AtomMask const& mask) const
{
  Darray mass;
  mass.reserve( mask.Nselected() );
  for (AtomMask::const_iterator atom = mask.begin(); atom != mask.end(); ++atom) 
    mass.push_back( currentParm[ *atom ].Mass() );
  return mass;
}

/** Setup PME calculation. */
int Ewald_ParticleMesh::Setup(Topology const& topIn, AtomMask const& maskIn) {
  CalculateCharges(topIn, maskIn);
  // NOTE: These dont need to actually be calculated if the lj ewald coeff
  //       is 0.0, but do it here anyway to avoid segfaults.
  CalculateC6params( topIn, maskIn );
  coordsD_.clear();
  coordsD_.reserve( maskIn.Nselected() * 3);
  SetupExcluded(topIn, maskIn);
  return 0;
}

/** Set up exclusion lists for selected atoms. */
void Ewald::SetupExcluded(Topology const& topIn, AtomMask const& maskIn)
{
  Excluded_.clear();
  Excluded_.resize( maskIn.Nselected() );
  // Create a character mask so we can see if atoms in excluded lists are
  // also selected.
  CharMask Cmask(maskIn.ConvertToCharMask(), maskIn.Nselected());
  // Create a map of atom number to maskIn index.
  int selectedIdx = 0;
  Iarray atToIdx( Cmask.Natom(), -1 );
  for (int cidx = 0; cidx != Cmask.Natom(); cidx++)
    if (Cmask.AtomInCharMask(cidx))
      atToIdx[cidx] = selectedIdx++;
  // Loop over selected atoms
  for (int idx = 0; idx != maskIn.Nselected(); idx++)
  {
    // Always exclude self
    Excluded_[idx].insert( idx );
    int at = maskIn[idx];
    for (Atom::excluded_iterator excluded_atom = topIn[at].excludedbegin();
                                 excluded_atom != topIn[at].excludedend();
                               ++excluded_atom)
    {
      if (Cmask.AtomInCharMask(*excluded_atom))
      {
        // Find excluded atoms index in maskIn
        int excluded_idx = atToIdx[*excluded_atom];
        Excluded_[idx         ].insert( excluded_idx );
        Excluded_[excluded_idx].insert( idx          );
      }
    }
  }
  unsigned int ex_size = 0;
  for (Iarray2D::const_iterator it = Excluded_.begin(); it != Excluded_.end(); ++it)
    ex_size += it->size();
  mprintf("\tMemory used by full exclusion list: %s\n",
          ByteString(ex_size * sizeof(int), BYTE_DECIMAL).c_str());
}

void Ewald::CalculateC6params(Topology const& topIn, AtomMask const& maskIn) {
  Cparam_.clear();
  if (lw_coeff_ > 0.0) {
    for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom)
    {
      double rmin = topIn.GetVDWradius( *atom );
      double eps  = topIn.GetVDWdepth( *atom );
      Cparam_.push_back( 8.0 * (rmin*rmin*rmin) * sqrt(2 * eps) );
      if (debug_ > 0)
        mprintf("DEBUG: C6 param atom %8i = %16.8f\n", *atom+1, Cparam_.back());
    }
  } else
    Cparam_.assign(maskIn.Nselected(), 0.0);
}

/** Place selected atoms into grid cells. Convert to fractional coords, wrap
  * into primary cell, then determine grid cell.
  */
void PairList::GridUnitCell(Frame const& frmIn, Matrix_3x3 const& ucell,
                             Matrix_3x3 const& recip, AtomMask const& maskIn)
{
  // Clear any existing atoms in cells.
  for (Carray::iterator cell = cells_.begin(); cell != cells_.end(); ++cell)
    cell->ClearAtoms();
  Frac_.clear();
  Frac_.reserve( maskIn.Nselected() );
  if (frmIn.BoxCrd().Type() == Box::ORTHO) {
    // Orthogonal imaging
    for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom)
    {
      const double* XYZ = frmIn.XYZ(*atom);
      Vec3 fc( XYZ[0]*recip[0],    XYZ[1]*recip[4],    XYZ[2]*recip[8]   );
      Vec3 fcw(fc[0]-floor(fc[0]), fc[1]-floor(fc[1]), fc[2]-floor(fc[2]));
      Vec3 ccw(fcw[0]*ucell[0],    fcw[1]*ucell[4],    fcw[2]*ucell[8]   );
#     ifdef DEBUG_PAIRLIST
      mprintf("DBG: o %6i fc=%7.3f%7.3f%7.3f  fcw=%7.3f%7.3f%7.3f  ccw=%7.3f%7.3f%7.3f\n",
              *atom+1, fc[0], fc[1], fc[2], fcw[0], fcw[1], fcw[2], ccw[0], ccw[1], ccw[2]);
#     endif
      GridAtom( atom-maskIn.begin(), fcw, ccw );
    }
  } else {
    // Non-orthogonal imaging
    for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom)
    {
      Vec3 fc = recip * Vec3(frmIn.XYZ(*atom));
      Vec3 fcw(fc[0]-floor(fc[0]), fc[1]-floor(fc[1]), fc[2]-floor(fc[2]));
      Vec3 ccw = ucell.TransposeMult( fcw );
#     ifdef DEBUG_PAIRLIST
      mprintf("DBG: n %6i fc=%7.3f%7.3f%7.3f  fcw=%7.3f%7.3f%7.3f  ccw=%7.3f%7.3f%7.3f\n",
              *atom+1, fc[0], fc[1], fc[2], fcw[0], fcw[1], fcw[2], ccw[0], ccw[1], ccw[2]);
#     endif
      GridAtom( atom-maskIn.begin(), fcw, ccw );
    }
  }
}

static inline Exec::RetType MaskError(AtomMask const& mask) {
  mprinterr("Error: Mask '%s' selects %i atoms, expected 1.\n",
            mask.MaskString(), mask.Nselected());
  return CpptrajState::ERR;
}

/** Set up bond parameters for bonds for which both atoms present in mask. */
void Action_CheckStructure::SetupBondList(AtomMask const& iMask, Topology const& top) {
  CharMask cMask( iMask.ConvertToCharMask(), iMask.Nselected() );
 
  ProcessBondArray(top.Bonds(),  top.BondParm(), cMask);
  ProcessBondArray(top.BondsH(), top.BondParm(), cMask);
}