diff options
Diffstat (limited to 't/t.jl')
| -rw-r--r-- | t/t.jl | 755 |
1 files changed, 0 insertions, 755 deletions
@@ -1,755 +0,0 @@ -# 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_0() |