# Reference:
# Granular packings with sliding, rolling, and twisting friction

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)

   num = lmp.get_natoms()
   nonRat = num - int(lmp.extract_compute('rat',0, 0))
   if nonRat == 0: return 1.
   print(nonRat)

   force = lmp.extract_compute('f', 2, 1) # pair force
   aforce = lmp.numpy.extract_compute('af', 1, 2) # atom force
   npair = lmp.extract_compute('f', 2, 0)

   if npair == 0: return 1.
   print(npair)

   isRat = lmp.extract_compute('rat', 1, 1)

   # pair force
   sum_force = 0.
   for i in range(npair):
     sum_force +=  force[i]

   total_force = 0.
   fnorm = 0.
   for i in range(num):
     #print(isRat[i])
     if isRat[i] == 1: continue
     #fnorm += aforce[i][0]**2+aforce[i][1]**2+aforce[i][2]**2
     total_force += np.linalg.norm(aforce[i])
   
   sum_force /= npair
   total_force /= nonRat
   print('fnorm:', math.sqrt(fnorm), 'average pair force:', sum_force, 'average total force:', total_force)
   return total_force/sum_force
"""

variable        tmr timer
variable        seednum equal ceil(${tmr}*100000000)
#variable        seednum equal round(random(1,100000000,12345))
variable        num_particles equal 100
variable        sigma equal 1
variable        vol_particles equal ${num_particles}*PI*${sigma}^3/6

variable        init_phi equal 0.005
variable        boxL equal (${vol_particles}/${init_phi})^(1/3)
region          box block 0 ${boxL} 0 ${boxL} 0 ${boxL}
create_box      1 box
create_atoms    1 random ${num_particles} ${seednum} box overlap 1

set             atom 1 diameter ${sigma}
set             atom 1 mass 1

# https://docs.lammps.org/pair_granular.html
pair_style      granular
# limit_damping or modify source code
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

compute         af nonrattler property/atom fx fy fz
compute         f nonrattler pair/local force

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

velocity        all set 0.0 0.0 0.0 units box

variable        Pa equal 1e-6
variable        Pdamp equal 2.25

#change_box      all triclinic

#fix             1 all nph/sphere x ${Pa} ${Pa} ${Pdamp} y ${Pa} ${Pa} ${Pdamp} z ${Pa} ${Pa} ${Pdamp} xy 0.0 0.0 ${Pdamp} yz 0.0 0.0 ${Pdamp} xz 0.0 0.0 ${Pdamp} nreset 1 pchain 0 ptemp 1
fix             1 all nph/sphere iso ${Pa} ${Pa} ${Pdamp} nreset 1 pchain 0 ptemp 1

compute         press all pressure NULL virial


variable        num_rattlers equal c_rat
variable        phi equal ${vol_particles}/lx/ly/lz
thermo_style    custom step c_press ke fnorm fmax v_avg_f c_rat v_phi v_num_rattlers
thermo_modify   format 8 %20.16f
thermo          10000
thermo_modify   lost warn

variable        ave_ke equal ke/${num_particles}
#dump            1 all xyz 10000 dump.xyz
dump            1 all custom 10000 mono_dump.xyz id radius x y z

variable        stop equal “(v_ave_ke < 1e-16) && (v_avg_f < 1e-4)”
fix             2 all halt 10000 v_stop > 0 error continue

#print "${ave_contact}" append analysis.txt

#fix             extra all print 10000 "${ave_contact}" append analysis.txt
run             100000000

unfix           1
unfix           2

velocity        all set 0.0 0.0 0.0 units box

run             0
write_dump      all custom ${num_particles}_silbert.dump id type radius x y z modify append no format line "%d %d %20.16e %20.16e %20.16e %20.16e"

variable        spress equal c_press
variable        sfmax equal fmax 
variable        sfnorm equal fnorm 

print           "${spress} ${sfmax} ${sfnorm} ${avg_f} ${phi}"
print           "Simulation Complete."