Trying to form Si 001 surface which shows surface reconstruction

LAMMPS LAMMPS Mailing List Mirror
1

I have fixed the bottom few layers, set the z BC non periodic, used tersoff potential, minimized energy. I even tried to relax the atoms in x and y directions using fix box/relax command.

But I don’t see why am I not getting 1x2 reconstruction of Si 001. There is something very basic I am missing. I am a beginner in MD simulations. Can anyone help me with that?

---------- Initialize Simulation ---------------------

clear
units metal
dimension 3
boundary p p s
atom_style atomic
atom_modify map array

---------- Create Atoms ---------------------

diamond lattice, lattice constant 5.5 AA

lattice diamond 5.5
region box block 0 5 0 5 0 10 units lattice
region boxfix block 0 5 0 5 0 2 units lattice
create_box 1 box

lattice diamond 5.5 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
mass 1 28.0855
create_atoms 1 box
#replicate 2 2 2
group fixed region boxfix

---------- Define Interatomic Potential ---------------------

pair_style tersoff
pair_coeff * * Si.tersoff Si
neighbor 4.0 bin
neigh_modify delay 10 check yes

#------------ dump xyz file with atoms and bonds --------

read with JMol

dump 1 all xyz 10000 sisurf.xyz
dump_modify 1 element Si

---------- Define Settings ---------------------

compute eng all pe/atom
compute eatoms all reduce sum c_eng

---------- Run Minimization ---------------------

reset_timestep 0
#fix 1 all box/relax x 0.0 y 0.0 vmax 0.001

fix 1 fixed setforce 0.0 0.0 0.0

thermo 10

thermo_style custom step pe lx ly lz c_eatoms
min_style cg
minimize 1e-25 1e-25 5000 10000

variable natoms equal “count(all)”
variable teng equal “c_eatoms”
variable length equal “lx”
variable ecoh equal “v_teng/v_natoms”

print “Total energy (eV) = {teng};" print "Number of atoms = {natoms};”
print “Lattice constant (Angstoms) = {length};" print "Cohesive energy (eV) = {ecoh};”

print “All done!”

2

Perhaps you should ask yourself the question of at what conditions
does surface reconstructions happen? Do Si surface reconstructions
occur at 0 K?

Regarding the input script, is it necessary to fix the bottom layers?
Is relaxing atoms only in the x and y directions physical?

Ray

3

Increase the box size in z direction (ie create empty space above your Si) and anneal it instead of minimize it.

Thats how I did it (I’m not sure whether anneal vs minimization matters or not)

Alexandre

4

Thanks Ray and Alex.
I took the pointers on increasing temperature and creating a buffer empty space on top.
I tried the simulation w/wo fixing the bottom layers. The problem I think is that the position of atoms is not getting updated with every iteration and are the same as initialized even though I used fix move or fix nve commands. I am finding it particularly hard to relax the atoms for a non periodic BC.

I have gone through the previous threads on this topic, but it was not of much help. Can someone tell me how to relax a non periodic surface? My goal is to relax the surface to 1x2 for Silicon 100.

Find minimum energy Si configuration

Jim Hannon 2013, based on Mark Tschopp, 2010

---------- Initialize Simulation ---------------------

clear
units metal
dimension 3
boundary p p s
atom_style atomic
atom_modify map array

---------- Create Atoms ---------------------

diamond lattice, lattice constant 5.5 AA

lattice diamond 5.4
region boxfull block 0 5 0 5 0 20 units lattice
region boxm block 0 5 0 5 2 12 units lattice
region boxf block 0 5 0 5 0 2 units lattice
region box union 2 boxf boxm
create_box 1 boxfull

lattice diamond 5.4 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
mass 1 28.0855
create_atoms 1 box
group afix region boxf

group amov region boxm

#---- velocity---------------------
velocity amov create 500 123

---------- Define Interatomic Potential ---------------------

pair_style tersoff
pair_coeff * * Si.tersoff Si
neighbor 4.0 bin
neigh_modify delay 10 check yes

#------------ dump xyz file with atoms and bonds --------

read with JMol

