SetDirectory[NotebookDirectory[]];

Mathematica notebook for the kinetic Monte Carlo concerning polypseudorotaxane assembly between α-cyclodextrin and end functionalized PEG in dilute solutions. This notebook concerns the assembly of a PEG chain with variable threading barriers on the chain end groups.

For description of the theory, please refer to:
C. D.; Ke, C.; Zhang, W. Soft Matter 2024, xx, xx-xx.

NR[]:=nR-Total[chains]; (*Total number of rings in bulk*)
volM=boxVol/1000;(*factor of 1000 accounts for L to m^3*)
totalAggs[data_]:=Module[{splits,aggSizes},
splits=Split[data];
aggSizes=DeleteCases[Table[Total[splits[[n]]],{n,Length[splits]}],0];
Table[Count[aggSizes,n],{n,1,spots}]
]
KMC[]:=Module[{},
rBulk=NR[];
middle=chains[[2;;-2]];
edges=Union[SequencePosition[chains[[2;;-2]],{0,1}][[;;,2]]+1,SequencePosition[chains[[2;;-2]],{1,0}][[;;,1]]+1];
rm=Length[edges];
countAggs=totalAggs[middle];
blockE1=chains[[2]]-chains[[1]];
blockE2=chains[[-2]]-chains[[-1]];
diffuse=(6*(diffaCD+diffPEG)*l0)/volM;
diffuseOut=(6*(diffaCD+diffPEG))/l0^2;
r1Din=diffuse*Abs[chains[[1]]-1]*rBulk;
r1Dout=diffuseOut*chains[[1]];
r1=If[blockE1==-1,k1,0];
r1m=If[blockE1==1,k1m ,0];
r2Din=diffuse*Abs[chains[[-1]]-1]*rBulk;
r2Dout=diffuseOut*chains[[-1]];
r2=If[blockE2==-1,k2,0];
r2m=If[blockE2==1,k2m,0];
r3=If[rm<chainLattice,khop*rm,0];
rTot=r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3;
p=RandomReal[];
notdone=True;
	Which[
p<= r1Din/rTot,(*Add CD to threading point on end group 1*)
		While[notdone,
		b=1;
			If[chains[[b]]==0,
			chains[[b]]=1;
			notdone=False;]],
r1Din/rTot<p<= (r1Din+r1Dout)/rTot,(*Remove CD from threading point on end group 1*)
		While[notdone,
		b=1;
			If[chains[[b]]==1,
			chains[[b]]=0;
			notdone=False;]],
(r1Din+r1Dout)/rTot<p<= (r1Din+r1Dout+r1)/rTot,(*Thread a CD on end group 1*)
		While[notdone,
		b=1;
			If[chains[[b]]==1&&chains[[2]]==0,
			chains[[b]]=0;
			chains[[2]]=1;
			notdone=False;]],
(r1Din+r1Dout+r1)/rTot<p<=  (r1Din+r1Dout+r1+r1m)/rTot,(*Dethread from an end group 1, place in threading point*)
		While[notdone,
		b=2;
			If[chains[[b]]==1&&chains[[1]]==0,
			chains[[b]]=0;
			chains[[1]]=1;
			notdone=False;]],
(r1Din+r1Dout+r1+r1m)/rTot<p<=(r1Din+r1Dout+r1+r1m+r2Din)/rTot,
		While[notdone,
		b=-1;
			If[chains[[b]]==0,
			chains[[b]]=1;
			notdone=False;,
			notdone=True;]],
(r1Din+r1Dout+r1+r1m+r2Din)/rTot<p<= (r1Din+r1Dout+r1+r1m+r2Din+r2Dout)/rTot,
		While[notdone,
		b=-1;
			If[chains[[b]]==1,
			chains[[b]]=0;
			notdone=False;]],
(r1Din+r1Dout+r1+r1m+r2Din+r2Dout)/rTot<p<= (r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2)/rTot,
		While[notdone,
		b=-1;
			If[chains[[b]]==1&&chains[[-2]]==0,
			chains[[b]]=0;
			chains[[-2]]=1;
			notdone=False;]],
(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2)/rTot<p<= (r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m)/rTot,
		While[notdone,
		b=-2;
			If[chains[[b]]==1&&chains[[-1]]==0,
			chains[[b]]=0;
			chains[[-1]]=1;
			notdone=False;]],
 (r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m)/rTot<p,
		While[notdone,
		b=RandomChoice[edges];
			Which[
			b==2&&chains[[b]]==1&&chains[[b+1]]==1,
			notdone=True;,

			b==2&&chains[[b]]==1&&chains[[b+1]]==0,
			newloc=3;
			chains[[newloc]]=1;
			chains[[b]]=0;
			notdone=False;,

			b==NPEO-1&&chains[[b]]==1&&chains[[b-1]]==1,
			notdone=True;,
			b==NPEO-1&&chains[[b]]==1,
			newloc=NPEO-2;
			chains[[newloc]]=1;
			chains[[b]]=0;
			notdone=False;,

			chains[[b-1]]==1&&chains[[b]]==1&&chains[[b+1]]==1,
			notdone=True;,

			chains[[b]]==1&&chains[[b+1]]==1,
			chains[[b-1]]=1;
			chains[[b]]=0;
			notdone=False;,

			chains[[b]]==1&&chains[[b-1]]==1,
			chains[[b+1]]=1;
			chains[[b]]=0;
			notdone=False;,

			chains[[b]]==1,
			newloc=b+RandomChoice[{-1,1}];
			chains[[newloc]]=1;
			chains[[b]]=0;
			notdone=False;,

			chains[[b]]==0,
			notdone=True]]];

rr=RandomReal[];
time=time-Log[rr]/rTot;
If[countAggs[[10]]>=1&&Length[aggTime[[1]]]<=1,AppendTo[aggTime[[1]],time]];
If[i>=5000000000,Break[]]]

