Questaal Home

First steps in lmf

Table of Contents

This tutorial describes the use of the full-potential code lmf, which is an accurate density functional implementation and the underlying band code providing single-particle states used in the QSGW and +DMFT methods.

The full-potential single-particle/DFT code in Questaal is called lmf:

  • full-potential linear muffin-tin-orbital method (FP-LMTO)
  • all-electron augmented method
  • basis of atom centred spherical waves (smoothed Hankel functions)
  • precise yet fast

Three cases are worked through, the cubic insulator SrTiO3, BCC iron and the f-electron compound ErAs in the rocksalt structure.

  1. command line flags
  2. plbnds utility
  3. symmetry line notes

Preparation for the exercises

Please connect to the departmental cluster using your SSH “identity” (private SSH key part), replacing tmpXX with the username assigned to you (see slack announcements):

    ssh -i ~/.ssh/id_ed25519_questaal_school

… and immediately request an interactive slot:


This sets up an interactive session for you with 12 cores (one Intel Xeon socket).

Copy the .cif files to where you are working:

    cp /work3/training/questaal/data/*.cif .

Example 1: SrTiO3

  1. interpret cif file using cif2cell
  2. reformat output to “init” file format.

     cif2cell SrTiO3.cif > cif2cell.txt
     cif2init cif2cell.txt init.srtio3
     cat init.srtio3
  • Questaal filenames usually type.extn, extn is a project name

Setup a control file with blm

blm is a utility to setup ctrl files for the various Questaal codes

blm reads the basic init file and produces main control files “ctrl” suitable for different purposes


    blm --help


    blm --simple init.srtio3 |tee lblm

The FP-LMTO basis requires:

  1. radii of the non-overlapping atomic spheres. In LMTO these should generally be as large as possible because the augmentation (within the sphere) better describes the wavefunction than the MTO “tails” do in the interstitial.
  2. maximum l for LMTO basis (determines size of the basis): lmx
  3. the maximum l used for expanding tails of MTO from other sites: lmxa
  4. the plane-wave cutoff for representing quantities on the smooth mesh: gmax

These, along with the crystal structure are setup by blm with good defaults: the muffin-tin construction is obtained by finding the maximum in V(r) resulting from overlapped atomic charges. The basis lmx and lmxa are read from internal tables; gmax is determined by the requirement of representing the basis to a sufficient tolerance.

The output of blm is the file actrl.ext. This should be copied to ctrl.srtio3:

    less actrl.srtio3
    mv actrl.srtio3 ctrl.srtio3

For a simple calculation, it doesn’t need much.

To ensure that the system is as you expect, it is sensible to check:

    lmchk ctrl.srtio3 

Tokens are arranged in groups (indented for each). To see all the tokens understood by lmf with a brief description, try:

    lmf --input ctrl.srtio3 | less

Only very few of these are commonly needed.

Running lmf

blm can’t tell us what k-mesh we need. Edit the ctrl file to provide a mesh. We need to change “null” to a suitable value:

  • 1 number (which is used in all 3 directions in reciprocal space)
  • 3 numbers (one for each direction)
  • 1 negative number (an mesh with approximately so many points is formed such that the mesh density in each direction is roughly constant)

Either use your editor of choice, or execute the sed command:

    sed -i 's/null/6 6 6/' ctrl.srtio3

lmf uses overlapped atomic charges (the Mattheis construction) as the initial charge density. We provide this by calling the atom solver lmfa:

    lmfa ctrl.srtio3 | tee llmfa


    mpirun -n 12 lmf ctrl.srtio3 | tee llmf

it doesn’t converge in 10 iterations, continue:

    mpirun -n 12 lmf ctrl.srtio3 | tee -a llmf
    grep gap llmf
    less llmf

The output of lmf is described in detail here.

Plotting energy bands

  1. setup a path in the BZ

     lmchk --syml ctrl.srtio3
  2. examine the path

     less syml.srtio3
  3. get eigenvalues (specifying the filename to be syml.srtio3):

     mpirun -n 12 lmf --band~fn=syml ctrl.srtio3
  4. produce a plot (-8 to 8 eV, set E_Fermi at zero, figure size 12cm x 8cm)

     echo -8 8 12 8 | plbnds -fplot -ef=0 -scl=13.6057 -lbl=G,M,X,G,R,XR,M bnds.srtio3
     fplot -f plot.plbnds
  5. copy back to your machine to view

     scp[directory shown]/ .

Plotting weighted energy bands

    mpirun -n 12 lmf --band~fn=syml~scol@z=22@l=2~scol2@z=8@l=1 ctrl.srtio3
    echo -8 8 12 8 | plbnds -fplot -ef=0 -scl=13.6057 -lbl=G,M,X,G,R,XR,M -lt=1,col=0,0,0,colw=1,0,0,colw2=0,0,1 bnds.srtio3
    fplot -f plot.plbnds

Three curves: “col” with colour 0,0,0 (black) (everything). “col1,” red for Ti-d and “col2,” blue for O-p. The Ti states, which are unoccupied here, form a set of localised states nicely separate from other bands.

(See here for help using plbnds)

SrTiO3 bands

Example 2: BCC Fe

Prepare an init file for Fe:

    cif2cell Fe.cif > cif2cell.txt
    cif2init cif2cell.txt init.fe

Cause blm to setup a magnetic problem (nspin=2) using Libxc’s PBE-GGA:

    blm --simple --mag --pbe init.fe
    mv actrl.fe ctrl.fe

Edit the ctrl file:

  • give Fe a starting moment of 2 mu_B in l=2
  • increase the mixing parameter b to 0.8
  • include a preprocessor directive for bz_nkabc

      # Autogenerated from init.fe using:
      # blm --simple --mag --pbe init.fe
      vers  lm:7 fp:7
      #symgrp i*r3d r2(1,1,0)
      #symgrp spglib
            gmax=   8.7
            autobas[pnu=1 loc=1 mto=4]
            nspin=  2
            so=     0
            xcfun="xc_gga_x_pbe xc_gga_c_pbe"
            nit=    10
            mix= b2,b=0.8,k=7
            conv=   1e-5
            convc=  3e-5
      % const nk=16
            nbas=1  nspec=1
            alat=   5.41689994
            plat= -0.5  0.5  0.5  0.5  -0.5  0.5  0.5  0.5  -0.5
        atom=Fe       z= 26  r= 2.345586  lmx=2 lmxa=4 a=0.015 mmom=0 0 2
      site                                  # Positions in Cartesian coordinates
        atom=Fe         pos=  0.0000000   0.0000000   0.0000000

To see the ctrl file after preprocessing (i.e. what is used), try:

    lmf --showp ctrl.fe
    lmf --showp -vnk=20 ctrl.fe

Setup starting (atomic) charges:

    lmfa ctrl.fe |tee llmfa

Investigate the convergence of moment and total energy with k-mesh using tetrahedron integration:

    for i in 16 20 24 28 32; do  mpirun -n 12 lmf -vnk=$i ctrl.fe |tee -a llmf; done

    grep ^c save.fe


    c nk=16 mmom=2.2329777 ehf=-2545.5633846 ehk=-2545.5633896
    c nk=20 mmom=2.2275563 ehf=-2545.5633772 ehk=-2545.5633816
    c nk=24 mmom=2.2242604 ehf=-2545.5633796 ehk=-2545.5633836
    c nk=28 mmom=2.2221787 ehf=-2545.5633752 ehk=-2545.5633778
    c nk=32 mmom=2.2199566 ehf=-2545.5633789 ehk=-2545.5633806

The moment is printed at each iteration. It is almost entirely contained in the muffin-tin.

     BZWTS : --- Tetrahedron Integration ---
             Est E_f           Window        Tolerance  DOS(E_f)
             0.023774   0.022206   0.024522   0.002316  14.153716
             0.023774   0.023757   0.023781   0.000023  14.154847
             0.023774   0.023773   0.023774   0.000000  14.154872
     BZINTS: Fermi energy:      0.023774;  14.000000 electrons;  D(Ef):   14.155
             Sum occ. bands:  -24.1801023  incl. Bloechl correction:  -0.0003573
             Mag. moment:       2.227841

     Saved qp weights ...

     mkrout:  Qtrue      sm,loc       local        true mm   smooth mm    local mm
       1   13.023780    3.755799    9.267981      2.286754    0.429379    1.857375 

More usefully, we can invoke lmf with different grids, eg when plotting densities of states:

  1. n(E) resolved by l, save weights generated by lmf

     mpirun -n 12 lmf -vnk=40 --pdos~mode=1 --quit=rho ctrl.fe
  2. sum weights into different spins/channels

     mpirun -n 12 lmdos -vnk=40 --pdos~mode=1 --quit=rho --dos:npts=1001:window=-0.5,0.5 ctrl.fe

     Channels in dos file generated by LMDOS:
     site class label   spin-1                       spin-2
        1    1   Fe     1:9:2                        2:10:2
  1. plot

     echo .6 20 -6 6 | pldos -fplot -ef=0 -esclxy=13.6057 -lst="1" -lst2 dos.fe
     fplot -f plot.dos
     echo .6 20 -6 6 | pldos -fplot -ef=0 -esclxy=13.6057 -lst="3" -lst2 dos.fe
     fplot -f plot.dos
     echo 6 20 -6 6 | pldos -fplot -ef=0 -esclxy=13.6057 -lst="5" -lst2 dos.fe
     fplot -f plot.dos

Fe d dos

What does the density of states look like for non-magnetic Fe? Paramagnetic Fe? How do we know that the ground state is ferromagnetic?

Example 3: ErAs (ferromagnetic)

Erbium arsenide has the rock salt structure; the 11 Er electrons take up localised, atomic-like states and form a moment of (7 “up”, 4 “down”). The LDA description, failing to distinguish between the occupied and unoccupied minority f-states, places these minority states at , giving a metallic description. LDA+U correctly splits the occupied and unoccupied states, moving the f-states away from the Fermi-level and restoring the observed semi-metallic behaviour. The +U description is a reasonable starting point for QSGW.


    cif2cell ErAs.cif > cif2cell.txt
    cif2init cif2cell.txt init.eras

    blm --simple --mag --nk=8 init.eras | tee lblm
    cp actrl.eras ctrl.eras
  • add starting moment for Er f-electrons
  • decrease mixing from 0.3 to 0.1, increase number of iterations

      # Autogenerated from init.eras using:
      # blm --simple --mag --nk=8 init.eras
      vers  lm:7 fp:7
      #symgrp i*r3(-1,1,1) r4z
      #symgrp spglib
            gmax=   9.4
            autobas[pnu=1 loc=1 mto=4 pfloat=1,0]
            nspin=  2
            so=     0
            nit=    60
            mix= b2,b=0.1,k=7
            conv=   1e-5
            convc=  3e-5
            nkabc= 8
            nbas=2  nspec=2
            alat=   10.83191015
            plat= 0.5  0.5  0  0.5  0  0.5  0  0.5  0.5
        atom=Er       z= 68  r= 2.846544  lmx=3 lmxa=6 mmom=0 0 0 3
        atom=As       z= 33  r= 2.569411  lmx=2 lmxa=3
      site                                  # Positions in Cartesian coordinates
        atom=Er         pos=  0.0000000   0.0000000   0.0000000
        atom=As         pos=  0.0000000   0.5000000   0.0000000
  • starting density:

      lmfa ctrl.eras |tee llmfa

LDA description

    mpirun -n 12 lmf ctrl.eras > llmf
    mpirun -n 12 lmf --quit=dos --dos~npts=2001~window=-1,1 ctrl.eras
    echo 25 20 -12 12 | pldos -esclxy=13.6057 -ef=0 -fplot -lst=1 -lst2 dos.eras
    fplot -f plot.dos

ErAs LDA dos

     BZINTS: Fermi energy:     -0.047790;  25.000000 electrons;  D(Ef):  888.689
             Sum occ. bands:  -14.2924244  incl. Bloechl correction:  -0.0004008
             Mag. moment:       2.725389

LDA+U description


  • IDU=0,1,2 for no +U, AMF or FLL double counting corrections
  • IDU+=10 flags the inclusion of as well as
  • pick a Hubbard U, J.

        atom=Er       z= 68  r= 2.846544  lmx=3 lmxa=6 mmom=0 0 0 3
        idu= 0 0 0 12
        uh=  0 0 0 {8.0/13.6057}
        jh=  0 0 0 {0.5/13.6057}
        atom=As       z= 33  r= 2.569411  lmx=2 lmxa=3
      mpirun -n 12 lmf ctrl.eras > llmf-plusu

(did it converge yet?)

    mpirun -n 12 lmf ctrl.eras >> llmf-plusu
    mpirun -n 12 lmf --quit=dos --dos~npts=2001~window=-1,1 ctrl.eras
    echo 25 20 -12 12 | pldos -esclxy=13.6057 -ef=0 -fplot -lst=1 -lst2 dos.eras
    fplot -f plot.dos

The situation at the Fermi level:

     BZINTS: Fermi energy:     -0.016411;  25.000000 electrons;  D(Ef):    3.677
             Sum occ. bands:  -20.1499962  incl. Bloechl correction:  -0.0006934
             Mag. moment:       3.005527

Show bands:

  1. setup our own path, we prefer L-G-X-W-G, with 50 points max per segment:

     lmchk --syml~n=50~lblq:G=0,0,0,L=0.5,0.5,0.5,X=0,1,0,W=0.5,1,0~lbl=LGXWG ctrl.eras


     # generated from text ~n=50~lblq:G=0,0,0,L=0.5,0.5,0.5,X=0,1,0,W=0.5,1,0~lbl=LGXWG
       43    0.5000000   0.5000000   0.5000000     0.0000000   0.0000000   0.0000000  L  to  G
       50    0.0000000   0.0000000   0.0000000     0.0000000   1.0000000   0.0000000  G  to  X
       25    0.0000000   1.0000000   0.0000000     0.5000000   1.0000000   0.0000000  X  to  W
       56    0.5000000   1.0000000   0.0000000     0.0000000   0.0000000   0.0000000  W  to  G
  2. highlight the f states in red:

     lmf --band~scol@z=68@l=3~fn=syml ctrl.eras
  3. plot spin-up and spin-down bands

     echo -10 10 10 10 | plbnds -fplot -spin1 -ef=0 -scl=13.6057 -lbl=L,G,X,W,G -lt=1,col=0,0,0,colw=1,0,0 bnds.eras
     fplot -f plot.plbnds
     echo -10 10 10 10 | plbnds -fplot -spin2 -ef=0 -scl=13.6057 -lbl=L,G,X,W,G -lt=1,col=0,0,0,colw=1,0,0 bnds.eras
     fplot -f plot.plbnds

Corresponding to the occnum.eras, we arrived at the configuration with maximum orbital moment. Other configurations, with different energies, can be reached from different starting configurations (try!).

Both the total and Er local moment are close to 3, and is sensible.

Including spin-orbit

  • edit ctrl.eras to set so=1 in the “ham” category, then converge again.

  • spin are coupled now, but we can plot projections on up, down
  • position of f-states in Hamiltonian:

      lmf --pr50 --quit=ham ctrl.eras


     Orbital positions in hamiltonian, resolved by l:
     Site  Spec  Total    By l ...
       1   Er    1:35   1:1(s)   2:4(p)   5:9(d)   10:16(f) 17:17(s) 18:20(p) 21:25(d) 26:32(f) 33:35(p)
       2   As   36:48   36:36(s) 37:39(p) 40:44(d) 45:45(s) 46:48(p)

     Orbital positions in hamiltonian, resolved by kappa =l and kappa=-l-1:
     Site  Spec  Total    By kappa ...
       1   Er    1:70   1:2(S)         3:4,5:8(P)     9:12,13:18(D)  19:24,25:32(F) 
                        33:34(S)       35:36,37:40(P) 41:44,45:50(D) 51:56,57:64(F) 65:66,67:70(P)
       2   As   71:96   71:72(S)       73:74,75:78(P) 79:82,83:88(D) 89:90(S)       91:92,93:96(P)

picking out the (f):

    lmf --band~col=10:16,26:32~col2=48+10:48+16,48+26:48+32~fn=syml eras

    echo -10 10 10 10 | plbnds -fplot -ef0 -scl=13.6057 -lbl=L,G,X,W,G -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0 bnds.eras
    fplot -f plot.plbnds

ErAs LDA+U bands


  1. Calculate GaAs and plot the bands.
  2. Add empty spheres and vary lmx and lmxa; what difference do these make to the energies, the bands?
  3. What happens when gmax is increased?
  4. Should Ga 3p be treated as semi-core or core? Arsenic?
  5. How do the bands change when spin-orbit is included?