dump 1 all xyz 10000 Si.xyz
dump_modify 1 element Si

---------- Define Settings ---------------------

compute eng all pe/atom
compute eatoms all reduce sum c_eng

---------- Run Minimization ---------------------

reset_timestep 0
fix myFreeze afix setforce 0.0 0.0 0.0

fix 1 amov move linear -0.3 NULL NULL units lattice
thermo 10

thermo_style custom step temp pe lx ly lz c_eatoms
min_style cg
minimize 1e-25 1e-25 5000 10000

variable natoms equal “count(all)”
variable teng equal “c_eatoms”
variable length equal “lx”
variable ecoh equal “v_teng/v_natoms”

print “Total energy (eV) = {teng};" print "Number of atoms = {natoms};”
print “Lattice constant (Angstoms) = {length};" print "Cohesive energy (eV) = {ecoh};”

print “All done!”

5

Surface reconstruction is a thermally activated process. You only
need normal MD run with appropriate thermalization, not minimization.
Minimize a system with initial velocity created make little sense, and
fix move makes even less. It is verified that Tersoff potential is
able to obtain a Si surface reconstruction.

Ray

6

Hi, so there are a couple of changes that I did to my script after reading more about lammps.
The surface is periodic in x and y and shrink wrapped in z.
I created 3 regions (bottom one fixed, middle one attached to a thermostat at 500K using fix temp/rescale, and the top surface layer are let to relax in nve ensemble using fix nve. I have also initialized the temperature of surface atoms to 500K). There is a buffer vacuum region above the surface ( Do I need the buffer region on the sides too?).

And now I am letting the atoms to thermalize (run) for 5000 timesteps ( Is this time too short?)

No desired results yet!

Form 2x1 SI 100 surface reconstruction

#------------------Initialize Simulation------------
clear
units metal
dimension 3
boundary p p s
atom_style atomic
atom_modify map array

#------------------Create atoms--------------------

---------- Create Atoms ---------------------

diamond lattice, lattice constant 5.43 AA

lattice diamond 5.43
region boxsim block 0 5 0 5 0 55 units lattice
region boxf block 0 5 0 5 0 2 units lattice
region boxm block 0 5 0 5 2 35 units lattice
region boxs block 0 5 0 5 35 50 units lattice
region box union 3 boxf boxm boxs
region boxd block 0 5 0 5 49.5 50 units lattice
region free union 2 boxm boxs
create_box 1 boxsim

lattice diamond 5.43 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
mass 1 28.0855
create_atoms 1 box
group afix region boxf
group amov region boxm
group asur region boxs
group afre region free
group adum region boxd

---------- Define Interatomic Potential ---------------------

pair_style tersoff
pair_coeff * * Si.tersoff Si
neighbor 4.0 bin
neigh_modify delay 10 check yes

#------------ dump xyz file with atoms and bonds --------

dump 1 all xyz 100 Si0823.xyz
dump_modify 1 element Si

dump 2 adum xyz 100 Sisurf0823.xyz
dump_modify 1 element Si

---------- Define Settings ---------------------

compute eng all pe/atom
compute eatoms all reduce sum c_eng

---------- Run Minimization ---------------------

velocity afre create 500 102486
reset_timestep 0
fix 1 afix setforce 0.0 0.0 0.0
fix 2 amov temp/rescale 10 500 500.0 100.0 1
fix 3 asur nve

thermo 10
thermo_style custom step temp pe lx ly lz c_eatoms

run 5000

variable natoms equal “count(all)”
variable teng equal “c_eatoms”
variable length equal “lx”
variable ecoh equal “v_teng/v_natoms”

print “Total energy (eV) = {teng};" print "Number of atoms = {natoms};”
print “Lattice constant (Angstoms) = {length};" print "Cohesive energy (eV) = {ecoh};”

print “All done!”

7

At this point it is probably best you talk to somebody local, that knows about MD and can fill you in on all the details that you seem to have little or no experience in. Trying to resolve this over email is very ineffective and time consuming.

Axel.

8

Why do you think 500K is sufficient?

If 5000 steps is too short, why not increase it?