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);