C++ Matrix_3x3类代码示例

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

// -----------------------------------------------------------------------------
void Image::WrapToCell0(std::vector<double>& CoordsIn, Frame const& frmIn,
                        AtomMask const& maskIn,
                        Matrix_3x3 const& ucell, Matrix_3x3 const& recip)
{
  double* uFrac = &CoordsIn[0];
  int nUatoms = maskIn.Nselected();
  int idx;
  double* result;
  const double* XYZ;
# ifdef _OPENMP
# pragma omp parallel private(idx, result, XYZ)
  {
# pragma omp for
# endif
  for (idx = 0; idx < nUatoms; idx++)
  {
    result = uFrac + idx*3;
    XYZ = frmIn.XYZ( maskIn[idx] );
    // Convert to frac coords
    recip.TimesVec( result, XYZ );
    // Wrap to primary unit cell
    result[0] = result[0] - floor(result[0]);
    result[1] = result[1] - floor(result[1]);
    result[2] = result[2] - floor(result[2]);
    // Convert back to Cartesian
    ucell.TransposeMult( result, result );
  }
# ifdef _OPENMP
  } // END pragma omp parallel
# endif
}

/** \param a1 First set of XYZ coordinates.
  * \param a2 Second set of XYZ coordinates.
  * \param ucell Unit cell vectors.
  * \param recip Fractional cell vectors.
  * \return the shortest vector from a1 to a2. 
  */
Vec3 MinImagedVec(Vec3 const& a1, Vec3 const& a2,
                  Matrix_3x3 const& ucell, Matrix_3x3 const& recip)
{
  Vec3 f1 = recip * a1;
  Vec3 f2 = recip * a2;
//  mprintf("DEBUG: a1= %g %g %g, f1= %g %g %g\n", a1[0], a1[1], a1[2], f1[0], f1[1], f1[2]);
//  mprintf("DEBUG: a2= %g %g %g, f2= %g %g %g\n", a2[0], a2[1], a2[2], f2[0], f2[1], f2[2]);
  for (unsigned int i = 0; i < 3; i++) {
    f1[i] = f1[i] - floor(f1[i]);
    f2[i] = f2[i] - floor(f2[i]);
  }
  // Self
  Vec3 c1 = ucell.TransposeMult( f1 );
  Vec3 c2 = ucell.TransposeMult( f2 );
  Vec3 minVec = c2 - c1;
  double minDist2 = minVec.Magnitude2();
  // Images
  for (int ix = -1; ix < 2; ix++) {
    for (int iy = -1; iy < 2; iy++) {
      for (int iz = -1; iz < 2; iz++) {
        if (ix != 0 || iy != 0 || iz != 0) { // Ignore a2 self
          Vec3 ixyz(ix, iy, iz);
          c2 = ucell.TransposeMult(f2 + ixyz); // a2 image back in Cartesian space
          Vec3 dxyz = c2 - c1;
          double dist2 = dxyz.Magnitude2();
          if (dist2 < minDist2) {
            minDist2 = dist2;
            minVec = dxyz;
          }
        }
      }
    }
  }
  return minVec;
}

/** \param Coord Coordinate to image.
  * \param truncoct If true, image in truncated octahedral shape.
  * \param origin If true, image w.r.t. coordinate origin.
  * \param ucell Unit cell matrix.
  * \param recip Reciprocal coordinates matrix.
  * \param fcom If truncoct, image translated coordinate w.r.t. this coord.
  * \return Vector containing image translation.
  */
Vec3 Image::Nonortho(Vec3 const& Coord, bool truncoct, 
                     bool origin, Matrix_3x3 const& ucell, Matrix_3x3 const& recip, 
                     Vec3 const& fcom, double min)
{
  int ixyz[3];

  Vec3 fc = recip * Coord;

  if ( origin )
    fc += 0.5; 

  Vec3 boxTransOut = ucell.TransposeMult( Vec3(floor(fc[0]), floor(fc[1]), floor(fc[2])) );
  boxTransOut.Neg();

  // Put into familiar trunc. oct. shape
  if (truncoct) {
    Vec3 TransCoord = recip * (Coord + boxTransOut);
    Vec3 f2 = recip * fcom;

    if (origin) {
      TransCoord += 0.5;
      f2 += 0.5;
    }

    DIST2_ImageNonOrthoRecip(TransCoord, f2, min, ixyz, ucell);
    if (ixyz[0] != 0 || ixyz[1] != 0 || ixyz[2] != 0) {
      boxTransOut += ucell.TransposeMult( ixyz );
      //if (debug > 2)
      //  mprintf( "  IMAGING, FAMILIAR OFFSETS ARE %i %i %i\n", 
      //          ixyz[0], ixyz[1], ixyz[2]);
    }
  }
  return boxTransOut;
}

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

