units           lj
dimension       3
boundary        p p p 
#processors      * * *
atom_style      sphere

python source here """
from lammps import lammps
import math
import numpy as np

def get_force(lmpptr):
   # initialize python lammps class from C++ pointer
   lmp = lammps(ptr=lmpptr)

   # conversion factor
   num = lmp.get_natoms()

   force = lmp.numpy.extract_compute('f', 2, 1) # pair force
   aforce = lmp.numpy.extract_compute('af', 1, 2) # atom force
   npair = lmp.extract_compute('f', 2, 0)
   nonRat = lmp.get_natoms() - int(lmp.extract_compute('rat',0, 0))
   isRat = lmp.extract_compute('rat', 1, 1)

   if npair == 0: return 0.

   # for computing the average
   total_force = np.mean(force)
   print(len(aforce))
   print(aforce)

  #  lst = [np.linalg.norm(f) for f in aforce]
  #  sumf = np.mean(np.array(lst))
   
   sumf = 0.
   fnorm = 0.
   for i in range(lmp.get_natoms()):
     if isRat[i] == 1: 
        print(i, np.linalg.norm(aforce[i]))
        continue
     fnorm+=aforce[i][0]**2+aforce[i][1]**2+aforce[i][2]**2
     sumf +=  np.linalg.norm(aforce[i])
     #if i<10:
      #  print(math.sqrt(aforce[i][0]**2+aforce[i][1]**2+aforce[i][2]**2))
   
   print(nonRat, npair)
   #total_force /= npair
   sumf /= nonRat
   #print('fnorm:', math.sqrt(fnorm), 'average pair force:', total_force, 'average force:', sumf)
   return sumf/total_force
"""

# Note: the minimum paticle diameter should be 1
read_data       test.data
compute         pro all property/atom mass
#dump            4a all custom 1 tmp.dump x y z fx fy fz radius

variable        num_particles equal atoms

# https://docs.lammps.org/pair_granular.html
pair_style      granular
pair_coeff      * * hooke 1 0.5 damping mass_velocity tangential linear_history 1 1 0 limit_damping
comm_modify     vel yes

neighbor        0.3 bin 
neigh_modify    every 1 delay 0 check yes
timestep        0.02 # 0.02*sqrt(mass/kn)

# Iterativly remove all rattlers
compute         rat all rattlers/atom radius 4 10000
run             0

variable        pp atom "c_rat == 0"
group           nonrattler variable pp
dump 2 nonrattler custom 1 test.dump id 

compute         af all property/atom fx fy fz
compute         f all pair/local force
#dump            4a all local 10000 tmp.dump c_f


velocity        all set 0.0 0.0 0.0 units box

variable        avg_f python get_force
python          get_force input 1 SELF return v_avg_f format pf

#compute         1 all contact/atom
compute         2 all pressure NULL virial

variable        num_rattlers equal c_rat

thermo_style    custom step press c_2 fnorm fmax v_avg_f ke density c_rat v_num_rattlers
thermo          1
thermo_modify   lost warn

run             0

variable        spress equal press
variable        sfmax equal fmax 
variable        sfnorm equal fnorm

print           "${spress} ${sfmax} ${sfnorm} ${avg_f}"
print           "Simulation Completes."