MultiStepKMC[PEGconc_,PEGlength_,aCDconc_,barrier1_,barrier2_,latticeSizeNM_,outputf_]:=Module[{},
spots=Round[PEGlength/2];
NPEO=spots+2; (*chain Length*)
chainLattice=NPEO-2;
l0=latticeSizeNM;
boxVol=(1/(PEGconc*6.022*10^23));
nR=Round[boxVol*6.022*10^23*aCDconc]; (*total # of CD*)
diffaCD=0.353*10^-9 (1-0.75(aCDconc));(*Ribeiro et al. J. Chem. Eng. Data 2007*)
diffPEG=((1.38*10^-23)(300))/((0.851*10^-3)(10^-9)(PEGlength/3)^0.6);
k1=9.0*10^9*Exp[-barrier1];(*Barrier for threading end 1*)
k1m=9.0*10^9 *Exp[-barrier1-8]; (*Barrier for dethreading end 1*)
k2=9.0*10^9*Exp[-barrier2];(*Barrier for threading end 2*)
k2m=9.0*10^9*Exp[-barrier2-8]; (*Barrier for threading end 2*)
khop=9.0*10^9*Exp[-4.4]; (*hopping along backbone*)
chains=Table[0,NPEO];
countAggs=totalAggs[chains[[2;;-2]]];
ϕList={};
timeList={};
confList={};
aggList={};
aggTime={{10}};
r=nR;
time=0;
SeedRandom[];
For[i=1,countAggs[[10]]==0,i++,
KMC[];If[Mod[i,outputf]==0,
AppendTo[ϕList,Total[middle]/chainLattice];
AppendTo[aggList,countAggs];
AppendTo[timeList,time];
AppendTo[confList,chains];]]]

MultiStepKMC[PEGconc_, PEGlength_, aCDconc_, barrier1_, barrier2_, latticeSizeNM_, outputf_]

Monitor[MultiStepKMC[0.001,22,0.05,5,5,(10^(-9)),20],time]

threaded=Transpose@{timeList,ϕList};
aggs=Transpose@{timeList,aggList};
confs=Table[confList[[n]],{n,1,Length[confList],5}];

ListLinePlot[threaded,Frame->True,FrameLabel->{"Time (s)","PSR coverage fraction"},FrameStyle->Directive[Black,Thickness[0.0025]],ImageSize->400,LabelStyle->Directive[Black,14]]

positions=Transpose@{{1,2,7,12,13},{"TP","1","6","11","TP"}};

graphic=Table[MatrixPlot[{confList[[n]]},PlotLabel->"Time = "<>ToString[Round[timeList[[n]]*(10^6),0.1]]<>" μs",FrameStyle->Directive[Black,Thickness[0.0025]],LabelStyle->Directive[Black,14],FrameTicks->{{None,None},{positions,positions}},ImageSize->400],{n,1,Length[timeList],10}];

ListAnimate[graphic]

MatrixPlot[confList[[;;Length[confList];;10,2;;12]],FrameLabel->{"← Increasing time","Backbone position"},LabelStyle->Directive[Black,14],FrameStyle->Directive[Black,Thickness[0.0025]],FrameTicks->{{None,None},{{1,6,11},{1,6,11}}},ImageSize->200]