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

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

//.........这里部分代码省略.........
           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];
    rtemp2 = con / pmom[1];
    rtemp3 = con / pmom[2];
    if (rtemp1 < 0.2) {
      outfile.Printf("\n Warning-- assumption of classical behavior for rotation\n");
      outfile.Printf(  "           may cause significant error\n");
    }