C++ Matrix_3x3::CalcRotationMatrix方法代码示例

// Action_DihedralScan::IntervalAngles()
void Action_DihedralScan::IntervalAngles(Frame& currentFrame) {
  Matrix_3x3 rotationMatrix;
  double theta_in_radians = interval_ * Constants::DEGRAD;
  // Write original frame
  if (!outfilename_.empty())
    outtraj_.WriteSingle(outframe_++, currentFrame);
  for (std::vector<DihedralScanType>::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",
               CurrentParm_->TruncResAtomName( (*dih).atom1 ).c_str(), 
               CurrentParm_->TruncResAtomName( (*dih).atom2 ).c_str(), interval_, maxVal_); 
    }
    for (int rot = 0; rot < maxVal_; ++rot) {
      // Rotate around axis
      currentFrame.Rotate(rotationMatrix, (*dih).Rmask);
      // Write output trajectory
      if (!outfilename_.empty())
        outtraj_.WriteSingle(outframe_++, currentFrame);
    }
  }
}

// Action_MakeStructure::DoAction()
Action::RetType Action_MakeStructure::DoAction(int frameNum, Frame* currentFrame, 
                                               Frame** frameAddress) 
{
  Matrix_3x3 rotationMatrix;
  for (std::vector<SecStructHolder>::iterator ss = secstruct_.begin();
                                              ss != secstruct_.end(); ++ss)
  {
    std::vector<float>::iterator theta = ss->thetas_.begin();
    std::vector<AtomMask>::iterator Rmask = ss->Rmasks_.begin();
    for (DihedralSearch::mask_it dih = ss->dihSearch_.begin();
                                 dih != ss->dihSearch_.end(); ++dih, ++theta, ++Rmask)
    {
      double theta_in_radians = (double)*theta;
      // Calculate current value of dihedral
      double torsion = Torsion( currentFrame->XYZ( (*dih).A0() ),
                                currentFrame->XYZ( (*dih).A1() ),
                                currentFrame->XYZ( (*dih).A2() ),
                                currentFrame->XYZ( (*dih).A3() ) );
      // Calculate delta needed to get to theta
      double delta = theta_in_radians - torsion;
      // Set axis of rotation
      Vec3 axisOfRotation = currentFrame->SetAxisOfRotation((*dih).A1(), (*dih).A2());
      // Calculate rotation matrix for delta 
      rotationMatrix.CalcRotationMatrix(axisOfRotation, delta);
      if (debug_ > 0) {
        std::string a0name = CurrentParm_->TruncResAtomName( (*dih).A0() );
        std::string a1name = CurrentParm_->TruncResAtomName( (*dih).A1() );
        std::string a2name = CurrentParm_->TruncResAtomName( (*dih).A2() );
        std::string a3name = CurrentParm_->TruncResAtomName( (*dih).A3() );
          mprintf("\tRotating Dih %i:%s (%i-%i-%i-%i) (@%.2f) by %.2f deg to get to %.2f.\n",
                  (*dih).ResNum()+1, (*dih).Name().c_str(),
                  (*dih).A0() + 1, (*dih).A1() + 1, (*dih).A2() + 1, (*dih).A3() + 1, 
                  torsion*Constants::RADDEG, delta*Constants::RADDEG, theta_in_radians*Constants::RADDEG);
      }
      // Rotate around axis
      currentFrame->Rotate(rotationMatrix, *Rmask);
    }
  }
  return Action::OK;
}

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

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

