From e8009267f9a463f72c572dc28df2040f04a3bfa7 Mon Sep 17 00:00:00 2001 From: Valentin Churavy Date: Wed, 15 Mar 2023 13:08:35 -0400 Subject: [PATCH 1/6] Add support for fix external --- examples/fix_external.jl | 38 ++++++++++++++++++++++++++++++++++++ src/LAMMPS.jl | 9 +++++++-- src/external.jl | 42 ++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 24 +++++++++++++++++++++++ 4 files changed, 111 insertions(+), 2 deletions(-) create mode 100644 examples/fix_external.jl create mode 100644 src/external.jl diff --git a/examples/fix_external.jl b/examples/fix_external.jl new file mode 100644 index 0000000..dd849d0 --- /dev/null +++ b/examples/fix_external.jl @@ -0,0 +1,38 @@ +using LAMMPS + +lmp = LMP() + +command(lmp, "units lj") +command(lmp, "atom_style atomic") +command(lmp, "atom_modify map array sort 0 0") +command(lmp, "box tilt large") + +# Setup box +x_hi = 10.0 +y_hi = 10.0 +z_hi = 10.0 +command(lmp, "boundary p p p") +command(lmp, "region cell block 0 $x_hi 0 $y_hi 0 $z_hi units box") +command(lmp, "create_box 1 cell") + +command(lmp, "fix julia_lj all external pf/callback 1 1") + +# Register external fix +lj = LAMMPS.FixExternal(lmp, "julia_lj") do fix, timestep, nlocal, ids, x, fexternal + @info "julia lj called" timestep nlocal + LAMMPS.energy_global!(fix, 0.0) +end + +# Setup atoms +natoms = 10 +command(lmp, "create_atoms 1 random $natoms 1 NULL") +command(lmp, "mass 1 1.0") + +# (x,y,z), natoms +positions = rand(3, 10) .* 5 +LAMMPS.API.lammps_scatter_atoms(lmp, "x", 1, 3, positions) + +command(lmp, "run 0") + +# extract forces +forces = extract_atom(lmp, "f") \ No newline at end of file diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index a9620d5..7a78a0d 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -49,6 +49,7 @@ end mutable struct LMP @atomic handle::Ptr{Cvoid} + external_fixes::Dict{String, Any} function LMP(args::Vector{String}=String[], comm::Union{Nothing, MPI.Comm}=nothing) if !isempty(args) @@ -67,7 +68,7 @@ mutable struct LMP end end - this = new(handle) + this = new(handle, Dict{String, Any}()) finalizer(close!, this) return this end @@ -82,8 +83,10 @@ Shutdown an LMP instance. function close!(lmp::LMP) handle = @atomicswap lmp.handle = C_NULL if handle !== C_NULL - API.lammps_close(handle) + empty!(lmp.external_fixes) + API.lammps_close(handle) end + return nothing end function LMP(f::Function, args=String[], comm=nothing) @@ -376,4 +379,6 @@ function gather_atoms(lmp::LMP, name, T, count) return data end +include("external.jl") + end # module diff --git a/src/external.jl b/src/external.jl new file mode 100644 index 0000000..1499c07 --- /dev/null +++ b/src/external.jl @@ -0,0 +1,42 @@ +function fix_external_callback end + +mutable struct FixExternal + lmp::LMP + name::String + callback + + function FixExternal(lmp::LMP, name::String, callback) + if haskey(lmp.external_fixes, name) + error("FixExternal has already been registered with $name") + end + + this = new(lmp, name, callback) + lmp.external_fixes[name] = this # preserves pair globally + + ctx = Base.pointer_from_objref(this) + callback = @cfunction(fix_external_callback, Cvoid, (Ptr{Cvoid}, Int64, Cint, Ptr{Cint}, Ptr{Ptr{Float64}}, Ptr{Ptr{Float64}})) + API.lammps_set_fix_external_callback(lmp, name, callback, ctx) + + return this + end +end + +FixExternal(callback, lmp::LMP, name::String) = FixExternal(lmp::LMP, name::String, callback) + +function fix_external_callback(ctx::Ptr{Cvoid}, timestep::Int64, nlocal::Cint, ids::Ptr{Cint}, x::Ptr{Ptr{Float64}}, fexternal::Ptr{Ptr{Float64}}) + fix = Base.unsafe_pointer_to_objref(ctx)::FixExternal + + @debug "Calling fix_external_callback on" fix timestep ids x fexternal + # TODO wrap as arrays + fix.callback(fix, timestep, nlocal, ids, x, fexternal) + return nothing +end + +function energy_global!(fix::FixExternal, energy) + API.lammps_fix_external_set_energy_global(fix.lmp, fix.name, energy) +end + +# TODO +# virial_global! +# function virial_global!(fix::FixExternal, ) + diff --git a/test/runtests.jl b/test/runtests.jl index 73b319a..ea7696f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,6 +13,16 @@ LMP(["-screen", "none"]) do lmp @test_throws ErrorException command(lmp, "nonsense") end + +function f() + lmp = LMP(["-screen", "none"]) + @test LAMMPS.version(lmp) >= 0 + command(lmp, "clear") + @test_throws ErrorException command(lmp, "nonsense") + LAMMPS.close!(lmp) +end + + @testset "Variables" begin LMP(["-screen", "none"]) do lmp command(lmp, "box tilt large") @@ -42,3 +52,17 @@ end MPI.mpiexec() do mpiexec @test success(pipeline(`$mpiexec -n 2 $(Base.julia_cmd()) mpitest.jl`, stderr=stderr, stdout=stdout)) end + +LMP(["-screen", "none"]) do lmp + called = Ref(false) + command(lmp, "boundary p p p") + command(lmp, "region cell block 0 1 0 1 0 1 units box") + command(lmp, "create_box 1 cell") + command(lmp, "fix julia all external pf/callback 1 1") + LAMMPS.FixExternal(lmp, "julia") do fix, timestep, nlocal, ids, x, fexternal + called[] = true + end + command(lmp, "mass 1 1.0") + command(lmp, "run 0") + @test called[] == true +end From 4d6c5fe848c86462225de227ad4a388280889ce2 Mon Sep 17 00:00:00 2001 From: Valentin Churavy Date: Wed, 15 Mar 2023 14:27:25 -0400 Subject: [PATCH 2/6] Update test/runtests.jl --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index ea7696f..ab3fd4d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -60,6 +60,7 @@ LMP(["-screen", "none"]) do lmp command(lmp, "create_box 1 cell") command(lmp, "fix julia all external pf/callback 1 1") LAMMPS.FixExternal(lmp, "julia") do fix, timestep, nlocal, ids, x, fexternal + LAMMPS.energy_global!(fix, 0.0) called[] = true end command(lmp, "mass 1 1.0") From 69d9ac6d1a28fa2687fa2350db3a00a7c7ba6094 Mon Sep 17 00:00:00 2001 From: Valentin Churavy Date: Wed, 15 Mar 2023 16:08:41 -0400 Subject: [PATCH 3/6] finish example --- examples/fix_external.jl | 33 +++++++++++++++++++--- src/LAMMPS.jl | 33 ++++++++++++++++++++++ src/external.jl | 59 ++++++++++++++++++++++++++++++++++++++-- 3 files changed, 119 insertions(+), 6 deletions(-) diff --git a/examples/fix_external.jl b/examples/fix_external.jl index dd849d0..44e0e21 100644 --- a/examples/fix_external.jl +++ b/examples/fix_external.jl @@ -15,14 +15,39 @@ command(lmp, "boundary p p p") command(lmp, "region cell block 0 $x_hi 0 $y_hi 0 $z_hi units box") command(lmp, "create_box 1 cell") +# Use `pair_style zero` to create neighbor list for `julia_lj` +cutoff = 2.5 +command(lmp, "pair_style zero $cutoff") +command(lmp, "pair_coeff * *") command(lmp, "fix julia_lj all external pf/callback 1 1") -# Register external fix -lj = LAMMPS.FixExternal(lmp, "julia_lj") do fix, timestep, nlocal, ids, x, fexternal - @info "julia lj called" timestep nlocal - LAMMPS.energy_global!(fix, 0.0) +const coefficients = Dict( + 1 => Dict( + 1 => [48.0, 24.0, 4.0,4.0] + ) +) + +function compute_force(rsq, itype, jtype) + coeff = coefficients[itype][jtype] + r2inv = 1.0/rsq + r6inv = r2inv^3 + lj1 = coeff[1] + lj2 = coeff[2] + return (r6inv * (lj1*r6inv - lj2))*r2inv +end + +function compute_energy(rsq, itype, jtype) + coeff = coefficients[itype][jtype] + r2inv = 1.0/rsq + r6inv = r2inv^3 + lj3 = coeff[3] + lj4 = coeff[4] + return (r6inv * (lj3*r6inv - lj4)) end +# Register external fix +lj = LAMMPS.PairExternal(lmp, "julia_lj", "zero", compute_force, compute_energy, cutoff) + # Setup atoms natoms = 10 command(lmp, "create_atoms 1 random $natoms 1 NULL") diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index 7a78a0d..cb0114e 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -379,6 +379,39 @@ function gather_atoms(lmp::LMP, name, T, count) return data end +function pair_neighbor_list(lmp, name, exact, nsub, request) + idx = API.lammps_find_pair_neighlist(lmp, name, exact, nsub, request) + if idx == -1 + error("Could not find neighbor list for pair $(name)") + end + return idx +end + +function fix_neighbor_list(lmp, name, request) + idx = API.lammps_find_fix_neighlist(lmp, name, request) + if idx == -1 + error("Could not find neighbor list for fix $(name)") + end + return idx +end + + +""" + neighbors(lmb::LMP, idx, element) + +Given a neighbor list `idx` and the element therein, +return the atom index, and it's neigbors. +""" +function neighbors(lmp, idx, element) + r_iatom = Ref{Cint}() + r_numneigh = Ref{Cint}() + r_neighbors = Ref{Ptr{Cint}}(0) + + API.lammps_neighlist_element_neighbors(lmp, idx, element - 1, r_iatom, r_numneigh, r_neighbors) + + return r_iatom[], Base.unsafe_wrap(Array, r_neighbors[], r_numneigh[]; own = false) +end + include("external.jl") end # module diff --git a/src/external.jl b/src/external.jl index 1499c07..57db4fc 100644 --- a/src/external.jl +++ b/src/external.jl @@ -25,9 +25,14 @@ FixExternal(callback, lmp::LMP, name::String) = FixExternal(lmp::LMP, name::Stri function fix_external_callback(ctx::Ptr{Cvoid}, timestep::Int64, nlocal::Cint, ids::Ptr{Cint}, x::Ptr{Ptr{Float64}}, fexternal::Ptr{Ptr{Float64}}) fix = Base.unsafe_pointer_to_objref(ctx)::FixExternal + nlocal = Int(nlocal) + + @debug "Calling fix_external_callback on" fix timestep nlocal + shape = (nlocal, 3) + x = unsafe_wrap(x, shape) + fexternal = unsafe_wrap(fexternal, shape) + ids = unsafe_wrap(ids, (nlocal,)) - @debug "Calling fix_external_callback on" fix timestep ids x fexternal - # TODO wrap as arrays fix.callback(fix, timestep, nlocal, ids, x, fexternal) return nothing end @@ -36,7 +41,57 @@ function energy_global!(fix::FixExternal, energy) API.lammps_fix_external_set_energy_global(fix.lmp, fix.name, energy) end +function neighbor_list(fix::FixExternal, request) + fix_neighbor_list(fix.lmp, fix.name, request) +end + # TODO # virial_global! # function virial_global!(fix::FixExternal, ) +function PairExternal(lmp, name, neigh_name, compute_force, compute_energy, cut_global) + cutsq = cut_global^2 + FixExternal(lmp, name) do fix, timestep, nlocal, ids, x, fexternal + # Full neighbor list + idx = pair_neighbor_list(fix.lmp, neigh_name, 1, 0, 0) + nelements = API.lammps_neighlist_num_elements(fix.lmp, idx) + + type = LAMMPS.extract_atom(lmp, "type") + + # zero-out fexternal (noticed some undef memory) + fexternal .= 0 + + for ii in 1:nelements + iatom, neigh = LAMMPS.neighbors(lmp, idx, ii) + iatom += 1 # 1-based indexing + xtmp, ytmp, ztmp = view(x, :, iatom) # TODO SArray? + itype = type[iatom] + for jj in 1:length(neigh) + jatom = neigh[jj] + 1 + delx = xtmp - x[1, jatom] + dely = ytmp - x[2, jatom] + delz = ztmp = x[3, jatom] + jtype = type[jatom] + + rsq = delx*delx + dely*dely + delz*delz; + + if rsq < cutsq + fpair = compute_force(rsq, itype, jtype) + + fexternal[1, iatom] += delx*fpair + fexternal[2, iatom] += dely*fpair + fexternal[3, iatom] += delz*fpair + if jatom <= nlocal # newton_pair + fexternal[1, jatom] -= delx*fpair + fexternal[2, jatom] -= dely*fpair + fexternal[3, jatom] -= delz*fpair + end + + # todo call compute_energy + # TODO eflag + # TODO evflag + end + end + end + end +end From 5099cdf71d572b29d3c02d84f75ea5748b0c17e8 Mon Sep 17 00:00:00 2001 From: Valentin Churavy Date: Wed, 15 Mar 2023 17:05:57 -0400 Subject: [PATCH 4/6] support newton_pair --- src/LAMMPS.jl | 7 +++++++ src/external.jl | 19 +++++++++++++++++-- 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index cb0114e..e64987d 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -157,6 +157,7 @@ end function extract_global(lmp::LMP, name, dtype=nothing) if dtype === nothing dtype = API.lammps_extract_global_datatype(lmp, name) + dtype == -1 && error("Could not find dataype for global $name") end dtype = API._LMP_DATATYPE_CONST(dtype) type = dtype2type(dtype) @@ -364,6 +365,12 @@ function extract_variable(lmp::LMP, name::String, group=nothing) end end +function extract_setting(lmp, name) + val = API.lammps_extract_setting(lmp, name) + val == -1 && error("Could not find setting $name") + return val +end + function gather_atoms(lmp::LMP, name, T, count) if T === Int32 dtype = 0 diff --git a/src/external.jl b/src/external.jl index 57db4fc..85a7afa 100644 --- a/src/external.jl +++ b/src/external.jl @@ -49,6 +49,10 @@ end # virial_global! # function virial_global!(fix::FixExternal, ) +const NEIGHMASK = 0x3FFFFFFF +const SBBITS = 30 +sbmask(atom) = (atom >> SBBITS) & 3 + function PairExternal(lmp, name, neigh_name, compute_force, compute_energy, cut_global) cutsq = cut_global^2 FixExternal(lmp, name) do fix, timestep, nlocal, ids, x, fexternal @@ -56,6 +60,13 @@ function PairExternal(lmp, name, neigh_name, compute_force, compute_energy, cut_ idx = pair_neighbor_list(fix.lmp, neigh_name, 1, 0, 0) nelements = API.lammps_neighlist_num_elements(fix.lmp, idx) + # TODO how to obtain in fix + eflag = false + evflag = false + + # How to get special_lj. + newton_pair = extract_setting(fix.lmp, "newton_pair") == 1 + type = LAMMPS.extract_atom(lmp, "type") # zero-out fexternal (noticed some undef memory) @@ -67,7 +78,11 @@ function PairExternal(lmp, name, neigh_name, compute_force, compute_energy, cut_ xtmp, ytmp, ztmp = view(x, :, iatom) # TODO SArray? itype = type[iatom] for jj in 1:length(neigh) - jatom = neigh[jj] + 1 + jatom = neigh[jj] + factor_lj = 1.0 + jatom &= NEIGHMASK + jatom += 1 # 1-based indexing + delx = xtmp - x[1, jatom] dely = ytmp - x[2, jatom] delz = ztmp = x[3, jatom] @@ -81,7 +96,7 @@ function PairExternal(lmp, name, neigh_name, compute_force, compute_energy, cut_ fexternal[1, iatom] += delx*fpair fexternal[2, iatom] += dely*fpair fexternal[3, iatom] += delz*fpair - if jatom <= nlocal # newton_pair + if jatom <= nlocal || newton_pair fexternal[1, jatom] -= delx*fpair fexternal[2, jatom] -= dely*fpair fexternal[3, jatom] -= delz*fpair From e47f567f18149324b1fa8b1293b475e32d167da8 Mon Sep 17 00:00:00 2001 From: Valentin Churavy Date: Wed, 15 Mar 2023 20:33:55 -0400 Subject: [PATCH 5/6] Full test and fix silly bug --- examples/fix_external.jl | 6 ++-- src/LAMMPS.jl | 1 + src/api.jl | 4 +++ src/external.jl | 38 +++++++++++---------- test/external_pair.jl | 72 ++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 ++ 6 files changed, 104 insertions(+), 19 deletions(-) create mode 100644 test/external_pair.jl diff --git a/examples/fix_external.jl b/examples/fix_external.jl index 44e0e21..27210e5 100644 --- a/examples/fix_external.jl +++ b/examples/fix_external.jl @@ -54,10 +54,12 @@ command(lmp, "create_atoms 1 random $natoms 1 NULL") command(lmp, "mass 1 1.0") # (x,y,z), natoms -positions = rand(3, 10) .* 5 +# positions = rand(3, 10) .* 5 +positions = [4.4955289268519625 3.3999909266656836 4.420245465344918 2.3923580632470216 1.9933183377321746 2.3367019702697096 0.014668174434679937 4.5978923623562356 2.9389893820585025 4.800351333939365; 4.523573662784505 3.1582899538900304 2.5562765646443 3.199496583966941 4.891026316235915 4.689641854106464 2.7591724192198575 0.7491156338926308 1.258994308308421 2.0419941687773937; 2.256261603545908 0.694847945108647 4.058244561946366 3.044596885569421 2.60225212714946 4.0030490608195555 0.9941423774290642 1.8076536961230087 1.9712395260164222 1.2705916409499818] + LAMMPS.API.lammps_scatter_atoms(lmp, "x", 1, 3, positions) command(lmp, "run 0") # extract forces -forces = extract_atom(lmp, "f") \ No newline at end of file +forces = extract_atom(lmp, "f") diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl index e64987d..49c89ad 100644 --- a/src/LAMMPS.jl +++ b/src/LAMMPS.jl @@ -202,6 +202,7 @@ function extract_atom(lmp::LMP, name, if dtype === nothing dtype = API.lammps_extract_atom_datatype(lmp, name) + dtype == -1 && error("Could not find dataype for atom $name") dtype = API._LMP_DATATYPE_CONST(dtype) end diff --git a/src/api.jl b/src/api.jl index 578439b..41dd51d 100644 --- a/src/api.jl +++ b/src/api.jl @@ -389,6 +389,10 @@ function lammps_flush_buffers(ptr) ccall((:lammps_flush_buffers, liblammps), Cvoid, (Ptr{Cvoid},), ptr) end +function lammps_fix_external_set_energy_peratom(handle, id, eng) + ccall((:lammps_fix_external_set_energy_peratom, liblammps), Cvoid, (Ptr{Cvoid}, Ptr{Cchar}, Ptr{Cdouble}), handle, id, eng) +end + function lammps_free(ptr) ccall((:lammps_free, liblammps), Cvoid, (Ptr{Cvoid},), ptr) end diff --git a/src/external.jl b/src/external.jl index 85a7afa..58bbd5c 100644 --- a/src/external.jl +++ b/src/external.jl @@ -52,6 +52,7 @@ end const NEIGHMASK = 0x3FFFFFFF const SBBITS = 30 sbmask(atom) = (atom >> SBBITS) & 3 +const special_lj = [1.0, 0.0, 0.0 ,0.0] function PairExternal(lmp, name, neigh_name, compute_force, compute_energy, cut_global) cutsq = cut_global^2 @@ -64,49 +65,52 @@ function PairExternal(lmp, name, neigh_name, compute_force, compute_energy, cut_ eflag = false evflag = false - # How to get special_lj. newton_pair = extract_setting(fix.lmp, "newton_pair") == 1 + # special_lj = extract_global(fix.lmp, "special_lj") type = LAMMPS.extract_atom(lmp, "type") # zero-out fexternal (noticed some undef memory) fexternal .= 0 + energies = zeros(nlocal) + for ii in 1:nelements + # local atom index (i.e. in the range [0, nlocal + nghost) iatom, neigh = LAMMPS.neighbors(lmp, idx, ii) iatom += 1 # 1-based indexing xtmp, ytmp, ztmp = view(x, :, iatom) # TODO SArray? itype = type[iatom] for jj in 1:length(neigh) jatom = neigh[jj] - factor_lj = 1.0 + factor_lj = special_lj[sbmask(jatom) + 1] jatom &= NEIGHMASK jatom += 1 # 1-based indexing delx = xtmp - x[1, jatom] dely = ytmp - x[2, jatom] - delz = ztmp = x[3, jatom] + delz = ztmp - x[3, jatom] jtype = type[jatom] rsq = delx*delx + dely*dely + delz*delz; - if rsq < cutsq - fpair = compute_force(rsq, itype, jtype) - - fexternal[1, iatom] += delx*fpair - fexternal[2, iatom] += dely*fpair - fexternal[3, iatom] += delz*fpair - if jatom <= nlocal || newton_pair - fexternal[1, jatom] -= delx*fpair - fexternal[2, jatom] -= dely*fpair - fexternal[3, jatom] -= delz*fpair + fpair = factor_lj * compute_force(rsq, itype, jtype) + + if iatom <= nlocal + fexternal[1, iatom] += delx*fpair + fexternal[2, iatom] += dely*fpair + fexternal[3, iatom] += delz*fpair + if jatom <= nlocal || newton_pair + fexternal[1, jatom] -= delx*fpair + fexternal[2, jatom] -= dely*fpair + fexternal[3, jatom] -= delz*fpair + end + energies[iatom] += compute_energy(rsq, itype, jtype) end - - # todo call compute_energy - # TODO eflag - # TODO evflag end end end + API.lammps_fix_external_set_energy_peratom(fix.lmp, fix.name, energies) + energy_global!(fix, sum(energies)) end end diff --git a/test/external_pair.jl b/test/external_pair.jl new file mode 100644 index 0000000..9fb18ce --- /dev/null +++ b/test/external_pair.jl @@ -0,0 +1,72 @@ +using LAMMPS +using Test + +lmp_native = LMP() +lmp_julia = LMP() + +for lmp in (lmp_native, lmp_julia) + command(lmp, "units lj") + command(lmp, "atom_style atomic") + command(lmp, "atom_modify map array sort 0 0") + command(lmp, "box tilt large") + + # Setup box + command(lmp, "boundary p p p") + command(lmp, "region cell block 0 10.0 0 10.0 0 10.0 units box") + command(lmp, "create_box 1 cell") +end + +cutoff = 2.5 +command(lmp_julia, "pair_style zero $cutoff") +command(lmp_julia, "pair_coeff * *") +command(lmp_julia, "fix julia_lj all external pf/callback 1 1") + +command(lmp_native, "pair_style lj/cut $cutoff") +command(lmp_native, "pair_coeff * * 1 1") + +const coefficients = Dict( + 1 => Dict( + 1 => [48.0, 24.0, 4.0,4.0] + ) +) + +function compute_force(rsq, itype, jtype) + coeff = coefficients[itype][jtype] + r2inv = 1.0/rsq + r6inv = r2inv^3 + lj1 = coeff[1] + lj2 = coeff[2] + return (r6inv * (lj1*r6inv - lj2))*r2inv +end + +function compute_energy(rsq, itype, jtype) + coeff = coefficients[itype][jtype] + r2inv = 1.0/rsq + r6inv = r2inv^3 + lj3 = coeff[3] + lj4 = coeff[4] + return (r6inv * (lj3*r6inv - lj4)) +end + +# Register external fix +lj = LAMMPS.PairExternal(lmp_julia, "julia_lj", "zero", compute_force, compute_energy, cutoff) + +# Setup atoms +natoms = 10 +positions = rand(3, 10) .* 5 +for lmp in (lmp_native, lmp_julia) + command(lmp, "create_atoms 1 random $natoms 1 NULL") + command(lmp, "mass 1 1.0") + + LAMMPS.API.lammps_scatter_atoms(lmp, "x", 1, 3, positions) + + command(lmp, "run 0") +end + +# extract forces +forces_native = extract_atom(lmp_native, "f") +forces_julia = extract_atom(lmp_julia, "f") + +@testset "External Pair" begin + @test forces_native == forces_julia +end diff --git a/test/runtests.jl b/test/runtests.jl index ab3fd4d..96a76c9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -67,3 +67,5 @@ LMP(["-screen", "none"]) do lmp command(lmp, "run 0") @test called[] == true end + +include("external_pair.jl") From 6aa7cf052b4c16200dc8d055715b134bdc3fc4bd Mon Sep 17 00:00:00 2001 From: Valentin Churavy Date: Mon, 3 Apr 2023 15:38:37 -0400 Subject: [PATCH 6/6] Add CUDA example for fit_external --- examples/cuda/Manifest.toml | 478 ++++++++++++++++++++++++++++++++++ examples/cuda/Project.toml | 4 + examples/fix_external_cuda.jl | 161 ++++++++++++ 3 files changed, 643 insertions(+) create mode 100644 examples/cuda/Manifest.toml create mode 100644 examples/cuda/Project.toml create mode 100644 examples/fix_external_cuda.jl diff --git a/examples/cuda/Manifest.toml b/examples/cuda/Manifest.toml new file mode 100644 index 0000000..6067763 --- /dev/null +++ b/examples/cuda/Manifest.toml @@ -0,0 +1,478 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.9.0-beta4" +manifest_format = "2.0" +project_hash = "3bbc7dc23aed25f92593de810445734885740cee" + +[[deps.AbstractFFTs]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "16b6dbc4cf7caee4e1e75c49485ec67b667098a0" +uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" +version = "1.3.1" +weakdeps = ["ChainRulesCore"] + + [deps.AbstractFFTs.extensions] + AbstractFFTsChainRulesCoreExt = "ChainRulesCore" + +[[deps.Adapt]] +deps = ["LinearAlgebra", "Requires"] +git-tree-sha1 = "cc37d689f599e8df4f464b2fa3870ff7db7492ef" +uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +version = "3.6.1" +weakdeps = ["StaticArrays"] + + [deps.Adapt.extensions] + AdaptStaticArraysExt = "StaticArrays" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.1" + +[[deps.ArraysOfArrays]] +deps = ["Adapt", "ChainRulesCore", "StaticArraysCore", "Statistics"] +git-tree-sha1 = "c59b725b0aadf7df93fb3de05b5e1b14029af2da" +uuid = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018" +version = "0.6.1" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.Atomix]] +deps = ["UnsafeAtomics"] +git-tree-sha1 = "c06a868224ecba914baa6942988e2f2aade419be" +uuid = "a9b6321e-bd34-4604-b9c9-b65b8de01458" +version = "0.1.0" + +[[deps.BFloat16s]] +deps = ["LinearAlgebra", "Printf", "Random", "Test"] +git-tree-sha1 = "dbf84058d0a8cbbadee18d25cf606934b22d7c66" +uuid = "ab4f0b2a-ad5b-11e8-123f-65d77653426b" +version = "0.4.2" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[deps.CEnum]] +git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90" +uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" +version = "0.4.2" + +[[deps.CUDA]] +deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "CompilerSupportLibraries_jll", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "Preferences", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "SpecialFunctions", "UnsafeAtomicsLLVM"] +git-tree-sha1 = "6591ddc73adb429b9d97145c8197a0ac81664ab4" +uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" +version = "4.1.3" + +[[deps.CUDA_Driver_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] +git-tree-sha1 = "10ca2b63b496edc09258b3de5d1aa64094b18b1d" +uuid = "4ee394cb-3365-5eb0-8335-949819d2adfc" +version = "0.5.0+0" + +[[deps.CUDA_Runtime_Discovery]] +deps = ["Libdl"] +git-tree-sha1 = "6c8fceaaa6850dea627288ac3bb86fdcdf05e326" +uuid = "1af6417a-86b4-443c-805f-a4643ffb695f" +version = "0.2.0" + +[[deps.CUDA_Runtime_jll]] +deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] +git-tree-sha1 = "802b1f2220fd43251d343219adf478e6b7992bd4" +uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2" +version = "0.5.0+0" + +[[deps.ChainRulesCore]] +deps = ["Compat", "LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "c6d890a52d2c4d55d326439580c3b8d0875a77d9" +uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +version = "1.15.7" + +[[deps.Compat]] +deps = ["UUIDs"] +git-tree-sha1 = "7a60c856b9fa189eb34f5f8a6f6b5529b7942957" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.6.1" +weakdeps = ["Dates", "LinearAlgebra"] + + [deps.Compat.extensions] + CompatLinearAlgebraExt = "LinearAlgebra" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.0.2+0" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[deps.Distributed]] +deps = ["Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[deps.DocStringExtensions]] +deps = ["LibGit2"] +git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.9.3" + +[[deps.Downloads]] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" + +[[deps.ExprTools]] +git-tree-sha1 = "c1d06d129da9f55715c6c212866f5b1bddc5fa00" +uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" +version = "0.1.9" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" + +[[deps.GPUArrays]] +deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"] +git-tree-sha1 = "9ade6983c3dbbd492cf5729f865fe030d1541463" +uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" +version = "8.6.6" + +[[deps.GPUArraysCore]] +deps = ["Adapt"] +git-tree-sha1 = "1cd7f0af1aa58abc02ea1d872953a97359cb87fa" +uuid = "46192b85-c4d5-4398-a991-12ede77f4527" +version = "0.1.4" + +[[deps.GPUCompiler]] +deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "TimerOutputs", "UUIDs"] +git-tree-sha1 = "590d394bad1055b798b2f9b308327ba871b7badf" +uuid = "61eb1bfa-7361-4325-ad38-22787b887f55" +version = "0.19.0" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[deps.IrrationalConstants]] +git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.2.2" + +[[deps.JLLWrappers]] +deps = ["Preferences"] +git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.4.1" + +[[deps.KernelAbstractions]] +deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "SnoopPrecompile", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] +git-tree-sha1 = "350a880e80004f4d5d82a17f737d8fcdc56c3462" +uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c" +version = "0.9.1" + +[[deps.LAMMPS]] +deps = ["CEnum", "LAMMPS_jll", "Libdl", "MPI", "Preferences"] +path = "../.." +uuid = "ee2e13b9-eee9-4449-aafa-cfa6a2dbe14d" +version = "0.2.0" + +[[deps.LAMMPS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "TOML"] +git-tree-sha1 = "c845a182c828051bf429174ee7b614c643e4221f" +uuid = "5b3ab26d-9607-527c-88ea-8fe5ba57cafe" +version = "2.3.2+0" + +[[deps.LLVM]] +deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"] +git-tree-sha1 = "a8960cae30b42b66dd41808beb76490519f6f9e2" +uuid = "929cbde3-209d-540e-8aea-75f648917ca0" +version = "5.0.0" + +[[deps.LLVMExtra_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] +git-tree-sha1 = "09b7505cc0b1cee87e5d4a26eea61d2e1b0dcd35" +uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" +version = "0.0.21+0" + +[[deps.LazyArtifacts]] +deps = ["Artifacts", "Pkg"] +uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.3" + +[[deps.LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" +version = "7.84.0+0" + +[[deps.LibGit2]] +deps = ["Base64", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.10.2+0" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.LogExpFunctions]] +deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "0a1b7c2863e44523180fdb3146534e265a91870b" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.3.23" + + [deps.LogExpFunctions.extensions] + LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" + LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" + LogExpFunctionsInverseFunctionsExt = "InverseFunctions" + + [deps.LogExpFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[deps.MPI]] +deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "Requires", "Serialization", "Sockets"] +git-tree-sha1 = "6d72bafd3960f9c119ceb8f034fef28346490fe5" +uuid = "da04e1cc-30fd-572f-bb4f-1f8673147195" +version = "0.20.8" + +[[deps.MPICH_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "d790fbd913f85e8865c55bf4725aff197c5155c8" +uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4" +version = "4.1.1+1" + +[[deps.MPIPreferences]] +deps = ["Libdl", "Preferences"] +git-tree-sha1 = "71f937129731a29eabe6969db2c90368a4408933" +uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" +version = "0.1.7" + +[[deps.MPItrampoline_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "ad88f863a5a16b3e26d14446afd3cd746266281b" +uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748" +version = "5.2.1+3" + +[[deps.MacroTools]] +deps = ["Markdown", "Random"] +git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.10" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.0+0" + +[[deps.MicrosoftMPI_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "a16aa086d335ed7e0170c5265247db29172af2f9" +uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf" +version = "10.1.3+2" + +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2022.10.11" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.21+0" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.1+0" + +[[deps.OpenMPI_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "f3080f4212a8ba2ceb10a34b938601b862094314" +uuid = "fe0851c0-eecd-5654-98d4-656369965a5c" +version = "4.1.5+0" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.5+0" + +[[deps.Pkg]] +deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +version = "1.9.0" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.3.0" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[deps.Random]] +deps = ["SHA", "Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.Random123]] +deps = ["Random", "RandomNumbers"] +git-tree-sha1 = "7a1a306b72cfa60634f03a911405f4e64d1b718b" +uuid = "74087812-796a-5b5d-8853-05524746bad3" +version = "1.6.0" + +[[deps.RandomNumbers]] +deps = ["Random", "Requires"] +git-tree-sha1 = "043da614cc7e95c703498a491e2c21f58a2b8111" +uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143" +version = "1.5.3" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.0" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.Scratch]] +deps = ["Dates"] +git-tree-sha1 = "30449ee12237627992a99d5e30ae63e4d78cd24a" +uuid = "6c6a2e73-6563-6170-7368-637461726353" +version = "1.2.0" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.SnoopPrecompile]] +deps = ["Preferences"] +git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c" +uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c" +version = "1.0.3" + +[[deps.Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[deps.SparseArrays]] +deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.SpecialFunctions]] +deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "ef28127915f4229c971eb43f3fc075dd3fe91880" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.2.0" +weakdeps = ["ChainRulesCore"] + + [deps.SpecialFunctions.extensions] + SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"] +git-tree-sha1 = "b8d897fe7fa688e93aef573711cb207c08c9e11e" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.5.19" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.0" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +version = "1.9.0" + +[[deps.SuiteSparse_jll]] +deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"] +uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" +version = "5.10.1+6" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.0" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.TimerOutputs]] +deps = ["ExprTools", "Printf"] +git-tree-sha1 = "f2fd3f288dfc6f507b0c3a2eb3bac009251e548b" +uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +version = "0.5.22" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[deps.UnsafeAtomics]] +git-tree-sha1 = "6331ac3440856ea1988316b46045303bef658278" +uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f" +version = "0.2.1" + +[[deps.UnsafeAtomicsLLVM]] +deps = ["LLVM", "UnsafeAtomics"] +git-tree-sha1 = "ea37e6066bf194ab78f4e747f5245261f17a7175" +uuid = "d80eeb9a-aca5-4d75-85e5-170c8b632249" +version = "0.1.2" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.13+0" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.4.0+0" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.48.0+0" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+0" diff --git a/examples/cuda/Project.toml b/examples/cuda/Project.toml new file mode 100644 index 0000000..b2be155 --- /dev/null +++ b/examples/cuda/Project.toml @@ -0,0 +1,4 @@ +[deps] +ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018" +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +LAMMPS = "ee2e13b9-eee9-4449-aafa-cfa6a2dbe14d" diff --git a/examples/fix_external_cuda.jl b/examples/fix_external_cuda.jl new file mode 100644 index 0000000..da06082 --- /dev/null +++ b/examples/fix_external_cuda.jl @@ -0,0 +1,161 @@ +using LAMMPS +using CUDA +using Adapt +using ArraysOfArrays + +lmp = LMP() + +command(lmp, "units lj") +command(lmp, "atom_style atomic") +command(lmp, "atom_modify map array sort 0 0") +command(lmp, "box tilt large") + +# Setup box +x_hi = 10.0 +y_hi = 10.0 +z_hi = 10.0 +command(lmp, "boundary p p p") +command(lmp, "region cell block 0 $x_hi 0 $y_hi 0 $z_hi units box") +command(lmp, "create_box 1 cell") + +# Use `pair_style zero` to create neighbor list for `julia_lj` +cutoff = 2.5 +command(lmp, "pair_style zero $cutoff") +command(lmp, "pair_coeff * *") +command(lmp, "fix julia_lj_cuda all external pf/callback 1 1") + +struct Potential{Coeff} + coefficients::Coeff +end + +function compute_force(p::Potential, rsq, itype, jtype) + coeff = p.coefficients[itype][jtype] + r2inv = 1.0/rsq + r6inv = r2inv^3 + lj1 = coeff[1] + lj2 = coeff[2] + return (r6inv * (lj1*r6inv - lj2))*r2inv +end + +function compute_energy(p::Potential, rsq, itype, jtype) + coeff = p.coefficients[itype][jtype] + r2inv = 1.0/rsq + r6inv = r2inv^3 + lj3 = coeff[3] + lj4 = coeff[4] + return (r6inv * (lj3*r6inv - lj4)) +end + +const NEIGHMASK = 0x3FFFFFFF +const SBBITS = 30 +sbmask(atom) = (atom >> SBBITS) & 3 +const special_lj = CuArray([1.0, 0.0, 0.0 ,0.0]) + +LAMMPS.FixExternal(lmp, "julia_lj_cuda") do fix, timestep, nlocal, ids, x, fexternal + # Full neighbor list + idx = LAMMPS.pair_neighbor_list(fix.lmp, "zero", 1, 0, 0) + nelements = LAMMPS.API.lammps_neighlist_num_elements(fix.lmp, idx) + + potential = Potential((((48.0, 24.0, 4.0, 4.0),),)) + cutsq::Float64 = cutoff^2 + + # TODO how to obtain in fix + eflag = false + evflag = false + + newton_pair = LAMMPS.extract_setting(fix.lmp, "newton_pair") == 1 + # special_lj = extract_global(fix.lmp, "special_lj") + + type = LAMMPS.extract_atom(lmp, "type") + + # zero-out fexternal (noticed some undef memory) + fexternal .= 0 + + energies = CUDA.zeros(nlocal) + cu_fexternal = adapt(CuArray, fexternal) + cu_x = adapt(CuArray, x) + cu_type = adapt(CuArray, type) + + neighbors_array = VectorOfArrays{Int64, 2}() + ilist = Vector{Int64}() + + # Copy neighbor_list to Julia datastructure + for ii in 1:nelements + # local atom index (i.e. in the range [0, nlocal + nghost) + iatom, neigh = LAMMPS.neighbors(lmp, idx, ii) + push!(neighbors_array, reshape(neigh, (1, length(neigh)))) + push!(ilist, iatom) + end + neighbors_array = adapt(CuArray, neighbors_array) + ilist = adapt(CuArray, ilist) + + function kernel(potential, x, fexternal, energies, ilist, neighbors, cutsq, nlocal, type, special_lj) + ii = threadIdx().x + iatom = ilist[ii] + neighs = neighbors[ii] + + iatom += 1 # 1-based indexing + xtmp = x[1, iatom] + ytmp = x[2, iatom] + ztmp = x[3, iatom] + + itype = type[iatom] + + for jatom in neighs + factor_lj = special_lj[sbmask(jatom) + 1] + jatom &= NEIGHMASK + jatom += 1 # 1-based indexing + + delx = xtmp - x[1, jatom] + dely = ytmp - x[2, jatom] + delz = ztmp - x[3, jatom] + + jtype = type[jatom] + + rsq = delx*delx + dely*dely + delz*delz; + if rsq < cutsq + fpair = factor_lj * compute_force(potential, rsq, itype, jtype) + + if iatom <= nlocal + CUDA.@atomic fexternal[1, iatom] += delx*fpair + CUDA.@atomic fexternal[2, iatom] += dely*fpair + CUDA.@atomic fexternal[3, iatom] += delz*fpair + if jatom <= nlocal || newton_pair + CUDA.@atomic fexternal[1, jatom] -= delx*fpair + CUDA.@atomic fexternal[2, jatom] -= dely*fpair + CUDA.@atomic fexternal[3, jatom] -= delz*fpair + end + energies[iatom] += compute_energy(potential, rsq, itype, jtype) + end + end + end + return nothing + end + + @cuda threads=nlocal kernel(potential, cu_x, cu_fexternal, energies, + ilist, neighbors_array, + cutsq, nlocal, cu_type, special_lj) + + copyto!(fexternal, cu_fexternal) # TODO async + total_energy = sum(energies) + energies = Array(energies) + + LAMMPS.API.lammps_fix_external_set_energy_peratom(fix.lmp, fix.name, energies) + LAMMPS.energy_global!(fix, total_energy) +end + +# Setup atoms +natoms = 10 +command(lmp, "create_atoms 1 random $natoms 1 NULL") +command(lmp, "mass 1 1.0") + +# (x,y,z), natoms +# positions = rand(3, 10) .* 5 +positions = [4.4955289268519625 3.3999909266656836 4.420245465344918 2.3923580632470216 1.9933183377321746 2.3367019702697096 0.014668174434679937 4.5978923623562356 2.9389893820585025 4.800351333939365; 4.523573662784505 3.1582899538900304 2.5562765646443 3.199496583966941 4.891026316235915 4.689641854106464 2.7591724192198575 0.7491156338926308 1.258994308308421 2.0419941687773937; 2.256261603545908 0.694847945108647 4.058244561946366 3.044596885569421 2.60225212714946 4.0030490608195555 0.9941423774290642 1.8076536961230087 1.9712395260164222 1.2705916409499818] + +LAMMPS.API.lammps_scatter_atoms(lmp, "x", 1, 3, positions) + +command(lmp, "run 0") + +# extract forces +forces = extract_atom(lmp, "f")