# Demonstrate bispectrum per-atom and grid computes

# Invoked from grid.py
# pass in values for ngridx, ngridy, ngridz, twojmax

# Initialize simulation

variable nsteps equal 0
variable nrep equal 1
variable a equal 3.316

units		metal

# generate the box and atom positions using a BCC lattice

variable nx equal ${nrep}
variable ny equal ${nrep}
variable nz equal ${nrep}

boundary	p p p

lattice		custom ${a} a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0 0 0 basis 0.5 0.5 0.5

region		box block 0 ${nx} 0 ${ny} 0 ${nz}
create_box	1 box
create_atoms	1 box

mass 1 180.88

# define grid compute and atom compute

group 		snapgroup type 1
variable 	rcutfac equal 4.67637
variable 	rfac0 equal 0.99363
variable 	rmin0 equal 0
variable 	wj equal 1
variable 	radelem equal 0.5
variable 	bzero equal 0
variable 	quad equal 0
variable 	switch equal 1

compute b all sna/atom ${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} bzeroflag ${bzero} quadraticflag ${quad} switchflag ${switch}

compute bgrid all sna/grid grid ${ngridx} ${ngridy} ${ngridz} ${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} bzeroflag ${bzero}	quadraticflag ${quad} switchflag ${switch}

compute bgridlocal all sna/grid/local grid ${ngridx} ${ngridy} ${ngridz} ${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} bzeroflag ${bzero}	quadraticflag ${quad} switchflag ${switch}

# create dummy potential for neighbor list

variable rcutneigh equal 2.0*${rcutfac}*${radelem}
pair_style zero ${rcutneigh}
pair_coeff * *

# define output

thermo_style	custom step temp ke pe vol c_bgrid[1][1]
thermo_modify norm yes

dump mydump_b all custom 1000 dump_b id c_b[*]

dump mydump_bgridlocal all local 1000 dump_bgridlocal index c_bgridlocal[1]

# run

run 0