Merge pull request #5 from JuliaDiffEq/improving_tests

Improved testing

Author: ChrisRackauckas

examples/liquid_argon.jl

@@ -1,25 +1,4 @@
-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
+using NBodySimulator, StaticArrays
 
 function generate_bodies_in_line(n::Int, m::Real, v_dev::Real, L::Real)
     dL = L / (ceil(n^(1 / 3)))
@@ -47,26 +26,27 @@ units = :real
 units = :reduced
 
 const T = 120.0 # °K
+const T0 = 90.0
 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 N = 216#floor(Int, ρ * L^3 / m)
 const L = (m*N/ρ)^(1/3)#10.229σ
-const R = 2.25σ   
+const R = 0.5*L   
 const v_dev = sqrt(kb * T / m)
-const τ = 1e-14 # σ/v
+const τ = 0.5e-15 # σ/v
 const t1 = 0τ
-const t2 = 2000τ
+const t2 = 30000τ
 #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)
+#thermostat = AndersenThermostat(0.02, T, kb)
 lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => jl_parameters));
-simulation = NBodySimulation(lj_system, (t1, t2), pbc, thermostat);
+simulation = NBodySimulation(lj_system, (t1, t2), pbc, kb);
 #result = run_simulation(simulation, Tsit5())
 result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)
 
@@ -93,4 +73,26 @@ timesteps = round(length(result.solution.t))
 @time animate(result, "D:/$Nactual liquid argon particles with $timesteps timesteps $time_now.gif")
 
 #plot(simulation)
+=#
+
+#=
+(rs, grf) = rdf(result)
+(ts, dr2) = msd(result)
+
+t = t1:τ:result.solution.t[end-1]
+temper = temperature.(result, t)
+
+time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
+Nactual = length(bodies)
+timesteps = round(length(result.solution.t))
+
+using MAT
+file = matopen("d:/liquid argon SI units and $timesteps timesteps_$time_now.mat", "w")
+write(file, "t", collect(t))
+write(file, "temper", temper)
+write(file, "rs", rs)
+write(file, "ts", ts)
+write(file, "grf", grf)
+write(file, "dr2", dr2)
+close(file)
 =#

examples/liquid_argon_omm_units.jl

@@ -1,36 +1,40 @@
 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/τ)
+const T = 120.0 # °K
+const T0 = 90.0 # °K
+const kb = 8.3144598e-3 # kJ/(K*mol)
+const ϵ = T * kb
+const σ = 0.34 # nm
+const ρ = 1374/1.6747# Da/nm^3
+const m = 39.95# Da
+const N = 8
+const L = (m*N/ρ)^(1/3)#10.229σ
+const R = 0.5*L   
+const v_dev = sqrt(kb * T / m)
+const bodies = generate_bodies_in_cell_nodes(N, m, v_dev, L)
+
+const τ = 0.5e-4 # ps or 1e-12 s
+const t1 = 0.0
+const t2 = 3τ
+
+const parameters = LennardJonesParameters(ϵ, σ, R)
+const lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => parameters));
+#const thermostat = AndersenThermostat(90, 0.01/τ)
 #thermostat = BerendsenThermostat(90, 2000τ)
-#thermostat = LangevinThermostat(90, 0.00)
+const thermostat = LangevinThermostat(70, 10.0)
 #thermostat = NoseHooverThermostat(T0, 200τ)
-pbc = CubicPeriodicBoundaryConditions(L)
-simulation = NBodySimulation(lj_system, (t1, t2), pbc, thermostat, kb);
-result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)
+const pbc = CubicPeriodicBoundaryConditions(L)
+const 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))
+result = @time run_simulation(simulation, EM(), dt=τ)
+
+#(rs, grf) = rdf(result)
+#(ts, dr2) = msd(result)
 
-#t = t1:τ:result.solution.t[end-1]
-#temper = temperature.(result, t)
+t = t1:τ:result.solution.t[end-1]
+temper = @time temperature.(result, t)
 
 
 #using Plots
@@ -41,13 +45,33 @@ result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)
 #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_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.τ)
 
 
+#=
+using MAT
+
+time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
+Nactual = length(bodies)
+timesteps = round(length(result.solution.t))
+
+file = matopen("d:/liquid argon $(length(bodies)) omm units at temperature $T0 andgamma $(thermostat.γ)  $timesteps timesteps $T0 _$time_now.mat", "w")
+write(file, "t", collect(t))
+write(file, "temper", temper)
+write(file, "T0", T0)
+if thermostat isa LangevinThermostat
+    write(file, "gamma", thermostat.γ)
+end
+#write(file, "rs", rs)
+#write(file, "ts", ts)
+#write(file, "grf", grf)
+#write(file, "dr2", dr2)
+close(file)
+=#

