# 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")