C++ AtomMask类代码示例

/** Print atoms for which the cumulative energy satisfies the given
  * cutoffs. Also create MOL2 files containing those atoms.
  */
int Action_Pairwise::PrintCutAtoms(Frame const& frame, int frameNum, EoutType ctype,
                                   Darray const& Earray, double cutIn)
{
  AtomMask CutMask;  // Hold atoms that satisfy the cutoff
  Darray CutCharges; // Hold evdw/eelec corresponding to CutMask atoms.

  if (Eout_ != 0) {
    if (nb_calcType_==COMPARE_REF)
      Eout_->Printf("\tPAIRWISE: Cumulative d%s:", CalcString[ctype]);
    else
      Eout_->Printf("\tPAIRWISE: Cumulative %s:", CalcString[ctype]);
    Eout_->Printf(" %s < %.4f, %s > %.4f\n", CalcString[ctype], -cutIn,
                 CalcString[ctype], cutIn);
  }
  for (AtomMask::const_iterator atom = Mask0_.begin(); atom != Mask0_.end(); ++atom)
  {
    if (fabs(Earray[*atom]) > cutIn)
    {
      if (Eout_ != 0) 
        Eout_->Printf("\t\t%[email protected]%s: %12.4f\n", *atom+1,
                    (*CurrentParm_)[*atom].c_str(), Earray[*atom]);
      CutMask.AddAtom(*atom);
      CutCharges.push_back(Earray[*atom]);
    }
  }
  // Write mol2 with atoms satisfying cutoff
  if (!mol2Prefix_.empty() && !CutMask.None()) {
    if (WriteCutFrame(frameNum, *CurrentParm_, CutMask, CutCharges, 
                      frame, mol2Prefix_ + CutName[ctype])) 
      return 1;
  }

  return 0;
}

/** Determine VDW long range correction prefactor. */
void Ewald::Setup_VDW_Correction(Topology const& topIn, AtomMask const& maskIn) {
  Vdw_Recip_term_ = 0.0;
  NB_ = static_cast<NonbondParmType const*>( &(topIn.Nonbond()) );
  if (!NB_->HasNonbond()) {
    mprintf("Warning: '%s' has no nonbonded parameters. Cannot calculate VDW correction.\n",
            topIn.c_str());
    return;
  }
  // Count the number of each unique nonbonded type.
  Iarray N_vdw_type( NB_->Ntypes(), 0 );
  for (AtomMask::const_iterator atm = maskIn.begin(); atm != maskIn.end(); ++atm)
    N_vdw_type[ topIn[*atm].TypeIndex() ]++;
  if (debug_ > 0) {
    mprintf("DEBUG: %zu VDW types.\n", N_vdw_type.size());
    for (Iarray::const_iterator it = N_vdw_type.begin(); it != N_vdw_type.end(); ++it)
      mprintf("\tType %li = %i\n", it-N_vdw_type.begin(), *it);
  }
  // Determine correction term from types and LJ B parameters
  for (unsigned int itype = 0; itype != N_vdw_type.size(); itype++)
  {
    unsigned int offset = N_vdw_type.size() * itype;
    for (unsigned int jtype = 0; jtype != N_vdw_type.size(); jtype++)
    {
      unsigned int idx = offset + jtype;
      int nbidx = NB_->NBindex()[ idx ];
      if (nbidx > -1)
        Vdw_Recip_term_ += N_vdw_type[itype] * N_vdw_type[jtype] * NB_->NBarray()[ nbidx ].B();
    }
  }
}

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

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

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

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

