fixup! Add CUDA example for fit_external

examples/fix_external_cuda.jl

@@ -77,21 +77,24 @@ LAMMPS.FixExternal(lmp, "julia_lj_cuda") do fix, timestep, nlocal, ids, x, fexte
     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)
-        _, neigh = LAMMPS.neighbors(lmp, idx, ii)
+        iatom, neigh = LAMMPS.neighbors(lmp, idx, ii)
         push!(neighbors_array, reshape(neigh, (1, length(neigh))))
+        push!(ilist, iatom)
     end
-    @show neighbors_array
-
     neighbors_array = adapt(CuArray, neighbors_array)
+    ilist = adapt(CuArray, ilist)
 
-    function kernel(potential, x, fexternal, energies, neighbors, cutsq, nlocal, type, special_lj)
-        iatom = threadIdx().x
-        neighs = neighbors[iatom]
+    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]
@@ -129,7 +132,8 @@ LAMMPS.FixExternal(lmp, "julia_lj_cuda") do fix, timestep, nlocal, ids, x, fexte
         return nothing
     end
 
-    @cuda threads=nlocal kernel(potential, cu_x, cu_fexternal, energies, neighbors_array,
+    @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