C++ Frame::BoxCrd方法代码示例

/** \param frameIn Frame to image.
  * \param bp Box + boundary.
  * \param bm Box - boundary.
  * \param center If true image w.r.t. center of atoms, otherwise first atom.
  * \param useMass If true calc center of mass, otherwise geometric center.
  */
void Image::Ortho(Frame& frameIn, Vec3 const& bp, Vec3 const& bm, Vec3 const& offIn,
                  bool center, bool useMass, PairType const& AtomPairs)
{
  Vec3 Coord;
  Vec3 offset(offIn[0] * frameIn.BoxCrd()[0],
              offIn[1] * frameIn.BoxCrd()[1],
              offIn[2] * frameIn.BoxCrd()[2]);
  // Loop over atom pairs
  for (PairType::const_iterator atom = AtomPairs.begin();
                                atom != AtomPairs.end(); atom++)
  {
    int firstAtom = *atom;
    ++atom;
    int lastAtom = *atom;
    //if (debug>2)
    //  mprintf( "  IMAGE processing atoms %i to %i\n", firstAtom+1, lastAtom);
    // Set up Coord with position to check for imaging based on first atom or 
    // center of mass of atoms first to last.
    if (center) {
      if (useMass)
        Coord = frameIn.VCenterOfMass(firstAtom,lastAtom);
      else
        Coord = frameIn.VGeometricCenter(firstAtom,lastAtom);
    } else 
      Coord = frameIn.XYZ( firstAtom );

    // boxTrans will hold calculated translation needed to move atoms back into box
    Vec3 boxTrans = Ortho(Coord, bp, bm, frameIn.BoxCrd()) + offset;

    // Translate atoms according to Coord
    frameIn.Translate(boxTrans, firstAtom, lastAtom);
  } // END loop over atom pairs
}

void Action_Closest::Action_ImageOrtho(Frame& frmIn, double maxD)
{
  double Dist;
  int solventMol;
  AtomMask::const_iterator solute_atom;
  Iarray::const_iterator solvent_atom;
    // Loop over all solvent molecules in original frame
  if (useMaskCenter_) {
    Vec3 maskCenter = frmIn.VGeometricCenter( distanceMask_ );
    for (solventMol=0; solventMol < NsolventMolecules_; solventMol++) {
      SolventMols_[solventMol].D = maxD;
      for (solvent_atom = SolventMols_[solventMol].solventAtoms.begin();
           solvent_atom != SolventMols_[solventMol].solventAtoms.end(); ++solvent_atom)
      {

        Dist = DIST2_ImageOrtho( maskCenter.Dptr(),
                      frmIn.XYZ(*solvent_atom),frmIn.BoxCrd());
        if (Dist < SolventMols_[solventMol].D) 
          SolventMols_[solventMol].D = Dist;
      }
    }
  } else {
    for (solventMol=0; solventMol < NsolventMolecules_; solventMol++) {
      if (debug_ > 1)
        mprintf("DEBUG: Calculating distance for molecule %i\n", solventMol);
      // Set the initial minimum distance for this solvent mol to be the
      // max possible distance.
      SolventMols_[solventMol].D = maxD;
      // Calculate distance between each atom in distanceMask and atoms in solvent Mask
      for (solvent_atom = SolventMols_[solventMol].solventAtoms.begin();
           solvent_atom != SolventMols_[solventMol].solventAtoms.end(); ++solvent_atom)
      {
        for (solute_atom = distanceMask_.begin(); 
             solute_atom != distanceMask_.end(); ++solute_atom)
        {
          Dist = DIST2_ImageOrtho(frmIn.XYZ(*solute_atom),
                       frmIn.XYZ(*solvent_atom), frmIn.BoxCrd());
          if (Dist < SolventMols_[solventMol].D) 
            SolventMols_[solventMol].D = Dist;
          if (debug_ > 2)
            mprintf("DEBUG: SolvMol %i, soluteAtom %i, solventAtom %i, D= %f, minD= %f\n",
                    solventMol, *solute_atom, *solvent_atom, Dist,
                    sqrt(SolventMols_[solventMol].D));
        }
      }
      if (debug_ > 1) mprintf("DEBUG:\tMol %8i minD= %lf\n",solventMol, SolventMols_[solventMol].D);
    } // END for loop over solventMol
  }

}

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