// Action_Center::Init()
Action::RetType Action_Center::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  // Get keywords
  useMass_ = actionArgs.hasKey("mass");
  ReferenceFrame refFrm = init.DSL().GetReferenceFrame( actionArgs );
  if (refFrm.error()) return Action::ERR;
  // Determine center mode.
  if (!refFrm.empty())
    centerMode_ = REF;
  else if (actionArgs.hasKey("origin"))
    centerMode_ = ORIGIN;
  else if (actionArgs.hasKey("point")) {
    centerMode_ = POINT;
    refCenter_[0] = actionArgs.getNextDouble(0.0);
    refCenter_[1] = actionArgs.getNextDouble(0.0);
    refCenter_[2] = actionArgs.getNextDouble(0.0);
  } else
    centerMode_ = BOXCTR;

  // Get Masks
  Mask_.SetMaskString( actionArgs.GetMaskNext() );
  // Get reference mask if reference specified.
  AtomMask refMask;
  if (centerMode_ == REF) {
    std::string rMaskExpr = actionArgs.GetMaskNext();
    if (rMaskExpr.empty())
      rMaskExpr = Mask_.MaskExpression();
    refMask.SetMaskString( rMaskExpr );
    if (refFrm.Parm().SetupIntegerMask( refMask, refFrm.Coord() ))
      return Action::ERR;
    // Get center of mask in reference
    if (useMass_)
      refCenter_ = refFrm.Coord().VCenterOfMass( refMask );
    else
      refCenter_ = refFrm.Coord().VGeometricCenter( refMask );
  }

  mprintf("    CENTER: Centering coordinates using");
  if (useMass_)
    mprintf(" center of mass");
  else
    mprintf(" geometric center");
  mprintf(" of atoms in mask (%s) to\n", Mask_.MaskString());
  switch (centerMode_) {
    case ORIGIN: mprintf("\tcoordinate origin.\n"); break;
    case BOXCTR: mprintf("\tbox center.\n"); break;
    case REF:
      mprintf("\tcenter of mask (%s) in reference '%s'.\n", refMask.MaskString(),
              refFrm.refName());
      break;
    case POINT:  mprintf("\tpoint (%g, %g, %g).\n",
                         refCenter_[0], refCenter_[1], refCenter_[2]);
      break;
  }
  return Action::OK;
}

/** Remove any selected solvent atoms from mask. */
static void removeSelectedSolvent( Topology const& parmIn, AtomMask& mask ) {
  AtomMask newMask = mask;
  newMask.ClearSelected();
  for (AtomMask::const_iterator atom = mask.begin(); atom != mask.end(); ++atom) {
    int molnum = parmIn[*atom].MolNum();
    if (!parmIn.Mol(molnum).IsSolvent())
      newMask.AddSelectedAtom( *atom );
  }
  mask = newMask;
}

/** Set up atom/residue indices corresponding to atoms selected in mask.
  * This is done to make creating an atom/residue contact map easier.
  */
Action_NativeContacts::Iarray Action_NativeContacts::SetupContactIndices(
                                AtomMask const& mask, Topology const& parmIn)
{
  Iarray contactIdx;
  for (AtomMask::const_iterator atom = mask.begin(); atom != mask.end(); ++atom)
    if (byResidue_)
      contactIdx.push_back( parmIn[*atom].ResNum() );
    else
      contactIdx.push_back( *atom );
  return contactIdx;
}

/** Check that all atoms in mask belong to same residue. */
int Action_NMRrst::CheckSameResidue(Topology const& top, AtomMask const& mask) const {
  if (mask.None()) return -1;
  int resnum = top[mask[0]].ResNum();
  for (AtomMask::const_iterator at = mask.begin(); at != mask.end(); ++at) {
    int r = top[*at].ResNum();
    if (r != resnum) {
      mprintf("Warning: Mask atom %i %s not in same residue as %i %s\n",
              *at + 1, top.AtomMaskName(*at).c_str(),
              mask[0] + 1, top.AtomMaskName(mask[0]).c_str());
    }
  }
  return resnum;
}

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);
}

