Molecular Statics in Se

This tutorial demonstrates molecular statics mode in lmf : relaxation of site positions to their equilibrium point where energy is minimum and forces vanish. Elemental Se is taken as a representative material.

For the study and relaxation of a complex interface, see this tutorial.

Command summary
Introduction
Tutorial
Additional Exercises (optional)

Command summary

Start from file init.se.

blm --brief --ctrl=ctrl --wsitex --nk~met --molstat init.se
lmchk se --shell --angles
lmfa se
lmf se > out.lmfsc
lmf se -vminx=10 --wpos=pos > out.relax
lmchk se --rpos=pos --angles

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 90^o 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 spⁿ bonded semiconductor, with n→∞.

The angle between Se atoms is not given by symmetry; it is approximately 90^o, 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.

# Se and Te structures. From Wyckoff pg. 36.  Space group D_3^4 (C3(_1)21).
% const se=t
% ifdef se # For Se
% const a=8.234 cbya=1.136 u=.217
% char nam=Se
% else     # For Te
% const a=8.406 cbya=1.330 u=.269
% char nam=Te
% endif

LATTICE
   ALAT={a} PLAT=  1 0 0  -0.5 sqrt(3)/2  0  0 0 {cbya}
# This one is for Bi2Te3
SITE
    ATOM={nam} X=  {u}   0.000   0
    ATOM={nam} X= -{u}   -{u}   1/3
    ATOM={nam} X=  0.000  {u}   2/3

Note:: the basis positions are expressed in multiples of the lattice vector PLAT. (X= tells blm to interpret them that way, instead of in Cartesian coordinates.) But most of the Questaal output (e.g. positions in the relaxation step) print out basis positions in Cartesian coordinates, as multiples of alat.

We will use the PBE functional. Enter the following command:

blm --brief --ctrl=ctrl --wsitex --nk~met --molstat init.se

Switches have the following meaning:

--brief Tells blm make an input file of moderate verbosity, and autogenerate the basis set
--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
--nk~met Tells blm to use a default setting for the k mesh needed to integrate over the Brillouin zone.
--molstat Set up a category for molecular statics; it has the effect of adding a DYN category to the input file:

% ifdef minx
dyn     mstat[mode=6 hess=t xtol=.001 gtol=0 step=.010 nkill=0] nit={minx}
% endif

The first and third line are preprocessor instructions. The middle line will become part of the input and parsed by the preprocessor only if variable minx is defined with a nonzero value. minx should be an integer; it will be the argument to DYN_NIT, which is the maximum number of relaxation steps.

1.1 Bond lengths and angles

To confirm the that bond angles are approximately 90^o, run lmchk with the following switches

lmchk se --shell --angles

You should see this table for bond pairs:

 Shell decomposition for class Se        class   1  z=34
 shell   d     nsh csiz  class ...
   1  0.000000   1    1  1:Se
   2  0.533531   2    3  1:Se(2)
   3  0.796025   4    7  1:Se(4)
   4  0.845471   2    9  1:Se(2)
   5  1.000000   6   15  1:Se(6)
   ...

This establishes that there are three NN bonds. For bond angles, this table should appear

 Bond angles for site 1, species 1:Se
 shl1    d1    shl2    d2     cl1      cl2       angle(s) ...
   1  0.533531   1  0.533531  Se       Se        104.81
   1  0.533531   2  0.796025  Se       Se         95.52   95.52   99.86  159.67
   ...

Angles between the NN Se atoms is somwhat larger at 105^o, somewhat larger than the simple nearest-neighbor tight-binding model would suggest, indicating that the Se s states are not totally inert.

2. Self-consistency

Make the system self-consistent with:

lmfa se
lmf se > out.lmfsc

If you have compiled it in parallel mode you can do it a bit faster, e.g.

mpirun -n 6 lmf se

lmf should converge to self-consistency in 7 iterations.

You should see this table at the end:

 Forces, with eigenvalue correction
  ib              estatic                     eigval                      total
   1    866.59     0.00     0.00   -821.63     0.00     0.00     44.96     0.00     0.00
   2   -433.29  -750.49     0.00    410.82   711.55     0.00    -22.48   -38.93     0.00
   3   -433.30   750.49     0.00    410.82  -711.55     0.00    -22.48    38.93     0.00
 Maximum Harris force = 45 mRy/au (site 2)  Max eval correction = 0.1

The forces are in units mRy/au. They are all equal in magnitude, but rotated by 120^o angles, as follows from crystal symmetry.

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 $O(n-n^*)^2$ , where $n$ and $n^*$ are the current and self-consistent densities. Rather it converges as $O(n{-}n^*)$ . 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

 Harris correction to forces: shift in free-atom density
  ib         delta-n dVes             delta-n dVxc               total
   1                                                    875.50   -0.00   -0.00
   2                                                   -437.75 -758.21    0.00
   3                                                   -437.75  758.21   -0.00
 shift forces to make zero average correction:            0.00    0.00    0.00

