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

// -----------------------------------------------------------------------------
// Exec_PermuteDihedrals::IntervalAngles()
void Exec_PermuteDihedrals::IntervalAngles(Frame const& frameIn, Topology const& topIn,
                                           double interval_in_deg)
{
  Matrix_3x3 rotationMatrix;
  double theta_in_radians = interval_in_deg * Constants::DEGRAD;
  int maxVal = (int) (360.0 / interval_in_deg);
  if (maxVal < 0) maxVal = -maxVal;
  // Write original frame
  if (outtraj_.IsInitialized())
    outtraj_.WriteSingle(outframe_++, frameIn);
  if (crdout_ != 0)
    crdout_->AddFrame( frameIn );
  Frame currentFrame = frameIn;
  for (std::vector<PermuteDihedralsType>::const_iterator dih = BB_dihedrals_.begin();
                                                     dih != BB_dihedrals_.end();
                                                     ++dih)
  {
    // Set axis of rotation
    Vec3 axisOfRotation = currentFrame.SetAxisOfRotation(dih->atom1, dih->atom2);
    // Calculate rotation matrix for interval 
    rotationMatrix.CalcRotationMatrix(axisOfRotation, theta_in_radians);
    if (debug_ > 0) {
      mprintf("\tRotating Dih %s-%s by %.2f deg %i times.\n",
               topIn.TruncResAtomName( dih->atom1 ).c_str(), 
               topIn.TruncResAtomName( dih->atom2 ).c_str(), interval_in_deg, maxVal); 
    }
    for (int rot = 0; rot != maxVal; ++rot) {
      // Rotate around axis
      currentFrame.Rotate(rotationMatrix, dih->Rmask);
      // Write output trajectory
      if (outtraj_.IsInitialized())
        outtraj_.WriteSingle(outframe_++, currentFrame);
      if (crdout_ != 0)
        crdout_->AddFrame( currentFrame );
    }
  }
}

/** Check for bad bond lengths. */
int Action_CheckStructure::CheckBonds(int frameNum, Frame const& currentFrame, Topology const& top)
{
  double D2;
  int idx;
  int Nproblems = 0;

  int bond_max = (int)bondList_.size();
# ifdef _OPENMP
# pragma omp parallel private(idx,D2) reduction(+: Nproblems)
  {
  //mprintf("OPENMP: %i threads\n",omp_get_num_threads());
  //mythread = omp_get_thread_num();
# pragma omp for
# endif
  for (idx = 0; idx < bond_max; idx++) {
    D2 = DIST2_NoImage( currentFrame.XYZ(bondList_[idx].a1_),
                        currentFrame.XYZ(bondList_[idx].a2_) );
    if (D2 > bondList_[idx].Req_off2_) {
      ++Nproblems;
      if (outfile_ != 0) {
#       ifdef _OPENMP
#       pragma omp critical
#       endif
        outfile_->Printf(
                    "%i\t Warning: Unusual bond length %i:%s to %i:%s (%.2lf)\n", frameNum,
                    bondList_[idx].a1_+1, top.TruncResAtomName(bondList_[idx].a1_).c_str(),
                    bondList_[idx].a2_+1, top.TruncResAtomName(bondList_[idx].a2_).c_str(),
                    sqrt(D2));
      }
    }
  } // END loop over bonds
# ifdef _OPENMP
  } // END pragma omp parallel
# endif
  return Nproblems;   
}

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

/** Check for bad overlaps. */
int Action_CheckStructure::CheckOverlap(int frameNum, Frame const& currentFrame, Topology const& top)
{
  double D2;
  Matrix_3x3 ucell, recip; // ToFrac, ToCart
  int nmask1, nmask2;
  int atom1, atom2;
  int Nproblems = 0;

  // Get imaging info for non-orthogonal box // TODO Check volume
  if (image_.ImageType()==NONORTHO)
    currentFrame.BoxCrd().ToRecip(ucell, recip);
  if ( Mask2_.MaskStringSet() ) {
    // Calculation of all atoms in Mask1 to all atoms in Mask2
    int outer_max = OuterMask_.Nselected();
    int inner_max = InnerMask_.Nselected();
#   ifdef _OPENMP
#   pragma omp parallel private(nmask1,nmask2,atom1,atom2,D2) reduction(+: Nproblems)
    {
    //mprintf("OPENMP: %i threads\n",omp_get_num_threads());
    //mythread = omp_get_thread_num();
#   pragma omp for
#   endif
    for (nmask1 = 0; nmask1 < outer_max; nmask1++) {
      atom1 = OuterMask_[nmask1];
      for (nmask2 = 0; nmask2 < inner_max; nmask2++) {
        atom2 = InnerMask_[nmask2];
        if (atom1 != atom2) {
          D2 = DIST2( currentFrame.XYZ(atom1), currentFrame.XYZ(atom2),
                      image_.ImageType(), currentFrame.BoxCrd(), ucell, recip);
          if (D2 < nonbondcut2_) {
            ++Nproblems;
            if (outfile_ != 0) {
#             ifdef _OPENMP
#             pragma omp critical
#             endif
              outfile_->Printf(
                    "%i\t Warning: Atoms %i:%s and %i:%s are close (%.2lf)\n", frameNum,
                    atom1+1, top.TruncResAtomName(atom1).c_str(),
                    atom2+1, top.TruncResAtomName(atom2).c_str(), sqrt(D2));
            }
          }
        }
      } // END loop over inner mask
    } // END loop over outer mask
#   ifdef _OPENMP
    } // END pragma omp parallel
#   endif
  } else {
    // Calculation of atoms in Mask1 to all other atoms in Mask1
    int mask1_max = Mask1_.Nselected();
#   ifdef _OPENMP
#   pragma omp parallel private(nmask1,nmask2,atom1,atom2,D2) reduction(+: Nproblems)
    {
    //mprintf("OPENMP: %i threads\n",omp_get_num_threads());
    //mythread = omp_get_thread_num();
#   pragma omp for schedule(dynamic)
#   endif
    for (nmask1 = 0; nmask1 < mask1_max; nmask1++) {
      atom1 = Mask1_[nmask1];
      for (nmask2 = nmask1 + 1; nmask2 < mask1_max; nmask2++) {
        atom2 = Mask1_[nmask2];
        D2 = DIST2( currentFrame.XYZ(atom1), currentFrame.XYZ(atom2),
                    image_.ImageType(), currentFrame.BoxCrd(), ucell, recip);
        if (D2 < nonbondcut2_) {
          ++Nproblems;
          if (outfile_ != 0) {
#           ifdef _OPENMP
#           pragma omp critical
#           endif
            outfile_->Printf(
                  "%i\t Warning: Atoms %i:%s and %i:%s are close (%.2lf)\n", frameNum,
                  atom1+1, top.TruncResAtomName(atom1).c_str(),
                  atom2+1, top.TruncResAtomName(atom2).c_str(), sqrt(D2));
          }
        }
      } // END inner loop over Mask1
    } // END outer loop over Mask1
#   ifdef _OPENMP
    } // END pragma omp parallel
#   endif
  }

  return Nproblems;
}