The problem with sputtering

[artemE4]

hi!
I’m trying to simulate sputtering of silicon with ~300 argon atoms.
I set the NVE relaxation for the area of the Ar atom deposit and the NVT relaxation for the rest of the material. The bottom layer is fixed.
As the Ar count increases, the energy in the NVE ensemble area increases rapidly, and the atoms fly apart. What could be the problem?

# ========================================================================================================

# Regions and groups
region          silicon1  block 5 100 5 100 10.01 ${zbox_b} units box
group           silicon1 region silicon1

# =============================================
region             nosehover_left block 0 4.9 0 INF 10.01 ${zbox_b} units box
region             nosehover_right block 99.9 120 0 INF 10.01 ${zbox_b} units box

region             nosehover_front block 0 INF 0 4.9 10.01 ${zbox_b} units box
region             nosehover_back block 0 INF 100.1 120 10.01 ${zbox_b} units box

region             nosehover_bottom block 0 INF 0 INF 5.1 10 units box

# =============================================
group            thermo_left region nosehover_left
group            thermo_right region nosehover_right
group            thermo_front region nosehover_front
group            thermo_back region nosehover_back
group            thermo_bottom region nosehover_bottom

# =============================================
fix             thermo_left_group thermo_left nvt temp 300 300 0.01
fix             thermo_right_group thermo_right nvt temp 300 300 0.01
fix             thermo_front_group thermo_front nvt temp 300 300 0.01
fix             thermo_back_group thermo_back nvt temp 300 300 0.01
fix             thermo_bottom_group thermo_bottom nvt temp 300 300 0.01

fix             1 silicon1 nve
fix             3 layer nvt temp 0.1 0.1 0.01 
velocity        layer set 0 0 0

[Germain]

Hi @artemE4,

Without proper description of your deposition procedure, model and objective, it is impossible to provide meaningful answer. There are 7 different integrators in your system, some of them using extended Hamiltonians. Also, it is experimentally known that some adatoms can be kicked out of the deposition surface depending on the deposition procedure (routine PVD for example).

There are several sputtering simulation procedures in the literature. And I’ve never seen one looking like yours.

[akohlmey]

@artemE4

What could be the problem?

How should anybody know without seeing your entire input deck and being able to reproduce?

@artemE4

in the NVE ensemble area

There is no such thing. You may be referring to a part of your system that is using fix nve, but a statistical mechanical ensemble can only refer to your entire system and yours is for obvious reasons representing none of the well know ensembles. Please also note that groups do not change their members, even if the original conditions are no longer met unless they are defined as dynamic groups, which has a whole bunch of complications. For your kind of system, this should not be a problem, since you would only want to apply a (dissipative) thermostat to mimic the exchange of kinetic energy with the bulk.

@artemE4

As the Ar count increases, the energy […] increases rapidly,

That may be the correct behavior. I would depend on many factors like how well your system has been equilibrated, how rapidly to thrust Ar atoms at the system, how much kinetic energy those atoms have and how you set up thermostatting in general. What you show in the few input lines you quote makes no sense to me.

There have been several discussions about how to simulate a sputtering experiment with LAMMPS, thus you should search the forum for them and learn from them.

[artemE4]

units           metal
dimension       3
boundary        p p m
atom_style      atomic
atom_modify     map hash

variable        zbox equal 20 # * 5.44Ang
variable        zbox_b equal ${zbox}*5.5  
 
variable        enAr equal 1000 #eV

variable        Ar_dep1 equal ${zbox_b}+5
variable        Ar_dep2 equal ${zbox_b}+10

variable        Sput_count equal ${zbox_b}+10
variable        Sput_name string yield_sput_${enAr}

variable        Ar_N equal 300
variable        Ar_timestep equal 3000

variable        Dump_name string Sput_${enAr}eV_therm
variable        Dump_vac_name string vac_${enAr}eV_therm

variable        tstep equal 0.001		# picosecond
variable        rn equal ${Ar_N}*${Ar_timestep}

# regions and groups ================================================================================
lattice         diamond 5.44
region          simbox block 0 20 0 20 0 25
create_box      2 simbox
region          silicon block 0 20 0 20 0 ${zbox}
create_atoms    1 region silicon
group           silicon region silicon

region          layer block INF INF INF INF -2 5 units box  
group           layer region layer

# Potentials ======================================================================================
pair_style      hybrid/overlay sw lj/cut 5.0 zbl 0.1 5.0
pair_coeff      * * sw Si.sw Si NULL
pair_coeff      1 2 zbl 14.0 18.0
pair_coeff      2 2 lj/cut 0.0103 3.405
newton          on

# masses ==========================================================================================
mass            1 28.086
mass            2 39.948


# Minimization and  initial heat ==================================================================
min_style cg
minimize 1.0e-4 1.0e-6 1000 10000

velocity        silicon create 100 12345 
fix             1 silicon nvt temp 300.0 300.0 0.05
fix             2 layer nvt temp 1 0.01 0.01 
velocity        layer set 0 0 0

thermo 50
thermo_style    custom step temp etotal ke 

timestep        0.001
run 50
unfix 1
unfix 2

reset_timestep  0