/** Calculate full nonbonded energy with PME */
double Ewald_ParticleMesh::CalcEnergy(Frame const& frameIn, AtomMask const& maskIn, double& e_vdw)
{
  t_total_.Start();
  Matrix_3x3 ucell, recip;
  double volume = frameIn.BoxCrd().ToRecip(ucell, recip);
  double e_self = Self( volume );
  double e_vdw_lr_correction;

  pairList_.CreatePairList(frameIn, ucell, recip, maskIn);

  // TODO make more efficient
  int idx = 0;
  coordsD_.clear();
  for (AtomMask::const_iterator atm = maskIn.begin(); atm != maskIn.end(); ++atm, ++idx) {
    const double* XYZ = frameIn.XYZ( *atm );
    coordsD_.push_back( XYZ[0] );
    coordsD_.push_back( XYZ[1] );
    coordsD_.push_back( XYZ[2] );
  }

//  MapCoords(frameIn, ucell, recip, maskIn);
  double e_recip = Recip_ParticleMesh( frameIn.BoxCrd() );

  // TODO branch
  double e_vdw6self, e_vdw6recip;
  if (lw_coeff_ > 0.0) {
    e_vdw6self = Self6();
    e_vdw6recip = LJ_Recip_ParticleMesh( frameIn.BoxCrd() );
    if (debug_ > 0) {
      mprintf("DEBUG: e_vdw6self = %16.8f\n", e_vdw6self);
      mprintf("DEBUG: Evdwrecip = %16.8f\n", e_vdw6recip);
    }
    e_vdw_lr_correction = 0.0;
  } else {
    e_vdw6self = 0.0;
    e_vdw6recip = 0.0;
    e_vdw_lr_correction = Vdw_Correction( volume );
  }

  e_vdw = 0.0;
  double e_direct = Direct( pairList_, e_vdw );
  if (debug_ > 0)
    mprintf("DEBUG: Eself= %20.10f   Erecip= %20.10f   Edirect= %20.10f  Evdw= %20.10f\n",
            e_self, e_recip, e_direct, e_vdw);
  e_vdw += (e_vdw_lr_correction + e_vdw6self + e_vdw6recip);
  t_total_.Stop();
  return e_self + e_recip + e_direct;
}

/** Convert charges to Amber units. Calculate sum of charges and squared charges. */
void Ewald::CalculateCharges(Topology const& topIn, AtomMask const& maskIn) {
  sumq_ = 0.0;
  sumq2_ = 0.0;
  Charge_.clear();
  TypeIndices_.clear();
  for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom) {
    double qi = topIn[*atom].Charge() * Constants::ELECTOAMBER;
    Charge_.push_back(qi);
    sumq_ += qi;
    sumq2_ += (qi * qi);
    // Store atom type indices for selected atoms.
    TypeIndices_.push_back( topIn[*atom].TypeIndex() );
  }
  //mprintf("DEBUG: sumq= %20.10f   sumq2= %20.10f\n", sumq_, sumq2_);
  Setup_VDW_Correction( topIn, maskIn );
}

// -----------------------------------------------------------------------------
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_Center::Init()
Action::RetType Action_Center::Init(ArgList& actionArgs, TopologyList* PFL, DataSetList* DSL, DataFileList* DFL, int debugIn)
{
  // Get keywords
  if (actionArgs.hasKey("origin"))
    centerMode_ = ORIGIN;
  else
    centerMode_ = BOXCTR;
  useMass_ = actionArgs.hasKey("mass");
  ReferenceFrame refFrm = DSL->GetReferenceFrame( actionArgs );
  if (refFrm.error()) return Action::ERR;

  // Get Masks
  Mask_.SetMaskString( actionArgs.GetMaskNext() );
  // Get reference mask if reference specified.
  AtomMask refMask;
  if (!refFrm.empty()) {
    std::string rMaskExpr = actionArgs.GetMaskNext();
    if (rMaskExpr.empty())
      rMaskExpr = Mask_.MaskExpression();
    refMask.SetMaskString( rMaskExpr );
    if (refFrm.Parm().SetupIntegerMask( refMask, refFrm.Coord() ))
      return Action::ERR;
    // Get center of mask in reference
    if (useMass_)
      refCenter_ = refFrm.Coord().VCenterOfMass( refMask );
    else
      refCenter_ = refFrm.Coord().VGeometricCenter( refMask );
    centerMode_ = POINT; 
  }

  mprintf("    CENTER: Centering coordinates using");
  if (useMass_)
    mprintf(" center of mass");
  else
    mprintf(" geometric center");
  mprintf(" of atoms in mask (%s) to\n", Mask_.MaskString());
  if (centerMode_ == POINT)
    mprintf("\tcenter of mask (%s) in reference '%s'.\n", refMask.MaskString(),
            refFrm.refName());
  else if (centerMode_ == ORIGIN)
    mprintf("\tcoordinate origin.\n");
  else
    mprintf("\tbox center.\n");

  return Action::OK;
}