I tried to make a workable example here. The problem is that it is very very slow. I must be doing something very stupid out here. Can somebody give me some clue how to evaluate this multiple sum faster. Right now it takes around 2 mins. However, for real problem I need the parameter Npol to be of the order of ~10.
Clear["Global`*"];
SetSystemOptions["CompileOptions"->"CompileReportExternal"->False];
DRSH=N[(94688879/10000000)*10^(-7)];
Npol=2;
nk=40;
PN=Evaluate[Table[0.0025 n,{n,1,Npol}]];
OmN=Evaluate[Table[(1.0/DRSH)^((-n+2)/Npol),{n,1,Npol}]];
fpol[n_, On_]:=(-0.38 Sqrt[4^-n])/(DRSH + 0.77 2.0^-n)/; On== OmN[[1]]
fpol[n_, On_]:=(-0.38 Sqrt[4^-n])/(DRSH + 0.77 2.0^-n) (0.77 2.0^-n-On)/(0.77 2.0^-n+On)
OverLaps1=Compile[{{k, _Integer},{n,_Real}, {m,_Real}},Exp[-0.5 Sum[If[q!= k,(fpol[q, n]-fpol[q, m])^2,0],{q,0,nk}]],CompilationTarget->"C"];
OverLaps2=Compile[{{k, _Integer},{n,_Real}, {m,_Real}},Exp[-0.5 Sum[If[q!= k,(fpol[q, n]+fpol[q, m])^2,0],{q,0,nk}]],CompilationTarget->"C"];
NpolSum=Compile[{{i, _Integer},{j, _Integer}, {k, _Integer}}, (Sum[ Sum[PN[[n]] PN[[m]] fpol[k, OmN[[n]]]^i fpol[k, OmN[[m]]]^j Exp[-0.5 (fpol[k, OmN[[n]]]^2+ fpol[k, OmN[[m]]]^2)] ( OverLaps1[k, OmN[[n]], OmN[[m]]]+ (-1.0)^j OverLaps2[k, OmN[[n]], OmN[[m]]]),{n,1,Npol}],{m,1,Npol}])^2,CompilationTarget->"C"];
S2Entrp=Compile[{{k, _Integer}},Sum[Sum[ ( 1/(i! j!) ) NpolSum[i, j, k],{i,0,10}],{j,0,10}]];
ListLogLinearPlot[Table[{k,S2Entrp[ k]}, {k,0, nk}],PlotStyle->Directive[PointSize[Medium],Red], PlotRange-> All]//AbsoluteTiming
I thank you for your time.
Sum[...,{}]
byTotal@Table[...,{}]
? $\endgroup$