# Molecular Statics in Se

This tutorial demonstrates how to carry out lattice relaxations of basis positions in **lmf**, by finding the equilibrium positions (molecular statics). It uses elemental Se as a representative material.

### Table of Contents

### Command summary

LDA self-consistency (starting from *init.fe*)

```
nano init.se
blm --ctrl=ctrl --pbe --gmax=5.4 --nk=6 --wsitex=site --molstat init.se
lmchk se --shell --angles
lmfa se --basfile=basp
lmf se
lmf se -vminx=t --wpos=pos
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 ^{n}* 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
```

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

```
$ blm --ctrl=ctrl --pbe --wsitex=site --molstat init.se
```

Switches have the following meaning:

**--ctrl=ctrl**Tells**blm**to write directly to the actual input file. By default it writes to*ctrl.se*so as not to overwrite the input file.*actrl.ext***--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.

```
$ blm --ctrl=ctrl --pbe --gmax=5.4 --nk=6 --wsitex=site --molstat init.se
$ lmfa se --basfile=basp
```

**--molstat** has the effect of adding the following category to the input file:

```
% ifdef minx
DYN MSTAT[MODE=6 HESS=t XTOL=.001 GTOL=0 STEP=.010 NKILL=0] NIT=10
% endif
```

This line will be parsed only if variable **minx** is defined with a nonzero value.

##### 1.1 bond lengths and angles

To confirm the that bond angles are approximately 90^{o} angles, 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 105^{o}, somewhat larger than the simple nearest-neighbor tight-binding model would suggest.

#### 2. Self-consistency

Make the system self-consistent with:

```
$ lmf se
```

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 9 iterations.

You should see this table at the end:

```
Forces, with eigenvalue correction
ib estatic eigval total
1 970.90 0.00 0.00 -917.24 0.00 0.00 53.67 0.00 0.00
2 -485.45 -840.83 0.00 458.62 794.35 0.00 -26.83 -46.47 0.00
3 -485.45 840.83 0.00 458.62 -794.35 0.00 -26.83 46.47 0.00
Maximum Harris force = 53.7 mRy/au (site 2) Max eval correction = 5.8
```

The forces are all equal in magnitude, but rotated by 120^{o} angles.

*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

```
Harris correction to forces: shift in free-atom density
ib delta-n dVes delta-n dVxc total
1 -846.44 0.00 0.00 -30.21 0.00 0.00 876.65 0.00 0.00
2 423.22 733.04 0.00 15.11 26.16 0.00 -438.32 -759.20 0.00
3 423.22 -733.04 0.00 15.11 -26.16 0.00 -438.32 759.20 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 152.09 0.00 0.00 -118.67 0.00 0.00 33.42 0.00 0.00
2 -76.05 -131.72 0.00 59.34 102.77 0.00 -16.71 -28.94 0.00
3 -76.05 131.72 0.00 59.34 -102.78 0.00 -16.71 28.94 0.00
Maximum Harris force = 33.4 mRy/au (site 1) Max eval correction = 876.6
```

*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

```
Harris correction to forces: shift in free-atom density
ib delta-n dVes delta-n dVxc total
1 -5.64 0.00 0.00 -0.14 0.00 0.00 5.78 0.00 0.00
2 2.82 4.88 0.00 0.07 0.12 0.00 -2.89 -5.00 0.00
3 2.82 -4.88 0.00 0.07 -0.12 0.00 -2.89 5.00 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 970.90 0.00 0.00 -917.24 0.00 0.00 53.67 0.00 0.00
2 -485.45 -840.83 0.00 458.62 794.35 0.00 -26.83 -46.47 0.00
3 -485.45 840.83 0.00 458.62 -794.35 0.00 -26.83 46.47 0.00
Maximum Harris force = 53.7 mRy/au (site 2) Max eval correction = 5.8
```

#### 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.

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.0931
brmin: max shift = 0.05376 is larger than dxmx. Scale by 0.1860
brmin: iter 1 max shift=0.01 |g|=0.09312 max g=0.05376
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.0931
Gradients:
-0.054 -0.000 -0.000 0.027 0.047 0.000 0.027 -0.047 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)
```

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 899.40 0.00 0.00 -879.01 0.00 0.00 20.39 0.00 0.00
2 -449.70 -778.90 0.00 439.50 761.24 0.00 -10.20 -17.66 0.00
3 -449.70 778.90 0.00 439.50 -761.24 0.00 -10.20 17.66 0.00
Maximum Harris force = 20.4 mRy/au (site 2) Max eval correction = 1.4
```

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 835.61 0.00 0.00 -835.12 0.00 0.00 0.48 0.00 0.00
2 -417.80 -723.65 0.00 417.56 723.24 0.00 -0.24 -0.42 0.00
3 -417.80 723.65 0.00 417.56 -723.24 0.00 -0.24 0.42 0.00
Maximum Harris force = 0.483 mRy/au (site 2) Max eval correction = 1.1
```

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:

```
$ grep ^c save.se
```

You should see

```
c minx=1 ehf=-14580.5132479 ehk=-14580.5132302
c minx=1 ehf=-14580.5221027 ehk=-14580.5220784
c minx=1 ehf=-14580.523951 ehk=-14580.5239234
c minx=1 ehf=-14580.524057 ehk=-14580.5240287
```

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

```
VBmax = 0.136522 CBmin = 0.201860 gap = 0.065338 Ry = 0.88860 eV
```

The experimental gap is about 1.8 eV.

### 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*.