/** Process an NOEtypeArray */
void Action_NMRrst::ProcessNoeArray(NOEtypeArray& Narray, Frame const& frameIn, int frameNum) const
{
  for (NOEtypeArray::iterator my_noe = Narray.begin();
                              my_noe != Narray.end(); ++my_noe)
  {
    unsigned int shortest_idx1 = 0, shortest_idx2 = 0;
    double shortest_dist2 = -1.0;
    for (unsigned int idx1 = 0; idx1 != my_noe->Site1().Nindices(); ++idx1)
    {
      for (unsigned int idx2 = 0; idx2 != my_noe->Site2().Nindices(); ++idx2)
      {
        double dist2 = DIST2(frameIn.XYZ(my_noe->Site1().Idx(idx1)),
                             frameIn.XYZ(my_noe->Site2().Idx(idx2)),
                             Image_.ImageType(), frameIn.BoxCrd(),
                             ucell_, recip_);
        if (shortest_dist2 < 0.0 || dist2 < shortest_dist2) {
          shortest_dist2 = dist2;
          shortest_idx1 = idx1;
          shortest_idx2 = idx2;
        }
      }
    }
    // NOTE: Saving d^2
    my_noe->UpdateNOE(frameNum, shortest_dist2, shortest_idx1, shortest_idx2);
  }
}

double Action_LIE::Calculate_LJ(Frame const& frameIn, Topology const& parmIn) const {
  double result = 0;
  // Loop over ligand atoms
  AtomMask::const_iterator mask1_end = Mask1_.end();
  AtomMask::const_iterator mask2_end = Mask2_.end();
  for (AtomMask::const_iterator maskatom1 = Mask1_.begin();
       maskatom1 != mask1_end; maskatom1++) {

    int crdidx1 = (*maskatom1) * 3; // index into coordinate array
    Vec3 atm1 = Vec3(frameIn.CRD(crdidx1));

    for (AtomMask::const_iterator maskatom2 = Mask2_.begin();
         maskatom2 != mask2_end; maskatom2++) {

      int crdidx2 = (*maskatom2) * 3; // index into coordinate array
      Vec3 atm2 = Vec3(frameIn.CRD(crdidx2));
      
      double dist2;
      // Get imaged distance
      Matrix_3x3 ucell, recip;
      switch( ImageType() ) {
        case NONORTHO:
          frameIn.BoxCrd().ToRecip(ucell, recip);
          dist2 = DIST2_ImageNonOrtho(atm1, atm2, ucell, recip);
          break;
        case ORTHO:
          dist2 = DIST2_ImageOrtho(atm1, atm2, frameIn.BoxCrd());
          break;
        default:
          dist2 = DIST2_NoImage(atm1, atm2);
      }

      if (dist2 > cut2vdw_) continue;
      // Here we add to our nonbonded (VDW) energy
      NonbondType const& LJ = parmIn.GetLJparam(*maskatom1, *maskatom2);
      double r2 = 1 / dist2;
      double r6 = r2 * r2 * r2;
      result += LJ.A() * r6 * r6 - LJ.B() * r6;
    }
  }

  return result;
}

double Action_LIE::Calculate_Elec(Frame const& frameIn) const {
  double result = 0;
  // Loop over ligand atoms
  AtomMask::const_iterator mask1_end = Mask1_.end();
  AtomMask::const_iterator mask2_end = Mask2_.end();
  for (AtomMask::const_iterator maskatom1 = Mask1_.begin();
       maskatom1 != mask1_end; maskatom1++) {

    int crdidx1 = (*maskatom1) * 3; // index into coordinate array
    Vec3 atm1 = Vec3(frameIn.CRD(crdidx1));

    for (AtomMask::const_iterator maskatom2 = Mask2_.begin();
         maskatom2 != mask2_end; maskatom2++) {

      int crdidx2 = (*maskatom2) * 3; // index into coordinate array
      Vec3 atm2 = Vec3(frameIn.CRD(crdidx2));
      
      double dist2;
      // Get imaged distance
      Matrix_3x3 ucell, recip;
      switch( ImageType() ) {
        case NONORTHO:
          frameIn.BoxCrd().ToRecip(ucell, recip);
          dist2 = DIST2_ImageNonOrtho(atm1, atm2, ucell, recip);
          break;
        case ORTHO:
          dist2 = DIST2_ImageOrtho(atm1, atm2, frameIn.BoxCrd());
          break;
        default:
          dist2 = DIST2_NoImage(atm1, atm2);
      }

      if (dist2 > cut2elec_) continue;
      // Here we add to our electrostatic energy
      double qiqj = atom_charge_[*maskatom1] * atom_charge_[*maskatom2];
      double shift = (1 - dist2 * onecut2_);
      result += qiqj / sqrt(dist2) * shift * shift;
    }
  }

  return result;
}

