Introduction

Se and Te crystallize in the same hexagonal structure, which has atoms per cell. The Se (Te) s state is very deep and participates only weaking in the bonding. Each atom has three neighbors at approximately 90o angles. This allows each of the three p orbital to form strong covalent bonds with each of its neighbors. Thus this system may be considered as a covalently bonded spn bonded semiconductor, with n→∞.

The angle between Se atoms is not given by symmetry; it is approximately 90o, but not exactly so. Thus there a degree of freedom in the site positions not determined by symmetry. This tutorial uses molecular statics to find the minimum-energy value of this free parameter.

Tutorial

1. Building the input file

For the structure, copy in the contents of the box below into file init.se. Note the parameter u is the undetermined parameter the molecular statics algorithm will find by minimizing the energy.

--ctrl=ctrl Tells blm to write directly to the actual input file ctrl.se. By default it writes to actrl.ext so as not to overwrite the input file.

--pbe Use the PBE functional

--wsitex=site Write the site file with site positions expressed in units of PLAT

--molstat Set up a category for molecular statics; see below

There are two quantities you must supply yourself, GMAX and the k-point mesh. For the former, lmfa can determine GMAX for you once the basis set is known. Since lmfa also supplies a basis set provided HAM_AUTOBAS_LMTO is appropriately set, it will knows what the basis set is and tell you what GMAX should be. Do the following

$ lmfa se

At the end of the standard output you should find

FREEAT: estimate HAM_GMAX from RSMH: GMAX=5.4

This tells you what GMAX to use. For the k mesh, we will use 6 k divisions (216 points in the full Brillouin zone). Rerun both blm and lmfa, now specifying these numbers and also telling lmfa to write the basis set directly to the basp.se. By default it writes to basp0.ext so as not to overwrite whatever is already present.

Note: Unlike the total energy forces are not variational in the density; thus they are more difficult to converge.

That is, the forces do not converge to their self-consistent value as , where and are the current and self-consistent densities. Rather it converges as . However, an estimate for the correction can be performed. One way is to find a contribution to the force that would arise if the free-atomic density were dragged along with the nucleus. That is the ansatz used here (this ansatz is used when HAM_FORCE=1). An alternative is to use HAM_FORCE=12.

Inspect the output from the first iteration. You should find this table

Without the correction term, the would be 33.42−876.65 mRy/au. But with this ansatz, the force is only 33.42 mRy/a.u.. So the correction term is huge — much larger than the force itself. The corrected force, even from the initial [Mattheis construction]/tutorial/lmf/si_calculation/#self-consistency), is not so far from the self-consistent one (see below).

Now look at the equivalent output from the last iteration. You should find

3. Molecular statics

lmf can relax the site positions to their equilibrium point where the forces vanish (molecular statics).

The parameters that control how this relaxation proceeds are given in as tokens in the DYN_MSTAT tag. Notice the last lines of the input file. To see whether the proprocessor includes this line or not, compare the output of the following two commands.

$ lmf se --showp
$ lmf se -vminx=t --showp

In the second command the preprocessor includes the DYN category in the input file.

Note:: you should consider values for tokens in DYN_MSTAT. Some default values were chosen, but it is very difficult to make a general algorithm that relaxes atoms. In this case it uses a Broyden scheme (a kind of Newton-Raphson algorithm), with an initial step size of 0.01. This a good choice in the present context, but in other contexts it may not be.

To relax the structure, do

$ lmf se -vminx=t --wpos=pos

lmf iterates the density to self-consistency; when it is reached it shifts the nuclear positions using the Broyden scheme.

Rather than build the density from a Mattheis construction, it takes the self-consistent density, adds a Mattheis construction calculated from the new positions and subtracts the same construction built from the starting positions. This provides a reasonable ansatz for the density at the new positions.

When self-consistency is reached with these positions the forces read:

This is indeed a very small force. The final site positions are printed to stdout and also to file pos.se (because of --wpos=pos).

The x coordinate of the first atom in pos.se should be 0.2346. Note that the coordinates last output by lmf are slightly different: these correspond to what would be the updated coordinates if the molecular statics would proceed another step.

It is interesting to see how the total energy changes with each relaxation step. Do the following: