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

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

// Exec_PermuteDihedrals::RandomizeAngles()
void Exec_PermuteDihedrals::RandomizeAngles(Frame& currentFrame, Topology const& topIn) {
  Matrix_3x3 rotationMatrix;
# ifdef DEBUG_PERMUTEDIHEDRALS
  // DEBUG
  int debugframenum=0;
  Trajout_Single DebugTraj;
  DebugTraj.PrepareTrajWrite("debugtraj.nc",ArgList(),(Topology*)&topIn,
                             CoordinateInfo(), BB_dihedrals_.size()*max_factor_,
                             TrajectoryFile::AMBERNETCDF);
  DebugTraj.WriteSingle(debugframenum++,currentFrame);
# endif
  int next_resnum;
  int bestLoop = 0;
  int number_of_rotations = 0;
  // Set max number of rotations to try.
  int max_rotations = (int)BB_dihedrals_.size();
  max_rotations *= max_factor_;

  // Loop over all dihedrals
  std::vector<PermuteDihedralsType>::const_iterator next_dih = BB_dihedrals_.begin();
  next_dih++;
  for (std::vector<PermuteDihedralsType>::const_iterator dih = BB_dihedrals_.begin();
                                                     dih != BB_dihedrals_.end(); 
                                                     ++dih, ++next_dih)
  {
    ++number_of_rotations;
    // Get the residue atom of the next dihedral. Residues up to and
    // including this residue will be checked for bad clashes 
    if (next_dih != BB_dihedrals_.end()) 
      next_resnum = next_dih->resnum;
    else
      next_resnum = dih->resnum - 1;
    // Set axis of rotation
    Vec3 axisOfRotation = currentFrame.SetAxisOfRotation(dih->atom1, dih->atom2);
    // Generate random value to rotate by in radians
    // Guaranteed to rotate by at least 1 degree.
    // NOTE: could potentially rotate 360 - prevent?
    // FIXME: Just use 2PI and rn_gen, get everything in radians
    double theta_in_degrees = ((int)(RN_.rn_gen()*100000) % 360) + 1;
    double theta_in_radians = theta_in_degrees * Constants::DEGRAD;
    // Calculate rotation matrix for random theta
    rotationMatrix.CalcRotationMatrix(axisOfRotation, theta_in_radians);
    int loop_count = 0;
    double clash = 0;
    double bestClash = 0;
    if (debug_>0) mprintf("DEBUG: Rotating dihedral %zu res %8i:\n", dih - BB_dihedrals_.begin(),
                          dih->resnum+1);
    bool rotate_dihedral = true;
    while (rotate_dihedral) {
      if (debug_>0) {
        mprintf("\t%8i %12s %12s, +%.2lf degrees (%i).\n",dih->resnum+1,
                topIn.AtomMaskName(dih->atom1).c_str(),
                topIn.AtomMaskName(dih->atom2).c_str(),
                theta_in_degrees,loop_count);
      }
      // Rotate around axis
      currentFrame.Rotate(rotationMatrix, dih->Rmask);
#     ifdef DEBUG_PERMUTEDIHEDRALS
      // DEBUG
      DebugTraj.WriteSingle(debugframenum++,currentFrame);
#     endif
      // If we dont care about sterics exit here
      if (!check_for_clashes_) break;
      // Check resulting structure for issues
      int checkresidue;
      if (!checkAllResidues_)
        checkresidue = CheckResidue(currentFrame, topIn, *dih, next_resnum, clash);
      else
        checkresidue = CheckResidue(currentFrame, topIn, *dih, topIn.Nres(), clash);
      if (checkresidue==0)
        rotate_dihedral = false;
      else if (checkresidue==-1) {
        if (dih - BB_dihedrals_.begin() < 2) {
          mprinterr("Error: Cannot backtrack; initial structure already has clashes.\n");
          number_of_rotations = max_rotations + 1;
        } else {
          dih--; //  0
          dih--; // -1
          next_dih = dih;
          next_dih++;
          if (debug_>0)
            mprintf("\tCannot resolve clash with further rotations, trying previous again.\n");
        }
        break;
      }
      if (clash > bestClash) {bestClash = clash; bestLoop = loop_count;}
      //n_problems = CheckResidues( currentFrame, second_atom );
      //if (n_problems > -1) {
      //  mprintf("%i\tCheckResidues: %i problems.\n",frameNum,n_problems);
      //  rotate_dihedral = false;
      //} else if (loop_count==0) {
      if (loop_count==0 && rotate_dihedral) {
        if (debug_>0)
          mprintf("\tTrying dihedral increments of +%i\n",increment_);
        // Instead of a new random dihedral, try increments
        theta_in_degrees = (double)increment_;
        theta_in_radians = theta_in_degrees * Constants::DEGRAD;
        // Calculate rotation matrix for new theta
        rotationMatrix.CalcRotationMatrix(axisOfRotation, theta_in_radians);
//.........这里部分代码省略.........

/** \return 1 if a new dihedral should be tried, 0 if no clashes
  * \return -1 if further rotations will not help.
  */
int Exec_PermuteDihedrals::CheckResidue( Frame const& FrameIn, Topology const& topIn,
                                         PermuteDihedralsType const& dih, 
                                         int nextres, double& clash ) 
{
  int resnumIn = dih.resnum;
  int rstart = ResCheck_[ resnumIn ].start;
  int rstop = ResCheck_[ resnumIn ].stop;
  int rcheck = ResCheck_[ resnumIn ].checkatom;
  // Check for clashes with self
# ifdef DEBUG_PERMUTEDIHEDRALS
  mprintf("\tChecking residue %i\n",resnumIn+1);
  mprintf("\tATOMS %i to %i\n",rstart+1,rstop);
# endif
  for (int atom1 = rstart; atom1 < rstop - 1; atom1++) {
    for (int atom2 = atom1 + 1; atom2 < rstop; atom2++) {
      // Skip bonded atoms
      bool isBonded = false;
      for (Atom::bond_iterator bndatm = topIn[atom1].bondbegin();
                               bndatm != topIn[atom1].bondend(); ++bndatm)
        if (*bndatm == atom2) {
          isBonded = true;
          break;
        }
      if (!isBonded) {
        double atomD2 = DIST2_NoImage(FrameIn.XYZ(atom1), FrameIn.XYZ(atom2));
        if (atomD2 < cutoff_) {
#         ifdef DEBUG_PERMUTEDIHEDRALS 
          mprintf("\t\tCurrent Res %i Atoms %s and %s are close (%.3lf)\n", resnumIn+1, 
                  topIn.AtomMaskName(atom1).c_str(),
                  topIn.AtomMaskName(atom2).c_str(), sqrt(atomD2));
#         endif
          clash = atomD2;
          return 1;
        }
      }
    }
  }
  // Check for clashes with previous residues, as well as clashes up to and
  // including the next residue in which a dihedral will be rotated.
  for (int res = 0; res <= nextres; res++) {
    if (res == resnumIn) continue;
    int rstart2 = ResCheck_[ res ].start;
    int rstop2 = ResCheck_[ res ].stop;
    int rcheck2 = ResCheck_[ res ].checkatom;
    double resD2 = DIST2_NoImage(FrameIn.XYZ(rcheck), FrameIn.XYZ(rcheck2));
    // If residues are close enough check each atom
    if (resD2 < rescutoff_) { 
#     ifdef DEBUG_PERMUTEDIHEDRALS
      mprintf("\tRES %i ATOMS %i to %i\n",res+1,rstart2+2,rstop2);
#     endif
      for (int atom1 = rstart; atom1 < rstop; atom1++) {
        for (int atom2 = rstart2; atom2 < rstop2; atom2++) {
          double D2 = DIST2_NoImage(FrameIn.XYZ(atom1), FrameIn.XYZ(atom2));
          if (D2 < cutoff_) {
#           ifdef DEBUG_PERMUTEDIHEDRALS
            mprintf("\t\tResCheck %i Atoms %s and %s are close (%.3lf)\n", res+1,
                    topIn.TruncResAtomName(atom1).c_str(),
                    topIn.TruncResAtomName(atom2).c_str(), sqrt(D2));
#           endif
            clash = D2;
            // If the clash involves any atom that will not be moved by further
            // rotation, indicate it is not possible to resolve clash by
            // more rotation by returning -1.
            //if (atom1 == dih.atom2 || atom1 == dih.atom1) return -1;
            for (std::vector<int>::const_iterator ca = dih.checkAtoms.begin();
                                                  ca != dih.checkAtoms.end(); ca++) 
            {
              if (atom1 == *ca) return -1;
            }
            return 1;
          }
        }
      }
    }
  }
  return 0;
} 

/** Find potential symmetric atoms. All residues up to the last selected
  * residue are considered.
  */
int SymmetricRmsdCalc::SetupSymmRMSD(Topology const& topIn, AtomMask const& tgtMask, bool remapIn)
{
  // Allocate space for remapping selected atoms in target frame. This will
  // also put the correct masses in based on the mask.
  tgtRemap_.SetupFrameFromMask(tgtMask, topIn.Atoms());
  // Create map of original atom numbers to selected indices
  Iarray SelectedIdx( topIn.Natom(), -1 );
  int tgtIdx = 0;
  for (int originalAtom = 0; originalAtom != topIn.Natom(); ++originalAtom)
    if ( originalAtom == tgtMask[tgtIdx] )
      SelectedIdx[originalAtom] = tgtIdx++;
  if (debug_ > 0) {
    mprintf("DEBUG: Original atom -> Selected Index mapping:\n");
    for (int originalAtom = 0; originalAtom != topIn.Natom(); ++originalAtom)
      mprintf("\t%8i -> %8i\n", originalAtom + 1, SelectedIdx[originalAtom] + 1);
  }
  // Create initial 1 to 1 atom map for all selected atoms; indices in 
  // SymmetricAtomIndices will correspond to positions in AMap.
  AMap_.resize( tgtRemap_.Natom() );
  // Determine last selected residue.
  int last_res = topIn[tgtMask.back()].ResNum() + 1;
  mprintf("\tResidues up to %s will be considered for symmetry correction.\n",
          topIn.TruncResNameNum(last_res-1).c_str());
  // In each residue, determine which selected atoms are symmetric.
  SymmetricAtomIndices_.clear();
  AtomMap resmap;
  if (debug_ > 1) resmap.SetDebug(1);
  for (int res = 0; res < last_res; ++res) {
    AtomMap::AtomIndexArray residue_SymmetricGroups;
    if (resmap.SymmetricAtoms(topIn, residue_SymmetricGroups, res)) {
      mprinterr("Error: Finding symmetric atoms in residue '%s'\n",
                topIn.TruncResNameNum(res).c_str());
      return 1;
    }
    if (!residue_SymmetricGroups.empty()) {
      // Which atoms in symmetric groups are selected?
      bool resHasSelectedSymmAtoms = false;
      for (AtomMap::AtomIndexArray::const_iterator symmGroup = residue_SymmetricGroups.begin();
                                                   symmGroup != residue_SymmetricGroups.end();
                                                 ++symmGroup)
      {
        Iarray selectedAtomIndices;
        for (Iarray::const_iterator atnum = symmGroup->begin();
                                    atnum != symmGroup->end(); ++atnum)
        {
          if ( SelectedIdx[*atnum] != -1 )
            selectedAtomIndices.push_back( SelectedIdx[*atnum] ); // Store tgtMask indices
        }
        if (!selectedAtomIndices.empty()) {
          SymmetricAtomIndices_.push_back( selectedAtomIndices );
          resHasSelectedSymmAtoms = true;
        }
      }
      // If remapping and not all atoms in a residue are selected, warn user.
      // TODO: Should they just be considered even if not selected?
      if (remapIn && resHasSelectedSymmAtoms) {
        for (int atom = topIn.Res(res).FirstAtom(); atom != topIn.Res(res).LastAtom(); ++atom)
          if (SelectedIdx[atom] == -1) {
            mprintf("Warning: Not all atoms selected in residue '%s'. Re-mapped\n"
                    "Warning:   structures may appear distorted.\n", 
                    topIn.TruncResNameNum(res).c_str());
            break;
          }
      }
    }
  }
  if (debug_ > 0) {
    mprintf("DEBUG: Potential Symmetric Atom Groups:\n");
    for (AtomIndexArray::const_iterator symmatoms = SymmetricAtomIndices_.begin();
                                        symmatoms != SymmetricAtomIndices_.end();
                                        ++symmatoms)
    {
      mprintf("\t%8u) ", symmatoms - SymmetricAtomIndices_.begin());
      for (Iarray::const_iterator atom = symmatoms->begin();
                                  atom != symmatoms->end(); ++atom)
        mprintf(" %s(%i)", topIn.AtomMaskName(tgtMask[*atom]).c_str(), tgtMask[*atom] + 1);
      mprintf("\n");
    } 
  }
  return 0;
}

// DEBUG
static void DebugContactList(AtomMask const& mask, Topology const& parmIn)
{
  for (AtomMask::const_iterator atom = mask.begin(); atom != mask.end(); ++atom)
    mprintf("\tPotential Contact %li: %s\n", atom - mask.begin(),
            parmIn.AtomMaskName(*atom).c_str());
}

int Action_NativeContacts::DetermineNativeContacts(Topology const& parmIn, Frame const& fIn)
#endif
{
# ifdef MPI
  Frame fIn = frameIn;
  if (trajComm_.Size() > 1) {
    // Ensure all threads use same reference
    if (trajComm_.Master())
      for (int rank = 1; rank < trajComm_.Size(); rank++)
        fIn.SendFrame( rank, trajComm_ );
    else
      fIn.RecvFrame( 0, trajComm_ );
  }
# endif
  if (pfile_ != 0 || nfile_ != 0) {
    refFrame_ = fIn; // Save frame for later PDB output.
    refParm_ = &parmIn;  // Save parm for later PDB output.
  }
  if ( SetupContactLists(parmIn, fIn) ) return 1;
  // If specified, set up contacts maps; base size on atom masks.
  if (nativeMap_ != 0) {
    int matrix_max;
    if (Mask2_.MaskStringSet()) {
      matrix_min_ = std::min( Mask1_[0], Mask2_[0] );
      matrix_max = std::max( Mask1_.back(), Mask2_.back() );
    } else {
      matrix_min_ = Mask1_[0];
      matrix_max = Mask1_.back();
    }
    std::string label("Atom");
    if (byResidue_) {
      matrix_min_ = parmIn[matrix_min_].ResNum();
      matrix_max = parmIn[matrix_max].ResNum();
      label.assign("Residue");
    }
    int matrix_cols = matrix_max - matrix_min_ + 1;
    if (nativeMap_->AllocateHalf( matrix_cols )) return 1;
    if (nonnatMap_->AllocateHalf( matrix_cols )) return 1;
    Dimension matrix_dim( matrix_min_+1, 1, label );
    nativeMap_->SetDim(Dimension::X, matrix_dim);
    nativeMap_->SetDim(Dimension::Y, matrix_dim);
    nonnatMap_->SetDim(Dimension::X, matrix_dim);
    nonnatMap_->SetDim(Dimension::Y, matrix_dim);
  }
  double maxDist2 = 0.0;
  double minDist2 = DBL_MAX;
  nativeContacts_.clear();
  std::pair<contactListType::iterator, bool> ret;
  if (determineNativeContacts_) {
    if ( Mask2_.MaskStringSet() ) {
      for (AtomMask::const_iterator c1 = Mask1_.begin(); c1 != Mask1_.end(); ++c1)
        for (AtomMask::const_iterator c2 = Mask2_.begin(); c2 != Mask2_.end(); ++c2)
        {
          SetNativeContact();
        }
    } else {
      for (AtomMask::const_iterator c1 = Mask1_.begin(); c1 != Mask1_.end(); ++c1)
        for (AtomMask::const_iterator c2 = c1 + 1; c2 != Mask1_.end(); ++c2)
        {
          SetNativeContact();
        }
    }
  }
  //mprintf("\tMinimum observed distance= %f, maximum observed distance= %f\n",
  //        sqrt(minDist2), sqrt(maxDist2));
  // Print contacts
  mprintf("\tSetup %zu native contacts:\n", nativeContacts_.size());
  for (contactListType::const_iterator contact = nativeContacts_.begin();
                                       contact != nativeContacts_.end(); ++contact)
  {
    int a1 = contact->first.first;
    int a2 = contact->first.second;
    mprintf("\t\tAtom '%s' to '%s'\n", parmIn.AtomMaskName(a1).c_str(),
            parmIn.AtomMaskName(a2).c_str());
  }
  return 0;  
}

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