/** \param frameIn Frame to image.
  * \param origin If true image w.r.t. coordinate origin.
  * \param fcom If truncoct is true, calc distance w.r.t. this coordinate.
  * \param ucell Unit cell matrix.
  * \param recip Reciprocal coordinates matrix.
  * \param truncoct If true imaging will occur using truncated octahedron shape.
  * \param center If true image w.r.t. center coords, otherwise use first atom coords.
  * \param useMass If true use COM, otherwise geometric center.
  * \param AtomPairs Atom pairs to image.
  */
void Image::Nonortho(Frame& frameIn, bool origin, Vec3 const& fcom, Vec3 const& offIn, 
                     Matrix_3x3 const& ucell, Matrix_3x3 const& recip,
                     bool truncoct, bool center,
                     bool useMass, PairType const& AtomPairs)
{
  Vec3 Coord;
  Vec3 offset = ucell.TransposeMult( offIn );
  double min = -1.0;

  if (truncoct)
    min = 100.0 * (frameIn.BoxCrd().BoxX()*frameIn.BoxCrd().BoxX()+
                   frameIn.BoxCrd().BoxY()*frameIn.BoxCrd().BoxY()+
                   frameIn.BoxCrd().BoxZ()*frameIn.BoxCrd().BoxZ());

  // Loop over atom pairs
  for (PairType::const_iterator atom = AtomPairs.begin();
                                atom != AtomPairs.end(); ++atom)
  {
    int firstAtom = *atom;
    ++atom;
    int lastAtom = *atom;
    //if (debug>2)
    //  mprintf( "  IMAGE processing atoms %i to %i\n", firstAtom+1, lastAtom);
    // Set up Coord with position to check for imaging based on first atom or 
    // center of mass of atoms first to last.
    if (center) {
      if (useMass)
        Coord = frameIn.VCenterOfMass(firstAtom,lastAtom);
      else
        Coord = frameIn.VGeometricCenter(firstAtom,lastAtom);
    } else 
      Coord = frameIn.XYZ( firstAtom );

    // boxTrans will hold calculated translation needed to move atoms back into box
    Vec3 boxTrans = Nonortho(Coord, truncoct, origin, ucell, recip, fcom, min) + offset;

    frameIn.Translate(boxTrans, firstAtom, lastAtom);

  } // END loop over atom pairs
}

/** Set up centering if putting nonortho cell into familiar trunc. oct. shape.
  * \param frameIn Frame to set up for.
  * \param ComMask If not null center is calcd w.r.t. center of atoms in mask.
  * \param useMass If true calculate COM, otherwise calc geometric center.
  * \param origin If true and ComMask is null use origin, otherwise use box center.
  * \return Coordinates of center.
  */
Vec3 Image::SetupTruncoct( Frame const& frameIn, AtomMask* ComMask, bool useMass, bool origin)
{
  if (ComMask!=0) {
    // Use center of atoms in mask
    if (useMass)
      return frameIn.VCenterOfMass( *ComMask );
    else
      return frameIn.VGeometricCenter( *ComMask );
  } else if (!origin) {
    // Use box center
    return frameIn.BoxCrd().Center(); 
  }
  //fprintf(stdout,"DEBUG: fcom = %lf %lf %lf\n",fcom[0],fcom[1],fcom[2]);
  return Vec3(0.0, 0.0, 0.0); // Default is origin {0,0,0}
}

// PairList::CreatePairList()
int PairList::CreatePairList(Frame const& frmIn, Matrix_3x3 const& ucell,
                             Matrix_3x3 const& recip, AtomMask const& maskIn)
{
  t_total_.Start();
  // Calculate translation vectors based on current unit cell.
  FillTranslateVec(ucell);
  // If box size has changed a lot this will reallocate grid
  t_gridpointers_.Start();
  if (SetupGrids(frmIn.BoxCrd().RecipLengths(recip))) return 1;
  t_gridpointers_.Stop();
  // Place atoms in grid cells
  t_map_.Start();
  GridUnitCell(frmIn, ucell, recip, maskIn);
  t_map_.Stop();
  t_total_.Stop();
  return 0;
}

