First steps in lmf
Table of Contents
- Links to docs
- Preparation for the exercises
- Example 1: SrTiO3
- Example 2: BCC Fe
- Example 3: ErAs (ferromagnetic)
- Exercise
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.
Links to docs
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 tmpXX@ui2.scarf.rl.ac.uk -i ~/.ssh/id_ed25519_questaal_school
… and immediately request an interactive slot:
get-session
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 .
ls
Example 1: SrTiO3
- interpret cif file using cif2cell
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
Try:
blm --help
and:
blm --simple init.srtio3 |tee lblm
The FP-LMTO basis requires:
- 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.
- maximum l for LMTO basis (determines size of the basis): lmx
- the maximum l used for expanding tails of MTO from other sites: lmxa
- 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
then
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
setup a path in the BZ
lmchk --syml ctrl.srtio3examine the path
less syml.srtio3get eigenvalues (specifying the filename to be syml.srtio3):
mpirun -n 12 lmf --band~fn=syml ctrl.srtio3produce 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.plbndscopy back to your machine to view
pwd scp tmpXX@ui2.scarf.rl.ac.uk:[directory shown]/fplot.ps . okular fplot.ps
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)
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 ham gmax= 8.7 autobas[pnu=1 loc=1 mto=4] nspin= 2 so= 0 xcfun="xc_gga_x_pbe xc_gga_c_pbe" iter nit= 10 mix= b2,b=0.8,k=7 conv= 1e-5 convc= 3e-5 % const nk=16 bz nkabc={nk} struc 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 spec 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
Gives:
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:
n(E) resolved by l, save weights generated by lmf
mpirun -n 12 lmf -vnk=40 --pdos~mode=1 --quit=rho ctrl.fesum 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
plot
echo .6 20 -6 6 | pldos -fplot -ef=0 -esclxy=13.6057 -lst="1" -lst2 dos.fe fplot -f plot.dos mv fplot.ps fe-dos-s.ps echo .6 20 -6 6 | pldos -fplot -ef=0 -esclxy=13.6057 -lst="3" -lst2 dos.fe fplot -f plot.dos mv fplot.ps fe-dos-p.ps echo 6 20 -6 6 | pldos -fplot -ef=0 -esclxy=13.6057 -lst="5" -lst2 dos.fe fplot -f plot.dos mv fplot.ps fe-dos-d.ps
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.
Setup
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 ham gmax= 9.4 autobas[pnu=1 loc=1 mto=4 pfloat=1,0] nspin= 2 so= 0 xcfun="lm_lda_bh" iter nit= 60 mix= b2,b=0.1,k=7 conv= 1e-5 convc= 3e-5 bz nkabc= 8 struc nbas=2 nspec=2 alat= 10.83191015 plat= 0.5 0.5 0 0.5 0 0.5 0 0.5 0.5 spec 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.0000000starting 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
mv fplot.ps eras-dos-lsda.ps
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
Notes:
- 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.
spec 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
mv fplot.ps eras-dos-ldau.ps
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:
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.erasyields:
# 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 Ghighlight the f states in red:
lmf --band~scol@z=68@l=3~fn=syml ctrl.erasplot 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 mv fplot.ps eras-bnd-ldau-up.ps 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 mv fplot.ps eras-bnd-ldau-dn.ps
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.
so=1- spin are coupled now, but we can plot projections on up, down
position of f-states in Hamiltonian:
lmf --pr50 --quit=ham ctrl.eras
gives:
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
mv fplot.ps eras-bnd-ldau-so.ps
Exercise
- Calculate GaAs and plot the bands.
- Add empty spheres and vary lmx and lmxa; what difference do these make to the energies, the bands?
- What happens when gmax is increased?
- Should Ga 3p be treated as semi-core or core? Arsenic?
- How do the bands change when spin-orbit is included?