N-body simulations are moving here from DiffEqPhysics (#2)

transfer of n-body simulation files from DIffEqPhysics

Author: Mikhail-Vaganov

README.md

@@ -1,4 +1,8 @@
 # NbodySimulator
 
 [![Build Status](https://travis-ci.org/JuliaDiffEq/NBodySimulator.jl.svg?branch=master)](https://travis-ci.org/JuliaDiffEq/NBodySimulator.jl)
-[![Build status](https://ci.appveyor.com/api/projects/status/1ofg9ianvcciq26v?svg=true)](https://ci.appveyor.com/project/Mikhail-Vaganov/nbodysimulator-jl)
+[![Build status](https://ci.appveyor.com/api/projects/status/1ofg9ianvcciq26v?svg=true)](https://ci.appveyor.com/project/Mikhail-Vaganov/nbodysimulator-jl)
+
+This project is under development at the moment. The implementation of potential calculations is fairly experimental and has not been extensively verified yet.
+
+You can test simulation of different systems now but be aware of possible changes in future. 

REQUIRE

@@ -1 +1,9 @@
-julia 0.7-beta2
+julia 0.7-beta2
+DiffEqBase 3.0.0
+OrdinaryDiffEq 3.0.0
+Reexport
+StaticArrays
+RecursiveArrayTools
+RecipesBase
+DiffEqCallbacks
+FileIO

examples/liquid_argon.jl

@@ -0,0 +1,96 @@
+using NBodySimulator
+
+function generate_bodies_in_cell_nodes(n::Int, m::Real, v_dev::Real, L::Real)
+   
+    rng = MersenneTwister(n);
+    velocities = v_dev * randn(rng, Float64, (3, n))
+    bodies = MassBody[]
+
+    count = 1
+    dL = L / (ceil(n^(1 / 3)))
+    for x = dL/2:dL:L, y = dL/2:dL:L, z = dL/2:dL:L        
+        if count > n
+            break
+        end
+        r = SVector(x, y, z)
+        v = SVector{3}(velocities[:,count])
+        body = MassBody(r, v, m)
+        push!(bodies, body)
+        count += 1           
+    end
+    return bodies
+end
+
+function generate_bodies_in_line(n::Int, m::Real, v_dev::Real, L::Real)
+    dL = L / (ceil(n^(1 / 3)))
+    n_line = floor(Int, L / dL)
+    rng = MersenneTwister(n);
+    velocities = v_dev * randn(rng, Float64, (3, n_line))
+    bodies = MassBody[]
+    x = y = L / 2
+    for i = 1:n_line      
+        r = SVector(x, y, i * dL)
+        v = SVector{3}(velocities[:,i])
+        body = MassBody(r, v, m)
+        push!(bodies, body)  
+    end
+    return bodies
+end
+
+function generate_random_directions(n::Int)
+    theta = acos.(1 - 2 * rand(n));
+    phi = 2 * pi * rand(n);
+    directions = [@SVector [sin(theta[i]) .* cos(phi[i]), sin(theta[i]) .* sin(phi[i]), cos(theta[i])] for i = 1:n]
+end
+
+units = :real
+units = :reduced
+
+const T = 120.0 # °K
+const kb = 1.38e-23 # J/K
+const ϵ = T * kb
+const σ = 3.4e-10 # m
+const ρ = 1374 # kg/m^3
+const m = 39.95 * 1.6747 * 1e-27 # kg
+const N = 125#floor(Int, ρ * L^3 / m)
+const L = (m*N/ρ)^(1/3)#10.229σ
+const R = 2.25σ   
+const v_dev = sqrt(kb * T / m)
+const τ = 1e-14 # σ/v
+const t1 = 0τ
+const t2 = 2000τ
+#bodies = generate_bodies_randomly(N, m, v_dev, L)
+bodies = generate_bodies_in_cell_nodes(N, m, v_dev, L)
+#bodies = generate_bodies_in_line(N, m, v_dev, L)
+jl_parameters = LennardJonesParameters(ϵ, σ, R)
+pbc = CubicPeriodicBoundaryConditions(L)
+thermostat = AndersenThermostat(0.02, T, kb)
+lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => jl_parameters));
+simulation = NBodySimulation(lj_system, (t1, t2), pbc, thermostat);
+#result = run_simulation(simulation, Tsit5())
+result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)
+
+#=
+using Plots
+import GR
+(rs, grf) = rdf(result)
+(ts, dr2) = msd(result)
+plot(rs/σ, grf, xlim=[0, 0.4999L/σ], label=["Radial distribution function"],ylabel="g(r)", xlabel="r/sigma")
+
+using JLD
+time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
+Nactual = length(bodies)
+timesteps = round(length(result.solution.t))
+#save("D:/water $Nactual molecules $timesteps steps.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2)
+save("D:/liquid argon $Nactual molecules $timesteps steps $time_now.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2, "e_tot", e_tot, "e_kin", e_kin, "e_pot", e_pot)
+=#
+#=
+using Plots
+import GR
+time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
+Nactual = length(bodies)
+timesteps = round(length(result.solution.t))
+@time animate(result, "D:/$Nactual liquid argon particles with $timesteps timesteps $time_now.gif")
+
+#plot(simulation)
+=#

examples/liquid_argon_omm_units.jl

@@ -0,0 +1,53 @@
+using NBodySimulator
+using StochasticDiffEq
+
+T = 120.0 # °K
+T0 = 90.0 # °K
+kb = 8.3144598e-3 # kJ/(K*mol)
+ϵ = T * kb
+σ = 0.34 # nm
+ρ = 1374/1.6747# Da/nm^3
+m = 39.95# Da
+N = 216
+L = (m*N/ρ)^(1/3)#10.229σ
+R = 0.5*L   
+v_dev = sqrt(kb * T / m)
+bodies = generate_bodies_in_cell_nodes(N, m, v_dev, L)
+
+τ = 0.5e-3 # ps or 1e-12 s
+t1 = 0.0
+t2 = 2000τ
+
+parameters = LennardJonesParameters(ϵ, σ, R)
+lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => parameters));
+thermostat = AndersenThermostat(90, 0.01/τ)
+#thermostat = BerendsenThermostat(90, 2000τ)
+#thermostat = LangevinThermostat(90, 0.00)
+#thermostat = NoseHooverThermostat(T0, 200τ)
+pbc = CubicPeriodicBoundaryConditions(L)
+simulation = NBodySimulation(lj_system, (t1, t2), pbc, thermostat, kb);
+result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)
+#result = @time run_simulation_sde(simulation, ISSEM(symplectic=true,theta=0.5))
+
+#t = t1:τ:result.solution.t[end-1]
+#temper = temperature.(result, t)
+
+
+#using Plots
+#pl=plot(t, temper, ylim=[0,200], xlabel="t, ps", ylabel = "T, °K", label="Temperature, °K", linewidth=2)
+#plot!(pl, t, T0*ones(length(t)), label = "90 °K", linewidth=2)
+#plot!(pl, title="Andersen thermostat at dt*v=$(thermostat.ν*τ) ")
+#plot!(pl, title="Berendsen thermostat at tau=$(thermostat.τ) ps")
+#plot!(pl, title="Nose-Hoover thermostat at tau=$(thermostat.τ) ps")
+
+
+#time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
+#Nactual = length(bodies)
+#timesteps = round(length(result.solution.t))
+
+#@time save_to_pdb(result, "D:/liquid argon simulation $Nactual molecules and $timesteps steps $time_now.pdb" )
+
+#using JLD
+#save("d:/nosehoover thermostat for water liquid argon $T0 and $(thermostat.τ) _$time_now.jld", "t",t, "temper", temper, "T0", T0, "τ", thermostat.τ)
+
+

examples/liquid_argon_reduced.jl

@@ -0,0 +1,66 @@
+using NBodySimulator
+
+function generate_bodies_in_cell_nodes(n::Int, m::Real, v_dev::Real, L::Real)
+   
+    rng = MersenneTwister(n);
+    velocities = v_dev * randn(rng, Float64, (3, n))
+    bodies = MassBody[]
+
+    count = 1
+    dL = L / (ceil(n^(1 / 3)))
+    for x = dL / 2:dL:L, y = dL / 2:dL:L, z = dL / 2:dL:L        
+        if count > n
+            break
+        end
+        r = SVector(x, y, z)
+        v = SVector{3}(velocities[:,count])
+        body = MassBody(r, v, m)
+        push!(bodies, body)
+        count += 1           
+    end
+    return bodies
+end
+
+const T = 120.0 # °K
+const kb = 1.38e-23 # J/K
+const ϵ = T * kb
+const σ = 3.4e-10 # m
+const ρ = 1374 # kg/m^3
+const m = 39.95 * 1.6747 * 1e-27 # kg
+const N = 216#floor(Int, ρ * L^3 / m)
+const L = (m*N/ρ)^(1/3)#10.229σ
+const R = 2.25σ   
+const v_dev = sqrt(kb * T / m)
+const τ = 0.5e-15 # σ/v, fs
+const t1 = 0.0
+const t2 = 3000τ
+
+const _L = L / σ
+const _σ = 1.0
+const _ϵ = 1.0
+const _m = 1.0
+const _v = v_dev / sqrt(ϵ / m)
+const _R = R / σ
+const _τ = τ * sqrt(ϵ / m) / σ
+const _t1 = t1 * sqrt(ϵ / m) / σ
+const _t2 = t2 * sqrt(ϵ / m) / σ
+bodies = generate_bodies_in_cell_nodes(N, _m, _v, _L)
+parameters = LennardJonesParameters(_ϵ, _σ, _R)
+lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => parameters));
+simulation = NBodySimulation(lj_system, (_t1, _t2), CubicPeriodicBoundaryConditions(_L), _ϵ/T);
+#result = run_simulation(simulation, Tsit5())
+result = @time run_simulation(simulation, VelocityVerlet(), dt=_τ)
+
+#(rs, grf) = rdf(result)
+#(ts, dr2) = msd(result)
+
+#using Plots
+#plot(rs, grf, xlim=[0, 0.4999_L], label=["Radial distribution function"],ylabel="g(r)", xlabel="r")
+#plot(rs/_σ, grf, xlim=[0, 0.4999_L/_σ], label=["Radial distribution function"],ylabel="g(r)", xlabel="r/sigma")
+
+
+#time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
+#Nactual = length(bodies)
+#timesteps = round(length(result.solution.t))
+
+#@time save_to_pdb(result, "D:/liquid argon simulation $Nactual molecules and $timesteps steps $time_now.pdb" )