/** Calculate Ewald energy.
  * Slow version that does not use pair list. Note that pair list is still
  * called since we require the fractional and imaged coords.
  * FIXME No Eadjust calc.
  */
double Ewald::CalcEnergy_NoPairList(Frame const& frameIn, Topology const& topIn,
                                    AtomMask const& maskIn)
{
  t_total_.Start();
  Matrix_3x3 ucell, recip;
  double volume = frameIn.BoxCrd().ToRecip(ucell, recip);
  double e_self = Self( volume );
  // Place atoms in pairlist. This calcs frac/imaged coords.
  pairList_.CreatePairList(frameIn, ucell, recip, maskIn);

  double e_recip = Recip_Regular( recip, volume );

  double e_direct = Direct( ucell, topIn, maskIn );

  //mprintf("DEBUG: Eself= %20.10f   Erecip= %20.10f   Edirect= %20.10f\n",
  //        e_self, e_recip, e_direct);
  t_total_.Stop();
  return e_self + e_recip + e_direct;
}

/** Place selected atoms into grid cells. Convert to fractional coords, wrap
  * into primary cell, then determine grid cell.
  */
void PairList::GridUnitCell(Frame const& frmIn, Matrix_3x3 const& ucell,
                             Matrix_3x3 const& recip, AtomMask const& maskIn)
{
  // Clear any existing atoms in cells.
  for (Carray::iterator cell = cells_.begin(); cell != cells_.end(); ++cell)
    cell->ClearAtoms();
  Frac_.clear();
  Frac_.reserve( maskIn.Nselected() );
  if (frmIn.BoxCrd().Type() == Box::ORTHO) {
    // Orthogonal imaging
    for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom)
    {
      const double* XYZ = frmIn.XYZ(*atom);
      Vec3 fc( XYZ[0]*recip[0],    XYZ[1]*recip[4],    XYZ[2]*recip[8]   );
      Vec3 fcw(fc[0]-floor(fc[0]), fc[1]-floor(fc[1]), fc[2]-floor(fc[2]));
      Vec3 ccw(fcw[0]*ucell[0],    fcw[1]*ucell[4],    fcw[2]*ucell[8]   );
#     ifdef DEBUG_PAIRLIST
      mprintf("DBG: o %6i fc=%7.3f%7.3f%7.3f  fcw=%7.3f%7.3f%7.3f  ccw=%7.3f%7.3f%7.3f\n",
              *atom+1, fc[0], fc[1], fc[2], fcw[0], fcw[1], fcw[2], ccw[0], ccw[1], ccw[2]);
#     endif
      GridAtom( atom-maskIn.begin(), fcw, ccw );
    }
  } else {
    // Non-orthogonal imaging
    for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom)
    {
      Vec3 fc = recip * Vec3(frmIn.XYZ(*atom));
      Vec3 fcw(fc[0]-floor(fc[0]), fc[1]-floor(fc[1]), fc[2]-floor(fc[2]));
      Vec3 ccw = ucell.TransposeMult( fcw );
#     ifdef DEBUG_PAIRLIST
      mprintf("DBG: n %6i fc=%7.3f%7.3f%7.3f  fcw=%7.3f%7.3f%7.3f  ccw=%7.3f%7.3f%7.3f\n",
              *atom+1, fc[0], fc[1], fc[2], fcw[0], fcw[1], fcw[2], ccw[0], ccw[1], ccw[2]);
#     endif
      GridAtom( atom-maskIn.begin(), fcw, ccw );
    }
  }
}

