Part 1

Number of Spins

nsp=3;

Input of Parameters

freqs= {372.2,364.4,342,6.083,5.8};
JJ=(({{0, .91, 17.9, 1, 1, 1}, {0, 0, 11.75, 1, 1, 1}, {0, 0, 0, 1, 1, 1}, {0, 0, 0, 0, 1, 1}, {0, 0, 0, 0, 0, 1}, {0, 0, 0, 0, 0, 0}}));

Do[ν[i]=freqs[[i]],{i,1,nsp}];
Do[J[i,k]=JJ[[i,k]],{i,1,nsp-1},{k,i+1,nsp}];

Basevectors

base=Apply[φ,IntegerDigits[#,2,nsp]+1&/@Table[j,{j,0,(2^nsp)-1}]/. 2->-1,{1}]

mT - Values

mT=Table[-nsp/2+j,{j,0,nsp}]

Number of lines in the spectrum  (  =  Sum[Binomial[nsp, k] Binomial[nsp, k + 1], {k, 0, nsp - 1}],  because only transitions between states of different mT's give lines on non-zero intensity )

numberoflines=Binomial[2 nsp,nsp-1]

Rule for scalar products

rscp={φ[x__]^2->1,φ[x__]φ[y__]->0}

SpinOperators

cs[a_,j_]:=Module[{t},t=a;b=t[[j]];t[[j]]=-b;t]
Ix[v_,j_]:=v/2/.φ[x___]:>φ@@cs[{x},j]
Iy[v_,j_]:=I v/2/.φ[x___]:>{x}[[j]]φ@@cs[{x},j]
Iz[v_,j_]:=v/2/.φ[x___]:>{x}[[j]]φ[x]

Example

{Ix[φ[1],1],Ix[φ[-1],1],Iy[φ[1],1],Iy[φ[-1],1],Iz[φ[1],1],Iz[φ[-1],1]}

Hamiltonian  - Matrixelements  (Subscript[H, i,j])  for (sub-)base b

HH[b_,i_,j_]:=(Sum[ν[m]b[[j]]Iz[b[[i]],m],{m,1,nsp}]+Sum[J[m,k]b[[j]](Ix[Ix[b[[i]],k],m]+Iy[Iy[b[[i]],k],m]+Iz[Iz[b[[i]],k],m]),{m,1,nsp-1},{k,m+1,nsp}]//Expand)/.rscp

Example

HH[base,1,1]

HH[base,1,2]

HH[base,3,5]

Spinfunctions according to mT - value

wf=Function[x,Select[base,Total[List@@#]==2 x&]]/@mT

Hamilton - Submatrices

HSM[j_]:=(nn=Length[wf[[j]]];Table[HH[wf[[j]],m,n],{m,1,nn},{n,1,nn}])

HSM[2];
%//MatrixForm

HSM[2]/.{ν[1]->-ν,ν[2]->0,ν[3]->ν,J[1,2]->J,J[1,3]->J,J[2,3]->J};

%//MatrixForm

Get Eigenstates  ≡ all sets { freq,  eigenvector  } for different spin-states (mT values)

frevec[n_]:=Module[{es},
es=Eigensystem[HSM[n]];
{#[[1]],#[[2]].wf[[n]]}&/@Transpose[es]
]

eigenstates=frevec/@Range[nsp+1]

Operator for calculating relative intensities

IOP[x_]:=Sum[Ix[x,n],{n,1,nsp}]

Calculating a spectral line  = difference of eigenvalues and intensity

line[a_,b_]:=Module[{},
freq=Abs[a[[1]]-b[[1]]];
m2=(Expand[a[[2]] IOP[b[[2]]]]/.rscp)^2;
{freq,(Sqrt[m2])}]

Lorentzfunction

LF[x_,x0_,a_,h_]:=Module[{},
If[h==0,h=1];
1/(Sqrt[Pi]) ((a h/2)/(h^2/4 +Pi(x-x0)^2))]

Calculating the spectrum

spec=Table[0,{numberoflines}];
nL=0;
Do[
lk=Length[eigenstates[[i]]];
lk1=Length[eigenstates[[i+1]]];
Do[
Do[
nL=nL+1;
spec[[nL]]=line[eigenstates[[i,m]],eigenstates[[i+1,n ]]],
{n,1,lk1}
],
{m,1,lk}
],
{i,1,Length[eigenstates]-1}];
normalizer=Max[Transpose[spec][[2]]];
bb={.95 Min[Transpose[spec][[1]]],1.05 Max[Transpose[spec][[1]]]};
spec={#[[1]],#[[2]]/normalizer}&/@spec;
spec
pl1=ListPlot[spec,Filling->Axis]

Show the spectrum with lines

pl2=Plot[Plus@@(LF[x,#[[1]],#[[2]],1.5]&/@spec),{x,bb[[1]],bb[[2]]},PlotRange->All,AxesOrigin->{320,0}]

For some physical reason spectra are recorded so that frequencies grow from right to left. So the plot is reversed and compared to the experimental spectrum ( see http://www.users.csbsju.edu/~frioux/nmr/Speclab4.htm ) which is given below the plot.

pl3=Plot[Plus@@(LF[-x,#[[1]],#[[2]],1.5]&/@spec),{x,-bb[[2]],-bb[[1]]},PlotRange->All]

;

(*********************************************)

(*********************************************)

Part 2

nsp=3;

Input of Parameters

freqs= {372.2,364.4,342,6.083,5.8};
JJ=(({{0, .91, 17.9, 1, 1, 1}, {0, 0, 11.75, 1, 1, 1}, {0, 0, 0, 1, 1, 1}, {0, 0, 0, 0, 1, 1}, {0, 0, 0, 0, 0, 1}, {0, 0, 0, 0, 0, 0}}));
Do[ν[i]=freqs[[i]],{i,1,nsp}];
Do[J[i,k]=JJ[[i,k]],{i,1,nsp-1},{k,i+1,nsp}];

number of lines to be expected and dimension ot (total) spin - space (at least for spin 1/2 )

numberoflines=Binomial[2 nsp,nsp-1]
dimspsp=(2^nsp)

spin operators for spin I = 1/2

ix=(({{0, 1}, {1, 0}}))/2;iy=(({{0, -I}, {I, 0}}))/2;iz=(({{1, 0}, {0, -1}}))/2;

this function constructs the spin operator of particle j as matrix in the spin space of n particles

Op[op_,n_,j_]:=Module[{x,m},
x=Join[Table[{{1,0},{0,1}},{j-1}],{op},Table[{{1,0},{0,1}},{n-j}]];
m=SparseArray[KroneckerProduct[Sequence@@x]]
]

oIx[j_]:=Op[ix,nsp,j]
oIy[j_]:=Op[iy,nsp,j]
oIz[j_]:=Op[iz,nsp,j]

Hamiltonian

HH=(∑_(i=1)^nsp)ν[i]oIz[i]+(∑_(i=1)^(nsp-1))(∑_(j=i+1)^nsp)J[i,j](oIx[i].oIx[j]+oIy[i].oIy[j]+oIz[i].oIz[j])

Eigensystem for the Hamilton - Operator

est=Transpose[Eigensystem[HH]]

Intensityoperator

IOP1=(∑_(j=1)^nsp)oIx[j]

line1[a_,b_]:=Module[{},
freq=Abs[a[[1]]-b[[1]]];
m2=(a[[2]].IOP1.b[[2]])^2;
{freq,(Sqrt[m2])}]

Calculate spectrum, show it and display result from part 1

spec1=Table[0,{Binomial[dimspsp,2]}];
nL=0;
Do[
Do[
nL=nL+1;
spec1[[nL]]=line1[est[[u]],est[[v]]],
{u,v+1,dimspsp}
],
{v,1,dimspsp-1}
]
spec1=Select[spec1,#[[2]]>0.&];
normalizer=Max[Transpose[spec1][[2]]];
bb={.95 Min[Transpose[spec][[1]]],1.05 Max[Transpose[spec][[1]]]};
spec1={#[[1]],#[[2]]/normalizer}&/@spec1;
spec1
ListPlot[spec1,Filling->Axis,FillingStyle->Directive[Red,Thick]]
Show[pl1]

Plot the spectrum and compare with the result of the 1 st part

pl4=Plot[Plus@@(LF[x,#[[1]],#[[2]],1.5]&/@spec1),{x,bb[[1]],bb[[2]]},PlotRange->All,PlotStyle->Red]
Show[pl2]