examples/liquid_argon_sde.jl

@@ -0,0 +1,47 @@
+using NBodySimulator
+using StochasticDiffEq
+
+T = 120.0 # °K
+T0 = 90
+kb = 8.3144598e-3 # kJ/(K*mol)
+ϵ = T * kb
+σ = 0.34 # nm
+ρ = 1374/1.6747# Da/nm^3
+m = 39.95# Da
+N = 216
+L = (m*N/ρ)^(1/3)#10.229σ
+R = 0.5*L   
+v_dev = sqrt(kb * T / m)
+bodies = generate_bodies_in_cell_nodes(N, m, v_dev, L)
+
+τ = 0.5e-3 # ps or 1e-12 s
+t1 = 0.0
+t2 = 10000τ
+
+parameters = LennardJonesParameters(ϵ, σ, R)
+lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => parameters));
+#thermostat = AndersenThermostat(T0, 1000)
+#thermostat = BerendsenThermostat(T0, 20τ)
+thermostat = LangevinThermostat(T0, 10.0)
+pbc = CubicPeriodicBoundaryConditions(L)
+simulation = NBodySimulation(lj_system, (t1, t2), pbc, thermostat, kb);
+#result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)
+result = @time run_simulation_sde(simulation, EM(),  dt=τ)
+
+#=
+time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
+Nactual = length(bodies)
+timesteps = round(length(result.solution.t))
+
+
+using Plots
+t = t1:τ:result.solution.t[end-1]
+temper = temperature.(result, t)
+pl=plot(t, temper, ylim=[0,500], xlabel="t, ps", ylabel = "T, °K", label="Temperature, °K", linewidth=2)
+plot!(pl, t, T0*ones(length(t)), label = "$T0 °K", linewidth=2)
+plot!(pl, title="Langevin thermostat at gamma=$(thermostat.γ)  1/ps")
+
+
+using JLD
+save("d:/$(typeof(thermostat)) thermostat for water $T0 _$time_now.jld", "t",t, "temper", temper, "T0", T0, "thermostat", thermostat)
+=#