From ab0f0003a481e1ce69b268efc0865a002632d454 Mon Sep 17 00:00:00 2001 From: Sikorski Date: Mon, 9 Oct 2023 16:28:19 +0200 Subject: [PATCH] move OverdampedLangevin, add export/test, remove potential --- src/Molly.jl | 1 - src/interactions/potential.jl | 22 ------ src/simulators.jl | 109 ++++++++++++++------------ test/simulation.jl | 143 +++++++++++++++++----------------- 4 files changed, 130 insertions(+), 145 deletions(-) delete mode 100644 src/interactions/potential.jl diff --git a/src/Molly.jl b/src/Molly.jl index 17a963df..46c0a4f4 100644 --- a/src/Molly.jl +++ b/src/Molly.jl @@ -59,7 +59,6 @@ include("interactions/periodic_torsion.jl") include("interactions/rb_torsion.jl") include("interactions/implicit_solvent.jl") include("interactions/muller_brown.jl") -include("interactions/potential.jl") include("energy.jl") include("constraints.jl") include("simulators.jl") diff --git a/src/interactions/potential.jl b/src/interactions/potential.jl deleted file mode 100644 index 9501bd25..00000000 --- a/src/interactions/potential.jl +++ /dev/null @@ -1,22 +0,0 @@ -struct PotentialFunction{F} - f::F - # fwdiff gradient caches -end - -function Molly.potential_energy(V::PotentialFunction, - sys, - neighbors=nothing; - n_threads=Threads.nthreads()) - - return V.f(sys.coords) -end - -function Molly.forces(V::PotentialFunction, - sys, - neighbors=nothing; - n_threads=Threads.nthreads()) - - return - Zygote.gradient(V.f, sys.coords) -end - - diff --git a/src/simulators.jl b/src/simulators.jl index dc4eedf5..670a8847 100644 --- a/src/simulators.jl +++ b/src/simulators.jl @@ -8,6 +8,7 @@ export StormerVerlet, Langevin, LangevinSplitting, + OverdampedLangevin, NoseHoover, TemperatureREMD, remd_exchange!, @@ -486,6 +487,63 @@ function B_step!(sys, dt_eff, acceleration_vector, compute_forces::Bool, return sys end +""" + OverdampedLangevin(; ) + +Simulates the overdamped Langevin equation using the Euler-Maruyama method. + +# Arguments +- `dt::S`: the time step of the simulation. +- `temperature::K`: the equilibrium temperature of the simulation. +- `friction::F`: the friction coefficient of the simulation. +- `remove_CM_motion::Bool=true`: whether to remove the center of mass motion + every time step. +""" +@kwdef struct OverdampedLangevin{S,K,F} + dt::S + temperature::K + friction::F + remove_CM_motion::Bool = true +end + +function simulate!(sys, + sim::OverdampedLangevin, + n_steps::Integer; + n_threads::Integer=Threads.nthreads(), + rng=Random.GLOBAL_RNG) + + sys.coords = wrap_coords.(sys.coords, (sys.boundary,)) + + sim.remove_CM_motion && remove_CM_motion!(sys) + neighbors = find_neighbors(sys, sys.neighbor_finder; n_threads=n_threads) + run_loggers!(sys, neighbors, 0; n_threads=n_threads) + + for step_n in 1:n_steps + a = accelerations(sys, neighbors; n_threads=n_threads) + # noise = kB * T / M + noise = random_velocities(sys, sim.temperature; rng=rng) + if !all(isfinite, reinterpret(Float64, a)) + @show step_n + @show reinterpret(Float64, a) + + end + sys.coords += a ./ sim.friction .* sim.dt .+ sqrt(2 / sim.friction * sim.dt) .* noise + + @assert all(isfinite, reinterpret(Float64, sys.coords)) + + sys.coords = wrap_coords.(sys.coords, (sys.boundary,)) + sim.remove_CM_motion && remove_CM_motion!(sys) + + run_loggers!(sys, neighbors, step_n; n_threads=n_threads) + + if step_n != n_steps + neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n; + n_threads=n_threads) + end + end + return sys +end + """ NoseHoover(; ) @@ -930,54 +988,3 @@ function random_unit_vector(float_type, dims) vec = randn(float_type, dims) return vec / norm(vec) end - - -""" - OverdampedLangevin(; ) - -Simulates the overdamped Langevin equation using the Euler-Maruyama method. - -# Arguments -- `dt::S`: the time step of the simulation. -- `temperature::K`: the equilibrium temperature of the simulation. -- `friction::F`: the friction coefficient of the simulation. -- `remove_CM_motion::Bool=true`: whether to remove the center of mass motion - every time step. -""" -@kwdef struct OverdampedLangevin{T,K,F} - dt::T - temperature::K - friction::F - remove_CM_motion::Bool = true -end - -function simulate!(sys, - sim::OverdampedLangevin, - n_steps::Integer; - n_threads::Integer=Threads.nthreads(), - rng=Random.GLOBAL_RNG) - - sys.coords = wrap_coords.(sys.coords, (sys.boundary,)) - - sim.remove_CM_motion && remove_CM_motion!(sys) - neighbors = find_neighbors(sys, sys.neighbor_finder; n_threads=n_threads) - run_loggers!(sys, neighbors, 0; n_threads=n_threads) - - for step_n in 1:n_steps - a = accelerations(sys, neighbors; n_threads=n_threads) - # noise = kB * T / M - noise = random_velocities(sys, sim.temperature; rng=rng) - sys.coords += a ./ sim.friction .* sim.dt .+ sqrt(2 / sim.friction * sim.dt) .* noise - - sys.coords = wrap_coords.(sys.coords, (sys.boundary,)) - sim.remove_CM_motion && remove_CM_motion!(sys) - - run_loggers!(sys, neighbors, step_n; n_threads=n_threads) - - if step_n != n_steps - neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n; - n_threads=n_threads) - end - end - return sys -end diff --git a/test/simulation.jl b/test/simulation.jl index 4a40e0b2..d15a4f90 100644 --- a/test/simulation.jl +++ b/test/simulation.jl @@ -21,7 +21,7 @@ coords=CoordinateLogger(100; dims=2), gen_temp=GeneralObservableLogger(gen_temp_wrapper, typeof(temp), 10), avg_temp=AverageObservableLogger(Molly.temperature_wrapper, - typeof(temp), 1; n_blocks=200), + typeof(temp), 1; n_blocks=200), ), ) random_velocities!(s, temp) @@ -207,6 +207,7 @@ end Verlet(dt=0.002u"ps", coupling=AndersenThermostat(temp, 10.0u"ps")), StormerVerlet(dt=0.002u"ps"), Langevin(dt=0.002u"ps", temperature=temp, friction=1.0u"ps^-1"), + OverdampedLangevin(dt=0.002u"ps", temperature=temp, friction=10.0u"ps^-1"), ] s = System( @@ -238,15 +239,15 @@ end push!(coords, coords[i] .+ [0.1, 0.0, 0.0]u"nm") end bonds = InteractionList2Atoms( - collect(1:(n_atoms ÷ 2)), - collect((1 + n_atoms ÷ 2):n_atoms), - [HarmonicBond(k=300_000.0u"kJ * mol^-1 * nm^-2", r0=0.1u"nm") for i in 1:(n_atoms ÷ 2)], + collect(1:(n_atoms÷2)), + collect((1+n_atoms÷2):n_atoms), + [HarmonicBond(k=300_000.0u"kJ * mol^-1 * nm^-2", r0=0.1u"nm") for i in 1:(n_atoms÷2)], fill("", n_atoms ÷ 2), ) eligible = trues(n_atoms, n_atoms) - for i in 1:(n_atoms ÷ 2) - eligible[i, i + (n_atoms ÷ 2)] = false - eligible[i + (n_atoms ÷ 2), i] = false + for i in 1:(n_atoms÷2) + eligible[i, i+(n_atoms÷2)] = false + eligible[i+(n_atoms÷2), i] = false end simulator = VelocityVerlet(dt=0.002u"ps", coupling=BerendsenThermostat(temp, 1.0u"ps")) @@ -275,8 +276,8 @@ end s.loggers.coords, boundary, temp_fp_viz; - connections=[(i, i + (n_atoms ÷ 2)) for i in 1:(n_atoms ÷ 2)], - trails=2, + connections=[(i, i + (n_atoms ÷ 2)) for i in 1:(n_atoms÷2)], + trails=2 ) end end @@ -305,14 +306,14 @@ end @testset "$inter" for inter in pairwise_inter_types if use_neighbors(inter) neighbor_finder = DistanceNeighborFinder(eligible=trues(n_atoms, n_atoms), n_steps=10, - dist_cutoff=1.5u"nm") + dist_cutoff=1.5u"nm") else neighbor_finder = NoNeighborFinder() end s = System( atoms=[Atom(charge=i % 2 == 0 ? -1.0 : 1.0, mass=10.0u"g/mol", σ=0.2u"nm", - ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms], + ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms], coords=place_atoms(n_atoms, boundary; min_dist=0.2u"nm"), boundary=boundary, velocities=[random_velocity(10.0u"g/mol", temp) .* 0.01 for i in 1:n_atoms], @@ -332,7 +333,7 @@ end @testset "Müller-Brown" begin atom_mass = 1.0u"g/mol" atoms = [Atom(mass=atom_mass, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1")] - boundary = RectangularBoundary(Inf*u"nm") + boundary = RectangularBoundary(Inf * u"nm") coords = [SVector(-0.5, 0.25)u"nm"] temp = 100.0u"K" velocities = [random_velocity(atom_mass, temp; dims=2)] @@ -392,9 +393,9 @@ end atoms=[Atom(charge=0.0, mass=10.0, σ=0.3, ϵ=0.2) for i in 1:n_atoms], coords=ustrip_vec.(coords), boundary=CubicBoundary(ustrip.(boundary)), - velocities=ustrip_vec.(u"nm/ps",velocities), + velocities=ustrip_vec.(u"nm/ps", velocities), pairwise_inters=(LennardJones(use_neighbors=true, force_units=NoUnits, - energy_units=NoUnits),), + energy_units=NoUnits),), neighbor_finder=DistanceNeighborFinder( eligible=trues(n_atoms, n_atoms), n_steps=10, @@ -449,13 +450,13 @@ end atom_selector(at, at_data) = at_data.atom_type == "A1" sys_res = add_position_restraints(sys, 100_000.0u"kJ * mol^-1 * nm^-2"; - atom_selector=atom_selector) + atom_selector=atom_selector) @time simulate!(sys_res, sim, n_steps) dists = norm.(vector.(starting_coords, Array(sys_res.coords), (boundary,))) @test maximum(dists[1:n_atoms_res]) < 0.1u"nm" - @test median(dists[(n_atoms_res + 1):end]) > 0.2u"nm" + @test median(dists[(n_atoms_res+1):end]) > 0.2u"nm" end end @@ -475,15 +476,15 @@ end velocities = [random_velocity(atom_mass, temp) for i in 1:n_atoms] eligible = trues(n_atoms, n_atoms) - for i in 1:(n_atoms ÷ 2) - eligible[i, i + (n_atoms ÷ 2)] = false - eligible[i + (n_atoms ÷ 2), i] = false + for i in 1:(n_atoms÷2) + eligible[i, i+(n_atoms÷2)] = false + eligible[i+(n_atoms÷2), i] = false end neighbor_finder = DistanceNeighborFinder(eligible=eligible, n_steps=10, dist_cutoff=1.5u"nm") - bond_lengths = [0.1u"nm" for i in 1:(n_atoms ÷ 2)] - sh = SHAKE(bond_lengths, collect(1:(n_atoms ÷ 2)), collect((1 + (n_atoms ÷ 2)):n_atoms)) + bond_lengths = [0.1u"nm" for i in 1:(n_atoms÷2)] + sh = SHAKE(bond_lengths, collect(1:(n_atoms÷2)), collect((1+(n_atoms÷2)):n_atoms)) constraints = (sh,) sys = System( @@ -498,7 +499,7 @@ end ) for i in eachindex(sys.coords) - sys.coords[i] += [rand()*0.01, rand()*0.01, rand()*0.01]u"nm" + sys.coords[i] += [rand() * 0.01, rand() * 0.01, rand() * 0.01]u"nm" end old_coords = sys.coords @@ -529,11 +530,11 @@ end coords = place_atoms(n_atoms ÷ 3, boundary, min_dist=0.3u"nm") - for i in 1:(n_atoms ÷ 3) + for i in 1:(n_atoms÷3) push!(coords, coords[i] .+ [0.13, 0.0, 0.0]u"nm") end - for i in 1:(n_atoms ÷ 3) + for i in 1:(n_atoms÷3) push!(coords, coords[i] .+ [0.26, 0.0, 0.0]u"nm") end @@ -541,20 +542,20 @@ end velocities = [random_velocity(atom_mass, temp) for i in 1:n_atoms] eligible = trues(n_atoms, n_atoms) - for i in 1:(n_atoms ÷ 3) - eligible[i, i + (n_atoms ÷ 3)] = false - eligible[i + (n_atoms ÷ 3), i] = false - eligible[i + (n_atoms ÷ 3), i + 2 * (n_atoms ÷ 3)] = false - eligible[i + 2 * (n_atoms ÷ 3), i + (n_atoms ÷ 3)] = false + for i in 1:(n_atoms÷3) + eligible[i, i+(n_atoms÷3)] = false + eligible[i+(n_atoms÷3), i] = false + eligible[i+(n_atoms÷3), i+2*(n_atoms÷3)] = false + eligible[i+2*(n_atoms÷3), i+(n_atoms÷3)] = false end neighbor_finder = DistanceNeighborFinder(eligible=eligible, n_steps=10, dist_cutoff=1.5u"nm") - bond_lengths = [0.1u"nm" for i in 1:(2 * (n_atoms ÷ 3))] + bond_lengths = [0.1u"nm" for i in 1:(2*(n_atoms÷3))] sh = SHAKE( bond_lengths, - collect(1:(2 * (n_atoms ÷ 3))), - collect(((n_atoms ÷ 3) + 1):n_atoms), + collect(1:(2*(n_atoms÷3))), + collect(((n_atoms÷3)+1):n_atoms), ) constraints = (sh,) @@ -570,7 +571,7 @@ end ) for i in eachindex(sys.coords) - sys.coords[i] += [rand()*0.01, rand()*0.01, rand()*0.01]u"nm" + sys.coords[i] += [rand() * 0.01, rand() * 0.01, rand() * 0.01]u"nm" end old_coords = sys.coords @@ -608,13 +609,13 @@ end rseed = 2022 simulator1 = Langevin(dt=0.002u"ps", temperature=temp, friction=1.0u"ps^-1") simulator2 = LangevinSplitting(dt=0.002u"ps", temperature=temp, - friction=10.0u"g * mol^-1 * ps^-1", splitting="BAOA") + friction=10.0u"g * mol^-1 * ps^-1", splitting="BAOA") @time simulate!(s1, simulator1, n_steps; rng=MersenneTwister(rseed)) - @test 280.0u"K" <= mean(s1.loggers.temp.history[(end - 100):end]) <= 320.0u"K" + @test 280.0u"K" <= mean(s1.loggers.temp.history[(end-100):end]) <= 320.0u"K" @time simulate!(s2, simulator2, n_steps; rng=MersenneTwister(rseed)) - @test 280.0u"K" <= mean(s2.loggers.temp.history[(end - 100):end]) <= 320.0u"K" + @test 280.0u"K" <= mean(s2.loggers.temp.history[(end-100):end]) <= 320.0u"K" @test maximum(maximum(abs.(v)) for v in (s1.coords .- s2.coords)) < 1e-5u"nm" end @@ -708,7 +709,7 @@ end exchange_time=2.5u"ps", ) - @time simulate!(repsys, simulator, n_steps; assign_velocities=true ) + @time simulate!(repsys, simulator, n_steps; assign_velocities=true) @time simulate!(repsys, simulator, n_steps; assign_velocities=false) efficiency = repsys.exchange_logger.n_exchanges / repsys.exchange_logger.n_attempts @@ -744,7 +745,7 @@ end replica_pairwise_inters = [(LennardJonesSoftCore(α=1, λ=λ_vals[i], p=2, use_neighbors=true),) for i in 1:n_replicas] - replica_loggers = [(temp=TemperatureLogger(10), ) for i in 1:n_replicas] + replica_loggers = [(temp=TemperatureLogger(10),) for i in 1:n_replicas] repsys = ReplicaSystem( atoms=atoms, @@ -770,7 +771,7 @@ end exchange_time=2.5u"ps", ) - @time simulate!(repsys, simulator, n_steps; assign_velocities=true ) + @time simulate!(repsys, simulator, n_steps; assign_velocities=true) @time simulate!(repsys, simulator, n_steps; assign_velocities=false) efficiency = repsys.exchange_logger.n_exchanges / repsys.exchange_logger.n_attempts @@ -802,10 +803,10 @@ end atoms=atoms, coords=coords, boundary=boundary, - pairwise_inters=(Coulomb(), ), + pairwise_inters=(Coulomb(),), neighbor_finder=neighbor_finder, loggers=( - coords=CoordinateLogger(10), + coords=CoordinateLogger(10), mcl=MonteCarloLogger(), avgpe=AverageObservableLogger(potential_energy, typeof(atoms[1].ϵ), 10), ), @@ -823,7 +824,7 @@ end trial_args=Dict(:shift_size => 0.1u"nm"), ) - @time simulate!(sys, simulator_uniform , n_steps) + @time simulate!(sys, simulator_uniform, n_steps) @time simulate!(sys, simulator_gaussian, n_steps) acceptance_rate = sys.loggers.mcl.n_accept / sys.loggers.mcl.n_trials @@ -837,7 +838,7 @@ end distance_sum = 0.0u"nm" for i in eachindex(sys) ci = sys.coords[i] - min_dist2 = Inf*u"nm^2" + min_dist2 = Inf * u"nm^2" for j in eachindex(sys) if i == j continue @@ -892,12 +893,12 @@ end @time simulate!(sys, lang, n_steps; n_threads=1) @test 280.0u"K" < mean(values(sys.loggers.temperature)) < 300.0u"K" - @test 50.0u"kJ * mol^-1" < mean(values(sys.loggers.total_energy )) < 120.0u"kJ * mol^-1" + @test 50.0u"kJ * mol^-1" < mean(values(sys.loggers.total_energy)) < 120.0u"kJ * mol^-1" @test 50.0u"kJ * mol^-1" < mean(values(sys.loggers.kinetic_energy)) < 120.0u"kJ * mol^-1" @test mean(values(sys.loggers.potential_energy)) < 0.0u"kJ * mol^-1" @test -5.0u"kJ * mol^-1" < mean(values(sys.loggers.virial)) < 5.0u"kJ * mol^-1" @test 1.7u"bar" < mean(values(sys.loggers.pressure)) < 2.2u"bar" - @test 0.1u"bar" < std( values(sys.loggers.pressure)) < 0.5u"bar" + @test 0.1u"bar" < std(values(sys.loggers.pressure)) < 0.5u"bar" @test all(values(sys.loggers.box_size) .== 8.0u"nm") @test sys.boundary == CubicBoundary(8.0u"nm") @@ -932,14 +933,14 @@ end @time simulate!(sys, sim, n_steps; n_threads=1) @test 260.0u"K" < mean(values(sys.loggers.temperature)) < 300.0u"K" - @test 50.0u"kJ * mol^-1" < mean(values(sys.loggers.total_energy )) < 120.0u"kJ * mol^-1" + @test 50.0u"kJ * mol^-1" < mean(values(sys.loggers.total_energy)) < 120.0u"kJ * mol^-1" @test 50.0u"kJ * mol^-1" < mean(values(sys.loggers.kinetic_energy)) < 120.0u"kJ * mol^-1" @test mean(values(sys.loggers.potential_energy)) < 0.0u"kJ * mol^-1" @test -5.0u"kJ * mol^-1" < mean(values(sys.loggers.virial)) < 5.0u"kJ * mol^-1" @test 0.8u"bar" < mean(values(sys.loggers.pressure)) < 1.2u"bar" - @test 0.1u"bar" < std( values(sys.loggers.pressure)) < 0.5u"bar" + @test 0.1u"bar" < std(values(sys.loggers.pressure)) < 0.5u"bar" @test 9.5u"nm" < mean(values(sys.loggers.box_size)) < 10.5u"nm" - @test 0.2u"nm" < std( values(sys.loggers.box_size)) < 1.0u"nm" + @test 0.2u"nm" < std(values(sys.loggers.box_size)) < 1.0u"nm" @test sys.boundary != CubicBoundary(8.0u"nm") end end @@ -968,8 +969,8 @@ end pressure_test_set = ( SVector(1.0u"bar", 1.0u"bar", 1.0u"bar"), # XYZ-axes coupled with the same pressure value SVector(1.5u"bar", 0.5u"bar", 1.0u"bar"), # XYZ-axes coupled with different pressure values - SVector(nothing , 1.0u"bar", nothing ), # Only Y-axis coupled - SVector(nothing , nothing , nothing ), # Uncoupled + SVector(nothing, 1.0u"bar", nothing), # Only Y-axis coupled + SVector(nothing, nothing, nothing), # Uncoupled ) for gpu in gpu_list @@ -1081,7 +1082,7 @@ end if !barostat.constant_volume && isnothing(barostat.pressure[3]) @test sys.boundary[3] == 8.0u"nm" @test 0.8u"bar" < mean(values(sys.loggers.pressure)) < 1.2u"bar" - @test 0.1u"bar" < std(values(sys.loggers.pressure)) < 0.5u"bar" + @test 0.1u"bar" < std(values(sys.loggers.pressure)) < 0.5u"bar" end end end @@ -1108,14 +1109,14 @@ end ),), loggers=(tot_eng=TotalEnergyLogger(100),), energy_units=u"kJ * mol^-1", - force_units=u"kJ * mol^-1 * nm^-1", + force_units=u"kJ * mol^-1 * nm^-1" ) sys_cp = System(sys) - @test sys_cp.atoms == sys.atoms + @test sys_cp.atoms == sys.atoms @test sys_cp.coords == sys.coords sys_mod = System(sys; coords=(sys.coords .* 0.5)) - @test sys_mod.atoms == sys.atoms + @test sys_mod.atoms == sys.atoms @test sys_mod.coords == sys.coords .* 0.5 σ = 0.34u"nm" @@ -1124,7 +1125,7 @@ end for i in eachindex(sys) push!(updated_atoms, Atom(index=sys.atoms[i].index, charge=sys.atoms[i].charge, - mass=sys.atoms[i].mass, σ=σ, ϵ=ϵ, solute=sys.atoms[i].solute)) + mass=sys.atoms[i].mass, σ=σ, ϵ=ϵ, solute=sys.atoms[i].solute)) end sys = System(sys; atoms=[updated_atoms...]) @@ -1137,7 +1138,7 @@ end @time simulate!(sys, simulator, 25_000; run_loggers=false) @time simulate!(sys, simulator, 25_000) - + @test length(values(sys.loggers.tot_eng)) == 251 @test -1850u"kJ * mol^-1" < mean(values(sys.loggers.tot_eng)) < -1650u"kJ * mol^-1" @@ -1177,7 +1178,7 @@ end simulator = VelocityVerlet(dt=f32 ? 0.02f0u"ps" : 0.02u"ps") k = f32 ? 10_000.0f0u"kJ * mol^-1 * nm^-2" : 10_000.0u"kJ * mol^-1 * nm^-2" r0 = f32 ? 0.2f0u"nm" : 0.2u"nm" - bonds = [HarmonicBond(k=k, r0=r0) for i in 1:(n_atoms ÷ 2)] + bonds = [HarmonicBond(k=k, r0=r0) for i in 1:(n_atoms÷2)] specific_inter_lists = (InteractionList2Atoms( gpu ? CuArray(Int32.(collect(1:2:n_atoms))) : Int32.(collect(1:2:n_atoms)), gpu ? CuArray(Int32.(collect(2:2:n_atoms))) : Int32.(collect(2:2:n_atoms)), @@ -1201,12 +1202,12 @@ end coords = CuArray(deepcopy(f32 ? starting_coords_f32 : starting_coords)) velocities = CuArray(deepcopy(f32 ? starting_velocities_f32 : starting_velocities)) atoms = CuArray([Atom(charge=f32 ? 0.0f0 : 0.0, mass=atom_mass, σ=f32 ? 0.2f0u"nm" : 0.2u"nm", - ϵ=f32 ? 0.2f0u"kJ * mol^-1" : 0.2u"kJ * mol^-1") for i in 1:n_atoms]) + ϵ=f32 ? 0.2f0u"kJ * mol^-1" : 0.2u"kJ * mol^-1") for i in 1:n_atoms]) else coords = deepcopy(f32 ? starting_coords_f32 : starting_coords) velocities = deepcopy(f32 ? starting_velocities_f32 : starting_velocities) atoms = [Atom(charge=f32 ? 0.0f0 : 0.0, mass=atom_mass, σ=f32 ? 0.2f0u"nm" : 0.2u"nm", - ϵ=f32 ? 0.2f0u"kJ * mol^-1" : 0.2u"kJ * mol^-1") for i in 1:n_atoms] + ϵ=f32 ? 0.2f0u"kJ * mol^-1" : 0.2u"kJ * mol^-1") for i in 1:n_atoms] end s = System( @@ -1231,22 +1232,22 @@ end end runs = [ - ("CPU" , [false, false, false, false]), - ("CPU f32" , [false, false, true , false]), - ("CPU NL" , [true , false, false, false]), - ("CPU f32 NL", [true , false, true , false]), + ("CPU", [false, false, false, false]), + ("CPU f32", [false, false, true, false]), + ("CPU NL", [true, false, false, false]), + ("CPU f32 NL", [true, false, true, false]), ] if run_parallel_tests - push!(runs, ("CPU parallel" , [false, true , false, false])) - push!(runs, ("CPU parallel f32" , [false, true , true , false])) - push!(runs, ("CPU parallel NL" , [true , true , false, false])) - push!(runs, ("CPU parallel f32 NL", [true , true , true , false])) + push!(runs, ("CPU parallel", [false, true, false, false])) + push!(runs, ("CPU parallel f32", [false, true, true, false])) + push!(runs, ("CPU parallel NL", [true, true, false, false])) + push!(runs, ("CPU parallel f32 NL", [true, true, true, false])) end if run_gpu_tests - push!(runs, ("GPU" , [false, false, false, true])) - push!(runs, ("GPU f32" , [false, false, true , true])) - push!(runs, ("GPU NL" , [true , false, false, true])) - push!(runs, ("GPU f32 NL", [true , false, true , true])) + push!(runs, ("GPU", [false, false, false, true])) + push!(runs, ("GPU f32", [false, false, true, true])) + push!(runs, ("GPU NL", [true, false, false, true])) + push!(runs, ("GPU f32 NL", [true, false, true, true])) end final_coords_ref, E_start_ref = test_sim(runs[1][2]...)