# Demonstrate bispectrum computes

# initialize simulation

variable 	nsteps index 0
variable 	nrep equal 1
variable 	a equal 2.0
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

atom_modify	map hash
lattice         bcc $a
region		box block 0 ${nx} 0 ${ny} 0 ${nz}
create_box	2 box
create_atoms	2 box

mass 		* 180.88

displace_atoms 	all random 0.1 0.1 0.1 123456

# choose SNA parameters

variable 	twojmax equal 2
variable 	rcutfac equal 1.0
variable 	rfac0 equal 0.99363
variable 	rmin0 equal 0
variable 	radelem1 equal 2.3
variable 	radelem2 equal 2.0
variable	wj1 equal 1.0
variable	wj2 equal 0.96
variable	quadratic equal 1
variable	bzero equal 0
variable	switch equal 0
variable 	snap_options string &
"${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"

# set up dummy potential to satisfy cutoff

pair_style 	zero ${rcutfac}
pair_coeff 	* *

# set up reference potential

variable 	zblcutinner equal 4
variable 	zblcutouter equal 4.8
variable 	zblz equal 73
pair_style 	zbl ${zblcutinner} ${zblcutouter}
pair_coeff 	* * ${zblz} ${zblz}

# set up per-atom computes

compute 	b all sna/atom ${snap_options}
compute 	vb all snav/atom ${snap_options}
compute 	db all snad/atom ${snap_options}

# perform sums over atoms

group 		snapgroup1 type 1
group 		snapgroup2 type 2
compute         bsum1 snapgroup1 reduce sum c_b[*]
compute         bsum2 snapgroup2 reduce sum c_b[*]
# fix 		bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
# fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
compute		vbsum all reduce sum c_vb[*]
# fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
variable	db_2_120 equal c_db[2][120]

# set up compute snap generating global array

compute 	snap all snap ${snap_options}
fix 		snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector

thermo 		100

# test output:   1: total potential energy
#                2: xy component of stress tensor
#                3: Sum(0.5*(B_{222}^i)^2, all i of type 2) 
#                4: xy component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j) 
#                5: z component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
#
#                followed by 5 counterparts from compute snap

thermo_style	custom &
		pe            pxy            c_bsum2[20]   c_vbsum[240]   v_db_2_120 & 
		c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40]
thermo_modify 	norm no

# dump 		mydump_db all custom 1000 dump_db id c_db[*]
# dump_modify 	mydump_db sort id

# Run MD

run             ${nsteps}