It says that the correction term to the force is about 880 mRy/au. This is a very large term!

The correction table is soon followed by

 Forces, with eigenvalue correction
  ib              estatic                     eigval                      total
   1    151.58     0.00     0.00    -94.11     0.00     0.00     57.47     0.00     0.00
   2    -75.79  -131.28     0.00     47.06    81.51     0.00    -28.73   -49.77     0.00
   3    -75.79   131.28     0.00     47.06   -81.51     0.00    -28.73    49.77     0.00
 Maximum Harris force = 57.5 mRy/au (site 1)  Max eval correction = 875.5

Without the correction term, the would be 57.47−875.5 mRy/au. But with this ansatz, the force is only 57 mRy/a.u.. So the correction term is huge — much larger than the force itself. The corrected force, even from the initial Mattheis construction, is not so far from the self-consistent one (see below).

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

 Harris correction to forces: shift in free-atom density
  ib         delta-n dVes             delta-n dVxc               total
   1                                                      0.12   -0.00   -0.00
   2                                                     -0.06   -0.10   -0.00
   3                                                     -0.06    0.10   -0.00
 shift forces to make zero average correction:            0.00    0.00    0.00

The correction term has become very small; indeed the tighter the tolerance the density is converged, the smaller this term becomes.

Look for the table shown below in the output of the last iteration, when self-consistency has been reached. It should read:

 Forces, with eigenvalue correction
  ib              estatic                     eigval                      total
   1    866.59     0.00     0.00   -821.63     0.00     0.00     44.96     0.00     0.00
   2   -433.29  -750.49     0.00    410.82   711.55     0.00    -22.48   -38.93     0.00
   3   -433.30   750.49     0.00    410.82  -711.55     0.00    -22.48    38.93     0.00
 Maximum Harris force = 45 mRy/au (site 2)  Max eval correction = 0.1

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=10 --showp

In the second command the preprocessor includes the DYN category in the input file, and sets up a molecular statics calculation to run for a maximum of 10 relaxation steps

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 number is a dimensionless step length in Cartesian coordinates (it is scaled by 1/ALAT). For step size in atomic units, multiply by ALAT. This a good choice in the present context, but in other contexts it may not be.

To relax the structure, do

mpirun -n 6 lmf se -vminx=10 --wpos=pos > out.relax

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

At the end of the first self-consistency cycle you should see

 RELAX: (warning) failed to read Hessian; set to unity
 gradzr: begin Broyden xtol=1e-3  dxmx=0.01  isw=211
 brmin: start  xtol=|1e-3|  isw=15  |g|=0.0779
 brmin: max shift =  0.04496 is larger than dxmx.  Scale by  0.2224
 brmin: iter 1  max shift=0.01  |g|=0.07787  max gi=0.04496
 Updated x: 
   0.227000  0.000000  0.000000 -0.113500 -0.196588  0.378667 -0.113500  0.196588
  -0.378667
 RELAX:  completed Broyden iter 1;  max shift=0.01000  |g|=0.0779
        Gradients:
  -0.045  -0.000  -0.000   0.022   0.039   0.000   0.022  -0.039   0.000
      Diagonal inverse Hessian:
   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000

 Updated atom positions:
  Site   Class                      Position(relaxed)
   1      Se         0.22700000(T)    0.00000000(T)    0.00000000(T)
   2      Se        -0.11350000(T)   -0.19658777(T)    0.37866663(T)
   3      Se        -0.11350000(T)    0.19658777(T)   -0.37866663(T)

|g| is a scalar average of the gradients over all the sites; when |g|=0 all the forces vanish.

The x coordinate of the first Se atom changed by 0.01 (this was controlled by STEP=.010, which limits the largest allowed change in any component).

Next follows this table:

 smshft:  add shifted free atom densities
   site                old pos                      new pos              shift
   1:Se       0.21700  0.00000  0.00000    0.22700  0.00000  0.00000   0.010000
   2:Se      -0.10850 -0.18793  0.37867   -0.11350 -0.19659  0.37867   0.010000
   3:Se      -0.10850  0.18793 -0.37867   -0.11350  0.19659 -0.37867   0.010000

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:

 Forces, with eigenvalue correction
  ib              estatic                     eigval                      total
   1    991.32     0.00     0.00   -971.52     0.00     0.00     19.80     0.00     0.00
   2   -495.66  -858.51     0.00    485.76   841.36     0.00     -9.90   -17.14     0.00
   3   -495.66   858.51     0.00    485.76  -841.36     0.00     -9.90    17.14     0.00
 Maximum Harris force = 19.8 mRy/au (site 2)  Max eval correction = 173.6

