SetDirectory[NotebookDirectory[]];

Mathematica notebook for the kinetic Monte Carlo concerning polypseudorotaxane assembly between α-cyclodextrin and intrachain functionalized PEG in dilute solutions. This notebook concerns the assembly of a PEG chain with four (4) variable intrachain diffusion barriers, as well as 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.

(*Check coverage in each region*)
checkBlock1[]:=Total[chains[[2;;innerPos1]]]
checkBlock2[]:=Total[chains[[innerPos1+1;;innerPos2]]]
checkBlock3[]:=Total[chains[[innerPos2+1;;innerPos3]]]
checkBlock4[]:=Total[chains[[innerPos3+1;;innerPos4]]]
checkBlock5[]:=Total[chains[[innerPos4+1;;-2]]]

(*Rings in bulk solution*)
NR[]:=nR-Total[chains];

(*kMC implementation*)
KMC[]:=Module[{},
rBulk=NR[];

(*Count threaded rings with available neighboring positions - excluding intrachain barriers.
Separating this section into "domains" would allow for the introduction of variable backbone diffusion constants for block copolymers*)
edgesRight=Count[Table[chains[[x;;x+1]],{x,{innerPos1,innerPos2,innerPos3,innerPos4,-3}}],{0,1}];
edgesLeft=Count[Table[chains[[x+2;;x+3]],{x,{0,innerPos1,innerPos2,innerPos3,innerPos4}}],{1,0}];
middle=
Total[Flatten[{Table[SequenceCount[chains[[x+2;;x+spacing+1]],{0,1,0}],{x,{0,innerPos1,innerPos2,innerPos3,innerPos4}}],
Table[SequenceCount[chains[[x+2;;x+spacing+1]],{0,1,1}],{x,{0,innerPos1,innerPos2,innerPos3,innerPos4}}],
Table[SequenceCount[chains[[x+2;;x+spacing+1]],{1,1,0}],{x,{0,innerPos1,innerPos2,innerPos3,innerPos4}}]}]];
rm=edgesRight+edgesLeft+middle;

(*Count rings at intrachain barrier positions*)
fLeft1=chains[[innerPos1]];
fLeft2=chains[[innerPos2]];
fLeft3=chains[[innerPos3]];
fLeft4=chains[[innerPos4]];
fRight1=chains[[innerPos1+1]];
fRight2=chains[[innerPos2+1]];
fRight3=chains[[innerPos3+1]];
fRight4=chains[[innerPos4+1]];

(*Find reactions that are blocked by occupied sites at barriers*)
blockE1=chains[[2]]-chains[[1]];
blockE2=chains[[-2]]-chains[[-1]];
blockI1=chains[[innerPos1+1]]-chains[[innerPos1]];
blockI2=chains[[innerPos2+1]]-chains[[innerPos2]];
blockI3=chains[[innerPos3+1]]-chains[[innerPos3]];
blockI4=chains[[innerPos4+1]]-chains[[innerPos4]];
checkBlocks=ContainsOnly[{checkBlock1[],checkBlock2[],checkBlock3[],checkBlock4[],checkBlock5[]},{0,spacing}];

(*Solution diffusion constants*)
diffuse=(6*(diffaCD+diffPEG)*l0)/volM;
diffuseOut=(6*(diffaCD+diffPEG))/l0^2;

(*Determine all rates*)
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[checkBlocks,0,khop*rm];
r4l2r=If[blockI1==-1,k3*fLeft1,0];
r4r2l=If[blockI1==1,k3*fRight1,0];
r5l2r=If[blockI2==-1,k4*fLeft2,0];
r5r2l=If[blockI2==1,k4*fRight2,0];
r6l2r=If[blockI3==-1,k4*fLeft3,0];
r6r2l=If[blockI3==1,k4*fRight3,0];
r7l2r=If[blockI4==-1,k3*fLeft4,0];
r7r2l=If[blockI4==1,k3*fRight4,0];
rTot=r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r+r5r2l+r6l2r+r6r2l+r7l2r+r7r2l;

(*Select and perform next step*)
p=RandomReal[];
notdone=True;
	Which[
(*Add CD to threading point on end group 1*)
p<= r1Din/rTot,
	While[notdone,
	b=1;
		If[chains[[b]]==0,
		chains[[b]]=1;
		notdone=False;]],

(*Remove CD from threading point on end group 1*)
r1Din/rTot<p<= (r1Din+r1Dout)/rTot,
	While[notdone,
	b=1;
		If[chains[[b]]==1,
		chains[[b]]=0;
		notdone=False;]],

(*Thread a CD on end group 1*)
(r1Din+r1Dout)/rTot<p<= (r1Din+r1Dout+r1)/rTot,
	While[notdone,
	b=1;
		If[chains[[b]]==1&&chains[[2]]==0,
		chains[[b]]=0;
		chains[[2]]=1;
		notdone=False;]],

(*Dethread from an end group 1, place in threading point*)
(r1Din+r1Dout+r1)/rTot<p<=  (r1Din+r1Dout+r1+r1m)/rTot,
	While[notdone,
	b=2;
		If[chains[[b]]==1&&chains[[1]]==0,
		chains[[b]]=0;
		chains[[1]]=1;
		notdone=False;]],

(*Add CD to threading point on end group 2*)
(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;]],

(*Remove CD from threading point on end group 2*)
(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;]],

