Creating Polycrystalline Samples via Simulated Sintering

LAMMPS input script

# Sintering Example, Work in Progress

# ENH ericnhahn@gmail.com

# updated 6/7/2020

# !!!WORK IN PROGRESS!!!

# This input script is still in its testing phase, may be incomplete, inaccurate, include redundant pieces of code, etc.

variable name string sinterfccpoly_Cu_WIP

units metal

boundary p p p

atom_style atomic

neighbor 2.0 bin

neigh_modify every 1 delay 2 check yes

##############################################

# Setting up Grains/Crystallites #

##############################################

variable latparam equal 3.615 #if you care more about a fully dense crystal structure than studying sintering a smaller latparam than equilibirium can be used, for ex 3.3

variable fcc equal 4

#dimensions can be altered from cubic by modifying below

variable dim equal 250

variable x equal ${dim}

variable y equal ${dim}

variable z equal ${dim}

variable 1z1 equal 0

variable 1z2 equal 0

variable 1z3 equal 1

variable 1y1 equal 0

variable 1y2 equal 1

variable 1y3 equal 0

variable 1x1 equal 1

variable 1x2 equal 0

variable 1x3 equal 0

variable xsize1 equal "sqrt(v_1x1^2 + v_1x2^2 + v_1x3^2)"

variable ysize1 equal "sqrt(v_1y1^2 + v_1y2^2 + v_1y3^2)"

variable zsize1 equal "sqrt(v_1z1^2 + v_1z2^2 + v_1z3^2)"

variable sx equal "ceil(v_x/v_latparam/v_xsize1/2)"

variable sy equal "ceil(v_y/v_latparam/v_ysize1/2)"

variable sz equal "ceil(v_z/v_latparam/v_zsize1/2)"

variable lsx equal "ceil(2*v_x/v_latparam/v_xsize1/2)"

variable lsy equal "ceil(2*v_y/v_latparam/v_ysize1/2)"

variable lsz equal "ceil(2*v_z/v_latparam/v_zsize1/2)"

lattice fcc ${latparam} orient x ${1x1} ${1x2} ${1x3} orient y ${1y1} ${1y2} ${1y3} orient z ${1z1} ${1z2} ${1z3} spacing ${xsize1} ${ysize1} ${zsize1}

region box block -${lsx} ${lsx} -${lsy} ${lsy} -${lsz} ${lsz} units lattice

create_box 1 box

pair_style eam/alloy

pair_coeff * * Cu_mishin1.eam.alloy Cu

### Grain size and filling parameters here ###

variable grainsize equal 40 #set approximate grain size

variable grainsig equal 10 #set some expected standard deviation of grain size

variable occupiedvolume equal 0.365 #what fraction of space will be occupied

variable expected equal v_occupiedvolume*v_sx*v_sy*v_sz*v_fcc*2*2*2

variable grains equal 0 #initialize grain counter

run 0

variable counter loop 5000000 #should be more than enough except for extreme circumstances. Also used to seed the next random grain.

label loop

variable randomr equal "normal(v_grainsize,v_grainsig,v_counter)"

variable holdr equal ${randomr}

if "${holdr} < 0" then "jump SELF negative" #else " "

variable randomx equal "random(-v_x/2+v_holdr,v_x/2-v_holdr,v_counter)"

variable randomy equal "random(-v_y/2+v_holdr,v_y/2-v_holdr,v_counter)"

variable randomz equal "random(-v_z/2+v_holdr,v_z/2-v_holdr,v_counter)"

variable rax1 equal "random(0,1,v_counter)"

variable ray1 equal "random(0,1,v_counter+1)"

variable raz1 equal "random(0,1,v_counter+2)"

variable rax2 equal "random(0,1,v_counter+3)"

variable ray2 equal "random(0,1,v_counter+4)"

variable raz2 equal "random(0,1,v_counter+5)"

variable drax equal "v_rax1-v_rax2"

variable dray equal "v_ray1-v_ray2"

variable draz equal "v_raz1-v_raz2"

