update naming

1 files changed, 0 insertions, 755 deletions

diff --git a/t/old/t.jl b/t/old/t.jl
deleted file mode 100644
index 6a06777..0000000
--- a/t/old/t.jl
+++ /dev/null
@@ -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()