C++ ActionFrame::Frm方法代码示例

Action::RetType Action_AtomicCorr::DoAction(int frameNum, ActionFrame& frm) {
  // On first pass through refframe will be empty and first frame will become ref.
  if (!previousFrame_.empty()) {
    ACvector::iterator atom_vector = atom_vectors_.begin();
    if (acorr_mode_ == ATOM) {
      // For each atom in mask, calc delta position.
      for (AtomMask::const_iterator atom = mask_.begin(); atom != mask_.end(); ++atom) 
      {
        const double* tgtxyz = frm.Frm().XYZ( *atom );
        const double* refxyz = previousFrame_.XYZ( *atom );
        atom_vector->push_back( (float)(tgtxyz[0] - refxyz[0]) );
        atom_vector->push_back( (float)(tgtxyz[1] - refxyz[1]) );
        atom_vector->push_back( (float)(tgtxyz[2] - refxyz[2]) );
        ++atom_vector;
      }
    } else {
      for (std::vector<AtomMask>::const_iterator rmask = resmasks_.begin();
                                                 rmask != resmasks_.end(); ++rmask)
      {
        Vec3 CXYZ = frm.Frm().VGeometricCenter( *rmask );
        Vec3 RXYZ = previousFrame_.VGeometricCenter( *rmask );
        atom_vector->push_back( (float)(CXYZ[0] - RXYZ[0]) );
        atom_vector->push_back( (float)(CXYZ[1] - RXYZ[1]) );
        atom_vector->push_back( (float)(CXYZ[2] - RXYZ[2]) );
        ++atom_vector;
      } 
    }
  }
  // Store this frame as new reference frame
  previousFrame_ = frm.Frm();
  return Action::OK;
}

// Action_Density::HistAction()
Action::RetType Action_Density::HistAction(int frameNum, ActionFrame& frm) {
  long int bin = 0;
  // Loop over masks
  for (unsigned int idx = 0; idx != masks_.size(); ++idx)
  {
    AtomMask const& mask = masks_[idx];
    HistType&       hist = histograms_[idx];
    std::map<long int, double> Sum;
    // Loop over mask atoms
    unsigned int midx = 0;
    for (AtomMask::const_iterator atm = mask.begin(); atm != mask.end(); ++atm, midx++)
    {
      const double* XYZ = frm.Frm().XYZ( *atm );
      if (XYZ[axis_] < 0) {
        // Coordinate is negative. Need to subtract off delta so that proper bin
        // is populated (-delta to 0.0).
        bin = (XYZ[axis_] - delta_) / delta_;
      } else {
        // Coordinate is positive.
        bin = XYZ[axis_] / delta_;
      }
      //mprintf("DEBUG: frm=%6i mask=%3u atm=%8i crd=%8.3f bin=%li\n", frameNum+1, idx, *atm+1, XYZ[axis_], bin);
      Sum[bin] += properties_[idx][midx];
    }
    // Accumulate sums
    hist.accumulate( Sum );
  }

  // Accumulate area
  Box const& box = frm.Frm().BoxCrd();
  area_.accumulate(box[area_coord_[0]] * box[area_coord_[1]]);

  return Action::OK;
}

// Action_Distance::DoAction()
Action::RetType Action_Distance::DoAction(int frameNum, ActionFrame& frm) {
  double Dist;
  Matrix_3x3 ucell, recip;
  Vec3 a1, a2;

  if (useMass_) {
    a1 = frm.Frm().VCenterOfMass( Mask1_ );
    a2 = frm.Frm().VCenterOfMass( Mask2_ );
  } else {
    a1 = frm.Frm().VGeometricCenter( Mask1_ );
    a2 = frm.Frm().VGeometricCenter( Mask2_ );
  }

  switch ( image_.ImageType() ) {
    case NONORTHO:
      frm.Frm().BoxCrd().ToRecip(ucell, recip);
      Dist = DIST2_ImageNonOrtho(a1, a2, ucell, recip);
      break;
    case ORTHO:
      Dist = DIST2_ImageOrtho(a1, a2, frm.Frm().BoxCrd());
      break;
    case NOIMAGE:
      Dist = DIST2_NoImage(a1, a2);
      break;
  }
  Dist = sqrt(Dist);

  dist_->Add(frameNum, &Dist);

  return Action::OK;
}

