periodic boundary conditions


I am trying to reproduce epitaxial growth of germanium on silicon as discussed in PRB vol. 34 pp 5621 (1986) by Chris Van de Walle.

the paper talks about superlattice structures of Si and Ge layers. In the limit when a thin layer of Ge is deposited on a thick layer of Silicon the “in-plane” lattice constant of Ge would be equal to lattice constant of Si (5.43A). this corresponds to a ~4% strain compared to its bulk lattice constant (a_ge = 5.65A)

As a consequence of this in-plane stress, the thin Ge layer distorts tetragonally, i.e. out of plane lattice constant become 5.82A. this is about 7.2% tensile strain referenced to silicon lattice spacing.

The paper gives analytic formulas to calculate the tetragonal distortion in terms of elastic constants of silicon and germanium.

I am trying to simulate this scenario with LAMMPS with molecular statics using Tersoff potential.

  1. I calculate elastic constant using Adian Thompson’s ELASTIC input script.

  2. I wrote a simple code which extracts strain tensor from LAMMPS xyz files using intial and final “minimized” geometry.

  3. I set up the problem with a two layer system. bottom layer is Silicon; top layer is Germanium. I choose periodic boundary conditions in the x, y directions and “shrink-wrapped” in the z-direction. additionally, since I do not want to include a lot of silicon atoms, I use setforce to zero forces in x-y-z directions for all silicon atoms (idea is to simulate to a very thick silicon which experience negligible distortion from top Germanium layer).

LAMMPS output from MS simulations indicate that the Germanium atoms do show negligible in-plane strain (i.e. they move very little from the intialized positions in the x-y plane) and a tetragonally distorted out-plane lattice spacing. the strain calculated from MD simulations is slightly lower than 7.2% value I quoted above.

I notice that the strain tensor elements using my method depend on the choice of atoms I “pick” for calculating strain. If my code selects atoms at the edge of x-y boundaries then the strain calculation give absurd answers. once I leave out the 1-2 atomic layers at the x-y edges of the simulation box, I get sensible answers. I have repeated simulations for different thickness of Ge layer and results consistently show the wisdom of leaving out the atoms located at the x-y edge of Germanium layer.

one may hypothesize possible reasons for disagreement between LAMMPS and analytic formulation (as discussed in the paper). Tersoff potential could be a source of problem including other sources errors. elastic constant calculated using Tersoff potential do not agree very well with expt.

however, I repeated above simulations with one modification – I used “setforce” command to zero out forces in the x-y plane (in-plane) for the Germanium atoms. what surprised me was that I get better agreement with the analytic formulas.

so I am left wondering if there is an issue with using periodic boundary conditions in this heterogenous system. I asked a friend who has expertise in continuum stress simulations and he told me that they have something called “reflecting boundary conditions” which imply that the displacement vectors at the x-y edges (if he were to setup analogous continnuum simulation) would be zero.

periodic boundary conditions with setforce 0.0 0.0 NULL seem to be emulating the “reflecting boundary conditions” situation (and hence better agreement with continuum like formulation). another interesting bit is that with this setup, I do not have to disregard atoms at the edge of x-y boundaries for my strain calculator.

But, I am perplexed why periodic boundary conditions by themselves are not sufficient ?

I know this is a very long read for anyone but if you have inputs please let me know.

thank you.


Maybe Aidan has some ideas. I am unclear if you
are periodic in xy, how you are making the box size
consistent with the natural periodicity of both the
Si and Ge lattice constants. If you are saying it
is consistent with one, and the other should go into
a strained mode to match the other, then that is fine.
You shouldn't see anomalies at the boundary if the
system has fully relaxed. You can verify that the
spacing of strained atoms in the middle of the
box is the same as near the boundaries by dumping
the coords.


I choose periodic boundary conditions in x, y and the box size is consistent with the “substrate layer” which is Silicon in my test case. z-direction is perpendicular to the plane of interface and the boundary condition in that direction are “shrink-wrapped”. the system relaxes in z-direction.


The "absurd" strain values may be caused by two different things, that I
can think of:

1. The germanium layering nucleates dislocations to relax the lateral
strain. Away from dislocations, the lateral strain will be close to zero.
2. You may be calculating strain incorrectly.

Either way, I recommend that you create some visualization of the
different cases in order to better understand what is going on. Atoms are
notorious for not obeying our idealized mental pictures.