From b2d0a79356136f4bece58097dffc035ed2e8bfc2 Mon Sep 17 00:00:00 2001
From: Matthew Grasinger <grasingerm@gmail.com>
Date: Mon, 15 May 2023 13:40:35 -0400
Subject: [PATCH 1/4] made a new phases script with less monomers per chain

---
 run/phases_2023-05-15.jl | 49 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 49 insertions(+)
 create mode 100644 run/phases_2023-05-15.jl

diff --git a/run/phases_2023-05-15.jl b/run/phases_2023-05-15.jl
new file mode 100644
index 0000000..10a6aa1
--- /dev/null
+++ b/run/phases_2023-05-15.jl
@@ -0,0 +1,49 @@
+println("Hello world!");
+
+@everywhere using Printf;
+
+@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]; [2.0; 1.0]; [1.0; 2.0]);
+E0s = vcat(0.05:0.05:1.0, 1.5:0.5:5.0, 10.0);
+ns = Int[15; 25];
+bs = [0.25; 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);

From 09ca9354ec9c71368227ce99c7c9ca7d5a25c358 Mon Sep 17 00:00:00 2001
From: Matthew Grasinger <grasingerm@gmail.com>
Date: Mon, 15 May 2023 15:23:45 -0400
Subject: [PATCH 2/4] script for randomly generated samples

---
 run/phases-random_2023-05-15.jl | 51 +++++++++++++++++++++++++++++++++
 1 file changed, 51 insertions(+)
 create mode 100644 run/phases-random_2023-05-15.jl

diff --git a/run/phases-random_2023-05-15.jl b/run/phases-random_2023-05-15.jl
new file mode 100644
index 0000000..e2539db
--- /dev/null
+++ b/run/phases-random_2023-05-15.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(35, 15);
+@everywhere bdist = [1e-3; 1e-2; 1e-1; 0.5; 1.0; 2.0; 1e2]; 
+@everywhere kTdist = [0.1; 1.0; 10.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] = rand(bdist);
+  case[:kT] = 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

From 7ad06e4baaa99169c62888b81213c9a345789f49 Mon Sep 17 00:00:00 2001
From: Matthew Grasinger <grasingerm@gmail.com>
Date: Thu, 18 May 2023 11:06:52 -0400
Subject: [PATCH 3/4] added functionality for interacting energies with a
 cutoff radius

---
 inc/eap_chain.jl             | 31 +++++++++++++++++++++++
 mcmc_clustering_eap_chain.jl |  6 ++++-
 run/phases-big_2023-05-18.jl | 49 ++++++++++++++++++++++++++++++++++++
 3 files changed, 85 insertions(+), 1 deletion(-)
 create mode 100644 run/phases-big_2023-05-18.jl

diff --git a/inc/eap_chain.jl b/inc/eap_chain.jl
index ba25799..5024520 100644
--- a/inc/eap_chain.jl
+++ b/inc/eap_chain.jl
@@ -82,6 +82,8 @@ function EAPChain(pargs::Dict)
                     InteractingEnergy();
                   elseif pargs["energy-type"] == "Ising"
                     IsingEnergy();
+                  elseif pargs["energy-type"] == "cutoff"
+                    UCutoff(pargs["cutoff-radius"]*pargs["mlen"])
                   else
                     error("energy-type is not understood.");
                   end,
@@ -144,6 +146,35 @@ function EAPChain(chain::EAPChain)
                  );
 end
 
+struct UCutoff <: Energy
+    cutoff_radius::Real;
+end
+
+# TODO: consider rewriting this to make better use of past calculations
+# this is probably by far the slowest calculation *shrug emoji*
+function (U_func::UCutoff)(chain::EAPChain)
+  crad2 = U_func.cutoff_radius * U_func.cutoff_radius;
+  U = 0.0;
+  @inbounds for i=1:n(chain)
+    @simd for j=i+1:n(chain) # once debugged, use inbounds
+      r = chain.xs[:, i] - chain.xs[:, j];
+      r2 = dot(r, r);
+      dU = if (r2 > crad2)
+          0.0
+      else
+          rmag = sqrt(r2);
+          r̂ = r / rmag;
+          r3 = r2*rmag;
+          μi = view(chain.μs, :, i);
+          μj = view(chain.μs, :, j);
+          (dot(μi, μj) - 3*dot(μi, r̂)*dot(μj, r̂)) / (4*π*r3);
+      end
+      U += dU;
+    end
+  end
+  return U;
+end
+
 # TODO: consider rewriting this to make better use of past calculations
 # this is probably by far the slowest calculation *shrug emoji*
 function U_interaction(chain::EAPChain)
diff --git a/mcmc_clustering_eap_chain.jl b/mcmc_clustering_eap_chain.jl
index f2b0e83..212088c 100644
--- a/mcmc_clustering_eap_chain.jl
+++ b/mcmc_clustering_eap_chain.jl
@@ -42,9 +42,13 @@ s = ArgParseSettings();
     default = 0.0
     help = "zero energy bond angle"
   "--energy-type", "-u"
-    help = "energy type (noninteracting|interacting|Ising)"
+    help = "energy type (interacting|cutoff|Ising|noninteracting)"
     arg_type = String
     default = "Ising"
+  "--cutoff-radius"
+    help = "cut off radius (units of monomer lengths)"
+    arg_type = Float64
+    default = 7.5
   "--kT", "-k"
     help = "dimensionless temperature"
     arg_type = Float64
diff --git a/run/phases-big_2023-05-18.jl b/run/phases-big_2023-05-18.jl
new file mode 100644
index 0000000..c7680b2
--- /dev/null
+++ b/run/phases-big_2023-05-18.jl
@@ -0,0 +1,49 @@
+println("Hello world!");
+
+@everywhere using Printf;
+
+@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.05:0.05:1.0, 1.5:0.5:5.0, 10.0);
+ns = Int[400];
+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 cutoff --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);

From 3f263e9f977011b06de20ec0b11b18122a2cf659 Mon Sep 17 00:00:00 2001
From: Matthew Grasinger <grasingerm@gmail.com>
Date: Thu, 18 May 2023 11:11:23 -0400
Subject: [PATCH 4/4] typo in study for long chains

---
 run/phases-big_2023-05-18.jl |  2 +-
 run/phases_2023-05-15b.jl    | 49 ++++++++++++++++++++++++++++++++++++
 2 files changed, 50 insertions(+), 1 deletion(-)
 create mode 100644 run/phases_2023-05-15b.jl

diff --git a/run/phases-big_2023-05-18.jl b/run/phases-big_2023-05-18.jl
index c7680b2..407eb2b 100644
--- a/run/phases-big_2023-05-18.jl
+++ b/run/phases-big_2023-05-18.jl
@@ -39,7 +39,7 @@ 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 cutoff --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)))`;
+    command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type cutoff -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
diff --git a/run/phases_2023-05-15b.jl b/run/phases_2023-05-15b.jl
new file mode 100644
index 0000000..e3fb9d6
--- /dev/null
+++ b/run/phases_2023-05-15b.jl
@@ -0,0 +1,49 @@
+println("Hello world!");
+
+@everywhere using Printf;
+
+@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], [2.0; 1.0], [1.0; 2.0]);
+E0s = vcat(0.05:0.05:1.0, 1.5:0.5:5.0, 10.0);
+ns = Int[40; 50];
+bs = [0.1; 0.25; 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);