# LAMMPS Input File for create simple steel 310S monocrystal 
# By Zbigniew Koziol, 2021. Based on code of Szymon Romanowski, Dec 2016

# These and other properties are from:
# https://www.azom.com/article.aspx?ArticleID=4392
# 310S Chemical composition (weight%)
# C	0.08 max
# Mn	2.00 max
# Si	1.50 max
# P	0.045 max
# S	0.030 max
# Cr	24.00-26.00
# Mo	0.75 max
# Ni	19.00-21.00
# Fe	53.095 min
#
# Here we assume Fe, Ni, Cr only. Hence, let us have (weight%): 
# Ni 20		# 58.6934 
# Cr 25		# 51.9961 
# Fe 55		# 55.8452 

# Here xFe, xNi, xCr is in percents, to be found from equations:
# GNUPLOT:
# mfe=55.8452;mni=58.6934;mcr=51.9961
# v=1./(0.55/mfe+0.2/mni+0.25/mcr) = 55.3579763872491
# xFe = 0.55 * v/mfe  = 0.54520150367421
# xNi = 0.20 * v/mni  = 0.188634416773433
# xCr = 0.25 * v/mcr  = 0.266164079552356

variable mfe equal 55.8452
variable mni equal 58.6934
variable mcr equal 51.9961

variable xfe equal 0.5452 # it is not used, though
variable xni equal 0.18863
# Take care here. We would like to have the final fraction 0.26616 for xcr.
# However, the order the fractions are created is that first xni fraction
# is created and after that xcr from what is left of iron.
# Therefore we must normalize 0.26616 by a factor 1/(1-xni)=1.23248 in order
# to obtain correct value of xcr=0.26616:
variable xcr equal 0.32804

# shell mkdir results310S
# shell cd results310S

# variable thisDir equal `pwd`
# print "xxxxxxxx" 
# print ${thisDir}


clear

# ---------- Define parametr of simulation --------------------- 
variable minsteps equal 150			# minimal number of minimization steps
variable latticeparam equal 3.56200762

clear
units metal 
dimension 3 
boundary p p p 
atom_style atomic 

# ---------- Create Box --------------------- 
lattice fcc ${latticeparam}  
region whole block 0 30 0 30 0 30 units lattice
create_box 3 whole 

# ---------- Define Interatomic Potential --------------------- 
pair_style eam/alloy
pair_coeff * * ../FeNiCr_Bonny_2013_ptDef.eam.alloy Fe Ni Cr
neighbor 2.0 bin 
neigh_modify delay 10 check yes 
 
mass 1 ${mfe} # Fe
mass 2 ${mni} # Ni
mass 3 ${mcr} # Cr

# ---------- Create Atomistic Structure --------------------- 
create_atoms 1 region whole  
set type 1 type/fraction 2 ${xni} 765556315
set type 1 type/fraction 3 ${xcr} 764688554

group iron type 1
group nickel type 2
group chromium type 3
# ---------- Displace atoms and delete overlapping atoms --------------------- 
delete_atoms overlap 1.4 all all

# ---------- Define Settings --------------------- 
compute csym all centro/atom fcc
compute eng all pe/atom 
compute eatoms all reduce sum c_eng 

#---------- Dump data into Data file ------------- 
reset_timestep 0 
#dump 		2 all custom 1 dump_*.cfg id type xs ys zs
#dump_modify     2 element Fe Ni Cr
#---------- Minimize  --------------------- 
reset_timestep 0 
thermo 10 
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms 
min_style cg 
minimize 1e-15 1e-15 ${minsteps} ${minsteps}

write_data steel310S.dat
print "$(count(all)) $(count(iron)) $(count(nickel)) $(count(chromium))"
print "All done"