support newton_pair

cesmix-mit · Mar 15, 2023 · 93bacd2 · 93bacd2
1 parent 879ba78
commit 93bacd2
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 2 deletions.
diff --git a/src/LAMMPS.jl b/src/LAMMPS.jl
@@ -141,6 +141,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)
@@ -285,6 +286,12 @@ function extract_compute(lmp::LMP, name, style, type)
     return nothing
 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
@@ -49,13 +49,24 @@ 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
         # Full neighbor list
         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