// MomentOfInertia()
static void MomentOfInertia(int natom, const double *X_, const double* Mass_, double* pmom)
{
  Matrix_3x3 IVEC;
  Vec3 eval;
  // Center of mass
  double cx = 0.0;
  double cy = 0.0;
  double cz = 0.0;
  double sumMass = 0.0;
  const double* mass = Mass_;
  int natom3 = natom * 3;
  for (int i = 0; i < natom3; i+=3) {
    sumMass += (*mass);
    cx += ( X_[i  ] * (*mass) );
    cy += ( X_[i+1] * (*mass) );
    cz += ( X_[i+2] * (*mass) );
    ++mass;
  }
  cx /= sumMass; 
  cy /= sumMass; 
  cz /= sumMass;

  // Moment of inertia
  double xx = 0.0; 
  double yy = 0.0; 
  double zz = 0.0; 
  double xy = 0.0; 
  double xz = 0.0; 
  double yz = 0.0; 
  mass = Mass_;
  for (int i = 0; i < natom3; i+=3) {
    double dx = X_[i  ] - cx;
    double dy = X_[i+1] - cy;
    double dz = X_[i+2] - cz;

    xx += *mass * ( dy * dy + dz * dz );
    yy += *mass * ( dx * dx + dz * dz );
    zz += *mass * ( dx * dx + dy * dy );
    xy -= *mass * dx * dy;
    xz -= *mass * dx * dz;
    yz -= *(mass++) * dy * dz;
  }
  IVEC[0] = xx;
  IVEC[1] = xy;
  IVEC[2] = xz;
  IVEC[3] = xy;
  IVEC[4] = yy;
  IVEC[5] = yz;
  IVEC[6] = xz;
  IVEC[7] = yz;
  IVEC[8] = zz;

  // NOTE: Diagonalize sorts evals/evecs in descending order, but
  //       thermo() expects ascending.
  IVEC.Diagonalize_Sort( eval );
  pmom[0] = eval[2];
  pmom[1] = eval[1];
  pmom[2] = eval[0];
}

/** Set box from unit cell matrix. */
void Box::SetBox(Matrix_3x3 const& ucell) {
  Vec3 x_axis = ucell.Row1();
  Vec3 y_axis = ucell.Row2();
  Vec3 z_axis = ucell.Row3();
  box_[0] = x_axis.Normalize(); // A
  box_[1] = y_axis.Normalize(); // B
  box_[2] = z_axis.Normalize(); // C
  box_[3] = y_axis.Angle( z_axis ) * Constants::RADDEG; // alpha
  box_[4] = x_axis.Angle( z_axis ) * Constants::RADDEG; // beta
  box_[5] = x_axis.Angle( y_axis ) * Constants::RADDEG; // gamma
  SetBoxType();
}

/** Use box coordinates to calculate unit cell and fractional matrix for use
  * with imaging routines. Return cell volume.
  */
double Box::ToRecip(Matrix_3x3& ucell, Matrix_3x3& recip) const {
  double u12x,u12y,u12z;
  double u23x,u23y,u23z;
  double u31x,u31y,u31z;
  double volume,onevolume;
  // If box lengths are zero no imaging possible
  if (box_[0]==0.0 || box_[1]==0.0 || box_[2]==0.0) {
    ucell.Zero();
    recip.Zero();
    return -1.0;
  }
  ucell[0] = box_[0]; // u(1,1)
  ucell[1] = 0.0;     // u(2,1)
  ucell[2] = 0.0;     // u(3,1)
  ucell[3] = box_[1]*cos(Constants::DEGRAD*box_[5]); // u(1,2)
  ucell[4] = box_[1]*sin(Constants::DEGRAD*box_[5]); // u(2,2)
  ucell[5] = 0.0;                                    // u(3,2)
  ucell[6] = box_[2]*cos(Constants::DEGRAD*box_[4]);
  ucell[7] = (box_[1]*box_[2]*cos(Constants::DEGRAD*box_[3]) - ucell[6]*ucell[3]) / ucell[4];
  ucell[8] = sqrt(box_[2]*box_[2] - ucell[6]*ucell[6] - ucell[7]*ucell[7]);

  // Get reciprocal vectors
  u23x = ucell[4]*ucell[8] - ucell[5]*ucell[7];
  u23y = ucell[5]*ucell[6] - ucell[3]*ucell[8];
  u23z = ucell[3]*ucell[7] - ucell[4]*ucell[6];
  u31x = ucell[7]*ucell[2] - ucell[8]*ucell[1];
  u31y = ucell[8]*ucell[0] - ucell[6]*ucell[2];
  u31z = ucell[6]*ucell[1] - ucell[7]*ucell[0];
  u12x = ucell[1]*ucell[5] - ucell[2]*ucell[4];
  u12y = ucell[2]*ucell[3] - ucell[0]*ucell[5];
  u12z = ucell[0]*ucell[4] - ucell[1]*ucell[3];
  volume=ucell[0]*u23x + ucell[1]*u23y + ucell[2]*u23z;
  onevolume = 1.0 / volume;

  recip[0] = u23x*onevolume;
  recip[1] = u23y*onevolume;
  recip[2] = u23z*onevolume;
  recip[3] = u31x*onevolume;
  recip[4] = u31y*onevolume;
  recip[5] = u31z*onevolume;
  recip[6] = u12x*onevolume;
  recip[7] = u12y*onevolume;
  recip[8] = u12z*onevolume;

  return volume;
}

