jgreener64
commented
2 years ago Thanks for doing this, updating the neighbour lists as required is a definite algorithmic improvement.
About the failing Float32 (not Int32) tests, adding the type to line 192 of neighbors.jl seems to fix it: Box([10*one(T),10*one(T),10*one(T)],one(T); T=T). Although this box is replaced, I guess later values are promoted to its type. Tests then pass up to the error shown below.
I ran the following benchmark script, which should be portable:
using Molly
using DelimitedFiles
data_dir = normpath(dirname(pathof(Molly)), "..", "data")
ff_dir = joinpath(data_dir, "force_fields")
openmm_dir = joinpath(data_dir, "openmm_6mrr")
ff = OpenMMForceField(joinpath.(ff_dir, ["ff99SBildn.xml", "tip3p_standard.xml", "his.xml"])...)
atoms, atoms_data, specific_inter_lists, general_inters, neighbor_finder, coords, box_size = setupsystem(
joinpath(data_dir, "6mrr_equil.pdb"), ff)
n_steps = 500
timestep = 0.0005u"ps"
temp = 300.0u"K"
velocities = SVector{3}.(eachrow(readdlm(joinpath(openmm_dir, "velocities_300K.txt"))))u"nm * ps^-1"
s = Simulation(
simulator=VelocityVerlet(),
atoms=atoms,
atoms_data=atoms_data,
specific_inter_lists=specific_inter_lists,
general_inters=general_inters,
coords=coords,
velocities=velocities,
temperature=temp,
box_size=box_size,
neighbor_finder=neighbor_finder,
thermostat=AndersenThermostat(1.0u"ps"),
timestep=timestep,
n_steps=n_steps,
)
parallel = false
simulate!(s, 5; parallel=parallel)
@time simulate!(s, n_steps; parallel=parallel)
parallel = true
simulate!(s, 5; parallel=parallel)
@time simulate!(s, n_steps; parallel=parallel)With the current Molly I get 120 s on 1 thread and 28 s on 16 threads. With the new change I get 18 mins on 1 thread and an error on 16 threads:
ERROR: LoadError: TaskFailedException
Stacktrace:
[1] wait
@ ./task.jl:322 [inlined]
[2] threading_run(func::Function)
@ Base.Threads ./threadingconstructs.jl:38
[3] macro expansion
@ ./threadingconstructs.jl:97 [inlined]
[4] UpdateCellList!(x::Vector{SVector{3, Float64}}, box::CellListMap.Box{CellListMap.OrthorhombicCell, 3, Float64, 9}, cl::CellListMap.CellList{3, Float64}, aux::CellListMap.AuxThreaded{3, Float64}; parallel::Bool)
@ CellListMap ~/.julia/packages/CellListMap/dfbQh/src/CellLists.jl:617
[5] find_neighbors!(s::Simulation{false, Quantity{Float64, π―, Unitful.FreeUnits{(K,), π―, nothing}}, Vector{Atom{Quantity{Float64, π π, Unitful.FreeUnits{(C,), π π, nothing}}, Quantity{Float64, π, Unitful.FreeUnits{(u,), π, nothing}}, Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}, Vector{AtomData}, Vector{SVector{3, Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}}}, Vector{SVector{3, Quantity{Float64, π π^-1, Unitful.FreeUnits{(nm, ps^-1), π π^-1, nothing}}}}, Tuple{LennardJones{false, DistanceCutoff{Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Quantity{Float64, π^2, Unitful.FreeUnits{(nm^2,), π^2, nothing}}, Quantity{Float64, π^-2, Unitful.FreeUnits{(nm^-2,), π^-2, nothing}}}, Float64, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), π π π^-1 π^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}, CoulombReactionField{Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Float64, Float64, Quantity{Float64, π^3 π π^-1 π^-2 π^-4, Unitful.FreeUnits{(C^-2, kJ, nm, mol^-1), π^3 π π^-1 π^-2 π^-4, nothing}}, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), π π π^-1 π^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}, Quantity{Float64, π^2, Unitful.FreeUnits{(nm^2,), π^2, nothing}}, Quantity{Float64, π^-3, Unitful.FreeUnits{(nm^-3,), π^-3, nothing}}, Quantity{Float64, π^-1, Unitful.FreeUnits{(nm^-1,), π^-1, nothing}}}}, Tuple{Vector{HarmonicBond{Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Quantity{Float64, π π^-1 π^-2, Unitful.FreeUnits{(kJ, nm^-2, mol^-1), π π^-1 π^-2, nothing}}}}, Vector{HarmonicAngle{Float64, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}, Vector{PeriodicTorsion{Float64, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}, Vector{PeriodicTorsion{Float64, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}}, SVector{3, Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}}, Quantity{Float64, π, Unitful.FreeUnits{(ps,), π, nothing}}, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), π π π^-1 π^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}, nf::CellListMapNeighborFinder{Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, 3, Float64}, step_n::Int64; parallel::Bool)
@ Molly ~/.julia/dev/Molly/src/neighbors.jl:257
[6] simulate!(s::Simulation{false, Quantity{Float64, π―, Unitful.FreeUnits{(K,), π―, nothing}}, Vector{Atom{Quantity{Float64, π π, Unitful.FreeUnits{(C,), π π, nothing}}, Quantity{Float64, π, Unitful.FreeUnits{(u,), π, nothing}}, Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}, Vector{AtomData}, Vector{SVector{3, Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}}}, Vector{SVector{3, Quantity{Float64, π π^-1, Unitful.FreeUnits{(nm, ps^-1), π π^-1, nothing}}}}, Tuple{LennardJones{false, DistanceCutoff{Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Quantity{Float64, π^2, Unitful.FreeUnits{(nm^2,), π^2, nothing}}, Quantity{Float64, π^-2, Unitful.FreeUnits{(nm^-2,), π^-2, nothing}}}, Float64, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), π π π^-1 π^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}, CoulombReactionField{Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Float64, Float64, Quantity{Float64, π^3 π π^-1 π^-2 π^-4, Unitful.FreeUnits{(C^-2, kJ, nm, mol^-1), π^3 π π^-1 π^-2 π^-4, nothing}}, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), π π π^-1 π^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}, Quantity{Float64, π^2, Unitful.FreeUnits{(nm^2,), π^2, nothing}}, Quantity{Float64, π^-3, Unitful.FreeUnits{(nm^-3,), π^-3, nothing}}, Quantity{Float64, π^-1, Unitful.FreeUnits{(nm^-1,), π^-1, nothing}}}}, Tuple{Vector{HarmonicBond{Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Quantity{Float64, π π^-1 π^-2, Unitful.FreeUnits{(kJ, nm^-2, mol^-1), π π^-1 π^-2, nothing}}}}, Vector{HarmonicAngle{Float64, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}, Vector{PeriodicTorsion{Float64, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}, Vector{PeriodicTorsion{Float64, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}}, SVector{3, Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}}, Quantity{Float64, π, Unitful.FreeUnits{(ps,), π, nothing}}, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), π π π^-1 π^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}, ::VelocityVerlet, n_steps::Int64; parallel::Bool)
@ Molly ~/.julia/dev/Molly/src/simulators.jl:41
[7] simulate!(s::Simulation{false, Quantity{Float64, π―, Unitful.FreeUnits{(K,), π―, nothing}}, Vector{Atom{Quantity{Float64, π π, Unitful.FreeUnits{(C,), π π, nothing}}, Quantity{Float64, π, Unitful.FreeUnits{(u,), π, nothing}}, Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}, Vector{AtomData}, Vector{SVector{3, Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}}}, Vector{SVector{3, Quantity{Float64, π π^-1, Unitful.FreeUnits{(nm, ps^-1), π π^-1, nothing}}}}, Tuple{LennardJones{false, DistanceCutoff{Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Quantity{Float64, π^2, Unitful.FreeUnits{(nm^2,), π^2, nothing}}, Quantity{Float64, π^-2, Unitful.FreeUnits{(nm^-2,), π^-2, nothing}}}, Float64, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), π π π^-1 π^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}, CoulombReactionField{Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Float64, Float64, Quantity{Float64, π^3 π π^-1 π^-2 π^-4, Unitful.FreeUnits{(C^-2, kJ, nm, mol^-1), π^3 π π^-1 π^-2 π^-4, nothing}}, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), π π π^-1 π^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}, Quantity{Float64, π^2, Unitful.FreeUnits{(nm^2,), π^2, nothing}}, Quantity{Float64, π^-3, Unitful.FreeUnits{(nm^-3,), π^-3, nothing}}, Quantity{Float64, π^-1, Unitful.FreeUnits{(nm^-1,), π^-1, nothing}}}}, Tuple{Vector{HarmonicBond{Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}, Quantity{Float64, π π^-1 π^-2, Unitful.FreeUnits{(kJ, nm^-2, mol^-1), π π^-1 π^-2, nothing}}}}, Vector{HarmonicAngle{Float64, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}, Vector{PeriodicTorsion{Float64, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}, Vector{PeriodicTorsion{Float64, Quantity{Float64, π^2 π π^-1 π^-2, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}}}}, SVector{3, Quantity{Float64, π, Unitful.FreeUnits{(nm,), π, nothing}}}, Quantity{Float64, π, Unitful.FreeUnits{(ps,), π, nothing}}, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), π π π^-1 π^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), π^2 π π^-1 π^-2, nothing}}, n_steps::Int64; parallel::Bool)
@ Molly ~/.julia/dev/Molly/src/simulators.jl:153
[8] top-level scope
@ ./timing.jl:210
nested task error: BoundsError: attempt to access 1-element Vector{CellListMap.CellList{3, Float64}} at index [2]
Stacktrace:
[1] getindex
@ ./array.jl:835 [inlined]
[2] macro expansion
@ ~/.julia/packages/CellListMap/dfbQh/src/CellLists.jl:618 [inlined]
[3] (::CellListMap.var"#45#threadsfor_fun#59"{Vector{SVector{3, Float64}}, CellListMap.Box{CellListMap.OrthorhombicCell, 3, Float64, 9}, CellListMap.AuxThreaded{3, Float64}, UnitRange{Int64}})(onethread::Bool)
@ CellListMap ./threadingconstructs.jl:85
[4] (::CellListMap.var"#45#threadsfor_fun#59"{Vector{SVector{3, Float64}}, CellListMap.Box{CellListMap.OrthorhombicCell, 3, Float64, 9}, CellListMap.AuxThreaded{3, Float64}, UnitRange{Int64}})()
@ CellListMap ./threadingconstructs.jl:52
in expression starting at /home/jgreener/.julia/dev/Molly/bench_cnl.jl:42Any idea what is going on there?
lmiq
codecov[bot]
@jgreener64
I am just opening the pull request so that the changes are easier to visualize. Needs polishing and fixing of the Int type for the tests to pass.
I have changed a little bit the implementation of the neighborlist structure to update in place the lists and cell lists from step to step in the simulations. This should bring some performance benefits (and much less memory usage).
Most tests are passing, but some are failing when you change from Int64 to Int32 integer types (I'm unsure where is this being defined, and it does not seem to be ready to be propagated to the type of neighbour list.
Otherwise the tests seem to be passing. If you want to take a look at this, and see if in the cases where it runs (Int64) the performance is improved, I can polish what is missing to make everything to work.
The current approach (in which NeighbourList is a structure that contains both the number of pairs and an array of pairs that does not need to be replaced each time) can be used for the other approaches of computing neighbor lists, with the same advantages.