diff options
Diffstat (limited to 'hw5/9-11.jl')
| -rw-r--r-- | hw5/9-11.jl | 713 |
1 files changed, 713 insertions, 0 deletions
diff --git a/hw5/9-11.jl b/hw5/9-11.jl new file mode 100644 index 0000000..802d9c6 --- /dev/null +++ b/hw5/9-11.jl @@ -0,0 +1,713 @@ +# molecular dynamics 2d. + +# usage: +# at the end of this script, under the header "DEMOS", +# you'll see some functions which implement demos from GN chapter 9. +# simply load the script in your development environment +# (I strongly recommend not using jupiter) +# and in the console/REPL run +# demo_0() +# etc. + +# demos 0,1,3 can optionally make an animated gif +# if you call it with the optional argument demo_3(gif=1) + +# lmk if this script is giving you grief or if you find any bugs +# kian@brown.edu + +using Statistics +using StatsPlots +using Plots.PlotMeasures + +mutable struct ParticleSystem + N::Int64 # number of particles + L::Float64 # square box side length + T₀::Float64 # initial temperature + t::Float64 # system time + dt::Float64 # time step + state::Vector{Float64} # state space array + steps::Int64 # number of steps + sampleInterval::Int64 # interval for sampling data + timeData::Vector{Float64} # array of sampled time points + energyData::Vector{Float64} # array of sampled energy values + tempData::Vector{Float64} # array of sampled temperature values + tempAccumulator::Float64 # temperature accumulator + squareTempAccumulator::Float64 # T^2 accumulator + virialAccumulator::Float64 # virial accumulator + xData::Vector{Vector{Float64}} # array of sampled position data + vData::Vector{Vector{Float64}} # array of sampled velocity data + forceType::String # force +end + +function ParticleSystem(N::Int64=64, L::Float64=10.0, T₀::Float64=1.0) + t = 0.0 + dt = 0.001 + state = zeros(4N) # state space array, [x1,y1,x2,y2,...,vx1,vy1,...] + steps = 0 + timeData = Float64[] + energyData = Float64[] + sampleInterval = 100 + tempData = Float64[] + tempAccumulator = 0.0 + squareTempAccumulator = 0.0 + virialAccumulator = 0.0 + xData = Vector{Float64}[] + vData = Vector{Float64}[] + forceType = "lennardJones" + + return ParticleSystem( + N, + L, + T₀, + t, + dt, + state, + steps, + sampleInterval, + timeData, + energyData, + tempData, + tempAccumulator, + squareTempAccumulator, + virialAccumulator, + xData, + vData, + forceType + ) +end + +# some useful "views" of the state array +# (read the performance tips chapter of the julia manual) +@views positions(state) = state[ 1:Int64(length(state)/2) ] +@views velocities(state) = state[ (Int64(length(state)/2)+1):end ] +@views xcomponent(vector) = vector[ 1:2:end ] +@views ycomponent(vector) = vector[ 2:2:end ] +@views particle(n, vector) = [ vector[2n-1], vector[2n] ] + +# INITIALIZATION +################################################################################ + +function set_random_positions!(sys::ParticleSystem) + println("\tposition configuration: random") + positions(sys.state) .= rand(2*sys.N) .* sys.L + cool!(sys) +end + +function set_square_lattice_positions!(sys::ParticleSystem) + println("\tposition configuration: square lattice") + + n = Int64(floor(sqrt(sys.N))) # num lattice points per axis + latticeSpacing = sys.L / n + + if sys.N != n^2 + println("\t\toops... your chosen N=$(sys.N) is not a square number") + println("\t\t-> resetting N to $(n^2).") + sys.N = n^2 + sys.state = zeros(4 * sys.N) + end + + for i in 0:(n-1) + for j in 0:(n-1) + sys.state[2*(i*n+j)+1] = (i + 0.5) * latticeSpacing + sys.state[2*(i*n+j)+2] = (j + 0.5) * latticeSpacing + end + end +end + +function set_triangular_lattice_positions!(sys::ParticleSystem) +end + +function add_position_jitter!(sys::ParticleSystem, jitter::Float64=0.5) + println("\tadding a wee bit of random jitter to particle positions...") + + for i = 1:length(positions(sys.state)) + sys.state[i] += rand() - jitter + end +end + +function set_random_velocities!(sys::ParticleSystem) + println("\tvelocity configuration: random") + + velocities(sys.state) .= rand(2*sys.N) .- 0.5 .+ velocities(sys.state) + zero_total_momentum!(sys) + velocities(sys.state) .*= sqrt(sys.T₀/temperature(sys)) +end + +function scale_velocities!(sys, scale) + velocities(sys.state) .*= scale +end + + +function zero_total_momentum!(sys::ParticleSystem) + xcomponent(velocities(sys.state)) .-= + mean(xcomponent(velocities(sys.state))) + ycomponent(velocities(sys.state)) .-= + mean(ycomponent(velocities(sys.state))) +end + + +# FORCES / POTENTIALS +################################################################################ + +function force(sys::ParticleSystem) + if sys.forceType == "lennardJones" + force, virial = lennard_jones_force(sys) + elseif sys.forceType == "powerLaw" + force, virial = power_law_force(sys) + end + + sys.virialAccumulator += virial + + return force +end + +# the minimum image approximation +# (periodic boundary conditions) +function minimum_image(xij::Float64, L::Float64) + if xij > (L/2) + xij -= L + elseif xij < -(L/2) + xij += L + end + return xij +end + +function lennard_jones_force(sys::ParticleSystem) + x = xcomponent(positions(sys.state)) + y = ycomponent(positions(sys.state)) + virial = 0.0 + force = zeros(2*sys.N) + + Threads.@threads for i = 1:(sys.N-1) + for j = (i+1):sys.N + dx = minimum_image(x[i] - x[j], sys.L) + dy = minimum_image(y[i] - y[j], sys.L) + + r2inv = 1.0 / (dx^2 + dy^2) + f = 48.0 * r2inv^7 - 24.0 * r2inv^4 + fx = dx * f + fy = dy * f + + force[2*i-1] += fx + force[2*i] += fy + force[2*j-1] -= fx + force[2*j] -= fy + + virial += fx * dx + fy * dy + end + end + + return force, 0.5 * virial +end + +function lennard_jones_potential(sys::ParticleSystem) + x = xcomponent(positions(sys.state)) + y = ycomponent(positions(sys.state)) + U = 0.0 + + Threads.@threads for i in 1:(sys.N-1) + for j in (i+1):sys.N + dx = minimum_image(x[i] - x[j], sys.L) + dy = minimum_image(y[i] - y[j], sys.L) + + r2inv = 1.0 / (dx^2 + dy^2) + U += r2inv^6 - r2inv^3 + end + end + return 4.0 * U +end + +function power_law_force(sys::ParticleSystem) +end + +function power_law_potential(sys::ParticleSystem) +end + +# TIME EVOLUTION +################################################################################ + +function keep_particles_in_box!(sys::ParticleSystem) + for i in 1:(2*sys.N) + if positions(sys.state)[i] > sys.L + positions(sys.state)[i] -= sys.L + elseif positions(sys.state)[i] < 0.0 + positions(sys.state)[i] += sys.L + end + end + +# # another way of doing this: use the ternary operator +# for i in 1:(2 * sys.N) +# positions(sys.state)[i] < 0.0 ? +# positions(sys.state)[i] % sys.L + sys.L : +# positions(sys.state)[i] % sys.L +# end +end + +function verlet_step!(sys::ParticleSystem) + # compute acceleration at current time + acceleration = force(sys) + + # compute positions at t + dt + positions(sys.state) .+= + velocities(sys.state) .* sys.dt .+ + 0.5 .* acceleration .* (sys.dt)^2 + + # enforce boundary conditions + # (basically check if any particles left the box and put them back) + # see function implementation for deets + keep_particles_in_box!(sys) + + # compute velocities at t + dt + velocities(sys.state) .+= + 0.5 * sys.dt .* (acceleration + force(sys)) +end + +function update_r_squared!(sys, r_squared) + # take the positions of the firs ttwo particles + p1 = particle(1, sys.state) + p2 = particle(2, sys.state) + + # calculate the distance between them + # compute the image difference + dx = minimum_image(p1[1] - p2[1], sys.L) + dy = minimum_image(p1[2] - p2[2], sys.L) + r = dx^2 + dy^2 + + # println("r_squared: ", r) + # println("t: ", sys.t) + + push!(r_squared, [sys.t, r]) +end + +function evolve!(sys::ParticleSystem, runtime::Float64=10.0, r_squared=[]) + numsteps = Int64(abs(runtime/sys.dt) + 1) + + print_evolution_message(runtime, numsteps) + + @time for step in 1:numsteps + # for problem 9.11 + update_r_squared!(sys, r_squared) + + verlet_step!(sys) + zero_total_momentum!(sys) + + if (step % sys.sampleInterval == 1) + push!(sys.timeData, sys.t) + push!(sys.energyData, energy(sys)) + push!(sys.xData, positions(sys.state)) + push!(sys.vData, velocities(sys.state)) + + T = temperature(sys) + push!(sys.tempData, T) + sys.tempAccumulator += T + sys.squareTempAccumulator += T^2 + end + + sys.t += sys.dt + sys.steps += 1 + end + println("done.") +end + +function reverse_time!(sys) + sys.dt *= -1 + println("\ntime reversed! dt = $(sys.dt)") +end + +function cool!(sys::ParticleSystem, cooltime::Float64=1.0) + numsteps = Int64(cooltime/sys.dt) + for step in 1:numsteps + verlet_step!(sys) + velocities(sys.state) .*= (1.0 - sys.dt) + end + reset_statistics!(sys) +end + +# MEASUREMENTS +################################################################################ + +function kinetic_energy(sys::ParticleSystem) + return 0.5 * sum(velocities(sys.state) .* velocities(sys.state)) +end + +function potential_energy(sys::ParticleSystem) + return lennard_jones_potential(sys) +end + +function temperature(sys::ParticleSystem) + return kinetic_energy(sys) / sys.N +end + +function energy(sys::ParticleSystem) + return potential_energy(sys) + kinetic_energy(sys) +end + +# STATISTICS +################################################################################ + +function reset_statistics!(sys::ParticleSystem) + sys.steps = 0 + sys.tempAccumulator = 0.0 + sys.squareTempAccumulator = 0.0 + sys.virialAccumulator = 0.0 + sys.xData = [] + sys.vData = [] +end + +function mean_temperature(sys::ParticleSystem) + return sys.tempAccumulator / sys.steps +end + +function mean_square_temperature(sys::ParticleSystem) + return sys.squareTempAccumulator / sys.steps +end + +function mean_pressure(sys::ParticleSystem) + # factor of half because force is calculated twice each step + meanVirial = 0.5 * sys.virialAccumulator / sys.steps + return 1.0 + 0.5 * meanVirial / (sys.N * mean_temperature(sys)) +end + +function heat_capacity(sys::ParticleSystem) + meanTemperature = mean_temperature(sys) + meanSquareTemperature = mean_square_temperature(sys) + σ2 = meanSquareTemperature - meanTemperature^2 + denom = 1.0 - σ2 * sys.N / meanTemperature^2 + return sys.N / denom +end + +function mean_energy(sys::ParticleSystem) + return mean(sys.energyData) +end + +function std_energy(sys::ParticleSystem) + return std(sys.energyData) +end + +# MATH / ADDITIONAL FUNCTIONS +################################################################################ + +function dot(v1::Vector{Float64}, v2::Vector{Float64}) + return sum(v1 .* v2) +end + +# GRAPHS +################################################################################ + +function initialize_plot() + plot( + size=(800,800), + titlefontsize=12, + guidefontsize=12, + ) +end + +function plot_positions_t(sys::ParticleSystem, t::Int64) + initialize_plot() + for n = 1:sys.N + scatter!( + [ sys.xData[t][2n-1] ], + [ sys.xData[t][2n] ], + markersize = 4.0, + markercolor = n, + markerstrokewidth = 0.4, + grid = true, + framestyle = :box, + legend = false, + ) + end +end + +function animate(sys::ParticleSystem, interval::Int64=1) + println("\ngenerating gif...") + + scatter!() + animation = @animate for t in 1:length(sys.xData) + scatter() + for n = 1:sys.N + scatter!( + [ sys.xData[t][2n-1] ], + [ sys.xData[t][2n] ], + #markersize = 4.0, + markercolor = n, + #markerstrokewidth = 0.4, + grid = true, + framestyle = :box, + legend = false, + ) + end + xlims!(0, sys.L) + ylims!(0, sys.L) + xlabel!("x") + ylabel!("y") + end every interval + + gif(animation, "./animation.gif") + println("done.") +end + +function plot_positions(sys::ParticleSystem) + initialize_plot() + for n = 1:sys.N + scatter!( + [ xcomponent(positions(sys.state))[n] ], + [ ycomponent(positions(sys.state))[n] ], + markersize = 4.0, + markercolor = n, + markerstrokewidth = 0.4, + grid = true, + framestyle = :box, + legend = false, + ) + end + xlims!(0, sys.L) + ylims!(0, sys.L) + xlabel!("x") + ylabel!("y") + title!("positions at t=$(round(sys.t, digits=4))") +end + +function plot_trajectories(sys::ParticleSystem, particles::Vector{Int64}=[ 1 ]) + initialize_plot() + for n = 1:sys.N + scatter!( + [ xcomponent(positions(sys.state))[n] ], + [ ycomponent(positions(sys.state))[n] ], + markersize = 4.0, + markercolor = n, + markerstrokewidth = 0.4, + grid = true, + framestyle = :box, + legend = false, + ) + end + + for n in collect(particles) + xdata = [ sys.xData[i][2n-1] for i in 1:length(sys.xData) ] + ydata = [ sys.xData[i][2n] for i in 1:length(sys.xData) ] + + # plot trajectory line for nth particle + scatter!( + xdata, + ydata, + color = n, + #markerstrokewidth = 0, + markerstrokecolor = n, + markersize = 0.7, + markeralpha = 0.5, + label = false, + widen = false, + ) + + # plot initial position for nth particle + scatter!( + [ sys.xData[1][2n-1] ], + [ sys.xData[1][2n] ], + markersize = 4.0, + markercolor = n, + markerstrokewidth = 0.4, + markeralpha = 0.3, + #label = "pcl. $n @t=t₀", + widen = false, + ) + + # plot final position for nth particle + scatter!( + [ sys.xData[end][2n-1] ], + [ sys.xData[end][2n] ], + markersize = 4.0, + markercolor = n, + markerstrokewidth = 0.4, + markeralpha = 1.0, + #label = "pcl $n @t=t", + widen = false, + ) + end + title!("positions & trajectories at time t=$(round(sys.t, digits=2))") + plot!() +end + +function plot_temperature(sys::ParticleSystem) + initialize_plot() + plot!( + sys.timeData, + sys.tempData, + #widen = true, + legend = false, + ) + ylims!( + mean(sys.tempData) - std(sys.tempData) * 1.5, + mean(sys.tempData) + std(sys.tempData) * 3, + ) + xlabel!("t") + ylabel!("T(t)") + title!("temperature vs time") +end + +function plot_energy(sys::ParticleSystem, ylimit::Float64=1.0) + initialize_plot() + plot!( + sys.timeData, + sys.energyData, + #widen = true, + legend = false, + ) + ylims!( + #ylimit * (mean(sys.energyData) - 1), + #ylimit * (mean(sys.energyData) + 1) + mean(sys.energyData) - std(sys.energyData) * 10, + mean(sys.energyData) + std(sys.energyData) * 10, + ) + xlabel!("t") + ylabel!("E(t)") + title!("energy vs time") +end + +function plot_speed_distribution(sys::ParticleSystem, numSamples::Int64=5) + initialize_plot() + + numDataPoints = Int64(length(sys.vData)) + interval = Int64(floor(numDataPoints / numSamples)) + + samples = collect(1:interval:numDataPoints) + for s in samples + speed = sqrt.( + xcomponent(sys.vData[s]).^2 .* + ycomponent(sys.vData[s]).^2 + ) + density!( + sys.vData[s], + normalize = :pdf, + label = "t = $(round(sys.timeData[s], digits=2))", + ) + end + xlabel!("speed") + title!("speed distribution") +end + +# CONSOLE PRINT DATA +################################################################################ + +function print_hello() + println("\nmolecular dynamics!") + println("number of threads: ", Threads.nthreads()) +end + +function print_bonjour() + println("\nbonjour") +end + +function print_system_parameters(sys::ParticleSystem) + println("\nsystem parameters:") + println("\tN = $(sys.N) (number of particles)") + println("\tL = $(sys.L) (side length of square box)") + println("\tDT = $(sys.dt) (time step)") +end + +function print_system_data(sys::ParticleSystem) + println("\nsystem data at time t=$(round(sys.t, digits=4))") + + if sys.steps == 0 + println("\ttemperature: $(temperature(sys))") + println("\tenergy: $(energy(sys))") + else + println("\tsteps evolved: $(sys.steps)") + println("\ttemperature: $(temperature(sys))") + println("\tenergy: $(energy(sys))") + println("\tmean energy: $(mean_energy(sys))") + println("\tstd energy: $(std_energy(sys))") + println("\theat capacity: $(heat_capacity(sys))") + println("\tPV/NkT: $(mean_pressure(sys))") + end +end + +function print_evolution_message(runtime, numsteps) + println("\nevolving...") +end + +# DEMOS +################################################################################ + +# DEMO 3: MELTING TRANSITION +function problem_9_11!(sys) + println("\nProblem 9.11") + println("----------------------------------------") + + # initialize system of particles on square lattice with zero velocity + set_square_lattice_positions!(sys) + # set_random_velocities!(sys) + print_system_data(sys) + + # evolve the system and watch them "crystallize" + # into a triangular lattice formation + r_squared = [] + evolve!(sys, 20.0, r_squared) + print_system_data(sys) + + # now, increase the temperature of the system by giving the particles + # some velocity. evolve the system and plot the trajectories. + scale_velocities!(sys, 1.5) + evolve!(sys, 10.0, r_squared) + print_system_data(sys) + + scale_velocities!(sys, 1.5) + evolve!(sys, 10.0, r_squared) + print_system_data(sys) + + scale_velocities!(sys, 1.5) + evolve!(sys, 10.0, r_squared) + print_system_data(sys) + + scale_velocities!(sys, 1.5) + evolve!(sys, 10.0, r_squared) + print_system_data(sys) + + # some more plots + p4 = plot_energy(sys, 0.0) + p5 = plot_temperature(sys) + + t = [] + r = [] + for i in 1:length(r_squared) + push!(t, r_squared[i][1]) + push!(r, r_squared[i][2]) + end + + # plot r_squared + p_r = plot(t, r, xlabel="t", ylabel="r^2", label="r^2(t)", title="Pair Separation for $(sys.N) particles, box length 4", legend=false, + top_margin=5mm, bottom_margin=5mm, left_margin=5mm) + + return p4, p5, p_r +end + +sys_9 = ParticleSystem(9, 4.0, 0.0) +p_e_9, p_t_9, p_r_9 = problem_9_11!(sys_9) + +sys_16 = ParticleSystem(16, 4.0, 0.0) +p_e_16, p_t_16, p_r_16 = problem_9_11!(sys_16) + +sys_25 = ParticleSystem(25, 4.0, 0.0) +p_e_25, p_t_25, p_r_25 = problem_9_11!(sys_25) + +sys_36 = ParticleSystem(36, 4.0, 0.0) +p_e_36, p_t_36, p_r_36 = problem_9_11!(sys_36) + +# plot( +# p_e_100, p_t_100, p_r_100, +# p_e_16, p_t_16, p_r_16, +# p_e_36, p_t_36, p_r_36, +# layout=(3,3), size=(1000,1000)) +# savefig("problem_9_11.png") + +# plot(p_r_16, p_t_16, layout=(1,2), size=(1000,500)) + +plot( + p_r_9, p_t_9, + p_r_16, p_t_16, + p_r_25, p_t_25, + p_r_36, p_t_36, + layout=(4,2), size=(1250, 1600) +) + +# only plot the 100 particle system +savefig("problem_9_11_t.png")
\ No newline at end of file |