/** Called every time a frame is read in. Calc distance RMSD.
  * If first is true, set the first frame read in as reference.
  */
Action::RetType Action_DistRmsd::DoAction(int frameNum, ActionFrame& frm) {
  // Perform any needed reference actions
  refHolder_.ActionRef( frm.Frm(), false, false );
  // Set selected frame atoms. Masses have already been set.
  SelectedTgt_.SetCoordinates(frm.Frm(), TgtMask_);
  double DR = SelectedTgt_.DISTRMSD( refHolder_.SelectedRef() );
  drmsd_->Add(frameNum, &DR);
  return Action::OK;
}

// Action_CheckStructure::DoAction()
Action::RetType Action_CheckStructure::DoAction(int frameNum, ActionFrame& frm) {
  int total_problems = CheckOverlap(frameNum+1, frm.Frm(), *CurrentParm_);
  if (bondcheck_)
    total_problems += CheckBonds(frameNum+1, frm.Frm(), *CurrentParm_);

  if (total_problems > 0 && skipBadFrames_)
    return Action::SUPPRESS_COORD_OUTPUT;
  return Action::OK;
}

// Action_Rmsd::DoAction()
Action::RetType Action_Rmsd::DoAction(int frameNum, ActionFrame& frm) {
  // Perform any needed reference actions
  REF_.ActionRef( frm.Frm(), fit_, useMass_ );
  // Calculate RMSD
  double rmsdval;
  Action::RetType err;
  // Set selected frame atoms. Masses have already been set.
  tgtFrame_.SetCoordinates(frm.Frm(), tgtMask_);
  if (!fit_) {
    rmsdval = tgtFrame_.RMSD_NoFit(REF_.SelectedRef(), useMass_);
    err = Action::OK;
  } else {
    rmsdval = tgtFrame_.RMSD_CenteredRef(REF_.SelectedRef(), rot_, tgtTrans_, useMass_);
    if (rmatrices_ != 0) rmatrices_->Add(frameNum, rot_.Dptr());
    if (rotate_)
      frm.ModifyFrm().Trans_Rot_Trans(tgtTrans_, rot_, REF_.RefTrans());
    else {
      tgtTrans_ += REF_.RefTrans();
      frm.ModifyFrm().Translate(tgtTrans_);
    }
    err = Action::MODIFY_COORDS;
  }
  rmsd_->Add(frameNum, &rmsdval);

  // ---=== Per Residue RMSD ===---
  // Set reference and selected frame for each residue using the previously
  // set-up masks in refResMask and tgtResMask. Use SetFrame instead
  // of SetCoordinates since each residue can be a different size.
  if (perres_) {
    for (perResArray::const_iterator PerRes = ResidueRMS_.begin();
                                     PerRes != ResidueRMS_.end(); ++PerRes)
    {
      if ( PerRes->isActive_ ) {
        ResRefFrame_.SetFrame(REF_.RefFrame(), PerRes->refResMask_);
        ResTgtFrame_.SetFrame(frm.Frm(),       PerRes->tgtResMask_);
        if (perrescenter_) {
          ResTgtFrame_.CenterOnOrigin( useMass_ );
          ResRefFrame_.CenterOnOrigin( useMass_ );
        }
        double R = ResTgtFrame_.RMSD_NoFit(ResRefFrame_, useMass_);
        PerRes->data_->Add(frameNum, &R);
      }
    }
  }

  if (REF_.Previous())
    REF_.SetRefStructure( frm.Frm(), fit_, useMass_ );

  return err;
}

// Action_Density::DensityAction()
Action::RetType Action_Density::DensityAction(int frameNum, ActionFrame& frm) {
  Matrix_3x3 ucell, recip;
  double volume = 0.0;
  if (image_.ImageType() == ORTHO)
    volume = frm.Frm().BoxCrd().BoxX() *
             frm.Frm().BoxCrd().BoxY() *
             frm.Frm().BoxCrd().BoxZ();
  else if (image_.ImageType() == NONORTHO)
    volume = frm.Frm().BoxCrd().ToRecip( ucell, recip );
  // Total mass is in delta_
  double density = delta_ / (volume * AMU_ANG_TO_G_CM3);
  density_->Add(frameNum, &density);
  return Action::OK;
}

