added burn in schedule capabilities to 3D version

2D/inc/dielectric_dipole.jl

@@ -1,4 +1,4 @@
 function μ(E0::Real, K1::Real, K2::Real, cϕ::Real, sϕ::Real)
   n̂i = n̂(cϕ, sϕ);
-  return ((K1 - K2) * E0 * cϕ * n̂i) + (K2 * [E0; 0.0]);
+  return ((K1 - K2) * E0 * sϕ * n̂i) + (K2 * [0.0; E0]);
 end

2D/inc/dipole_response.jl

@@ -6,7 +6,7 @@ abstract type DipoleResponse end
 
 function μ_dielectric(E0::Real, K1::Real, K2::Real, cϕ::Real, sϕ::Real)
   n̂i = n̂(cϕ, sϕ);
-  return ((K1 - K2) * E0 * cϕ * n̂i) + (K2 * [E0; 0.0]);
+  return ((K1 - K2) * E0 * sϕ * n̂i) + (K2 * [0.0; E0]);
 end
 
 struct DielectricResponse <: DipoleResponse

2D/inc/eap_chain.jl

@@ -61,7 +61,7 @@ function update_xs!(chain::EAPChain)
 end
 =#
 
-@inline u(E0::Real, μ::FdVector) = -1/2*E0*μ[1];
+@inline u(E0::Real, μ::FdVector) = -1/2*E0*μ[2];
 
 function EAPChain(pargs::Dict)
   ϕs = rand(ϕ_dist, pargs["num-monomers"]);
@@ -193,8 +193,7 @@ pflip_linear(ip::Real) = (1 + ip) / 2;
 
 function cluster_flip!(chain::EAPChain, idx::Int; 
                        pflip::Function = pflip_linear,
-                       ϵflip::Real = 1.0,
-                       ηerr::Real = 1e-10
+                       ϵflip::Real = 1.0
                       )
 
   # grow the cluster to the right
@@ -237,17 +236,6 @@ function cluster_flip!(chain::EAPChain, idx::Int;
       flip_n!(chain, i);
     end
 
-    if norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) > ηerr
-      println("prev");
-      display(prev_n̂s)
-      println();
-      println("curr");
-      display(chain.n̂s[:, lower_idx:upper_idx])
-      println();
-      @show(norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf))
-    end
-    @assert(norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) < ηerr);
-
     new_upper_p = if upper_idx < n(chain)
       pflip_linear(dot(n̂j(chain, upper_idx), n̂j(chain, upper_idx+1)));
     else

2D/run/Ising_2024-09-11.jl

@@ -0,0 +1,55 @@
+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]))_run-$(case[:run])";
+end
+
+cases = Any[];
+E0s = [1e-1;];
+K1s = [0.0]; 
+#K2s = [0.5; 1.0; 2.0]; 
+K2s = 1e-2*[1.0; 4.0]; 
+kTs = [1.0];
+Fxs = [0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 1.0; 2.0; 3.0; 4.0; 5.0; 7.5; 10.0; 12.5; 15.0; 20.0; 25.0; 30.0; 
+       35.0; 40.0; 45.0; 50.0; 60.0; 70.0; 80.0; 90.0; 100.0; 125.0; 150.0; 
+       175.0; 200.0; 250.0; 300.0; 400.0; 500.0];
+Fzs = [0.0];
+#ns = Int[100; 200; 500];
+ns = Int[25];
+#bs = [0.5; 1.0; 2.0];
+bs = [1.0; 4.0];
+runs = 1:10
+for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, K1 in K1s, K2 in K2s, run in runs
+  push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2,
+                    :kT => kT, :Fz => Fz,
+                    :Fx => Fx, :n => n, :b => b, :run => run));
+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 = `/home/grasinmj/julia-1.10.4/bin/julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 10000000 --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);

2D/run/Ising_2024-10-17.jl

@@ -0,0 +1,54 @@
+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]))_run-$(case[:run])";
+end
+
+cases = Any[];
+E0s = [1e-1;];
+K1s = [0.0]; 
+#K2s = [0.5; 1.0; 2.0]; 
+K2s = 1e-1*[0.1; 0.4; 1.0; 4.0]; 
+kTs = [1.0];
+Fxs = [0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 1.0; 2.0; 3.0; 4.0; 5.0; 7.5; 10.0; 12.5; 15.0; 20.0; 25.0; 30.0; 
+       35.0; 40.0; 45.0; 50.0];
+Fzs = [0.0];
+#ns = Int[100; 200; 500];
+ns = Int[25];
+#bs = [0.5; 1.0; 2.0];
+bs = [1.0];
+runs = 1:10
+for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, K1 in K1s, K2 in K2s, run in runs
+  push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2,
+                    :kT => kT, :Fz => Fz,
+                    :Fx => Fx, :n => n, :b => b, :run => run));
+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 = `/home/grasinmj/julia-1.10.4/bin/julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 10000000 --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);

2D/run/Ising_2024-10-21.jl

