initg[nn_,s_:1]:=Block[{},
n=nn;(*global for number of monomers*)
xx=s Transpose[{initg1[n],initg1[n],initg1[n]}];
];
initg1[nn_]:=Block[{x1},
(*set random values with variance 1*)
x1=Table[ RandomReal[NormalDistribution[0,1]],{nn}];
x1=Fourier[x1];(*ft*)
x1[[1]]=0;(*not wanted*)
Do[ x1[[i]]/= Sqrt[3.]2 Sin[2 Pi(i-1)/n/2], (*polymerise*)
{i,2,n}
];
x1=InverseFourier[x1]//Re
];
(*badinit*)
(*badinit*)
init[nn_,s_:1]:=Block[{},n=nn;
xx=Table[s(2  Random[] -1),{nn},{3}];
];
energy[xx_]:=Block[{xbf},xbf=xx-RotateRight[xx]//Flatten;1.5  xbf.xbf];
energy2=Compile[{{xx,_Real,2}},Block[{xbf},xbf=xx-RotateRight[xx]//Flatten;1.5  xbf.xbf] ];
sizesq[xx_,disp_:Floor[n/2] ]:=Block[{xxrf},
xxrf=xx-RotateRight[xx,disp]//Flatten;
xxrf.xxrf/n];
mcsimple[moves_,dxm_]:=Block[{p=0},
Do[xx1=xx+ Table[dxm(2  Random[] -1),{n},{3}];
If[Random[]<Exp[energy[xx]-energy[xx1] ],
p+=1;xx=xx1],
{moves}];
p/moves//N]

(*compiled version requires kmax*)
mcshearc=Compile[{{xx,_Real,2},moves,sm,kmax},
Module[{p=0.,k={0.,0.,0.},a={0.,0.,0.},phi=0.,xx1,xx2,phase,xbf,denergy },
xx1=xx;
Do[
(*find random spherical k on unit scale*)
While[k={2 Random[]-1,2 Random[]-1,2 Random[]-1};k.k>1,];
(*find random disp amplitudes on unit scale*)
While[a={2 Random[]-1,2 Random[]-1,2 Random[]-1};a.a>1,];
(*remove longitudinal part*)
a=a-(k.a/k.k)k;
(*scale k to kmax*)
k*=kmax;
(*scale a to sm/k*)
a=sm a/Sqrt[(a.a) (k.k)];
(*random phase*)
phi=2 Pi Random[];
(*construct new coordinates and finish off per simple mc*)
phase=Cos[xx1.k+phi];
xx2=xx1+Map[Times[a,#]&,phase];
(*xx2=xx1+Map[Times[a,#]&,Cos[xx1.k+phi]];*)
(*the above two lines reveal the actual timestep*)
xbf=xx1-RotateRight[xx1]//Flatten;denergy=-1.5  xbf.xbf;
xbf=xx2-RotateRight[xx2]//Flatten;denergy+=1.5  xbf.xbf;
If[Random[]<Exp[-denergy ],
p+=1;xx1=xx2],
{moves}];
xx1]
];
mcshear[moves_,sm_,kmax_:1]:=Block[{p=0,k,a,phi,xx1,phase(*,xxc*)},
Do[
(*find random spherical k on unit scale*)
While[k={2 Random[]-1,2 Random[]-1,2 Random[]-1};k.k>1,];
(*find random disp amplitudes on unit scale*)
While[a={2 Random[]-1,2 Random[]-1,2 Random[]-1};a.a>1,];
(*remove longitudinal part*)
a=a-(k.a/k.k)k;
(*scale k to kmax*)
k*=kmax;
(*scale a to sm/k*)
a=sm a/Sqrt[(a.a) (k.k)];
(*random phase*)
phi=2 Pi Random[];
(*construct new coordinates and finish off per simple mc*)
phase=Cos[xx.k+phi];
xx1=xx+Map[Times[a,#]&,phase];
(*xxc=xx+Table[a Cos[k.xx[[i]]+phi],{i,n}]-xx1//Flatten;Print[xxc.xxc];*)
(*the above two lines reveal the actual timestep*)
If[Random[]<Exp[energy[xx]-energy[xx1] ],
p+=1;xx=xx1],
{moves}];
p/moves//N]

(*utility to generate a time sequence*)
SetDirectory["H:\calcs\wavehydro"]
dynamics[nn_,steps_,outfile_,{mps_,sm_,kmax_}]:=Module[{m,outfilemore},
outfilemore=outfile<>"n="<>ToString[nn]<>"mvs="<>ToString[steps]<>"x"<>ToString[mps]<>"ep="<>ToString[Round[100sm]]<>"k="<>ToString[kmax];
str=OpenWrite[outfilemore];
initg[nn];
Write[str,{0,xx}];
Do[
xx=mcshearc[xx,mps,sm,kmax];
Write[str,{m,xx}],
{m,1,steps}];
Close[str];outfilemore
];
diffplot[infile_]:=Module[{timelist,xxlist,diff,diffr,diffc,count,nn,dx,dxc,dcom},
{time,xxlist}=ReadList[infile]//Transpose;
{name,n1,steps,mps,ep,k}=StringSplit[infile,{"n=","mvs=","x","ep=","k="}]//ToExpression;
tstep=mps (ep/100)^2;
count=diffc=diffr=diff=Table[0,{Length[xxlist]-1}];
nn=Length[xxlist[[1]]];
Do[
count[[j-i]]+=1;
dx=xxlist[[j]]-xxlist[[i]];
dcom=Apply[Plus,dx]/nn ;
dxc=dx-Table[dcom,{nn}] //Flatten;
dx=dx//Flatten;
diff[[j-i]]+=dx.dx/nn;
diffr[[j-i]]+=dxc.dxc/nn;
diffc[[j-i]]+=dcom.dcom;
,{i,Length[xxlist]-1},{j,i+1,Length[xxlist]}];
diff/=count;diffr/=count;diffc/=count;
count=Table[i tstep,{i,Length[diff]}];
ListPlot[Join[{count,diff}//Transpose,{count,diffr}//Transpose,{count,diffc}//Transpose]//Log]
];

dynamics[100,100,"hun2",{100,0.1,Pi}]//Timing

{1.3130000000001019`,"hun2n=100mvs=100x100ep=10k=Pi"}

diffplot[%[[2]]]

dynamics[100,100,"hun2",{100,0.1,2 Pi}]//Timing

diffplot[%[[2]]]

dynamics[100,100,"hun2",{500,0.1,Pi}]//Timing
diffplot[%[[2]]]

diffplot["hun2n=100mvs=100x500ep=10k=Pi"]

dynamics[100,100,"hun2",{500,0.1,2 Pi}]//Timing
diffplot[%[[2]]]

diffplot["hun2n=100mvs=100x500ep=10k=2 Pi"]

dynamics[200,100,"CC",{500,0.1,2 Pi}]//Timing
diffplot[%[[2]]]

dynamics[100,100,"hun2",{500,0.05,2 Pi}]//Timing
diffplot[%[[2]]]

Show[%130,%128,%133]

(*expect 40s*)dynamics[200,100,"CC",{2000,0.1,2 Pi}]//Timing

(*expect 160s*)dynamics[200,100,"CC",{8000,0.1,2 Pi}]//Timing

(*expect 10min*)dynamics[200,100,"CC",{32000,0.1,2 Pi}]//Timing

(*expect 40min*)dynamics[200,100,"CC",{128000,0.1,2 Pi}]//Timing

(*expect 2.5hrs*)dynamics[200,100,"CC",{512000,0.1,2 Pi}]//Timing
(*expect 10hrs*)dynamics[200,100,"CC",{2048000,0.1,2 Pi}]//Timing

Show[diffplot["CCn=200mvs=100x500ep=10k=2 Pi"],diffplot["CCn=200mvs=100x2000ep=10k=2 Pi"],diffplot["CCn=200mvs=100x8000ep=10k=2 Pi"],
diffplot["CCn=200mvs=100x32000ep=10k=2 Pi"],PlotRange->All]

Show[diffplot["CCn=200mvs=100x500ep=10k=2 Pi"],diffplot["CCn=200mvs=100x2000ep=10k=2 Pi"],diffplot["CCn=200mvs=100x8000ep=10k=2 Pi"],
diffplot["CCn=200mvs=100x32000ep=10k=2 Pi"],diffplot["CCn=200mvs=100x128000ep=10k=2 Pi"],PlotRange->All]

(*testread*)
diffplot["hundredn=100mps=100ep=10k=Pi"]

diffplot["hvvlongn=100mps=100000ep=10k=2 Pi"]

diffplot["hvvlongn=200mps=50000ep=5k=2 Pi"]

dynamics[100,100,"hundred",{100,0.1,Pi}]

diffplot["hundredn=100mps=100ep=10k=Pi"]

outfile=dynamics[100,100,"hundred",{100,0.1,2Pi}]
diffplot[outfile]

outfile=dynamics[100,200,"hlonger",{1000,0.1,2Pi}]
plot=diffplot[outfile]

diffplot[outfile]