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 b1e3de2..f597c5b 100644
--- a/src/LAMMPS.jl
+++ b/src/LAMMPS.jl
@@ -186,6 +186,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 c39d09f..3bdee1c 100644
--- a/src/api.jl
+++ b/src/api.jl
@@ -295,6 +295,10 @@ function lammps_fix_external_set_virial_global(arg1, arg2, arg3)
     ccall((:lammps_fix_external_set_virial_global, liblammps), Cvoid, (Ptr{Cvoid}, Ptr{Cchar}, Ptr{Cdouble}), arg1, arg2, arg3)
 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 52ef562..4634ebd 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -31,3 +31,5 @@ LMP(["-screen", "none"]) do lmp
     command(lmp, "run 0")
     @test called[] == true
 end
+
+include("external_pair.jl")