void Action_Vector::Principal(Frame const& currentFrame) {
  Matrix_3x3 Inertia;
  Vec3 Eval;

  // Origin is center of atoms in mask_ 
  Vec3 OXYZ = currentFrame.CalculateInertia( mask_, Inertia );
  // NOTE: Diagonalize_Sort_Chirality places sorted eigenvectors in rows.
  Inertia.Diagonalize_Sort_Chirality( Eval, 0 );
  // Eval.Print("PRINCIPAL EIGENVALUES");
  // Inertia.Print("PRINCIPAL EIGENVECTORS (Rows)");
  if ( mode_ == PRINCIPAL_X ) 
    Vec_->AddVxyz( Inertia.Row1(), OXYZ ); // First row = first eigenvector
  else if ( mode_ == PRINCIPAL_Y )
    Vec_->AddVxyz( Inertia.Row2(), OXYZ ); // Second row = second eigenvector
  else // PRINCIPAL_Z
    Vec_->AddVxyz( Inertia.Row3(), OXYZ ); // Third row = third eigenvector
}

/** Fill the translate vector array with offset values based on this
  * unit cell. Only need forward direction, so no -Z.
  */
void PairList::FillTranslateVec(Matrix_3x3 const& ucell) {
  int iv = 0;
  for (int i3 = 0; i3 < 2; i3++)
    for (int i2 = -1; i2 < 2; i2++)
      for (int i1 = -1; i1 < 2; i1++)
        translateVec_[iv++] = ucell.TransposeMult( Vec3(i1, i2, i3) );
  //for (int i = 0; i < 18; i++)
  //  mprintf("TRANVEC %3i%12.5f%12.5f%12.5f\n", i+1, translateVec_[i][0],
  //          translateVec_[i][1], translateVec_[i][2]);
}

/** Calculate direct space energy. This is the slow version that doesn't
  * use a pair list; for debug purposes only.
  */