variable theta equal "random(0,62.4,v_counter)" #~62.4 for cubic misorientations, more realistic distribution could be implemented here

variable holdx equal ${randomx}

variable holdy equal ${randomy}

variable holdz equal ${randomz}

region grain sphere ${holdx} ${holdy} ${holdz} ${holdr} units box

group grain region grain

variable overlap equal "count(grain)" #this variable will count atoms already in this space defined by the grain region

if "${overlap} > 0.5" then "jump SELF foundone" else " " #skip grain creation if an atom is found

create_atoms 1 region grain

variable grains equal ${grains}+1 #counts how many grains/crystallites have been created

group grain delete #clear grain group and region

group grain region grain

displace_atoms grain rotate ${holdx} ${holdy} ${holdz} ${drax} ${dray} ${draz} ${theta} units box

#write_data data.${name}.${counter} #data file can be written to view grain placement process

label foundone

variable crit equal atoms/v_expected*100

variable current equal atoms

print "Created ${grains} grains ; Current Atoms = ${current} ; % Complete = ${crit}" #The process should slow down considerably as %->100%

if "${crit} > 100" then "jump SELF after" else " "

region grain delete

group grain delete

label negative

print "${counter}"

print " "

next counter

jump SELF loop

label after

print "Done making initial grains"

#######################################################

# SINTERING PROCESS #

#######################################################

#fix 3 all box/relax aniso 0.0 vmax 1.0e-4 nreset 100

thermo 25

thermo_style custom atoms step temp pe etotal pxx pyy pzz lx ly lz

minimize 0.0 1.0e-6 1000 1000 #initial minimize to eliminate any overlapping atoms if they exist for some reason

#unfix 3

change_box all boundary s s s

run 0

change_box all boundary p p p

minimize 1e-15 1.0e-15 1000 10000 #make sure that creating periodic boundaries did not create any issues

compute 3 all pe/atom

compute 4 all stress/atom NULL ke pair

dump ini all custom 1000 dump.eve_${name}.* id type x y z vx vy vz c_3 c_4[1] c_4[2] c_4[3] #dump all atoms

dump sli all custom 500 dump.sli_${name}.* id type x y z vx vy vz c_3 c_4[1] c_4[2] c_4[3] #and thin middle slice for viz

dump_modify sli thresh y > -10

dump_modify sli thresh y < 10

#Current temperatures and pressures are initial guesses for a Cu (Tm ~ 1,358 K) nanocrystalline system.

#final anneal temperature ~0.8Tm (1100K). Mid anneal temp ~0.65Tm (875K).

#applied compressive pressures 0.5 GPa and 1 GPa

variable initT equal 300

variable midT equal 875

variable highT equal 1100

variable midP equal 5000

variable highP equal 10000

reset_timestep 0

velocity all create 600 12345

#ramp up the up the temperature

fix npt all npt temp ${initT} ${midT} 0.075 aniso 0 0 .075

run 20000

unfix npt

#hold temp and ramp up pressure

fix npt all npt temp ${midT} ${midT} 0.075 aniso 0 ${midP} .075

run 20000

unfix npt

#bring temperature up to 0.8 Tm

fix npt all npt temp ${midT} ${highT} 0.075 aniso ${midP} ${midP} .075

run 20000

unfix npt

#bring pressure up to 1 GPa

fix npt all npt temp ${highT} ${highT} 0.075 aniso ${midP} ${highP} .075

run 20000

unfix npt

#dwell at high pressure and temperature

fix npt all npt temp ${highT} ${highT} 0.075 aniso ${highP} ${highP} .075

run 60000

unfix npt

#ramp down both pressure and temp

fix npt all npt temp ${highT} ${midT} 0.075 aniso ${highP} ${midP} .075

run 20000

unfix npt

#pressure and temp brought to target values

fix npt all npt temp ${midT} ${initT} 0.075 aniso ${midP} 0 .075

run 20000

#final dwell at equlibrium target P & T

fix npt all npt temp ${initT} ${initT} 0.075 aniso 0 0 .075

run 20000

unfix npt

#output data file

write_data data.${name}_init

print "Jobs done"