I am working on a problem to investigate biaxial stress in a system where a thin film is grown pseudomorphically on top of a substrate. A simple way to describe pseudomorphic growth is that the atoms of the thin film material (say material A) align in the x-y plane with the atoms of the substrate (say material B) but since the lattice constants of A and B are dissimilar, atoms of the thin-film (material A) elongate in the z-dimension. So in a situation where both material A and B have cubic crystal structure, after pseudomorphic growth thin-film atoms will elongate in the z-direction to give a tetragonal shape. There is no stress in the normal direction (z-dimension) but a strain in z-direction. Biaxial stress exists in x and y-directions.
I have few question related to setting up this simulation in LAMMPS.
Since the atoms of the thin-film are allowed to move in the z-direction, what kind of boundary conditions should I specify ? I read the documentation and am thinking of periodic in x and y and probably “free” or “shrink-wrapped” in z-direction.
What is the correct ensemble to simulate this scenario ? I specified NPT in my input file with isobaric pressure of 0 but the code returns an error - “ERROR: Cannot use fix nvt/npt/nph on a non-periodic dimension” when I specify shrink-wrapped boundary conditions in the z-direction.
any input or experience with simulating biaxial stress will be very helpful for me.
I am working on a problem to investigate biaxial stress in a system where a
thin film is grown pseudomorphically on top of a substrate. A simple way to
describe pseudomorphic growth is that the atoms of the thin film material
(say material A) align in the x-y plane with the atoms of the substrate (say
material B) but since the lattice constants of A and B are dissimilar, atoms
of the thin-film (material A) elongate in the z-dimension. So in a
situation where both material A and B have cubic crystal structure, after
pseudomorphic growth thin-film atoms will elongate in the z-direction to
give a tetragonal shape. There is no stress in the normal direction
(z-dimension) but a strain in z-direction. Biaxial stress exists in x and
y-directions.
I have few question related to setting up this simulation in LAMMPS.
1. Since the atoms of the thin-film are allowed to move in the z-direction,
what kind of boundary conditions should I specify ? I read the
documentation and am thinking of periodic in x and y and probably "free" or
"shrink-wrapped" in z-direction.
This should be fine.
2. What is the correct ensemble to simulate this scenario ? I specified
NPT in my input file with isobaric pressure of 0 but the code returns an
error - "ERROR: Cannot use fix nvt/npt/nph on a non-periodic dimension" when
I specify shrink-wrapped boundary conditions in the z-direction.
You can use NPT, but only in the xy dimensions, which the fix npt command
lets you do. It doesn't make sense to control pressure
in z, since you effectively have a free surface. The system will
expand,/contract
in z freely on its own.
The materials in my simulation are Si (substrate) and SiGe (thin-film). I have access to SiCGe.tersoff file which comes bundled with LAMMPS (at least with the July 2 2010 version I have) and the documentation says that this can in principle be used to descibe SiGe potential parameters. Can you help me understand how to specify the parameters ?
I have the following in the input file but am not sure
pair_style tersoff
pair_coeff * * SiCGe.tersoff Si Ge
do we also need long range Coulomb interaction ?
since the substrate layer is in general much thicker than the thin-film, I thought of setting the force on atoms of the substrate layer to 0.0 using the following fix command. Does this make sense ?
1. The materials in my simulation are Si (substrate) and SiGe (thin-film).
I have access to SiCGe.tersoff file which comes bundled with LAMMPS (at
least with the July 2 2010 version I have) and the documentation says that
this can in principle be used to descibe SiGe potential parameters. Can you
help me understand how to specify the parameters ?
The params are in the SiCGe file - what do you need to specify?
I have the following in the input file but am not sure
pair_style tersoff
pair_coeff * * SiCGe.tersoff Si Ge
That is correct syntax.
do we also need long range Coulomb interaction ?
I doubt it for SiGe. But to be honest, you sound like
you are just guessing at how to model this system. I highly
recommend you find paper(s) in the literature that describe
how to model SiGe, read them, and then figure out how
to do what they describe in LAMMPS.
2. since the substrate layer is in general much thicker than the thin-film,
I thought of setting the force on atoms of the substrate layer to 0.0 using
the following fix command. Does this make sense ?
I am trying to simulate a situation of heteroepitaxial growth of a thin-film on top of a substrate. Thin film is grown pseudomorphically on top of a substrate which means that the atoms of the thin film material (say material A) align in the x-y plane with the atoms of the substrate (say material B) but since the lattice constants of A and B are dissimilar, atoms of the thin-film (material A) elongate in the z-dimension. So in a situation where both material A and B have cubic crystal structure, after pseudomorphic growth thin-film atoms will elongate in the z-direction to give a tetragonal shape. There is no stress in the normal direction (z-dimension) but a strain in z-direction. Biaxial stress exists in x and y-directions. This scenario is valid till a critical thickness of the thin-film is reached and then the dislocations begin to appear. I am concerned with a situation much below the cirtical thickness.
I had asked few questions related to this problem earlier but the results from LAMMPS have so far not been in agreement with what is observed experimentally.
So I will try to rephrase the questions and ask again.
My simulation box has periodic boundary conditions in x and y and “free” in z-direction. This is to simulate a situation the thin-film only moves in the z-direction. The box height in z-direction is about 2-3 lattice spacings more than maximum z-coordinate in the starting geometry. I start with a situation where layers A and B have the same lattice constants. Since the substrate atoms do not move I am using “fix setforce” command to zero the force on all atoms of layer B (substrate).
I specified “fix thin-film npt temp 300 300 10. x 0.0 0.0 10. y 0.0 0.0 10.” in my input file. There is no pressure control in the z-direction (direction of thin-film expansion). Moreover, the fix is applied to only atoms of the thin-film layer (I tried once will fix npt for all atoms also but the results were similar). However, LAMMPS shows me that the system expands in all directions (x, y and z). I do not know to avoid this behavior since I would like atoms to move in only the z-direction which is consistent with experiments.
another option which I explored without success was to use the " fix lineforce" command to set the forces on atoms in the x-y plane to be zero.
I would like to add that if I use the NVE ensemble (fix thin-film nve) then the results are in reasonable agreement with experiments. The atoms move only in the z-direction. However, I cannot unerstand the physcial basis of why this is happening.
any input or experience with simulating similar behavior will be very helpful for me.
If your system wants to expand in x and y when running NPT,
it is because it has a positive Pxx and Pyy pressure. This likely
has nothing to do with adding atoms and growing in the z
direction. If you have chosen a lattice constant (in x,y) or
a potential that is not at a minimum at that lattice
constant of your initial box, then you will get a positive
pressure and it will expand. Presumably it is exapnding
to a new equilibrium size where the Pxx and Pyy pressures
are 0.0.
I have a system composed of two regions – Si and Ge. I would like to place the Ge atoms on the Si lattice sites initially and then allow them to relax during a dynamics run. bounrday conditions are specified to be “f p p”.
the Ge region can expand in the x-direction
output from lammps log file looks like this
#initialization
units metal
boundary f p p
atom_style atomic
#atom definition
lattice diamond 5.432 orient x 1 0 0 orient y 0 1 1 orient z 0 -1 1
Lattice spacing in x,y,z = 5.432 7.68201 7.68201
lattice spacings
region mybox block -1.0 5.0 -2.0 2.0 -3.0 3.0
create_box 2 mybox
Created orthogonal box = (-5.432 -15.364 -23.046) to (27.16 15.364 23.046)
1 by 1 by 2 processor grid
region Ge_region block 0.1 4 -2.0 2.0 -3.0 3.0
region Si_region block -1 0.0 -2.0 2.0 -3.0 3.0
create_atoms 2 region Ge_region
Created 1476 atoms
create_atoms 1 region Si_region
Created 456 atoms
…
, inteaction is specified using Tersoff potential
…
then I specify “fix nve” followed by “fix langevin” to run for 5000 timesteps.
however, when I look at the atom coordinates dumped at the end of 5000 steps and check for the maximum y-coord of any atom it gives me
15.520
and minimum y-coord of any atom is
-15.450
I am not able to understand how can the y-coord of any atom be more than the box dimension in the y-direction (-15.364 to 15.364) ? it looks as if the atoms expanded outside the box.
however, when I look at the atom coordinates dumped at the end of 5000 steps
and check for the maximum y-coord of any atom it gives me
15.520
and minimum y-coord of any atom is
-15.450
I am not able to understand how can the y-coord of any atom be more than the
box dimension in the y-direction (-15.364 to 15.364) ? it looks as if the
atoms expanded outside the box.
this is a FAQ!
periodic boundaries are only enforced on
neighborlist rebuilds. mind you, the physics
of the system is totally unaffected by this.
what you see is normal and expected. if you
want to avoid this, you have to enforce rebuilding
the neighborlist in every step, but that will
impact your simulation performance. better to
wrap atoms back into the principal unitcell
in postprocessing.
Did you lose the atoms when they moved outside
the y boundary. This is what will happen on
the next reneighbor with "f" boundary conditions.
If you don't want that, then use "s". As Axel says
it is normal for atoms to be slightly outside the box
until the next reneighboring.
I was wondering if MD simulations are the correct approach to study critical thickness for epitaxial layers. I have not seen many papers in recent years on this type of study.
well-known theory of Matthews-Bakeslee underestimates critical thickness of modern day epitaxial layers (J. Appl. Physc, 97, 043519, 2005) and I did not come across any papers reporting MD simulations in the last few years.
if anyone can give pointers to work done in this area or share their experience it will be very useful.
What changes when you are less vs more than the critical
thickness and is the associated phenomena at the right
time and length scale to study with MD?
In case of heteroepitaxial layers, deposited films relax by forming dislocations once the film height exceeds cirtical thickness. However, I am not sure what is the time scale for dislocations to appear.