examples/liquid_argon_reduced.jl

@@ -1,26 +1,5 @@
 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

examples/water_spc_fw.jl

@@ -1,27 +1,5 @@
 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  
-    #for x = 1.5*dL:2dL:5.6*dL, y = 1.5*dL:2dL:5.6*dL, z = 1.5*dL:2dL:5.6*dL        
-        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)
@@ -41,10 +19,11 @@ end
 const qe = 1.6e-19
 const Na = 6.022e23
 const T = 298.16 # °K
+const T0 = 275 # °K
 const kb = 1.38e-23 # J/K
 const ϵOO = 0.1554253*4184/Na
 const σOO = 3.165492e-10 # m
-const ρ = 1000 # kg/m^3
+const ρ = 216 # kg/m^3
 const mO = 15.999 * 1.6747 * 1e-27 # kg
 const mH = 1.00794 * 1.6747 * 1e-27#
 const mH2O = mO+2*mH
@@ -55,32 +34,26 @@ const Rel = 0.49*L
 const v_dev = sqrt(kb * T /mH2O)
 const τ = 0.5e-15 # σ/v
 const t1 = 0τ
-const t2 = 6.5e-9 
+const t2 = 6000τ
 const k_bond = 1059.162*4184*1e20/Na # J/m^2
 const k_angle = 75.90*4184/Na # J/rad^2
 const rOH = 1.012e-10 # m
 const ∠HOH = 113.24*pi/180 # rad
 const qH = 0.41*qe
 const qO = -0.82*qe
 const k = 9e9 #
-bodies = generate_bodies_in_cell_nodes(N, mH2O, v_dev, L)
-#bodies = generate_bodies_in_line(N, mH2O, v_dev, L)
+#bodies = generate_bodies_in_cell_nodes(N, mH2O, v_dev, L)
+bodies = load_water_molecules_from_pdb("C:/Users/Michael/Desktop/GSoC/pdbs/output_4.pdb")
 jl_parameters = LennardJonesParameters(ϵOO, σOO, R)
 e_parameters = ElectrostaticParameters(k, Rel)
 spc_paramters = SPCFwParameters(rOH, ∠HOH, k_bond, k_angle)
 pbc = CubicPeriodicBoundaryConditions(L)
 water = WaterSPCFw(bodies, mH, mO, qH, qO,  jl_parameters, e_parameters, spc_paramters);
-simulation = NBodySimulation(water, (t1, t2), pbc);
+#thermostat = BerendsenThermostat(T0, 200τ)
+simulation = NBodySimulation(water, (t1, t2), pbc, kb);
 #result = run_simulation(simulation, Tsit5())
-result = @time run_simulation(simulation, VelocityVerlet(), dt=τ, saveat=0.5e-12)
-
-
-time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
-Nactual = length(bodies)
-timesteps = round(length(result.solution.t))
+result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)
 
-(rs, grf) = @time rdf(result)
-(ts, dr2) = @time msd(result)
 
 #using JLD
 #save("D:/water $Nactual molecules $timesteps steps $time_now.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2)
@@ -91,3 +64,27 @@ timesteps = round(length(result.solution.t))
 #plot(rs/1e-10, grf, xlim=[0, 0.4999L/1e-10], label=["Radial distribution function"],ylabel="g(r)", xlabel="r, Å")
 
 #@time animate(result, "D:/$Nactual H20 particles with $timesteps timesteps $time_now.gif")
+
+#=
+(rs, grf) = rdf(result)
+(ts, dr2) = msd(result)
+
+t = t1:τ:result.solution.t[end-1]
+temper = temperature.(result, t)
+
+time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
+Nactual = length(bodies)
+timesteps = round(length(result.solution.t))
+
+using MAT
+file = matopen("d:/water real units without thermostat and $timesteps timesteps_$time_now.mat", "w")
+write(file, "t", collect(t))
+write(file, "temper", temper)
+write(file, "rs", rs)
+write(file, "ts", ts)
+write(file, "grf", grf)
+write(file, "dr2", dr2)
+close(file)
+
+@time NBodySimulator.save_to_pdb(result, "d:/water in real units from PDB $timesteps timesteps $time_now.pdb" )
+=#