double Ewald::Direct(Matrix_3x3 const& ucell, Topology const& tIn, AtomMask const& mask)
{
  t_direct_.Start();
  double cut2 = cutoff_ * cutoff_;
  double Eelec = 0.0;
  Varray const& Image = pairList_.ImageCoords();
  Varray const& Frac = pairList_.FracCoords();
  unsigned int maxidx = Image.size();
  for (unsigned int idx1 = 0; idx1 != maxidx; idx1++)
  {
    // Set up coord for this atom
    Vec3 const& crd1 = Image[idx1];
    // Set up exclusion list for this atom
    int atom1 = mask[idx1];
    Atom::excluded_iterator excluded_atom = tIn[atom1].excludedbegin();
    for (unsigned int idx2 = idx1 + 1; idx2 != maxidx; idx2++)
    {
      int atom2 = mask[idx2];
      // If atom is excluded, just increment to next excluded atom.
      if (excluded_atom != tIn[atom1].excludedend() && atom2 == *excluded_atom) {
        ++excluded_atom;
        //mprintf("ATOM: Atom %4i to %4i excluded.\n", atom1+1, atom2+1);
      } else {
        // Only need to check nearest neighbors.
        Vec3 const& frac2 = Frac[idx2];
        for (Varray::const_iterator ixyz = Cells_.begin(); ixyz != Cells_.end(); ++ixyz)
        {
          Vec3 dxyz = ucell.TransposeMult(frac2 + *ixyz) - crd1;
          double rij2 = dxyz.Magnitude2();
          if ( rij2 < cut2 ) {
            double rij = sqrt( rij2 );
            // Coulomb
            double qiqj = Charge_[idx1] * Charge_[idx2];
            t_erfc_.Start();
            //double erfc = erfc_func(ew_coeff_ * rij);
            double erfc = ERFC(ew_coeff_ * rij);
            t_erfc_.Stop();
            double e_elec = qiqj * erfc / rij;
            Eelec += e_elec;
            //mprintf("EELEC %4i%4i%12.5f%12.5f%12.5f%3.0f%3.0f%3.0f\n",
//            mprintf("EELEC %6i%6i%12.5f%12.5f%12.5f\n", atom1, atom2, rij, erfc, e_elec);
            // TODO can we break here?
          } //else
            //mprintf("ATOM: Atom %4i to %4i outside cut, %6.2f > %6.2f %3.0f%3.0f%3.0f\n",
            //mprintf("ATOM: Atom %4i to %4i outside cut, %6.2f > %6.2f\n",
            //         atom1, atom2,sqrt(rij2),cutoff_);
        }
      }
    }
  }
  t_direct_.Stop();
  return Eelec;
}

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

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

//.........这里部分代码省略.........
  //        ground electronic state is taken to be the zero of
  //        electronic energy.

  //     for monatomics print and return.
  if (avgcrd_.size() <= 3){
    outfile.Printf("\n internal energy:   %10.3f joule/mol         %10.3f kcal/mol\n",
           etran, etran * tokcal);
    outfile.Printf(  " entropy:           %10.3f joule/k-mol       %10.3f cal/k-mol\n",
           stran, stran * tocal);
    outfile.Printf(  " heat capacity cv:  %10.3f joule/k-mol       %10.3f  cal/k-mol\n",
           ctran, ctran * tocal);
    return 0;
  }

  Frame AVG;
  AVG.SetupFrameXM( avgcrd_, mass_ );
  // Allocate workspace memory
  // vtemp   vibrational temperatures, in kelvin.
  // evibn   contribution to e from the vibration n.
  // cvibn   contribution to cv from the vibration n.
  // svibn   contribution to s from the vibration n.
  double* WorkSpace = new double[ 4 * nmodes_ ];
  double* vtemp = WorkSpace;
  double* evibn = WorkSpace + nmodes_;
  double* cvibn = WorkSpace + nmodes_*2;
  double* svibn = WorkSpace + nmodes_*3;

  //     compute contributions due to rotation.

  //     Compute the principal moments of inertia, get the rotational
  //     symmetry number, see if the molecule is linear, and compute
  //     the rotational temperatures.  Note the imbedded conversion
  //     of the moments to SI units.
  Matrix_3x3 Inertia;
  AVG.CalculateInertia( AtomMask(0, AVG.Natom()), Inertia );
  // NOTE: Diagonalize_Sort sorts evals/evecs in descending order, but
  //       thermo() expects ascending.
  // pmom      principal moments of inertia, in amu-bohr**2 and in ascending order.
  Vec3 pmom;
  Inertia.Diagonalize_Sort( pmom );
  rtemp = pmom[0];
  pmom[0] = pmom[2];
  pmom[2] = rtemp;
  outfile.Printf("\n principal moments of inertia (nuclei only) in amu-A**2:\n");
  outfile.Printf(  "      %12.2f%12.2f%12.2f\n", pmom[0], pmom[1], pmom[2]);

  bool linear = false;
  // Symmetry number: only for linear molecules. for others symmetry number is unity
  double sn = 1.0;
  if (AVG.Natom() <= 2) {
    linear = true;
    if (AVG.Mass(0) == AVG.Mass(1)) sn = 2.0;
  }
  outfile.Printf("\n rotational symmetry number %3.0f\n", sn);

  double con = planck / (boltz*8.0*pipi);
  con = (con / tokg)  *  (planck / (tomet*tomet));
  if (linear) {
    rtemp = con / pmom[2];
    if (rtemp < 0.2) {
      outfile.Printf("\n Warning-- assumption of classical behavior for rotation\n");
      outfile.Printf(  "           may cause significant error\n");
    }
    outfile.Printf("\n rotational temperature (kelvin) %12.5f\n", rtemp);
  } else {
    rtemp1 = con / pmom[0];