Random Polymer on a Grid

(* Define a function to take a random step along the grid in 3d *)
(* Mathematica's built in function RandomChoice[list] acts on a list and picks an item from the list at random with each item in the list being equally likely. *)
RandomGridStep[]:=
RandomChoice[{{1, 0, 0}, {-1, 0, 0}, {0, 1, 0}, {0, -1, 0}, {0, 0, 1}, {0, 0, -1}}]

(* Now you can call our new function to get the random step along the grid *)
RandomGridStep[]

(* Or you can make a table of M steps *)
M = 500; 
S = M-1; (* M monomers requires M-1 steps *)
MSteps = Table[RandomGridStep[], {i, 1, S}]

(* Mathematica has a function called Accumulate[] that adds up all the previous elements of a list. Let's see this in an example: *)
Accumulate[
{1, 1, 1, 1}
]

(* Or another example *)
Accumulate[
{1, 2, 3}
]

(* Notice that MSteps is a list of each independent step. If I act on that list with accumulate, I'll get the position of monomers number 2 through M. Monomer number 1's position is what we define to be our origin, so I should add that to the beginning of the list with the function prepend. *)
RandomPolymer = Prepend[
Accumulate[
MSteps
],
{0, 0, 0}
]

(* Let's plot our polymer *)
Show[
ListPointPlot3D[RandomPolymer],
Graphics3D[Line[RandomPolymer]]
]

(* Define a function to compute the end-to-end-distance from a polymer *)

(* Norm is a built in Mathematica function that computes the norm (length) of a vector *)
R[Polymer_]:=Norm[ 
RandomPolymer[[-1]]-RandomPolymer[[1]] (* This is the vector that connects the ends of the polymer *)
];

R[RandomPolymer]

(* Collect the end-to-end distances for a bunch of samples random polymers *)
NumberOfPolymersToSample = 50000;

SampledRandomPolymerDistances = Table[
MSteps = Table[RandomGridStep[], {i, 1, S}];
RandomPolymer = Prepend[
Accumulate[
MSteps
],
{0, 0, 0}
];
(* Make a list of the end to end distance for each random polymer *)
R[RandomPolymer],
{PolymerNumber, 1, NumberOfPolymersToSample}
];

(* Plot the histogram of the distances *)
Show[
Histogram[SampledRandomPolymerDistances, 30],
Frame->True,
AspectRatio->1,
ImageSize->300,
Axes->False,
FrameLabel->{Style["|R|", FontSize->18],
Style["P(R)", FontSize->18], Style["P(R) for polymers on grid in 3d", FontSize->18]}
]

(* Plot the histogram of the squared end-to-end distances *)
Show[
Histogram[(SampledRandomPolymerDistances^2), 30],
Frame->True,
AspectRatio->1,
ImageSize->300,
Axes->False,
FrameLabel->{Style["(R^2)", FontSize->18],
Style["P((R^2))", FontSize->18], Style["P((R^2)) for polymers on grid in 3d", FontSize->18]}
]

(* Find the sample mean of that distribution, which you should be able to compare with Problem 1 (vii) *)
N[Mean[(SampledRandomPolymerDistances^2)]]