// Action_LIE::action()
Action::RetType Action_LIE::DoAction(int frameNum, ActionFrame& frm) {
  
  if (doelec_) {
    double e = Calculate_Elec(frm.Frm());
    elec_->Add(frameNum, &e);
  }

  if (dovdw_) {
    double e = Calculate_LJ(frm.Frm(), *CurrentParm_);
    vdw_->Add(frameNum, &e);
  }

  return Action::OK;
}

// Action_CreateReservoir::DoAction()
Action::RetType Action_CreateReservoir::DoAction(int frameNum, ActionFrame& frm) {
  int bin = -1;
  if (bin_ != 0) bin = (int)bin_->Dval(frm.TrajoutNum());
  if (reservoir_.WriteReservoir(nframes_++, frm.Frm(), ene_->Dval(frm.TrajoutNum()), bin))
    return Action::ERR;
  return Action::OK;
}

// Action_Principal::DoAction()
Action::RetType Action_Principal::DoAction(int frameNum, ActionFrame& frm) {
  Matrix_3x3 Inertia;
  Vec3 Eval;

  frm.Frm().CalculateInertia( mask_, Inertia );

  // NOTE: Diagonalize_Sort_Chirality places sorted eigenvectors in rows.
  Inertia.Diagonalize_Sort_Chirality( Eval, debug_ );
  if (outfile_ != 0) {
    int fn = frameNum+1; 
    outfile_->Printf("%i EIGENVALUES: %f %f %f\n%i EIGENVECTOR 0: %f %f %f\n%i EIGENVECTOR 1: %f %f %f\n%i EIGENVECTOR 2: %f %f %f\n", 
      fn, Eval[0], Eval[1], Eval[2],
      fn, Inertia[0], Inertia[1], Inertia[2],
      fn, Inertia[3], Inertia[4], Inertia[5],
      fn, Inertia[6], Inertia[7], Inertia[8]);
    //Eval.Print("PRINCIPAL EIGENVALUES");
    //Inertia.Print("PRINCIPAL EIGENVECTORS (Rows)");
  }
  if (vecData_ != 0) {
    vecData_->AddMat3x3( Inertia );
    valData_->AddVxyz( Eval );
  }
  
  // Rotate - since Evec is already transposed (eigenvectors
  // are returned in rows) just do plain rotation to affect an
  // inverse rotation.
  if (doRotation_) {
    frm.ModifyFrm().Rotate( Inertia );
    return Action::MODIFY_COORDS;
  }
  return Action::OK;
}

// Action_ClusterDihedral::DoAction()
Action::RetType Action_ClusterDihedral::DoAction(int frameNum, ActionFrame& frm) {
  // For each dihedral, calculate which bin it should go into and store bin#
  int bidx = 0;
  for (std::vector<DCmask>::const_iterator dih = DCmasks_.begin(); 
                                           dih != DCmasks_.end(); ++dih)
  {
    double PHI = Torsion( frm.Frm().XYZ(dih->A1()),
                          frm.Frm().XYZ(dih->A2()),
                          frm.Frm().XYZ(dih->A3()),
                          frm.Frm().XYZ(dih->A4()) );
    // NOTE: Torsion is in radians; should bins be converted to rads as well?
    PHI *= Constants::RADDEG;
    //mprintf("[%6i]Dihedral=%8.3f", dih->A1(), PHI); // DEBUG
    PHI -= dih->Min();
    //mprintf(" Shifted=%8.3f", PHI); // DEBUG
    if (PHI < 0) PHI += 360;
    //mprintf(" Wrapped=%8.3f", PHI); // DEBUG
    PHI /= dih->Step();
    int phibin = (int)PHI;
    //mprintf(" Bin=%3i\n", phibin); // DEBUG
    Bins_[bidx++] = phibin;
  }
  // DEBUG - print bins
  //mprintf("[");
  //for (std::vector<int>::const_iterator bin = Bins_.begin(); bin != Bins_.end(); ++bin)
  //  mprintf("%3i,",*bin);
  //mprintf("]\n");
  // Now search for this bin combo in the DCarray
  std::vector<DCnode>::iterator DC = dcarray_.begin();
  for (; DC != dcarray_.end(); ++DC)
  {
    if ( DC->BinMatch( Bins_ ) ) break;
  }
  if (DC == dcarray_.end()) {
    // No match; create new bin combo and store frame num
    //mprintf("NEW DCNODE.\n");
    dcarray_.push_back( DCnode( Bins_, frameNum ) );
  } else {
    // Match; increment bin count and store frame num
    //mprintf("DCNODE ALREADY PRESENT.\n");
    DC->Increment();
    DC->AddFrame( frameNum );
  }
  // Store frame number
  lastframe_ = frameNum;
  return Action::OK;
}