(*Thread a CD on end group 2*)
(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;]],

(*Dethread from an end group 1, place in threading point*)
(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;]],

(*Backbone diffusion*)
 (r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m)/rTot<p<=(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3)/rTot,
	While[notdone,
	b=RandomInteger[{2,NPEO-1}];
			Which[
			chains[[b]]==0,
			notdone=True,

			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]]==1&&chains[[b]]==1&&b==innerPos1,
			notdone=True;,

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

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

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

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

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

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

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

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

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

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

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

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

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

			chains[[b+1]]==0&&chains[[b]]==1&&b==innerPos4+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&&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]],

(*Intrachain jumping*)
(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3)/rTot<p<=(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r)/rTot,
	While[notdone,
	b=innerPos1;
		If[chains[[b]]==1&&chains[[innerPos1+1]]==0,
		chains[[b]]=0;
		chains[[innerPos1+1]]=1;
		notdone=False;]],

(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r)/rTot<p<=(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l)/rTot,
	While[notdone,
	b=innerPos1+1;
		If[chains[[b]]==1&&chains[[innerPos1]]==0,
		chains[[b]]=0;
		chains[[innerPos1]]=1;
		notdone=False;]],

(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l)/rTot<p<= (r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r)/rTot,
	While[notdone,
	b=innerPos2;
		If[chains[[b]]==1&&chains[[innerPos2+1]]==0,
		chains[[b]]=0;
		chains[[innerPos2+1]]=1;
		notdone=False;]],

(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r)/rTot<p<=(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r+r5r2l)/rTot,
	While[notdone,
	b=innerPos2+1;
		If[chains[[b]]==1&&chains[[innerPos2]]==0,
		chains[[b]]=0;
		chains[[innerPos2]]=1;
		notdone=False;]],

(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r+r5r2l)/rTot<p<=(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r+r5r2l+r6l2r)/rTot,
	While[notdone,
	b=innerPos3;
		If[chains[[b]]==1&&chains[[innerPos3+1]]==0,
		chains[[b]]=0;
		chains[[innerPos3+1]]=1;
		notdone=False;]],

(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r+r5r2l+r6l2r)/rTot<p<=(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r+r5r2l+r6l2r+r6r2l)/rTot,
	While[notdone,
	b=innerPos3+1;
		If[chains[[b]]==1&&chains[[innerPos3]]==0,
		chains[[b]]=0;
		chains[[innerPos3]]=1;
		notdone=False;]],

