(* Additional translation-rotation to put center of mass to zero and axes of tensor of inertia to x,y,z *)
(* Transforms data that was parsed from Gaussian output *)

Do[qatoms[[natom]] = qatoms[[natom]] - rrcom, {natom, matoms}];
qatoms = matx.#&/@qatoms;
Do[qmodes[[nmode]] = matx.#&/@qmodes[[nmode]], {nmode, nmodes}];

(* Calculation of coordinates of center of mass, tensor of inertia, ... *)

tra = {d1, d2, d3} = 
      Transpose[Table[{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {matoms}]];
(* Center of mass *)
rrcom = 
      Sum [atmass[[natom]]qatoms[[natom]], {natom, matoms}]/
        Sum [atmass[[natom]], {natom, matoms}];
If[iprint === 1, Print["Coordinates of center of mass = ", rrcom]];
rcom = Table[qatoms[[natom]] - rrcom, {natom, matoms}];
(* Tensor of inertia *)
inert = Sum [{x, y, z} = rcom[[natom]];
      atmass[[natom]]{{y^2 + z^2, -x y, -x z},
          {-y x, x^2 + z^2, -y z}, {-z x, -z y, x^2 + y^2}}, {natom, 
        matoms}];
{eig, matx} = Eigensystem[inert];
rot = {d4, d5, d6} = Transpose[Table[
          {px, py, pz} = matx.rcom[[natom]];
          {x, y, z} = matx;
          {py z - pz y, pz x - px z, px y - py x}, {natom, matoms}]];
basis = (sqm#) & /@ Join[tra, rot, qmodes];
fbasis = GramSchmidt[Flatten /@ basis];
basis = unflatten/@fbasis;
qmodescorr = Table[a = basis[[n + 6]]/sqm;
      a/Sqrt[Flatten[a].Flatten[a]],
      {n, nmodes}];