From e650ed1e1e908e51c78c1b047bec0da7c4fea366 Mon Sep 17 00:00:00 2001 From: sotech117 Date: Sat, 27 Apr 2024 04:25:23 -0400 Subject: testing --- hw5/9-11.jl | 713 +++++++++++++++++++++++++++++++++++ hw5/9-14.jl | 846 +++++++++++++++++++++++++++++++++++++++++ hw5/9-14.png | Bin 0 -> 59247 bytes hw5/9-2.jl | 848 ++++++++++++++++++++++++++++++++++++++++++ hw5/9-2.png | Bin 0 -> 903130 bytes hw5/HarmonicOscillator.jl | 71 ++++ hw5/MolecularDynamics2.jl | 755 +++++++++++++++++++++++++++++++++++++ hw5/demo_2-f.png | Bin 0 -> 1842809 bytes hw5/demo_2.png | Bin 0 -> 2321892 bytes hw5/fpu-2.jl | 138 +++++++ hw5/fpu-3.jl | 78 ++++ hw5/fpu.gif | Bin 0 -> 1551508 bytes hw5/fpu.jl | 35 ++ hw5/fpu.png | Bin 0 -> 44884 bytes hw5/lfg-modes.png | Bin 0 -> 40240 bytes hw5/modes-beta15.png | Bin 0 -> 57359 bytes hw5/modes-beta3.png | Bin 0 -> 59535 bytes hw5/plot_175.svg | 105 ++++++ hw5/plz-modes.png | Bin 0 -> 49160 bytes hw5/plz.gif | Bin 0 -> 455391 bytes hw5/plz.jl | 148 ++++++++ hw5/plz.png | Bin 0 -> 49078 bytes hw5/problem_9_-test.png | Bin 0 -> 196293 bytes hw5/problem_9_11.png | Bin 0 -> 139775 bytes hw5/problem_9_11_100.png | Bin 0 -> 15223 bytes hw5/problem_9_11_16.png | Bin 0 -> 84632 bytes hw5/problem_9_11_25.png | Bin 0 -> 74954 bytes hw5/problem_9_11_9.png | Bin 0 -> 47089 bytes hw5/problem_9_11_old.png | Bin 0 -> 259026 bytes hw5/problem_9_11_t.png | Bin 0 -> 229515 bytes hw5/problem_9_2-fail.png | Bin 0 -> 185918 bytes hw5/problem_9_2-w.png | Bin 0 -> 175635 bytes hw5/problem_9_2.png | Bin 0 -> 171722 bytes hw5/problem_9_2_1000.png | Bin 0 -> 194071 bytes hw5/problem_9_2_L.png | Bin 0 -> 195201 bytes hw5/problem_9_2_final.png | Bin 0 -> 139872 bytes hw5/problem_9_2_long-100p.png | Bin 0 -> 179419 bytes hw5/problem_9_2_long-255p.png | Bin 0 -> 169825 bytes hw5/problem_9_2_long.png | Bin 0 -> 177946 bytes hw5/r_squared-1.png | Bin 0 -> 122244 bytes 40 files changed, 3737 insertions(+) create mode 100644 hw5/9-11.jl create mode 100644 hw5/9-14.jl create mode 100644 hw5/9-14.png create mode 100644 hw5/9-2.jl create mode 100644 hw5/9-2.png create mode 100644 hw5/HarmonicOscillator.jl create mode 100644 hw5/MolecularDynamics2.jl create mode 100644 hw5/demo_2-f.png create mode 100644 hw5/demo_2.png create mode 100644 hw5/fpu-2.jl create mode 100644 hw5/fpu-3.jl create mode 100644 hw5/fpu.gif create mode 100644 hw5/fpu.jl create mode 100644 hw5/fpu.png create mode 100644 hw5/lfg-modes.png create mode 100644 hw5/modes-beta15.png create mode 100644 hw5/modes-beta3.png create mode 100644 hw5/plot_175.svg create mode 100644 hw5/plz-modes.png create mode 100644 hw5/plz.gif create mode 100644 hw5/plz.jl create mode 100644 hw5/plz.png create mode 100644 hw5/problem_9_-test.png create mode 100644 hw5/problem_9_11.png create mode 100644 hw5/problem_9_11_100.png create mode 100644 hw5/problem_9_11_16.png create mode 100644 hw5/problem_9_11_25.png create mode 100644 hw5/problem_9_11_9.png create mode 100644 hw5/problem_9_11_old.png create mode 100644 hw5/problem_9_11_t.png create mode 100644 hw5/problem_9_2-fail.png create mode 100644 hw5/problem_9_2-w.png create mode 100644 hw5/problem_9_2.png create mode 100644 hw5/problem_9_2_1000.png create mode 100644 hw5/problem_9_2_L.png create mode 100644 hw5/problem_9_2_final.png create mode 100644 hw5/problem_9_2_long-100p.png create mode 100644 hw5/problem_9_2_long-255p.png create mode 100644 hw5/problem_9_2_long.png create mode 100644 hw5/r_squared-1.png (limited to 'hw5') 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 diff --git a/hw5/9-14.jl b/hw5/9-14.jl new file mode 100644 index 0000000..7829e40 --- /dev/null +++ b/hw5/9-14.jl @@ -0,0 +1,846 @@ +# 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 + +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 + zero_total_momentum!(sys) + velocities(sys.state) .*= sqrt(sys.T₀/temperature(sys)) +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_histogram!(histogram_data, bins, t_interval, sys) + # for problem 9-14 + # if sys.t % t_interval != 0 + # return # only update histogram at t = n * t_interval + # end + + # take the first particle as the origin + origin = particle(1, sys.state) + # println("origin: ", origin) + + # for each particle, calculate the distance from the origin + # and update the histogram + for i in 2:sys.N + pcl = particle(i, sys.state) + dx = minimum_image(pcl[1] - origin[1], sys.L) + dy = minimum_image(pcl[2] - origin[2], sys.L) + r = sqrt(dx^2 + dy^2) + + # find the bin that r belongs to + for j in 1:length(bins)-1 + if r >= bins[j] && r < bins[j+1] + histogram_data[j] += 1 + break + end + end + end +end + +function evolve!(sys, runtime, histogram_data, bins, t_interval) + numsteps = Int64(abs(runtime/sys.dt) + 1) + + print_evolution_message(runtime, numsteps) + + @time for step in 1:numsteps + verlet_step!(sys) + zero_total_momentum!(sys) + + if (step % sys.sampleInterval == 1) + # for problem, 9-14 update historgram here + update_histogram!(histogram_data, bins, t_interval, sys) + + 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, + ) + ylims!( + mean(sys.tempData) - std(sys.tempData) * 20, + mean(sys.tempData) + std(sys.tempData) * 20, + ) + 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, + ) + 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 0: APPROACH TO EQUILIBRIUM +function demo_0(;gif=0) + println("\nDEMO 0: APPROACH TO EQUILIBRIUM") + println("----------------------------------------") + + sys = ParticleSystem(64, 120.0, 1.0) + print_system_parameters(sys) + + set_square_lattice_positions!(sys) + set_random_velocities!(sys) + print_system_data(sys) + p1 = plot_positions(sys) + + evolve!(sys, 20.0) + print_system_data(sys) + + p2 = plot_trajectories(sys, collect(1:64)) + p3 = plot_energy(sys) + p4 = plot_temperature(sys) + + # make gif + if gif == 1 + animate(sys, 1) + end + + plot( + p1, p2, p3, p4, + layout = grid(2,2, heights=[0.7,0.3]), + size = (1280,720) + ) +end + +# DEMO 1: TIME REVERSAL TEST +function demo_1(;gif=0) + println("\nDEMO 1: TIME REVERSAL TEST") + println("----------------------------------------") + + sys = ParticleSystem(64, 120.0, 1.0) + print_system_parameters(sys) + + set_square_lattice_positions!(sys) + set_random_velocities!(sys) + print_system_data(sys) + p1 = plot_positions(sys) + + evolve!(sys, 50.0) + #p2 = plot_trajectories(sys, collect(1:64)) + p2 = plot_positions(sys) + + reverse_time!(sys) + evolve!(sys, 50.0) + print_system_data(sys) + #p3 = plot_trajectories(sys, collect(1:64)) + p3 = plot_positions(sys) + + # make gif + if gif == 1 + animate(sys, 4) + end + + plot( + p1, p2, p3, + layout = (1,3), + size = (1200,400) + ) +end + +# DEMO 2: SPEED DISTRIBUTION +function demo_2() + println("\nDEMO 2: SPEED DISTRIBUTION") + println("----------------------------------------") + + sys = ParticleSystem[] + + # array for speed distribution plots + ps = Plots.Plot{Plots.GRBackend}[] + + # array for trajectory plots + pt = Plots.Plot{Plots.GRBackend}[] + + # initialize three systems with different initial conditions + # but same KE and PE, evolve, and save plots + for i = 1:3 + push!(sys, ParticleSystem(64, 120.0, 1.0)) + + println("\nSYSTEM $i") + print_system_parameters(sys[i]) + + set_square_lattice_positions!(sys[i]) + add_position_jitter!(sys[i]) + set_random_velocities!(sys[i]) + print_system_data(sys[i]) + + evolve!(sys[i], 40.0) + print_system_data(sys[i]) + push!(ps, plot_speed_distribution(sys[i], 5)) + push!(pt, plot_trajectories(sys[i], collect(1:64)) ) + end + + + # plot speed distribution and trajectory plots + plot( + ps[1], ps[2], ps[3], + pt[1], pt[2], pt[3], + layout = (2,3), + size = (1920,1080) + ) +end + +# DEMO 3: MELTING TRANSITION +function demo_3(;gif=0) + println("\nDEMO 3: MELTING TRANSITION") + println("----------------------------------------") + + # initialize system of particles on square lattice with zero velocity + sys = ParticleSystem(100, 10.0, 5.0) + set_square_lattice_positions!(sys) + print_system_data(sys) + p1 = plot_positions(sys) + + # evolve the system and watch them "crystallize" + # into a triangular lattice formation + evolve!(sys, 20.0) + print_system_data(sys) + p2 = plot_trajectories(sys, collect(1:100)) + + # now, increase the temperature of the system by giving the particles + # some velocity. evolve the system and plot the trajectories. + set_random_velocities!(sys) + evolve!(sys, 60.0) + print_system_data(sys) + p3 = plot_trajectories(sys, collect(1:100)) + + # some more plots + p4 = plot_energy(sys, 0.0) + p5 = plot_temperature(sys) + p6 = plot_speed_distribution(sys, 20) + + # make gif + if gif == 1 + animate(sys, 1) + end + + plot( + p1, p2, p3, p4, p5, p6, + layout = (2,3), + size = (1280,720) + ) +end + +function normalize_histogram(histogram) + println("normalizing histogram...") + println(histogram) + return histogram ./ sum(histogram) +end + +# problem 9-14 +function problem_9_14(temp, bins, t_interval) + println("\nPROBLEM 9-14: RADIAL DISTRIBUTION FUNCTION") + println("----------------------------------------") + + # initialize system of particles on square lattice with zero velocity + sys = ParticleSystem(16, 4.0, temp) + set_square_lattice_positions!(sys) + set_random_velocities!(sys) + print_system_data(sys) + p1 = plot_positions(sys) + + # evolve the system and watch them "crystallize" + # into a triangular lattice formation + histogram_data = zeros(length(bins) - 1) + evolve!(sys, 500, histogram_data, bins, t_interval) + print_system_data(sys) + + # normalize the histogram + histogram_data = normalize_histogram(histogram_data) + # get the max y value for the plot + y_max = maximum(histogram_data) + + # center the bins + centered_bins = bins .+ (bins[2] - bins[1]) / 2 + + # plot the histogram + return plot( + centered_bins[2:end], histogram_data, + marker = :circle, + xlabel = "r", + ylabel = "g(r)", + title = "Pair Correlation (T=$temp)", + label="pdf(r)" + ), y_max +end + +bins = range(.5, stop=3, length=40) +p_05, y_max_1 = problem_9_14(.5, bins, 10.0) +# ont this plot, plot vertical lines at 1.1, 1.9, 2.8 +y_max_1 += 0.025 +plot!(p_05, [1.1, 1.1], [0, y_max_1], line = :dash, label = "r = 1.1", color = :green, ylim=(0, y_max_1), lw=1.5, legend=:topright) +plot!(p_05, [1.9, 1.9], [0, y_max_1], line = :dash, label = "r = 1.9", color = :orange, lw=1.5) +plot!(p_05, [2.8, 2.8], [0, y_max_1], line = :dash, label = "r = 2.8", color = :red, lw=1.5) + +p_5, y_max = problem_9_14(5.0, bins, 10.0) +plot!(p_5, [1.1, 1.1], [0, y_max_1], line = :dash, label = "r = 1.1", color = :green, ylim=(0, y_max_1), lw=1.5, legend=:topright) +plot!(p_5, [1.9, 1.9], [0, y_max_1], line = :dash, label = "r = 1.9", lw=1.5, color = :orange) +plot!(p_5, [2.8, 2.8], [0, y_max_1], line = :dash, label = "r = 2.8", lw=1.5, color = :red) + +using Plots.PlotMeasures # for padding + +plot(p_05, p_5, layout = (1,2), size = (1280,600), bottom_margin = 10mm, left_margin = 10mm, right_margin = 10mm, top_margin = 10mm) + +savefig("9-14.png") \ No newline at end of file diff --git a/hw5/9-14.png b/hw5/9-14.png new file mode 100644 index 0000000..2521999 Binary files /dev/null and b/hw5/9-14.png differ diff --git a/hw5/9-2.jl b/hw5/9-2.jl new file mode 100644 index 0000000..ef7af10 --- /dev/null +++ b/hw5/9-2.jl @@ -0,0 +1,848 @@ +# 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 # for margins + +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) + println("\tposition configuration: triangular 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 + + for i in 0:(n-1) + for j in 0:(n-1) + sys.state[2*(i*n+j)+1] += (j % 2) * latticeSpacing / 2 + end + end +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 + zero_total_momentum!(sys) + velocities(sys.state) .*= sqrt(sys.T₀/temperature(sys)) +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 evolve!(sys::ParticleSystem, runtime::Float64=10.0) + numsteps = Int64(abs(runtime/sys.dt) + 1) + + print_evolution_message(runtime, numsteps) + + @time for step in 1:numsteps + 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 ], t="trajectories") + 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, + #label = "pcl. $n @t=t₀", + widen = false, + legend = false, + top_margin=10mm, + bottom_margin=10mm, + left_margin=10mm, + right_margin=10mm + ) + + # 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!(t) + plot!() +end + +function plot_temperature(sys::ParticleSystem) + initialize_plot() + plot!( + sys.timeData, + sys.tempData, + #widen = true, + ) + ylims!( + mean(sys.tempData) - std(sys.tempData) * 20, + mean(sys.tempData) + std(sys.tempData) * 20, + ) + 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, + ) + 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, title="speed distribution") + 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 + ) + # assert that speed is not negative + speed_min = minimum(speed) + speed_max = maximum(speed) + density!( + speed, + normalize = :pdf, + xlims = (speed_min, speed_max), + label = "t = $(round(sys.timeData[s], digits=2))", + top_margin=10mm, + bottom_margin=10mm, + left_margin=10mm, + right_margin=10mm, + ) + end + + + xlabel!("speed") + ylabel!("pdf(speed)") + title!(title) +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 0: APPROACH TO EQUILIBRIUM +function demo_0(;gif=0) + println("\nDEMO 0: APPROACH TO EQUILIBRIUM") + println("----------------------------------------") + + sys = ParticleSystem(64, 120.0, 1.0) + print_system_parameters(sys) + + set_square_lattice_positions!(sys) + set_random_velocities!(sys) + print_system_data(sys) + p1 = plot_positions(sys) + + evolve!(sys, 40.0) + print_system_data(sys) + + p2 = plot_trajectories(sys, collect(1:64)) + p3 = plot_energy(sys) + p4 = plot_temperature(sys) + + # make gif + if gif == 1 + animate(sys, 1) + end + + plot( + p1, p2, p3, p4, + layout = grid(2,2, heights=[0.7,0.3]), + size = (1280,720) + ) +end + +# DEMO 1: TIME REVERSAL TEST +function demo_1(;gif=0) + println("\nDEMO 1: TIME REVERSAL TEST") + println("----------------------------------------") + + sys = ParticleSystem(64, 120.0, 1.0) + print_system_parameters(sys) + + set_square_lattice_positions!(sys) + set_random_velocities!(sys) + print_system_data(sys) + p1 = plot_positions(sys) + + evolve!(sys, 50.0) + #p2 = plot_trajectories(sys, collect(1:64)) + p2 = plot_positions(sys) + + reverse_time!(sys) + evolve!(sys, 50.0) + print_system_data(sys) + #p3 = plot_trajectories(sys, collect(1:64)) + p3 = plot_positions(sys) + + # make gif + if gif == 1 + animate(sys, 4) + end + + plot( + p1, p2, p3, + layout = (1,3), + size = (1200,400) + ) +end + +# DEMO 2: SPEED DISTRIBUTION +function demo_2() + println("\nDEMO 2: SPEED DISTRIBUTION") + println("----------------------------------------") + + sys = ParticleSystem[] + + # array for speed distribution plots + ps = Plots.Plot{Plots.GRBackend}[] + + # array for trajectory plots + pt = Plots.Plot{Plots.GRBackend}[] + + # initialize three systems with different initial conditions + # but same KE and PE, evolve, and save plots + for i = 1:3 + push!(sys, ParticleSystem(64, 120.0, 1.0)) + + println("\nSYSTEM $i") + print_system_parameters(sys[i]) + + set_square_lattice_positions!(sys[i]) + add_position_jitter!(sys[i]) + set_random_velocities!(sys[i]) + print_system_data(sys[i]) + + evolve!(sys[i], 500.0) + print_system_data(sys[i]) + push!(ps, plot_speed_distribution(sys[i], 5)) + push!(pt, plot_trajectories(sys[i], collect(1:64)) ) + end + + + # plot speed distribution and trajectory plots + plot( + ps[1], ps[2], ps[3], + pt[1], pt[2], pt[3], + layout = (2,3), + size = (1920,1080) + ) +end + +# DEMO 3: MELTING TRANSITION +function demo_3(;gif=0) + println("\nDEMO 3: MELTING TRANSITION") + println("----------------------------------------") + + # initialize system of particles on square lattice with zero velocity + sys = ParticleSystem(100, 10.0, 5.0) + set_square_lattice_positions!(sys) + print_system_data(sys) + p1 = plot_positions(sys) + + # evolve the system and watch them "crystallize" + # into a triangular lattice formation + evolve!(sys, 20.0) + print_system_data(sys) + p2 = plot_trajectories(sys, collect(1:100)) + + # now, increase the temperature of the system by giving the particles + # some velocity. evolve the system and plot the trajectories. + set_random_velocities!(sys) + evolve!(sys, 60.0) + print_system_data(sys) + p3 = plot_trajectories(sys, collect(1:100)) + + # some more plots + p4 = plot_energy(sys, 0.0) + p5 = plot_temperature(sys) + p6 = plot_speed_distribution(sys, 20) + + # make gif + if gif == 1 + animate(sys, 1) + end + + plot( + p1, p2, p3, p4, p5, p6, + layout = (2,3), + size = (1280,720) + ) +end + +function problem_9_2() + println("\nProblem 9.2") + println("----------------------------------------") + + sys = ParticleSystem[] + + # array for speed distribution plots + ps = Plots.Plot{Plots.GRBackend}[] + + # array for trajectory plots + pt = Plots.Plot{Plots.GRBackend}[] + + # initialize three systems with different initial conditions + # but same KE and PE, evolve, and save plots + for i = 1:3 + push!(sys, ParticleSystem(81, 1000.0, 1.0)) # large box so particles are far away + + println("\nSYSTEM $i") + print_system_parameters(sys[i]) + if i == 1 + set_square_lattice_positions!(sys[i]) + elseif i == 2 + set_triangular_lattice_positions!(sys[i]) + elseif i == 3 + set_random_positions!(sys[i]) + else + println("oops") + exit() + end + # add_position_jitter!(sys[i]) + set_random_velocities!(sys[i]) + print_system_data(sys[i]) + + evolve!(sys[i], 300.0) + print_system_data(sys[i]) + push!(ps, plot_speed_distribution(sys[i], 2, "speed distribution over time, system $i")) + push!(pt, plot_trajectories(sys[i], collect(1:64), "initial positions of particles, system $i")) + end + + + # plot speed distribution and trajectory plots + plot( + pt[1], ps[1], + pt[2], ps[2], + pt[3], ps[3], + layout = (3,2), + size = (1080,1200) + ) + + savefig("problem_9_2_long-255p.png") +end + +# demo_0() +problem_9_2() \ No newline at end of file diff --git a/hw5/9-2.png b/hw5/9-2.png new file mode 100644 index 0000000..a59182a Binary files /dev/null and b/hw5/9-2.png differ diff --git a/hw5/HarmonicOscillator.jl b/hw5/HarmonicOscillator.jl new file mode 100644 index 0000000..adfaff8 --- /dev/null +++ b/hw5/HarmonicOscillator.jl @@ -0,0 +1,71 @@ +#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia + +# Simulate anharmonic oscillator that may be damped and driven + +using Plots # for plotting trajectory +using DifferentialEquations # for solving ODEs + +ω0 = 1.0 # ω0^2 = k/m +β = 0.0 # β = b/m = friction +c = 10.0 # anharmonic term +f = 0.3 # forcing amplitude +ω = 2.0 # forcing frequency +param = (ω0, β, c, f, ω) # parameters of anharmonic oscillator + +function tendency!(dxp::Vector{Float64}, xp::Vector{Float64}, param, t::Float64) + + (x, p) = xp # 2d phase space + + (ω0, β, c, f, ω) = param + + a = -ω0^2 * x - β * p - c * x^3 + f * forcing(t, ω) # acceleration with m = 1 + + dxp[1] = p + dxp[2] = a + +end + +function forcing(t::Float64, ω::Float64) + + return cos(ω * t) + +end + +function energy(xp::Vector{Float64}, param) + + (x, p) = xp + + (ω0, β, c, f, ω) = param + + pe = 0.5 * ω0^2 * x^2 + (c/4.0) * x^4 + ke = 0.5 * p^2 + + return pe + ke + +end + +x0 = 0.0 # initial position in meters +p0 = 0.0 # initial velocity in m/s +xp0 = [x0, p0] # initial condition in phase space +t_final = 100.0 # final time of simulation + +tspan = (0.0, t_final) # span of time to simulate + +prob = ODEProblem(tendency!, xp0, tspan, param) # specify ODE +sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # solve using Tsit5 algorithm to specified accuracy + +sample_times = sol.t +println("\n\t Results") +println("final time = ", sample_times[end]) +println("Initial energy = ", energy(sol[:,1], param)) +println("Final energy = ", energy(sol[:, end], param)) + +(ω0, β, c, f, ω) = param + +# Plot of position vs. time +xt = plot(sample_times, [sol[1, :] f * forcing.(sample_times, ω)], xlabel = "t", ylabel = "x(t)", legend = false, title = "x vs. t") + +# Phase space plot +xp = plot(sol[1, :], sol[2, :], xlabel = "x", ylabel = "p", legend = false, title = "phase space") + +plot(xt, xp) \ No newline at end of file diff --git a/hw5/MolecularDynamics2.jl b/hw5/MolecularDynamics2.jl new file mode 100644 index 0000000..6580a40 --- /dev/null +++ b/hw5/MolecularDynamics2.jl @@ -0,0 +1,755 @@ +# 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 + +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 + zero_total_momentum!(sys) + velocities(sys.state) .*= sqrt(sys.T₀/temperature(sys)) +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 evolve!(sys::ParticleSystem, runtime::Float64=10.0) + numsteps = Int64(abs(runtime/sys.dt) + 1) + + print_evolution_message(runtime, numsteps) + + @time for step in 1:numsteps + 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, + ) + ylims!( + mean(sys.tempData) - std(sys.tempData) * 20, + mean(sys.tempData) + std(sys.tempData) * 20, + ) + 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, + ) + 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 0: APPROACH TO EQUILIBRIUM +function demo_0(;gif=0) + println("\nDEMO 0: APPROACH TO EQUILIBRIUM") + println("----------------------------------------") + + sys = ParticleSystem(64, 120.0, 1.0) + print_system_parameters(sys) + + set_square_lattice_positions!(sys) + set_random_velocities!(sys) + print_system_data(sys) + p1 = plot_positions(sys) + + evolve!(sys, 20.0) + print_system_data(sys) + + p2 = plot_trajectories(sys, collect(1:64)) + p3 = plot_energy(sys) + p4 = plot_temperature(sys) + + # make gif + if gif == 1 + animate(sys, 1) + end + + plot( + p1, p2, p3, p4, + layout = grid(2,2, heights=[0.7,0.3]), + size = (1280,720) + ) +end + +# DEMO 1: TIME REVERSAL TEST +function demo_1(;gif=0) + println("\nDEMO 1: TIME REVERSAL TEST") + println("----------------------------------------") + + sys = ParticleSystem(64, 120.0, 1.0) + print_system_parameters(sys) + + set_square_lattice_positions!(sys) + set_random_velocities!(sys) + print_system_data(sys) + p1 = plot_positions(sys) + + evolve!(sys, 50.0) + #p2 = plot_trajectories(sys, collect(1:64)) + p2 = plot_positions(sys) + + reverse_time!(sys) + evolve!(sys, 50.0) + print_system_data(sys) + #p3 = plot_trajectories(sys, collect(1:64)) + p3 = plot_positions(sys) + + # make gif + if gif == 1 + animate(sys, 4) + end + + plot( + p1, p2, p3, + layout = (1,3), + size = (1200,400) + ) +end + +# DEMO 2: SPEED DISTRIBUTION +function demo_2() + println("\nDEMO 2: SPEED DISTRIBUTION") + println("----------------------------------------") + + sys = ParticleSystem[] + + # array for speed distribution plots + ps = Plots.Plot{Plots.GRBackend}[] + + # array for trajectory plots + pt = Plots.Plot{Plots.GRBackend}[] + + # initialize three systems with different initial conditions + # but same KE and PE, evolve, and save plots + for i = 1:3 + push!(sys, ParticleSystem(64, 120.0, 1.0)) + + println("\nSYSTEM $i") + print_system_parameters(sys[i]) + + set_square_lattice_positions!(sys[i]) + add_position_jitter!(sys[i]) + set_random_velocities!(sys[i]) + print_system_data(sys[i]) + + evolve!(sys[i], 40.0) + print_system_data(sys[i]) + push!(ps, plot_speed_distribution(sys[i], 5)) + push!(pt, plot_trajectories(sys[i], collect(1:64)) ) + end + + + # plot speed distribution and trajectory plots + plot( + ps[1], ps[2], ps[3], + pt[1], pt[2], pt[3], + layout = (2,3), + size = (1920,1080) + ) +end + +# DEMO 3: MELTING TRANSITION +function demo_3(;gif=0) + println("\nDEMO 3: MELTING TRANSITION") + println("----------------------------------------") + + # initialize system of particles on square lattice with zero velocity + sys = ParticleSystem(100, 10.0, 5.0) + set_square_lattice_positions!(sys) + print_system_data(sys) + p1 = plot_positions(sys) + + # evolve the system and watch them "crystallize" + # into a triangular lattice formation + evolve!(sys, 20.0) + print_system_data(sys) + p2 = plot_trajectories(sys, collect(1:100)) + + # now, increase the temperature of the system by giving the particles + # some velocity. evolve the system and plot the trajectories. + set_random_velocities!(sys) + evolve!(sys, 60.0) + print_system_data(sys) + p3 = plot_trajectories(sys, collect(1:100)) + + # some more plots + p4 = plot_energy(sys, 0.0) + p5 = plot_temperature(sys) + p6 = plot_speed_distribution(sys, 20) + + # make gif + if gif == 1 + animate(sys, 1) + end + + plot( + p1, p2, p3, p4, p5, p6, + layout = (2,3), + size = (1280,720) + ) +end + +demo_2() \ No newline at end of file diff --git a/hw5/demo_2-f.png b/hw5/demo_2-f.png new file mode 100644 index 0000000..1595c9b Binary files /dev/null and b/hw5/demo_2-f.png differ diff --git a/hw5/demo_2.png b/hw5/demo_2.png new file mode 100644 index 0000000..c84dec6 Binary files /dev/null and b/hw5/demo_2.png differ diff --git a/hw5/fpu-2.jl b/hw5/fpu-2.jl new file mode 100644 index 0000000..3190f2b --- /dev/null +++ b/hw5/fpu-2.jl @@ -0,0 +1,138 @@ +using Plots + +N = 32 # number of masses +b =.3 # cubic string spring +A = 10 # amplitude +modes = 3 # number of modes to plot +final_time = 30 # seconds +dt = .05 # seconds + + +# get the intial positions +function calculate_x_i(mass_num, node_num, num_masses, amplitutude) + return A * + sqrt(2 / (num_masses + 1)) * + sin((mass_num * node_num * π) / (num_masses + 1)) +end + +# dynamics calculations +function dynamics!( + mass_displacement, # 2d array + params) + (beta, delta_t, num_masses, num_steps) = params + + for step in 1:num_steps-1 + if step == 1 + continue + end + for mass_num in 2:num_masses-1 + + if step == 1 + prev_mass_pos = 0 + else + prev_mass_pos = mass_displacement[mass_num, step - 1] + end + + right_mass_pos = mass_displacement[mass_num + 1, step] + left_mass_pos = mass_displacement[mass_num - 1, step] + + mass_pos = mass_displacement[mass_num, step] + + xs[mass_num, step + 1] = caluclate_next_displacement( + mass_pos, prev_mass_pos, + left_mass_pos, right_mass_pos, beta, delta_t) + end + + # update the end points + mass_displacement[1, step + 1] = caluclate_next_displacement( + 0, 0, + 0, mass_displacement[2, step], beta, delta_t) + mass_displacement[num_masses, step + 1] = caluclate_next_displacement( + 0, 0, + mass_displacement[num_masses - 1, step], 0, beta, delta_t) + end + +end + +function caluclate_next_displacement( + mass_pos, prev_mass_pos, + left_mass_pos, right_mass_pos, + beta, delta_t +) + # println(mass_pos, " " , prev_mass_pos, " ", left_mass_pos, " ", right_mass_pos, '\n') + return 2 * mass_pos - prev_mass_pos + + delta_t^(2) * (right_mass_pos + left_mass_pos - 2 * mass_pos) + + beta * delta_t^(2) * ((right_mass_pos - mass_pos)^3 + (left_mass_pos - mass_pos)^3) +end + + +# energy calcuations, after the dynamics +function calculate_a_k(k, num_masses, delta_t, xs_n, beta) + sum = 0 + for i in 1:num_masses + sum += xs_n[i] * sin((k * i * π) / (num_masses + 1)) + end + return sqrt(2 / (num_masses + 1)) * sum +end + +function calculate_energy(a_k, a_k_plus1, delta_t, mode, num_masses) + kinetic = .5 * ((a_k_plus1 - a_k) / delta_t)^2 + + # sum over the three modes + w_k = 2 * sin((mode * π) / (2 * (num_masses + 1))) + potential = .5 * (w_k^2 * a_k^2) + + return kinetic + potential +end + + +# do the simulation +num_steps = Int(final_time / dt) +params = (b, dt, N, num_steps) +# build the 2d array of mass displacements to time +xs = zeros(N, num_steps) +# fill in the initial displacements +for mass_num in 2:N-1 + xs[mass_num, 1] = calculate_x_i(mass_num, 1, N, A) +end +dynamics!(xs, params) + +# println(xs) + + +# plot these displacements over time as an animation +a = @animate for i in 1:num_steps + plot(xs[:, i], legend=false, + marker = :circle, xlabel="Mass Number", ylabel="Displacement", + yaxis = (-30, 30), title="Displacement Over Time" + ) +end + +gif(a, "fpu.gif", fps=15) + +# plot the first two timespets positions +# plot(xs[:, 1], label="t=0", marker=:circle, xlabel="Mass Number", ylabel="Displacement", title="First Two Time Steps") +# plot!(xs[:, 2], label="t=1", marker=:circle) +# plot!(xs[:, 3], label="t=2", marker=:circle) +# plot!(xs[:, 4], label="t=3", marker=:circle) +# plot!(xs[:, 5], label="t=4", marker=:circle) +# plot!(xs[:, 6], label="t=5", marker=:circle) +# plot!(xs[:, 7], label="t=6", marker=:circle) + + +# # calculate the energies for each mode from the displacements +# energies = zeros(modes, num_steps) +# for mode in 1:modes +# energies[mode, :] = calculate_energy_for_mode(mode, xs, N, dt, b) +# end + +# make a range time steps +# time = range(0, final_time, length=num_steps) + +# println("e:", length(energies[1, :])) +# println("t:", length(time)) + +# plot the energies for each mode +# scatter(time, energies[1, :], label="Mode 1", marker=:circle, xlabel="Time", ylabel="Energy", title="Energy for First Three Modes") +# scatter!(time, energies[2, :], label="Mode 2", marker=:circle) +# scatter!(time, energies[3, :], label="Mode 3", marker=:circle) \ No newline at end of file diff --git a/hw5/fpu-3.jl b/hw5/fpu-3.jl new file mode 100644 index 0000000..f2ceb38 --- /dev/null +++ b/hw5/fpu-3.jl @@ -0,0 +1,78 @@ +N = 32 # number of masses +b = 10 # cubic string spring +A = 1 # amplitude +modes = 3 # number of modes to plot +final_time = 10 # seconds +dt = .05 # seconds + +function kinetic_energy(a_k, a_k_plus1, delta_t) + return .5 * ((a_k_plus1 - a_k) / delta_t)^2 +end + +function potential_energy(a_k, a_k_plus1, mode, num_masses) + w_k = 2 * sin((mode * π) / (2 * (num_masses + 1))) + return .5 * (w_k^2 * a_k^2) +end + +function calculate_energy(a_k, a_k_plus1, delta_t, mode, num_masses) + k = kinetic_energy(a_k, a_k_plus1, delta_t) + p = potential_energy(a_k, a_k_plus1, mode, num_masses) + + return k + p +end + +function update_state!(state, prev_state, dt, beta) + x_prev = prev_state + x = copy(state) + + # update the left particle (set to rest) + state[1] = 0 + + # update the right particle (set to rest) + state[N] = 0 + + # update the middle particles + for i in 2:N-1 + state[i] = 2 * x[i] - x_prev[i] + + dt * dt * ((x[i + 1] - 2 * x[i] + x[i - 1]) + + dt * dt * beta * ((x[i + 1] - x[i])^3 - (x[i - 1] - x[i])^3) + ) + end +end + +function initial_state(N, modes, beta, A) + # make the range of masses + masses = 1:N + # make the range of modes + init_state = A * sin.((modes * π * masses) / (N + 1)) + init_state[1] = 0 + init_state[N] = 0 + return init_state +end + +function run_simulation(N, modes, beta, A, dt, num_steps) + data = [] + state = initial_state(N, modes, beta, A) + prev_state = copy(state) + println("Initial state: ", state) + for i in 1:num_steps + update_state!(state, prev_state, dt, beta) + push!(data, state) + prev_state = copy(state) + end + + return data +end + +steps = Int(final_time / dt) + +p = run_simulation(N, modes, b, A, dt, steps) + +using Plots +# plot the first and final position +plot(p[1], label="Initial position") +plot!(p[end], label="Final position") + +# plot displacement for the first particle over time +t_span = 0:dt:final_time +plot(t_span, p, label="Displacement for first particle") diff --git a/hw5/fpu.gif b/hw5/fpu.gif new file mode 100644 index 0000000..e1d1612 Binary files /dev/null and b/hw5/fpu.gif differ diff --git a/hw5/fpu.jl b/hw5/fpu.jl new file mode 100644 index 0000000..327cf3a --- /dev/null +++ b/hw5/fpu.jl @@ -0,0 +1,35 @@ +N = 32 # number of masses +beta =.3 # cubic string spring +A = 10 # amplitude +modes = 3 # number of modes to plot + + +function calculate_x_i(mass_num, node_num, num_masses, amplitutude) + return A * + sqrt(2 / (num_masses + 1)) * + sin((mass_num * node_num * π) / (num_masses + 1)) +end + +function calculate_x_s_for_mode(mass_num, node_num, num_masses, amplitutude) + x = 0 + for i in 1:num_masses + x += calculate_x_i(mass_num, node_num, num_masses, amplitutude) + end + return x +end + +# calculate the x_s for the first three modes + +x_s = zeros(modes, N+2) # add two for zeros on ends +# add the first and last 0s +for i in 1:modes + for j in 1:N + x_s[i, j+1] = calculate_x_s_for_mode(j, i, N, A) + end +end + +# plot the first three modes +using Plots +plot(x_s[1, :], label="Mode 1", marker=:circle, xlabel="Mass Number", ylabel="Displacement", title="First Three Modes") +plot!(x_s[2, :], label="Mode 2", marker=:circle) +plot!(x_s[3, :], label="Mode 3", marker=:circle) \ No newline at end of file diff --git a/hw5/fpu.png b/hw5/fpu.png new file mode 100644 index 0000000..52b6839 Binary files /dev/null and b/hw5/fpu.png differ diff --git a/hw5/lfg-modes.png b/hw5/lfg-modes.png new file mode 100644 index 0000000..ceffab8 Binary files /dev/null and b/hw5/lfg-modes.png differ diff --git a/hw5/modes-beta15.png b/hw5/modes-beta15.png new file mode 100644 index 0000000..6536741 Binary files /dev/null and b/hw5/modes-beta15.png differ diff --git a/hw5/modes-beta3.png b/hw5/modes-beta3.png new file mode 100644 index 0000000..1b2833f Binary files /dev/null and b/hw5/modes-beta3.png differ diff --git a/hw5/plot_175.svg b/hw5/plot_175.svg new file mode 100644 index 0000000..4733b20 --- /dev/null +++ b/hw5/plot_175.svg @@ -0,0 +1,105 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw5/plz-modes.png b/hw5/plz-modes.png new file mode 100644 index 0000000..465252b Binary files /dev/null and b/hw5/plz-modes.png differ diff --git a/hw5/plz.gif b/hw5/plz.gif new file mode 100644 index 0000000..febf678 Binary files /dev/null and b/hw5/plz.gif differ diff --git a/hw5/plz.jl b/hw5/plz.jl new file mode 100644 index 0000000..d4abcbd --- /dev/null +++ b/hw5/plz.jl @@ -0,0 +1,148 @@ +function dynamics!( + state, + prev_state, + params, + states +) + # Unpack the parameters + N, modes, beta, A, dt, num_steps = params + + for i in 1:num_steps + next_state = zeros(N) + + # update the left particle (set to rest) + state[1] = 0 + + # update the right particle (set to rest) + state[N] = 0 + + # update the middle particles + for i in 2:N-1 + a = 2 * state[i] - prev_state[i] + b = dt * dt * (state[i + 1] - 2 * state[i] + state[i - 1]) + c = dt * dt * beta * ((state[i + 1] - state[i])^3 + (state[i - 1] - state[i])^3) + next_state[i] = a + b + c + end + + push!(states, next_state) + + # update the previous state + for i in 1:N + prev_state[i] = state[i] + end + # update the current state + for i in 1:N + state[i] = next_state[i] + end + end +end + +function get_initial_state( + N, + modes, + beta, + A +) + state = zeros(N) + amp = A * sqrt(2 / (N - 1)) + for i in 2:N-1 + state[i] = amp * sin((modes * π * (i - 1)) / (N - 1)) + end + return state +end + +function run_simulation( + N, + modes, + beta, + A, + dt, + num_steps +) + states = [] + state = get_initial_state(N, 1, beta, A) + prev_state = zeros(N) + for i in 1:N + prev_state[i] = state[i] + end + params = (N, modes, beta, A, dt, num_steps) + dynamics!(state, prev_state, params, states) + return states +end + +# Run the simulation +N = 32 # number of masses +beta = 1.5 # cubic string spring +A = 10 # amplitude +modes = 3 # number of modes to plot +final_time = 50 # seconds +dt = .5 # seconds +num_steps = Int(final_time / dt) +params = (N, modes, beta, A, dt, num_steps) +s = run_simulation(N, modes, beta, A, dt, num_steps) +println("Final state: ", s[end]) + +# plot the inital positions and the final positions +using Plots +# plot(get_initial_state(N, 1, beta, A), label="Initial", marker=:circle, xlabel="Mass Number", ylabel="Displacement", title="First Three Modes") +# plot!(s[end], label="Final", marker=:circle) + +# animate the s array of positions +anim = @animate for i in 1:length(s) + plot(s[i], label="t = $(i * dt)", marker=:circle, xlabel="Mass Number", ylabel="Displacement", ylim=(-3, 3)) +end +gif(anim, "plz.gif", fps = 30) + +# function caluclate_energies_for_mode(states, mode) +# total = [] +# kinetic = [] +# potential = [] +# prev_a_k = 0 +# for i in 1:length(states) - 1 +# total_energy = 0 +# kinetic_energy = 0 +# potential_energy = 0 + +# sum = 0 +# for j in 1:length(states[i]) +# sum += states[i][j] * sin((mode * j * π) / (N - 1)) +# end +# amp = A * sqrt(2 / (N - 1)) +# a_k = amp * sum + +# k = .5 * ((a_k - prev_a_k) / dt)^2 +# if i == 1 +# k = 0 +# end +# p = (2 * sin((mode * π) / (2 * (N - 1)))^2 * a_k^2) + + +# total_energy += k + p +# kinetic_energy += k +# potential_energy += p +# push!(total, total_energy) +# push!(kinetic, kinetic_energy) +# push!(potential, potential_energy) + +# prev_a_k = a_k +# end +# return total, kinetic, potential +# end + +# # calculate the energies for each mode from the displacements +# total_1, kinetic_1, potential_2 = caluclate_energies_for_mode(s, 1) +# total_2, kinetic_2, potential_2 = caluclate_energies_for_mode(s, 2) +# total_3, kinetic_3, potential_3 = caluclate_energies_for_mode(s, 3) + +# # build a timestep space +# # timesteps = [i * dt for i in 1:num_steps - 1] +# # # plot the modes energies on the same y-axis +# plot(timesteps, total_1, label="Mode 1", xlabel="Time", ylabel="Energy of Modes", title="Energy Over Time (\$\\beta = 1.5\$)") +# plot!(timesteps, total_2, label="Mode 3") +# plot!(timesteps, total_3, label="Mode 5") + +# savefig("modes-beta15.png") + + + + diff --git a/hw5/plz.png b/hw5/plz.png new file mode 100644 index 0000000..d6e7209 Binary files /dev/null and b/hw5/plz.png differ diff --git a/hw5/problem_9_-test.png b/hw5/problem_9_-test.png new file mode 100644 index 0000000..894c5be Binary files /dev/null and b/hw5/problem_9_-test.png differ diff --git a/hw5/problem_9_11.png b/hw5/problem_9_11.png new file mode 100644 index 0000000..abe5054 Binary files /dev/null and b/hw5/problem_9_11.png differ diff --git a/hw5/problem_9_11_100.png b/hw5/problem_9_11_100.png new file mode 100644 index 0000000..e1025cb Binary files /dev/null and b/hw5/problem_9_11_100.png differ diff --git a/hw5/problem_9_11_16.png b/hw5/problem_9_11_16.png new file mode 100644 index 0000000..98b4a68 Binary files /dev/null and b/hw5/problem_9_11_16.png differ diff --git a/hw5/problem_9_11_25.png b/hw5/problem_9_11_25.png new file mode 100644 index 0000000..f4f6e94 Binary files /dev/null and b/hw5/problem_9_11_25.png differ diff --git a/hw5/problem_9_11_9.png b/hw5/problem_9_11_9.png new file mode 100644 index 0000000..1565697 Binary files /dev/null and b/hw5/problem_9_11_9.png differ diff --git a/hw5/problem_9_11_old.png b/hw5/problem_9_11_old.png new file mode 100644 index 0000000..a99ed6a Binary files /dev/null and b/hw5/problem_9_11_old.png differ diff --git a/hw5/problem_9_11_t.png b/hw5/problem_9_11_t.png new file mode 100644 index 0000000..aa8f3a6 Binary files /dev/null and b/hw5/problem_9_11_t.png differ diff --git a/hw5/problem_9_2-fail.png b/hw5/problem_9_2-fail.png new file mode 100644 index 0000000..95eb54f Binary files /dev/null and b/hw5/problem_9_2-fail.png differ diff --git a/hw5/problem_9_2-w.png b/hw5/problem_9_2-w.png new file mode 100644 index 0000000..e72e705 Binary files /dev/null and b/hw5/problem_9_2-w.png differ diff --git a/hw5/problem_9_2.png b/hw5/problem_9_2.png new file mode 100644 index 0000000..79c56f7 Binary files /dev/null and b/hw5/problem_9_2.png differ diff --git a/hw5/problem_9_2_1000.png b/hw5/problem_9_2_1000.png new file mode 100644 index 0000000..552725b Binary files /dev/null and b/hw5/problem_9_2_1000.png differ diff --git a/hw5/problem_9_2_L.png b/hw5/problem_9_2_L.png new file mode 100644 index 0000000..612553e Binary files /dev/null and b/hw5/problem_9_2_L.png differ diff --git a/hw5/problem_9_2_final.png b/hw5/problem_9_2_final.png new file mode 100644 index 0000000..98d5a91 Binary files /dev/null and b/hw5/problem_9_2_final.png differ diff --git a/hw5/problem_9_2_long-100p.png b/hw5/problem_9_2_long-100p.png new file mode 100644 index 0000000..be00a71 Binary files /dev/null and b/hw5/problem_9_2_long-100p.png differ diff --git a/hw5/problem_9_2_long-255p.png b/hw5/problem_9_2_long-255p.png new file mode 100644 index 0000000..11bb204 Binary files /dev/null and b/hw5/problem_9_2_long-255p.png differ diff --git a/hw5/problem_9_2_long.png b/hw5/problem_9_2_long.png new file mode 100644 index 0000000..353903b Binary files /dev/null and b/hw5/problem_9_2_long.png differ diff --git a/hw5/r_squared-1.png b/hw5/r_squared-1.png new file mode 100644 index 0000000..a91d331 Binary files /dev/null and b/hw5/r_squared-1.png differ -- cgit v1.2.3-70-g09d2