in.equilibrate

# Simple script to try and initialize a block of atoms at a certain temperature and zero pressure.
# Of course all the warnings about modeling reality with a limited number of atoms apply.
# Values for box size (Nx) and temperature should be increased in a "real" run. Below are small values
# that run quickly on my laptop.


# read in some parameters from command line
variable	Nx index 10  # half-width of cube in units of lattice constants
variable	T index 300  # temperature
variable	seed index 12345  # random number seed

# basic setup of cube
units		metal
atom_style	atomic
atom_modify	map yes  # for compute voronoi
lattice		bcc 2.855312
region		box block -${Nx} ${Nx} -${Nx} ${Nx} -${Nx} ${Nx}
create_box	1 box
create_atoms	1 box

# use this pair potential
pair_style	eam/fs
pair_coeff	* * potentials/Fe_2.eam.fs Fe

# get ready to output
variable	boxstring string "Nx-${Nx}-T-$T"
variable	outdir string "outputs/${boxstring}/equilibrate"
shell		mkdir outputs  # intermediate because 'mkdir -p' doesn't work
shell		mkdir outputs/${boxstring}  # intermediate because 'mkdir -p' doesn't work
shell		mkdir ${outdir}

# minimize the box
fix		box_relax all box/relax iso 0.0
minimize	1.0e-10 1.0e-10 1000 1000
unfix		box_relax
reset_timestep	0

# initialize velocity
variable	double_T equal 2.0*${T}
velocity	all create ${double_T} ${seed} dist gaussian

# output info
thermo		100
thermo_style	custom step time temp press
fix		output all ave/time 1 5 100 c_thermo_temp c_thermo_press file ${outdir}/time-history.txt

# try to get to equilibrate at zero pressure
print		"*** NPT equilibration ***"
fix		fix_npt all npt temp ${T} ${T} $(100.0*dt) iso 0.0 0.0 $(1000.0*dt) drag 1.0
run		5000
unfix		fix_npt

# little more NVE to check
print		"*** more NVE equilibration ***"
fix		fix_nve all nve
run             5000
unfix		fix_nve

# define a spherical shell region and group to pick a PKA from later
variable	shell_width equal 5.0
variable	outer_radius equal xhi/2.0
variable	inner_radius equal ${outer_radius}-${shell_width}
region		big_sphere sphere 0.0 0.0 0.0 ${outer_radius} side in units box
region		little_sphere sphere 0.0 0.0 0.0 ${inner_radius} side out units box
region		shell_region intersect 2 big_sphere little_sphere
group		shell_group region shell_region

# write to file
write_dump	shell_group custom ${outdir}/dump.shell.txt id type mass x y z

# write restart and atom coords so we can restart from here
write_restart	${outdir}/restart
write_dump	all custom ${outdir}/dump.end.txt id type mass x y z