(* Finding dispacements for exc. st. in respect to the ground state *)
 Clear["@"];
dir = "D:\\sergeev\\bgu\\math\\zilberg";
SetDirectory[dir];
filetab = "mail0807.txt";
nstate = 1;
<< readgaus.m;
<< trarot.m;
qatoms1 = qatoms; freqs1 = freqs; basis1 = basis;
nstate = 3;
<< readgaus.m;
<< trarot.m;
dq = qatoms - qatoms1;
dqn = (Flatten /@ basis1).Flatten[dq  sqm];
tab = Table[{n - 6, Chop[dqn[[n]]], 
        If[n <= 6, "", freqs[[n - 6]]/freqs1[[n - 6]]]}, {n, nmodes + 6}];
tab = Prepend[tab, {"Mode", "Displacement", "Ratio of frequencies"}];
Print[tab // TableForm];
scpr = (Flatten /@ basis1).Transpose[(Flatten /@ basis)];
tab = Table[Round[100 scpr[[m, n]]]/100., {m, nmodes + 6}, {n, nmodes + 6}];
tab = Prepend[tab, Table[n - 6, {n, nmodes + 6}]];
tab = Transpose[
      Prepend[Transpose[tab], Prepend[Table[n - 6, {n, nmodes + 6}], ""]]];
Print["Matrix of Duschinsky rotation:"];
Print[tab // TableForm];
Do[
    b = 0;
    Do[
      a1 = scpr[[m + 6, n + 6]];
      b1 = a1^2;
      If[b1 > b, b = b1; a = a1; nmax = n]
      , {n, nmodes}];
    Print["Mode ", m, " - ground -> Mode ", nmax,
      " - excited (", a, "). Ratio of frequencies: ", 
      freqs[[nmax]]/freqs1[[m]]];
    , {m, nmodes}];
(* Save results *)
file = "displace.mx";
freqs2 = freqs;
displ = Chop[Drop[dqn, 6]];
dush = Table[Chop[scpr[[i + 6, j + 6]]], {i, nmodes}, {j, nmodes}];
DeleteFile[file];
Save[file, {nmodes, freqs1, freqs2, displ, dush}];