# ========================================================================================================

# Regions and groups
region          silicon1  block 10 100 10 100 10.01 ${zbox_b} units box
group           silicon1 region silicon1

# =============================================
region             nosehover_left block 0 9.99 0 INF 20.01 ${zbox_b} units box
region             nosehover_right block 100.1 INF 0 INF 20.01 ${zbox_b} units box

region             nosehover_front block 0 INF 0 9.99 20.01 ${zbox_b} units box
region             nosehover_back block 0 INF 100.1 INF 20.01 ${zbox_b} units box

region             nosehover_bottom block 0 INF 0 INF 5.1 20 units box

# =============================================
group            thermo_left region nosehover_left
group            thermo_right region nosehover_right
group            thermo_front region nosehover_front
group            thermo_back region nosehover_back
group            thermo_bottom region nosehover_bottom

# =============================================
fix             thermo_left_group thermo_left nvt temp 300 300 0.01
fix             thermo_right_group thermo_right nvt temp 300 300 0.01
fix             thermo_front_group thermo_front nvt temp 300 300 0.01
fix             thermo_back_group thermo_back nvt temp 300 300 0.01
fix             thermo_bottom_group thermo_bottom nvt temp 300 300 0.01

fix             1 silicon1 nve
fix             3 layer nvt temp 0.1 0.1 0.01 
velocity        layer set 0 0 0


# Ar velocity ============================================= =============================================
variable        ar_vel equal sqrt(2*${enAr}*1.60219e-19/(39.948*1.66054e-27))/100
variable        ar_energy equal (39.948*(${ar_vel}/100)^2)/2
print           "Argon velocity = ${ar_vel} A/ps"
print           "Argon energy = ${ar_energy} eV"



# Ar deposit ============================================= ===============================================
region          ardeposit block 10 100 10 100 ${Ar_dep1} ${Ar_dep2} units box
group           argon region ardeposit
fix             argon argon deposit ${Ar_N} 2 ${Ar_timestep} 12345 region ardeposit vz -${ar_vel} -${ar_vel} units box
variable        deposit_argon equal count(argon)
fix             ar_nve argon nve



# Neighbors ============================================= ===============================================
neighbor         1 bin
neigh_modify     every 5 delay 0 check no
neigh_modify     binsize 0.0
neigh_modify     one 50000
neigh_modify     page 500000

# Sput count ============================================= ===============================================
region          out_of_surface block INF INF INF INF ${Sput_count} INF units box
group           sputtered dynamic silicon region out_of_surface
variable        sputtered_atoms equal count(sputtered)


# Sput out ============================================= =================================================

fix sputter_check all print ${Ar_timestep} "Step: $(step), Sputtered atoms: ${sputtered_atoms}, Deposit Argon: ${deposit_argon}" file ${Sput_name}.txt screen no

# Voronoi analysis and other computes ====================== =============================================
compute         1 argon ke  
compute         2 silicon voronoi/atom occupation
compute         3 silicon voronoi/atom
variable        vac_pos atom "(c_2[1] == 2 || c_2[1] == 0) && (c_3[1] > 21)"

compute         def_ke silicon ke/atom

dump            vac all custom 500 ${Dump_vac_name}.dump id type x y z c_2[1] c_2[2] c_3[1] c_def_ke
dump_modify     vac thresh v_vac_pos == 1.00

dump            sput   all custom 500 ${Dump_name}.dump id type x y z c_3[1]

# thermo style ============================================= =============================================
thermo          ${Ar_timestep}
thermo_style    custom step temp etotal c_1 v_sputtered_atoms v_deposit_argon
thermo_modify   lost ignore

# Run ============================================= =======================================================
timestep        ${tstep}
run             ${rn}

@Germain

There are several sputtering simulation procedures in the literature. And I’ve never seen one looking like yours.

hi! I have been studying the literature, looking for similar tasks and have seen similar techniques like mine…

@akohlmey

There is no such thing.

oh, thanks. i will be more precise in my wording.

@akohlmey

I would depend on many factors like how well your system has been equilibrated, how rapidly to thrust Ar atoms at the system, how much kinetic energy those atoms have and how you set up thermostatting in general.

I studied the relaxation time of the defective region after implantation. I set a deposit time greater than the relaxation time

[artemE4]

@akohlmey

There have been several discussions about how to simulate a sputtering experiment with LAMMPS, thus you should search the forum for them and learn from them.

Yes, I have studied them, however there are very specific problems there.
I initially started with the “sput” example in older versions of LAMMPS, then studied the literature on sputter modeling

[akohlmey]

@artemE4

have seen similar techniques like mine…

Sorry, but considering the purpose and functionality of Nose-Hoover thermostats, I can hardly believe that any serious publication would promote a use case like you are showing.

@artemE4

however there are very specific problems there.

Please note, that you are asking for advice without providing sufficient context and in this statement again you are not explaining crucial information. Yet you are dismissing any suggestions given based on the information provided. How can you expect meaningful suggestions, when you make it impossible to provide the specific help you are seeking? You will only alienate people that answer in good faith. Please keep in mind that this forum depends on people volunteering their time to answer questions and provide advice.