testing

1 files changed, 713 insertions, 0 deletions

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