diff --git a/mcmc_clustering_eap_chain.jl b/mcmc_clustering_eap_chain.jl
index 5e864f0..1eaab54 100644
--- a/mcmc_clustering_eap_chain.jl
+++ b/mcmc_clustering_eap_chain.jl
@@ -239,7 +239,10 @@ function mcmc(nsteps::Int, pargs, chain::EAPChain)
                       ]);
 
   outfile = open("$(pargs["prefix"])_trajectory.csv", "w");
-  writedlm(outfile, ["step" "r1" "r2" "r3" "p1" "p2" "p3" "U"], ',');
+  writedlm(outfile, hcat(["step" "r1" "r2" "r3" "p1" "p2" "p3" "U"], 
+                         reshape(vcat(reshape(["phi$i" for i=1:n(chain)], 1, :), reshape(["theta$i" for i=1:n(chain)], 1, :)), 1, :), 
+                         reshape(vcat(reshape(["mux$i" for i=1:n(chain)], 1, :), reshape(["muy$i" for i=1:n(chain)], 1, :), reshape(["muz$i" for i=1:n(chain)], 1, :)), 1, :)), 
+           ',');
   rollfile = open("$(pargs["prefix"])_rolling.csv", "w");
   writedlm(rollfile, ["step" "r1" "r2" "r3" "r1sq" "r2sq" "r3sq" "rsq" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psq" "U" "Usq" "Ealign" "psi"], ',');
 
@@ -297,7 +300,10 @@ function mcmc(nsteps::Int, pargs, chain::EAPChain)
     if step % pargs["stepout"] == 0
       writedlm(outfile, 
                hcat(step, transpose(chain.r), 
-                    transpose(chain_μ(chain)), chain.U), 
+                    transpose(chain_μ(chain)), chain.U, 
+                    reshape(vcat(transpose(chain.ϕs), transpose(chain.θs)), 1, :),
+                    reshape(chain.μs, 1, :)
+                   ), 
                ',');
       writedlm(rollfile, 
                hcat(
@@ -338,9 +344,11 @@ end
   mcmc(5, pargs); # run first to compile code
   #@profview mcmc(pargs["num-steps"], pargs);
   error("Not currently implemented...");
+  #=
   println("Press ENTER to continue...");
   readline(stdin);
   exit(1);
+  =#
 else
   burned_in_chain = mcmc(pargs["burn-in"], pargs)[1];
   mcmc(pargs["num-steps"], pargs, burned_in_chain);
diff --git a/run/phases-random_2023-07-11.jl b/run/phases-random_2023-07-11.jl
new file mode 100644
index 0000000..d0d900b
--- /dev/null
+++ b/run/phases-random_2023-07-11.jl
@@ -0,0 +1,51 @@
+println("Hello world!");
+
+@everywhere using Printf;
+@everywhere using Distributions;
+
+@show workdir = if length(ARGS) > 0
+  ARGS[1];
+else
+  ".";
+end
+
+@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x));
+@everywhere fmt_int(x) = @sprintf("%03d", x);
+
+@everywhere function prefix(case)
+  "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))_kappa-$(fmt(case[:kappa]))";
+end
+
+@everywhere Kdist = collect(zip([1.0; 0.0], [0.0; 1.0], [2.0; 1.0], [1.0; 2.0], [3.0; 2.0], [2.0; 3.0]));
+@everywhere Edist = Normal(2.5, 1.0);
+@everywhere ndist = Normal(50, 15);
+@everywhere bdist = Normal(1.0, 1.0);
+@everywhere kTdist = Uniform(0.0, 2.0);
+@everywhere Fdist = Normal(2.0, 1.0);
+
+mkpath(workdir);
+
+@sync @distributed for i=1:1000000
+
+  case = Dict();
+  case[:kappa] = 0.0;
+  case[:E0] = max(0.05, rand(Edist));
+  Ks = rand(Kdist)
+  case[:K1], case[:K2] = Ks[1], Ks[2]
+  case[:n] = max(10, round(Int, rand(ndist)));
+  case[:b] = max(0.01, rand(bdist));
+  case[:kT] = 10.0^(rand([-1; 1])*rand(kTdist));
+  case[:Fx] = max(0.0, rand(Fdist));
+  case[:Fz] = max(0.0, rand(Fdist));
+
+  outfile = joinpath(workdir, "$(prefix(case)).out");
+  if !isfile(outfile)
+    println("Running case: $case.");
+    command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type interacting -b $(case[:b]) --bend-mod $(case[:kappa]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 2000000 --burn-in 200000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`;
+    output = read(command, String);
+    write(outfile, output); 
+  else
+    println("Case: $case has already been run.");
+  end
+
+end
diff --git a/run/phases_random.jl b/run/phases_random.jl
new file mode 100644
index 0000000..28f449f
--- /dev/null
+++ b/run/phases_random.jl
@@ -0,0 +1,50 @@
+println("Hello world!");
+
+@everywhere using Printf;
+@everywhere using Distributions;
+
+@show workdir = if length(ARGS) > 0
+  ARGS[1];
+else
+  ".";
+end
+
+@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x));
+@everywhere fmt_int(x) = @sprintf("%03d", x);
+
+@everywhere function prefix(case)
+  "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))_kappa-$(fmt(case[:kappa]))";
+end
+
+cases = Any[];
+κs = [0.0];
+Ks = zip([1.0; 0.0], [0.0; 1.0]);
+E0s = vcat(0.0, 0.05:0.05:1.0, 1.5:0.5:5.0, 10.0);
+ns = Int[100; 200];
+bs = [0.5; 1.0];
+kTs = [0.1; 1.0];
+Fs = zip([0.0; 1.0; 0.0; 1.0/sqrt(2.0)], [0.0; 0.0; 1.0; 1.0/sqrt(2.0)]);
+for b in bs, κ in κs, n in ns, Kvec in Ks, E0 in E0s, kT in kTs, F in Fs
+  K1, K2 = Kvec;
+  push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2,
+                    :kT => kT, :Fz => F[2], :kappa => κ,
+                    :Fx => F[1], :n => n, :b => b));
+end
+
+@info "total number of cases to run: $(length(cases))";
+
+mkpath(workdir);
+
+pmap(case -> begin;
+
+  outfile = joinpath(workdir, "$(prefix(case)).out");
+  if !isfile(outfile)
+    println("Running case: $case.");
+    command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type interacting -b $(case[:b]) --bend-mod $(case[:kappa]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 2000000 --burn-in 200000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`;
+    output = read(command, String);
+    write(outfile, output); 
+  else
+    println("Case: $case has already been run.");
+  end
+
+end, cases);