// Image::UnwrapOrtho()
void Image::UnwrapOrtho( Frame& tgtIn, Frame& refIn, PairType const& AtomPairs,
                         bool center, bool useMass )
{
  Vec3 vtgt, vref, boxTrans;
  Vec3 boxVec = tgtIn.BoxCrd().Lengths();
  // Loop over atom pairs
  for (PairType::const_iterator atom = AtomPairs.begin();
                                atom != AtomPairs.end(); ++atom)
  {
    int firstAtom = *atom;
    ++atom;
    int lastAtom = *atom;
    if (center) {
      // Use center of coordinates between first and last atoms.
      if (useMass) {
        vtgt = tgtIn.VCenterOfMass(firstAtom, lastAtom);
        vref = refIn.VCenterOfMass(firstAtom, lastAtom);
      } else {
        vtgt = tgtIn.VGeometricCenter(firstAtom, lastAtom);
        vref = refIn.VGeometricCenter(firstAtom, lastAtom);
      }
    } else {
      // Use position first atom only.
      vtgt = tgtIn.XYZ( firstAtom );
      vref = refIn.XYZ( firstAtom );
    }
    Vec3 dxyz = vtgt - vref;
    boxTrans[0] = -floor( dxyz[0] / boxVec[0] + 0.5 ) * boxVec[0];
    boxTrans[1] = -floor( dxyz[1] / boxVec[1] + 0.5 ) * boxVec[1];
    boxTrans[2] = -floor( dxyz[2] / boxVec[2] + 0.5 ) * boxVec[2];
    // Translate atoms from first to last
    tgtIn.Translate(boxTrans, firstAtom, lastAtom);
    // Save new ref positions
    int i3 = firstAtom * 3;
    std::copy( tgtIn.xAddress()+i3, tgtIn.xAddress()+(lastAtom*3), refIn.xAddress()+i3 );
  } // END loop over atom pairs
}

/** Find the minimum distance between atoms in distanceMask and each 
  * solvent Mask.
  */