//.........这里部分代码省略.........
  // Determine whether we are setting or incrementing.
  enum ModeType { SET = 0, INCREMENT };
  ModeType mode = SET;
  if (argIn.Contains("value"))
    mode = SET;
  else if (argIn.Contains("increment"))
    mode = INCREMENT;
  else {
    mprinterr("Error: Specify 'value <value>' or 'increment <increment>'\n");
    return CpptrajState::ERR;
  }
  double value = argIn.getKeyDouble(ModeStr[mode], 0.0);
  switch (mode) {
    case SET: mprintf("\tDihedral will be set to %g degrees.\n", value); break;
    case INCREMENT: mprintf("\tDihedral will be incremented by %g degrees.\n", value); break;
  }
  // Convert to radians
  value *= Constants::DEGRAD;

  // Select dihedral atoms
  int A1, A2, A3, A4;
  if (argIn.Contains("type")) {
    // By type
    ArgList typeArg = argIn.GetStringKey("type");
    if (typeArg.empty()) {
      mprinterr("Error: No dihedral type specified after 'type'\n");
      return CpptrajState::ERR;
    }
    DihedralSearch dihSearch;
    dihSearch.SearchForArgs( typeArg );
    if (dihSearch.NoDihedralTokens()) {
      mprinterr("Error: Specified dihedral type not recognized.\n");
      return CpptrajState::ERR;
    }
    // Get residue
    int res = argIn.getKeyInt("res", -1);
    if (res <= 0) {
      mprinterr("Error: If 'type' specified 'res' must be specified and > 0.\n");
      return CpptrajState::ERR;
    }
    // Search for dihedrals. User residue #s start from 1.
    if (dihSearch.FindDihedrals(CRD->Top(), Range(res-1)))
      return CpptrajState::ERR;
    DihedralSearch::mask_it dih = dihSearch.begin();
    A1 = dih->A0();
    A2 = dih->A1();
    A3 = dih->A2();
    A4 = dih->A3();
  } else {
    // By masks
    AtomMask m1( argIn.GetMaskNext() );
    AtomMask m2( argIn.GetMaskNext() );
    AtomMask m3( argIn.GetMaskNext() );
    AtomMask m4( argIn.GetMaskNext() );
    if (CRD->Top().SetupIntegerMask( m1 )) return CpptrajState::ERR;
    if (CRD->Top().SetupIntegerMask( m2 )) return CpptrajState::ERR;
    if (CRD->Top().SetupIntegerMask( m3 )) return CpptrajState::ERR;
    if (CRD->Top().SetupIntegerMask( m4 )) return CpptrajState::ERR;
    if (m1.Nselected() != 1) return MaskError( m1 );
    if (m2.Nselected() != 1) return MaskError( m2 );
    if (m3.Nselected() != 1) return MaskError( m3 );
    if (m4.Nselected() != 1) return MaskError( m4 );
    A1 = m1[0];
    A2 = m2[0];
    A3 = m3[0];
    A4 = m4[0];
  }
  mprintf("\tRotating dihedral defined by atoms '%s'-'%s'-'%s'-'%s'\n",
          CRD->Top().AtomMaskName(A1).c_str(),
          CRD->Top().AtomMaskName(A2).c_str(),
          CRD->Top().AtomMaskName(A3).c_str(),
          CRD->Top().AtomMaskName(A4).c_str());
  // Set mask of atoms that will move during dihedral rotation
  AtomMask Rmask = DihedralSearch::MovingAtoms(CRD->Top(), A2, A3);
  // Calculate current value of dihedral
  double torsion = Torsion( FRM.XYZ(A1), FRM.XYZ(A2), FRM.XYZ(A3), FRM.XYZ(A4) );
  // Calculate delta needed to get to target value.
  double delta;
  switch (mode) {
    case SET:       delta = value - torsion; break;
    case INCREMENT: delta = value; break;
  }
  mprintf("\tOriginal torsion is %g, rotating by %g degrees.\n",
          torsion*Constants::RADDEG, delta*Constants::RADDEG);
  // Set axis of rotation
  Vec3 axisOfRotation = FRM.SetAxisOfRotation( A2, A3 );
  // Calculate rotation matrix for delta.
  Matrix_3x3 rotationMatrix;
  rotationMatrix.CalcRotationMatrix(axisOfRotation, delta);
  // Rotate around axis
  FRM.Rotate(rotationMatrix, Rmask);
  // RMS-fit the non-moving part of the coords back on original
  AtomMask refMask = Rmask;
  refMask.InvertMask();
  FRM.Align( refFrame, refMask );
  // Update coords
  OUT->SetCRD( outframe, FRM );

  return CpptrajState::OK;
}