From f4962deaecb87f89c32ebba27ba128f6ac432430 Mon Sep 17 00:00:00 2001 From: Joroks Date: Wed, 4 Sep 2024 23:32:52 +0200 Subject: [PATCH 1/9] implement neighbor list access --- src/LAMMPS.jl | 49 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index 820e294..2e5a658 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -5,6 +5,7 @@ include("api.jl") export LMP, command, create_atoms, get_natoms, extract_atom, extract_compute, extract_global, extract_setting, gather, gather_bonds, gather_angles, gather_dihedrals, gather_impropers, scatter!, group_to_atom_ids, get_category_ids, extract_variable, LAMMPSError, + find_compute_neighlist, find_fix_neighlist, find_pair_neighlist, # _LMP_DATATYPE LAMMPS_NONE, @@ -956,4 +957,52 @@ end _check_valid_category(category::String) = category in ("compute", "dump", "fix", "group", "molecule", "region", "variable") || error("$category is not a valid category name!") +struct NeighListElement <: AbstractVector{Int32} + iatom::Int32 + numneigh::Int + neighbors::Ptr{Int32} +end + +function Base.getindex(nle::NeighListElement, i::Integer) + 1 <= i <= nle.numneigh+1 || throw(BoundsError(nle, i)) + return i == 1 ? nle.iatom+1 : unsafe_load(nle.neighbors, i-1)+1 +end + +Base.size(nle::NeighListElement) = (nle.numneigh+1,) + +struct NeighList <: AbstractVector{NeighListElement} + lmp::LMP + idx::Int +end + +function Base.getindex(nl::NeighList, element::Integer) + iatom = Ref{Int32}() + numneigh = Ref{Int32}() + neighbors = Ref{Ptr{Int32}}() + API.lammps_neighlist_element_neighbors(nl.lmp, nl.idx, element-1, iatom, numneigh, neighbors) + iatom[] == -1 && throw(BoundsError(nl, element)) + return NeighListElement(iatom[], numneigh[], neighbors[]) +end + +Base.size(nl::NeighList) = (API.lammps_neighlist_num_elements(nl.lmp, nl.idx),) + +function find_compute_neighlist(lmp::LMP, id::String; request=0) + idx = API.lammps_find_compute_neighlist(lmp, id, request) + idx == -1 && throw(KeyError(id)) + return NeighList(lmp, idx) +end + +function find_fix_neighlist(lmp::LMP, id::String; request=0) + idx = API.lammps_find_compute_neighlist(lmp, id, request) + idx == -1 && throw(KeyError(id)) + return NeighList(lmp, idx) +end + +function find_pair_neighlist(lmp::LMP, style::String; exact=false, nsub=0, request=0) + idx = API.lammps_find_pair_neighlist(lmp, style, exact, nsub, request) + idx == -1 && throw(KeyError(style)) + return NeighList(lmp, idx) +end + + end # module From 40518fe126134cc431d642e711b6f7dacf84a036 Mon Sep 17 00:00:00 2001 From: Joroks Date: Wed, 4 Sep 2024 23:47:32 +0200 Subject: [PATCH 2/9] add documentation --- src/LAMMPS.jl | 38 +++++++++++++++++++++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index 2e5a658..828c9bf 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -965,6 +965,7 @@ end function Base.getindex(nle::NeighListElement, i::Integer) 1 <= i <= nle.numneigh+1 || throw(BoundsError(nle, i)) + # add one to account for 1-based indexing in julia. return i == 1 ? nle.iatom+1 : unsafe_load(nle.neighbors, i-1)+1 end @@ -979,25 +980,60 @@ function Base.getindex(nl::NeighList, element::Integer) iatom = Ref{Int32}() numneigh = Ref{Int32}() neighbors = Ref{Ptr{Int32}}() - API.lammps_neighlist_element_neighbors(nl.lmp, nl.idx, element-1, iatom, numneigh, neighbors) + API.lammps_neighlist_element_neighbors(nl.lmp, nl.idx, element-1 #= 0-based indexing =#, iatom, numneigh, neighbors) iatom[] == -1 && throw(BoundsError(nl, element)) return NeighListElement(iatom[], numneigh[], neighbors[]) end Base.size(nl::NeighList) = (API.lammps_neighlist_num_elements(nl.lmp, nl.idx),) +""" + find_compute_neighlist(lmp::LMP, id::String; request=0) + +Find index of a neighbor list requested by a compute + +The neighbor list request from a compute is identified by the compute ID and the request ID. +The request ID is typically 0, but will be > 0 in case a compute has multiple neighbor list requests. + +Each neighbor list contains vectors of local indices of neighboring atoms. +These can be used to index into Arrays returned form `extract_atom`. +""" function find_compute_neighlist(lmp::LMP, id::String; request=0) idx = API.lammps_find_compute_neighlist(lmp, id, request) idx == -1 && throw(KeyError(id)) return NeighList(lmp, idx) end +""" + find_fix_neighlist(lmp::LMP, id::String; request=0) + +Find index of a neighbor list requested by a fix + +The neighbor list request from a fix is identified by the fix ID and the request ID. +The request ID is typically 0, but will be > 0 in case a fix has multiple neighbor list requests. + +Each neighbor list contains vectors of local indices of neighboring atoms. +These can be used to index into Arrays returned form `extract_atom`. +""" function find_fix_neighlist(lmp::LMP, id::String; request=0) idx = API.lammps_find_compute_neighlist(lmp, id, request) idx == -1 && throw(KeyError(id)) return NeighList(lmp, idx) end +""" + find_pair_neighlist(lmp::LMP, style::String; exact=false, nsub=0, request=0) + +This function determines which of the available neighbor lists for pair styles matches the given conditions. It first matches the style name. If exact is true the name must match exactly, +if exact is false, a regular expression or sub-string match is done. If the pair style is hybrid or hybrid/overlay the style is matched against the sub styles instead. +If a the same pair style is used multiple times as a sub-style, the nsub argument must be > 0 and represents the nth instance of the sub-style (same as for the pair_coeff command, for example). +In that case nsub=0 will not produce a match and this function will Error. + +The final condition to be checked is the request ID (reqid). This will normally be 0, but some pair styles request multiple neighbor lists and set the request ID to a value > 0. + +Each neighbor list contains vectors of local indices of neighboring atoms. +These can be used to index into Arrays returned form `extract_atom`. +""" function find_pair_neighlist(lmp::LMP, style::String; exact=false, nsub=0, request=0) idx = API.lammps_find_pair_neighlist(lmp, style, exact, nsub, request) idx == -1 && throw(KeyError(style)) From 1ea87e7ad5af328c65069d30b09eb1863eef3ee9 Mon Sep 17 00:00:00 2001 From: Joroks Date: Thu, 5 Sep 2024 00:18:26 +0200 Subject: [PATCH 3/9] add tests --- test/runtests.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 1bb0f14..324b6b4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -351,4 +351,29 @@ end end end +@testset "Neighbor lists" begin + LMP(["-screen", "none"]) do lmp + command(lmp, """ + atom_modify map yes + region cell block 0 3 0 3 0 3 + create_box 1 cell + lattice sc 1 + create_atoms 1 region cell + mass 1 1 + + pair_style zero 1.0 + pair_coeff * * + + fix runfix all nve + + run 1 + """) + + neighlist = find_pair_neighlist(lmp, "zero") + @test length(neighlist) == 27 + @test length(neighlist[1]) == 4 + @test_throws KeyError find_pair_neighlist(lmp, "nonesense") + end +end + @test success(pipeline(`$(MPI.mpiexec()) -n 2 $(Base.julia_cmd()) mpitest.jl`, stderr=stderr, stdout=stdout)) From 62f8ba12ed65e727c50596fb3c5e5091fae61374 Mon Sep 17 00:00:00 2001 From: Joroks Date: Thu, 5 Sep 2024 00:24:46 +0200 Subject: [PATCH 4/9] fixup! typo --- src/LAMMPS.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index 828c9bf..c5a5e81 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -996,7 +996,7 @@ The neighbor list request from a compute is identified by the compute ID and the The request ID is typically 0, but will be > 0 in case a compute has multiple neighbor list requests. Each neighbor list contains vectors of local indices of neighboring atoms. -These can be used to index into Arrays returned form `extract_atom`. +These can be used to index into Arrays returned from `extract_atom`. """ function find_compute_neighlist(lmp::LMP, id::String; request=0) idx = API.lammps_find_compute_neighlist(lmp, id, request) @@ -1013,7 +1013,7 @@ The neighbor list request from a fix is identified by the fix ID and the request The request ID is typically 0, but will be > 0 in case a fix has multiple neighbor list requests. Each neighbor list contains vectors of local indices of neighboring atoms. -These can be used to index into Arrays returned form `extract_atom`. +These can be used to index into Arrays returned from `extract_atom`. """ function find_fix_neighlist(lmp::LMP, id::String; request=0) idx = API.lammps_find_compute_neighlist(lmp, id, request) @@ -1032,7 +1032,7 @@ In that case nsub=0 will not produce a match and this function will Error. The final condition to be checked is the request ID (reqid). This will normally be 0, but some pair styles request multiple neighbor lists and set the request ID to a value > 0. Each neighbor list contains vectors of local indices of neighboring atoms. -These can be used to index into Arrays returned form `extract_atom`. +These can be used to index into Arrays returned from `extract_atom`. """ function find_pair_neighlist(lmp::LMP, style::String; exact=false, nsub=0, request=0) idx = API.lammps_find_pair_neighlist(lmp, style, exact, nsub, request) From a9272659e15673e1f41c6d1e407fa2c050024e5b Mon Sep 17 00:00:00 2001 From: Joroks Date: Thu, 5 Sep 2024 00:34:34 +0200 Subject: [PATCH 5/9] use pair for neigh list entry --- src/LAMMPS.jl | 16 +++++++--------- test/runtests.jl | 4 +++- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index c5a5e81..06e1425 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -957,21 +957,19 @@ end _check_valid_category(category::String) = category in ("compute", "dump", "fix", "group", "molecule", "region", "variable") || error("$category is not a valid category name!") -struct NeighListElement <: AbstractVector{Int32} - iatom::Int32 +struct NeighListVec <: AbstractVector{Int32} numneigh::Int neighbors::Ptr{Int32} end -function Base.getindex(nle::NeighListElement, i::Integer) - 1 <= i <= nle.numneigh+1 || throw(BoundsError(nle, i)) - # add one to account for 1-based indexing in julia. - return i == 1 ? nle.iatom+1 : unsafe_load(nle.neighbors, i-1)+1 +function Base.getindex(nle::NeighListVec, i::Integer) + 1 <= i <= nle.numneigh || throw(BoundsError(nle, i)) + return unsafe_load(nle.neighbors, i)+1 end -Base.size(nle::NeighListElement) = (nle.numneigh+1,) +Base.size(nle::NeighListVec) = (nle.numneigh,) -struct NeighList <: AbstractVector{NeighListElement} +struct NeighList <: AbstractVector{Pair{Int32, NeighListVec}} lmp::LMP idx::Int end @@ -982,7 +980,7 @@ function Base.getindex(nl::NeighList, element::Integer) neighbors = Ref{Ptr{Int32}}() API.lammps_neighlist_element_neighbors(nl.lmp, nl.idx, element-1 #= 0-based indexing =#, iatom, numneigh, neighbors) iatom[] == -1 && throw(BoundsError(nl, element)) - return NeighListElement(iatom[], numneigh[], neighbors[]) + return iatom[]+1 => NeighListVec(numneigh[], neighbors[]) end Base.size(nl::NeighList) = (API.lammps_neighlist_num_elements(nl.lmp, nl.idx),) diff --git a/test/runtests.jl b/test/runtests.jl index 324b6b4..946b821 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -371,7 +371,9 @@ end neighlist = find_pair_neighlist(lmp, "zero") @test length(neighlist) == 27 - @test length(neighlist[1]) == 4 + iatom, neihgs = neighlist[1] + @test iatom == 1 # account for 1-based indexing + @test length(neihgs) == 3 @test_throws KeyError find_pair_neighlist(lmp, "nonesense") end end From 01e3677807188dcf8928d60ba9895d88a887a10c Mon Sep 17 00:00:00 2001 From: Joroks Date: Thu, 5 Sep 2024 00:47:17 +0200 Subject: [PATCH 6/9] add boundscheck annotation --- src/LAMMPS.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index 06e1425..a3118c9 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -963,7 +963,7 @@ struct NeighListVec <: AbstractVector{Int32} end function Base.getindex(nle::NeighListVec, i::Integer) - 1 <= i <= nle.numneigh || throw(BoundsError(nle, i)) + @boundscheck 1 <= i <= nle.numneigh || throw(BoundsError(nle, i)) return unsafe_load(nle.neighbors, i)+1 end From 832331e8de48cffb4e8fef6e2ed23ca665d962ed Mon Sep 17 00:00:00 2001 From: Joroks Date: Thu, 5 Sep 2024 13:12:32 +0200 Subject: [PATCH 7/9] add example --- src/LAMMPS.jl | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index a3118c9..f6feb43 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -1031,6 +1031,35 @@ The final condition to be checked is the request ID (reqid). This will normally Each neighbor list contains vectors of local indices of neighboring atoms. These can be used to index into Arrays returned from `extract_atom`. + +## Examples +```julia +lmp = LMP() + +command(lmp, \""" + region cell block 0 3 0 3 0 3 + create_box 1 cell + lattice sc 1 + create_atoms 1 region cell + mass 1 1 + + pair_style zero 1.0 + pair_coeff * * + + run 1 +\""") + +x = extract_atom(lmp, "x", LAMMPS_DOUBLE_2D; with_ghosts=true) + +for (iatom, neighs) in find_pair_neighlist(lmp, "zero") + for jatom in neighs + ix = @view x[:, iatom] + jx = @view x[:, jatom] + + println(ix => jx) + end +end +``` """ function find_pair_neighlist(lmp::LMP, style::String; exact=false, nsub=0, request=0) idx = API.lammps_find_pair_neighlist(lmp, style, exact, nsub, request) From 6c5a9326ebebc6433db1063f9bee566c03d5d5de Mon Sep 17 00:00:00 2001 From: Joroks <32484985+Joroks@users.noreply.github.com> Date: Thu, 5 Sep 2024 13:28:43 +0200 Subject: [PATCH 8/9] apply suggestion from review Co-authored-by: Valentin Churavy --- src/LAMMPS.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index 3fb95e2..7cc132d 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -986,9 +986,9 @@ struct NeighList <: AbstractVector{Pair{Int32, NeighListVec}} end function Base.getindex(nl::NeighList, element::Integer) - iatom = Ref{Int32}() - numneigh = Ref{Int32}() - neighbors = Ref{Ptr{Int32}}() + iatom = Ref{Cint}(-1) + numneigh = Ref{Cint}() + neighbors = Ref{Ptr{Cint}}() API.lammps_neighlist_element_neighbors(nl.lmp, nl.idx, element-1 #= 0-based indexing =#, iatom, numneigh, neighbors) iatom[] == -1 && throw(BoundsError(nl, element)) return iatom[]+1 => NeighListVec(numneigh[], neighbors[]) From 7249dc227a612a8b7106f84e59f0a31f2b5882ba Mon Sep 17 00:00:00 2001 From: Joroks Date: Thu, 5 Sep 2024 13:36:34 +0200 Subject: [PATCH 9/9] update naming --- src/LAMMPS.jl | 16 ++++++++-------- test/runtests.jl | 4 ++-- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index 7cc132d..30b1207 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -5,7 +5,7 @@ include("api.jl") export LMP, command, create_atoms, get_natoms, extract_atom, extract_compute, extract_global, extract_setting, gather, gather_bonds, gather_angles, gather_dihedrals, gather_impropers, scatter!, group_to_atom_ids, get_category_ids, extract_variable, LAMMPSError, - find_compute_neighlist, find_fix_neighlist, find_pair_neighlist, + compute_neighborlist, fix_neighborlist, pair_neighborlist, # _LMP_DATATYPE LAMMPS_NONE, @@ -997,7 +997,7 @@ end Base.size(nl::NeighList) = (API.lammps_neighlist_num_elements(nl.lmp, nl.idx),) """ - find_compute_neighlist(lmp::LMP, id::String; request=0) + compute_neighborlist(lmp::LMP, id::String; request=0) Find index of a neighbor list requested by a compute @@ -1007,14 +1007,14 @@ The request ID is typically 0, but will be > 0 in case a compute has multiple ne Each neighbor list contains vectors of local indices of neighboring atoms. These can be used to index into Arrays returned from `extract_atom`. """ -function find_compute_neighlist(lmp::LMP, id::String; request=0) +function compute_neighborlist(lmp::LMP, id::String; request=0) idx = API.lammps_find_compute_neighlist(lmp, id, request) idx == -1 && throw(KeyError(id)) return NeighList(lmp, idx) end """ - find_fix_neighlist(lmp::LMP, id::String; request=0) + fix_neighborlist(lmp::LMP, id::String; request=0) Find index of a neighbor list requested by a fix @@ -1024,14 +1024,14 @@ The request ID is typically 0, but will be > 0 in case a fix has multiple neighb Each neighbor list contains vectors of local indices of neighboring atoms. These can be used to index into Arrays returned from `extract_atom`. """ -function find_fix_neighlist(lmp::LMP, id::String; request=0) +function fix_neighborlist(lmp::LMP, id::String; request=0) idx = API.lammps_find_compute_neighlist(lmp, id, request) idx == -1 && throw(KeyError(id)) return NeighList(lmp, idx) end """ - find_pair_neighlist(lmp::LMP, style::String; exact=false, nsub=0, request=0) + pair_neighborlist(lmp::LMP, style::String; exact=false, nsub=0, request=0) This function determines which of the available neighbor lists for pair styles matches the given conditions. It first matches the style name. If exact is true the name must match exactly, if exact is false, a regular expression or sub-string match is done. If the pair style is hybrid or hybrid/overlay the style is matched against the sub styles instead. @@ -1062,7 +1062,7 @@ command(lmp, \""" x = extract_atom(lmp, "x", LAMMPS_DOUBLE_2D; with_ghosts=true) -for (iatom, neighs) in find_pair_neighlist(lmp, "zero") +for (iatom, neighs) in pair_neighborlist(lmp, "zero") for jatom in neighs ix = @view x[:, iatom] jx = @view x[:, jatom] @@ -1072,7 +1072,7 @@ for (iatom, neighs) in find_pair_neighlist(lmp, "zero") end ``` """ -function find_pair_neighlist(lmp::LMP, style::String; exact=false, nsub=0, request=0) +function pair_neighborlist(lmp::LMP, style::String; exact=false, nsub=0, request=0) idx = API.lammps_find_pair_neighlist(lmp, style, exact, nsub, request) idx == -1 && throw(KeyError(style)) return NeighList(lmp, idx) diff --git a/test/runtests.jl b/test/runtests.jl index f9c9f15..33afe0e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -369,12 +369,12 @@ end run 1 """) - neighlist = find_pair_neighlist(lmp, "zero") + neighlist = pair_neighborlist(lmp, "zero") @test length(neighlist) == 27 iatom, neihgs = neighlist[1] @test iatom == 1 # account for 1-based indexing @test length(neihgs) == 3 - @test_throws KeyError find_pair_neighlist(lmp, "nonesense") + @test_throws KeyError pair_neighborlist(lmp, "nonesense") end end