They are about half as large as the starting force. At the close of the fourth relaxation step the forces read

 Forces, with eigenvalue correction
  ib              estatic                     eigval                      total
   1    767.75     0.00     0.00   -767.07     0.00     0.00      0.69     0.00     0.00
   2   -383.88  -664.89     0.00    383.53   664.30     0.00     -0.34    -0.60     0.00
   3   -383.88   664.89     0.00    383.53  -664.30     0.00     -0.34     0.60     0.00
 Maximum Harris force = 0.688 mRy/au (site 2)  Max eval correction = 13.7

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 about 0.2354. 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:

grep ^c save.se

You should see

c ehf=-14561.1229831 ehk=-14561.122984
c minx=10 ehf=-14561.1229825 ehk=-14561.1229825
c minx=10 ehf=-14561.1305456 ehk=-14561.1305457
c minx=10 ehf=-14561.1316825 ehk=-14561.1316825
c minx=10 ehf=-14561.1317279 ehk=-14561.1317253

The energy drops monotonically, with a large change after the first relaxation step, and almost no change at the last step. This is because the energy stops changing as the forces go to zero.

You can see how the Se-Se bond length and angles change on relaxation:

lmchk se --rpos=pos --angles

You should find that the bond angle drops to 101^o, and the bond length increases a few percent.

Note also the bandgap (grep gap out.relax). In the last iteration it should read

 VBmax = 0.126549  CBmin = 0.191860  gap = 0.065311 Ry = 0.88860 eV

The experimental gap is about 1.8 eV. As is typical, DFT underestimates bandgaps.

3.1 Continuing incomplete relaxations

It may be that the relaxation is incomplete. In such a case, lmf would save the rst file not with the last site positions, but restore the point where the corresponding forces were smallest, before writing rst.se. (If you do only one iteration, this will always be the initial positions, since they are the only positions for which the forces are known!). However the pos file it writes (you must use --wpos=pos to accomplish that) uses positions from the last iteration.

What final positions are saved in the rst file

lmf keeps track of which collection of coordinates were found corresponding to the minimum value of |g|. Unless you specify otherwise, lmf will restore lattice coordinates to that point. (Note that the density was calculated from coordinates of the last cycle, and it will not be consistent with the restored coordinates, unless the minimum-|g| iteration is the final one).

You can also write the site positions --wpos=fnam as noted above. Positions estimated by the last call to the relaxation algorithm are written to file fnam.ext.

You can modify the default using DYN_MSTAT_LASTR. If not specified, this value defaults to −1, which precipates the default action. If you specify DYN_MSTAT_LASTR=0, it retains the coordinates of the final iteration instead. If you specify DYN_MSTAT_LASTR=1, it will save in the rst file the same positions as would be written to fnam.ext if --wpos=fnam is used.

The safest (but often not most efficient) is to use the default (DYN_MSTAT_LASTR=−1). If the relaxation is not troublesome, it is usually better to use DYN_MSTAT_LASTR=0.

Continuing relaxations

If you restart the calculation intending to resume relaxations, note that lmf reads the positions from the rst file by default. It may be that’s what you want to do. But you may wish to use some other positions, e.g. those generated by lmf estimated for next relaxation step. To do so, rather than the more conservative one in the rst file), you can do the following.

One way is to lmf without using the restart file, but with positions dumped by the relaxation algorithm

rm mixm.se
lmf se --rpos=pos --rs=0

Alternatively, run lmf using the rst file. You can to remake the rst with updated positions:

rm mixm.se
lmf se --rpos=pos --rs=11,1,1 

The switches do the following

--rpos=pos tells lmf to read site positions from pos.se. Note that positions may be read as many as three times, in the following order:
- when read from the ctrl or site files
- when read from pos.se, superseding the original positions
- when read from rst.se, superseding the current positions
You want lmf to read the second set, but not the third. This is the purpose of –rs=11,1,1:
- The 10’s digit of the first number (11) causes lmf to adjust the density read from rst.se, by adding a Mattheis construction of the density at the positions it will actually use, and subtract the corresponding construction of the density assuming the positions in rst.se. The 1’s digit tells lmf to read the rst file.
- The second number causes lmf to write the rst file after it makes a band pass.
- The third number tells lmf to use the site positions it already has, and suppresses using those contained in the rst file.

After running this step, lmf should generate a self-consistent density at the positions contained in the pos file, and they are also updated in the rst file. You can start a new relaxation with these positions.

Additional Exercises (optional)

Check out how the relaxed x coordinate depends on the choice of functional. The PBESOL functional (lxcf1=116 lxcf2=133) and the LDA (lxcf=2) both yield very similar values for u.

Table of Contents