// Action_Density::DoAction()
Action::RetType Action_Density::DoAction(int frameNum, ActionFrame& frm) {
  long slice;
  unsigned long i, j;
  Vec3 coord;
  Box box;

  i = 0;

  for (std::vector<AtomMask>::const_iterator mask = masks_.begin();
       mask != masks_.end();
       mask++)
  {
    j = 0;

    std::map<long,double> minus_histo, plus_histo;

    for (AtomMask::const_iterator idx = mask->begin();
         idx != mask->end();
         idx++)
    {
      coord = frm.Frm().XYZ(*idx);
      slice = (unsigned long) (coord[axis_] / delta_);

      if (coord[axis_] < 0) {
        minus_histo[slice] += properties_[i][j];
      } else {
        plus_histo[slice] += properties_[i][j];
      }

      j++;
    }

    if (minus_histo.size() > 0)
      minus_histograms_[i].accumulate(minus_histo);

    if (plus_histo.size() > 0)
      plus_histograms_[i].accumulate(plus_histo);

    i++;
  }

  box = frm.Frm().BoxCrd();
  area_.accumulate(box[area_coord_[0]] * box[area_coord_[1]]);

  return Action::OK;
}

// Action_Unwrap::DoAction()
Action::RetType Action_Unwrap::DoAction(int frameNum, ActionFrame& frm) {
  Matrix_3x3 ucell, recip;
  // Set reference structure if not already set
  if (RefFrame_.empty()) {
    RefFrame_ = frm.Frm();
    return Action::OK;
  }
 
  if (orthogonal_)
    Image::UnwrapOrtho( frm.ModifyFrm(), RefFrame_, imageList_, center_, true );
  else {
    frm.Frm().BoxCrd().ToRecip( ucell, recip );
    Image::UnwrapNonortho( frm.ModifyFrm(), RefFrame_, imageList_, ucell, recip, center_, true );
  }

  return Action::MODIFY_COORDS;
}

// Action_MultiVector::DoAction()
Action::RetType Action_MultiVector::DoAction(int frameNum, ActionFrame& frm) {
  for (unsigned int nv = 0; nv < CrdIdx1_.size(); ++nv) {
    Vec3 CXYZ( frm.Frm().CRD( CrdIdx1_[nv] ) );
    Vec3 VXYZ( frm.Frm().CRD( CrdIdx2_[nv] ) );
    VXYZ -= CXYZ;
    data_[nv]->AddVxyz(VXYZ, CXYZ);
  }

  return Action::OK;
}

// Action_Angle::action()
Action::RetType Action_Angle::DoAction(int frameNum, ActionFrame& frm) {
  Vec3 a1, a2, a3;
  if (useMass_) {
    a1 = frm.Frm().VCenterOfMass( Mask1_ );
    a2 = frm.Frm().VCenterOfMass( Mask2_ );
    a3 = frm.Frm().VCenterOfMass( Mask3_ );
  } else {
    a1 = frm.Frm().VGeometricCenter( Mask1_ );
    a2 = frm.Frm().VGeometricCenter( Mask2_ );
    a3 = frm.Frm().VGeometricCenter( Mask3_ );
  }
  double aval = CalcAngle( a1.Dptr(), a2.Dptr(), a3.Dptr() );

  aval *= Constants::RADDEG;

  ang_->Add(frameNum, &aval);

  return Action::OK;
}