(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r+r5r2l+r6l2r+r6r2l)/rTot<p<=(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r+r5r2l+r6l2r+r6r2l+r7l2r)/rTot,
	While[notdone,
	b=innerPos4;
		If[chains[[b]]==1&&chains[[innerPos4+1]]==0,
		chains[[b]]=0;
		chains[[innerPos4+1]]=1;
		notdone=False;]],

(r1Din+r1Dout+r1+r1m+r2Din+r2Dout+r2+r2m+r3+r4l2r+r4r2l+r5l2r+r5r2l+r6l2r+r6r2l+r7l2r)/rTot<p,
	While[notdone,
	b=innerPos4+1;
		If[chains[[b]]==1&&chains[[innerPos4]]==0,
		chains[[b]]=0;
		chains[[innerPos4]]=1;
		notdone=False;]]
	];
rr=RandomReal[];
time=time-Log[rr]/rTot;(*Update time*)
If[i>=5000000000,Break[]](*Prevent extensive job time*)
]

MultiStepKMC[PEGconc_,PEGlength_,aCDconc_,barrier1_,barrier2_,interBarrier1_,interBarrier2_,barrierSpacing_,latticeSizeNM_,outputf_]:=Module[{},
spots=Round[PEGlength/2];(*backbone positions*)
spacing=Round[barrierSpacing/2];(*Intrachain barrier positions*)
NPEO=spots+2; (*chain Length*)
chainLattice=NPEO-2;
barriers=Table[n+1,{n,spacing,spots,spacing}];
innerPos1=barriers[[1]];
innerPos2=barriers[[2]];
innerPos3=barriers[[3]];
innerPos4=barriers[[4]];
l0=latticeSizeNM;
boxVol=(1/(PEGconc*6.022*(10^23)));
volM=boxVol/1000;
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);(*Zimm model polymer diffusion*)
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*)
k3=9.0*10^9*Exp[-interBarrier1];(*Barrier for intra-chain jump - secondary regions*)
k4=9.0*10^9*Exp[-interBarrier2];(*Barrier for intra-chain jump - tertiary region*)
khop=9.0*10^9*Exp[-4.4]; (*hopping along backbone*)
chains=Table[0,{j,NPEO}];
ϕList={};
timeList={};
confList={};
aggList={};
r=nR;
rm=0;
time=0;
SeedRandom[];
For[i=1,
Total[chains[[2;;-2]]]/chainLattice<0.8(*cancel run at 80% coverage*),
i++,
KMC[];
If[Mod[i,outputf]==0,
AppendTo[ϕList,Total[chains[[2;;-2]]]/chainLattice];
AppendTo[timeList,time];
AppendTo[confList,chains];]]]

MultiStepKMC[PEGconc_, PEGlength_, aCDconc_, barrierEnd1_, barrierEnd2_, interBarrier1_, interBarrier2_, barrierSpacing_, latticeSizeNM_, outputf_]

MultiStepKMC[0.001,50,0.05,0,0,10,15,10,(10^(-9)),200]

positions=Transpose@{{1,2,7,12,16,21,26,27},{"TP","1","6","11","15","20","25","TP"}};

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

ListAnimate[graphic]

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

(*Batch run and export of data*)
Do[MultiStepKMC[0.001,50,0.05,0,0,15,15,10,(10^(-9)),200];
threaded=Transpose@{timeList,ϕList};
confs=Table[confList[[n]],{n,1,Length[confList],5}];
end=confList[[-1]];
Export["confs_barrier15_20_"<>ToString[n]<>".csv",AppendTo[confs,{time,chains}],"CSV"];
Export["phi_barrier15_20_"<>ToString[n]<>".csv",AppendTo[threaded,{time,Total[chains[[2;;-2]]]/chainLattice}],"CSV"];
Export["end_barrier15_20_"<>ToString[n]<>".csv",end,"CSV"];
,{n,1}]