@@ -0,0 +1,54 @@
+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]))_run-$(case[:run])";
+end
+
+cases = Any[];
+E0s = [1e-1;];
+K2s = [0.0]; 
+#K2s = [0.5; 1.0; 2.0]; 
+K1s = 1e-1*[0.1; 0.4; 1.0; 4.0]; 
+kTs = [1.0];
+Fxs = [0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 1.0; 2.0; 3.0; 4.0; 5.0; 7.5; 10.0; 12.5; 15.0; 20.0; 25.0; 30.0; 
+       35.0; 40.0; 45.0; 50.0];
+Fzs = [0.0];
+#ns = Int[100; 200; 500];
+ns = Int[25];
+#bs = [0.5; 1.0; 2.0];
+bs = [1.0];
+runs = 1:10
+for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, K1 in K1s, K2 in K2s, run in runs
+  push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2,
+                    :kT => kT, :Fz => Fz,
+                    :Fx => Fx, :n => n, :b => b, :run => run));
+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 = `/home/grasinmj/julia-1.10.4/bin/julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 10000000 --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);

2D/run/Ising_2024-11-06.jl

@@ -0,0 +1,54 @@
+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]))_run-$(case[:run])";
+end
+
+cases = Any[];
+E0s = [1e-1;];
+K2s = [0.0]; 
+#K2s = [0.5; 1.0; 2.0]; 
+K1s = 1e-1*[0.1; 0.4; 1.0; 4.0]; 
+kTs = [1.0];
+Fzs = [0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 1.0; 2.0; 3.0; 4.0; 5.0; 7.5; 10.0; 12.5; 15.0; 20.0; 25.0; 30.0; 
+       35.0; 40.0; 45.0; 50.0];
+Fxs = [0.0];
+#ns = Int[100; 200; 500];
+ns = Int[25];
+#bs = [0.5; 1.0; 2.0];
+bs = [1.0];
+runs = 1:10
+for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, K1 in K1s, K2 in K2s, run in runs
+  push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2,
+                    :kT => kT, :Fz => Fz,
+                    :Fx => Fx, :n => n, :b => b, :run => run));
+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 = `/home/grasinmj/julia-1.10.4/bin/julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 10000000 --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);

inc/eap_chain.jl

@@ -1,6 +1,5 @@
 using Distributions;
 using NLsolve;
-using ForwardDiff;
 
 include(joinpath(@__DIR__, "dipole_response.jl"));
 

mcmc_clustering_eap_chain.jl

@@ -135,6 +135,10 @@ s = ArgParseSettings();
     help = "steps for burn-in; i.e. steps before averaging"
     arg_type = Int
     default = 50000
+  "--burn-schedule"
+    help = "temperature schedule for burn-in"
+    arg_type = String
+    default = "[1000; 100; 10; 2; 1]"
   "--x0"
     help = "initial configuration"
     arg_type = String
@@ -357,7 +361,19 @@ end
   exit(1);
   =#
 else
-  burned_in_chain = mcmc(pargs["burn-in"], pargs)[1];
+  kT_multipliers = eval(Meta.parse(pargs["burn-schedule"]));
+  kT_base = pargs["kT"];
+  burnargs = copy(pargs);
+  burnargs["kT"] = kT_base * kT_multipliers[1];
+  burned_in_chain = mcmc(pargs["burn-in"], burnargs)[1];
+  if length(kT_multipliers) > 1
+      for kT_mult in kT_multipliers[2:end]
+          global burned_in_chain
+          local burnargs = copy(pargs);
+          burnargs["kT"] = kT_base * kT_mult;
+          burned_in_chain = mcmc(pargs["burn-in"], burnargs, burned_in_chain)[1];
+      end
+  end
   mcmc(pargs["num-steps"], pargs, burned_in_chain);
 end
 

run/Ising_2024-11-12.jl

@@ -0,0 +1,54 @@
+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]))_run-$(case[:run])";
+end
+
+cases = Any[];
+κs = [0.0];
+#Ks = zip(1e-2*[0.0; 1.0; 0.0; 10.0], 1e-2*[1.0; 0.0; 10.0; 0.0]);
+Ks = zip(1e-1*[0.0; 1.0], 1e-1*[1.0; 0.0]);
+#E0s = vcat(0.0, 0.005:0.005:1.0);
+E0s = [0.1; 1.0; 10.0];
+ns = Int[25];
+bs = [1.0];
+kTs = [0.1; 1.0];
+Fzs = [0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 1.0; 2.0; 3.0; 4.0; 5.0; 7.5; 10.0; 12.5; 15.0; 20.0; 25.0; 30.0; 
+       35.0; 40.0; 45.0; 50.0];
+Fxs = [0.0];
+runs = 1:10
+for b in bs, κ in κs, n in ns, Kvec in Ks, E0 in E0s, kT in kTs, Fx in Fxs, Fz in Fzs, run in runs
+  K1, K2 = Kvec;
+  push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2,
+                    :kT => kT, :Fz => Fz, :kappa => κ,
+                    :Fx => Fx, :n => n, :b => b, :run => run));
+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 -t 1 -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -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 2500000 --burn-in 100000 -v 2 --prefix $(joinpath(workdir, prefix(case))) --stepout 250`;
+    output = read(command, String);
+    write(outfile, output); 
+  else
+    println("Case: $case has already been run.");
+  end
+
+end, cases);

scripts/aggregate_by.jl

@@ -39,7 +39,7 @@ for datafile in datafiles
   @show pattern = (runflag) ? pattern[1:end-5] * "*" * pattern[end-3:end] : pattern
   @show value = datafile[value_start_index:value_end_index-1];
   @show outfile = joinpath(outdir, join(filtered_params, "_")*".csv");
-  @show command = `/home/grasingerm/julia-1.6.7/bin/julia scripts/aggregate_mcmc.jl $outfile $indir "$pattern" $chain_type $dims $kappaflag $runflag`;
+  @show command = `julia scripts/aggregate_mcmc.jl $outfile $indir "$pattern" $chain_type $dims $kappaflag $runflag`;
   @show process_output = readlines(command);
 end
 

scripts/install_dependencies.jl

@@ -11,6 +11,7 @@ import Pkg;
                   "Optim";
                   "NLopt";
                   "Cubature";
+                  "Glob"
                  ]);
 
 foreach(Pkg.add, depends);