Action_Closest::RetType Action_Closest::DoAction(int frameNum, Frame& frmIn) { 
  double maxD;
  Matrix_3x3 ucell, recip;
  AtomMask::const_iterator solute_atom;
  Iarray::const_iterator solvent_atom;

  if (image_.ImagingEnabled()) {
    frmIn.BoxCrd().ToRecip(ucell, recip);
    // Calculate max possible imaged distance
    maxD = frmIn.BoxCrd().BoxX() + frmIn.BoxCrd().BoxY() + 
           frmIn.BoxCrd().BoxZ();
    maxD *= maxD;
  } else {
    // If not imaging, set max distance to an arbitrarily large number
    maxD = DBL_MAX;
  }

  //subroutines to find the distance
  // if (image_.ImageType() == NOIMAGE)
  //   Action_NoImage(frmIn,maxD);
  // else if (image_.ImageType() == ORTHO)
  //   Action_ImageOrtho(frmIn,maxD);
  // else
  //   Action_ImageNonOrtho(frmIn,maxD, ucell,recip);

//remove this ..TODOi

useMaskCenter_ = false;
cudaEvent_t start_event, stop_event;
float elapsed_time_seq;

bool v[2] = { true, false };
int type = 2;   //keep it no imaging (as of now)
char* dict[3] = {"NONE", "ORTHO", "NON-ORTHO"};
for (int type  = 0; type < 3 ; type++){
  printf("-------------------------------------------------------\n");
for(int k =0 ; k < 2 ; k++)
{
  printf("Solute Center : %s\n", v[k] ? "YES" : "NO");
  printf("Imaging : %s\n", dict[type]); 
  useMaskCenter_ = v[k];

cudaEventCreate(&start_event);
cudaEventCreate(&stop_event);
cudaEventRecord(start_event, 0);

//serial section of the code
if (type == 0 )
  Action_NoImage(frmIn,maxD);
else if (type == 1)
  Action_ImageOrtho(frmIn,maxD);
else if  (type == 2)
  Action_ImageNonOrtho(frmIn, maxD, ucell, recip);
else{
  printf("Error, invalid imaging type\n");
  exit(1);
}



cudaThreadSynchronize();
cudaEventRecord(stop_event, 0);
cudaEventSynchronize(stop_event);
cudaEventElapsedTime(&elapsed_time_seq,start_event, stop_event );
printf("Done with kernel SEQ Kernel Time: %.2f\n", elapsed_time_seq);



bool result = true;
float elapsed_time_gpu;
if (useMaskCenter_)
  result = cuda_action_center(frmIn,maxD,ucell ,recip ,type ,elapsed_time_gpu); 
else
  result = cuda_action_no_center(frmIn,maxD,ucell ,recip ,type ,elapsed_time_gpu);//handling all the data formatting and copying etc
// we will only care about kernel time
//fixing the overhead will be later

  if(result){
    printf("CUDA PASS\n");
  }
  else{
    printf("CUDA FAIL!\n");
    exit(0);
  }

  printf("->#Seq Time:  = %0.2f\n", elapsed_time_seq);
  printf("->#CUDA Time: = %0.2f\n", elapsed_time_gpu);
  printf("->#Speedup =  %0.2f\n", elapsed_time_seq/elapsed_time_gpu);
}
}

  // Sort distances
  std::sort( SolventMols_.begin(), SolventMols_.end(), moldist_cmp() );
  // Add first closestWaters solvent atoms to stripMask
  std::vector<MolDist>::iterator solventend = SolventMols_.begin() + closestWaters_;
  for ( std::vector<MolDist>::const_iterator solvent = SolventMols_.begin();
                                             solvent != solventend;
//.........这里部分代码省略.........

int Traj_GmxTrX::writeFrame(int set, Frame const& frameOut) {
    int tsize;
    // Write header
    file_.Write( &Magic_, 4 );
    tsize = (int)Title().size() + 1;
    file_.Write( &tsize, 4);
    --tsize;
    file_.Write( &tsize, 4);
    file_.Write( Title().c_str(), Title().size() );
    file_.Write( &ir_size_, 4 );
    file_.Write( &e_size_, 4 );
    file_.Write( &box_size_, 4 );
    file_.Write( &vir_size_, 4 );
    file_.Write( &pres_size_, 4 );
    file_.Write( &top_size_, 4 );
    file_.Write( &sym_size_, 4 );
    file_.Write( &x_size_, 4 );
    file_.Write( &v_size_, 4 );
    file_.Write( &f_size_, 4 );
    file_.Write( &natoms_, 4 );
    file_.Write( &step_, 4 );
    file_.Write( &nre_, 4 );
    dt_ = (float)set;
    file_.Write( &dt_, 4 ); // TODO: Write actual time
    file_.Write( &lambda_, 4 );
    // Write box
    // NOTE: GROMACS units are nm
    if (box_size_ > 0) {
        double ucell[9];
        double by = frameOut.BoxCrd().BoxY() * 0.1;
        double bz = frameOut.BoxCrd().BoxZ() * 0.1;
        ucell[0] = frameOut.BoxCrd().BoxX() * 0.1;
        ucell[1] = 0.0;
        ucell[2] = 0.0;
        ucell[3] = by*cos(Constants::DEGRAD*frameOut.BoxCrd().Gamma());
        ucell[4] = by*sin(Constants::DEGRAD*frameOut.BoxCrd().Gamma());
        ucell[5] = 0.0;
        ucell[6] = bz*cos(Constants::DEGRAD*frameOut.BoxCrd().Beta());
        ucell[7] = (by*bz*cos(Constants::DEGRAD*frameOut.BoxCrd().Alpha()) - ucell[6]*ucell[3]) /
                   ucell[4];
        ucell[8] = sqrt(bz*bz - ucell[6]*ucell[6] - ucell[7]*ucell[7]);
        if (precision_ == sizeof(float)) {
            float f_ucell[9];
            for (int i = 0; i < 9; i++)
                f_ucell[i] = (float)ucell[i];
            file_.Write( f_ucell, box_size_ );
        } else // double
            file_.Write( ucell, box_size_ );
    }
    // Write coords/velo
    // NOTE: GROMACS units are nm
    const double* Xptr = frameOut.xAddress();
    const double* Vptr = frameOut.vAddress();
    int ix = 0;
    if (precision_ == sizeof(float)) {
        for (; ix < natom3_; ix++)
            farray_[ix] = (float)(Xptr[ix] * 0.1);
        if (v_size_ > 0)
            for (int iv = 0; iv < natom3_; iv++, ix++)
                farray_[ix] = (float)(Vptr[iv] * 0.1);
        file_.Write( farray_, x_size_ + v_size_ );
    } else { // double
        for (; ix < natom3_; ix++)
            darray_[ix] = (Xptr[ix] * 0.1);
        if (v_size_ > 0)
            for (int iv = 0; iv < natom3_; iv++, ix++)
                darray_[ix] = (Vptr[iv] * 0.1);
        file_.Write( darray_, x_size_ + v_size_ );
    }

    return 0;
}

void Action_Vector::MinImage(Frame const& frm) {
  Matrix_3x3 ucell, recip;
  frm.BoxCrd().ToRecip( ucell, recip );
  Vec3 com1 = frm.VCenterOfMass(mask_);
  Vec_->AddVxyz( MinImagedVec(com1, frm.VCenterOfMass(mask2_), ucell, recip), com1 );
}