Setting up the lmf ctrl file and running calculations
The behaviour of all the Questaal codes is governed by the ctrl file. This file has a common format that is understood by all the programs in the Questaal suite, although there are many keys/options that are specific to different methods. Filenames in Questaal are organised such that the root of the filename labels the kind of file and the extension is a project identifier – this makes finding files from the command line easier.
There are many tokens for the control file and they are arranged into sections: options relating to the Hamiltonian belong in the “HAM” section, structural details in “STRUC”, species data in “SPEC”, etc. The codes provide quick reference to the input tokens via the command line flag “–input”; other command line options are show by invoking the program with “–help.”
lmf --input | less
lmf --help
Minimal ctrl file for GaN
# wurtzite GaN
VERS LM=7 FP=7
HAM
AUTOBAS[MTO=2 LMTO=5]
XCFUN=1 # Ceperley-Alder LDA
GMAX=0.0
BZ
NKABC=1 1 1
ITER
NIT=20
STRUC
NSPEC=2
NBAS=4
ALAT=6.02822874 # expt. lattice constant 3.19A,c/a=1.6266458
PLAT=sqrt(3)/2 -1/2 0
0 1 0
0 0 1.6266458
SPEC
ATOM=Ga Z=31 R=2.0 LMX=2
ATOM=N Z=7 R=1.6 LMX=2
SITE
ATOM=Ga XPOS=2/3 1/3 0
ATOM=Ga XPOS=1/3 2/3 1/2
ATOM=N XPOS=2/3 1/3 0.377
ATOM=N XPOS=1/3 2/3 0.877
Notes:
- white space is unimportant (but subsections must be indented if on a new line)
- tokens are, eg: STRUC_NBAS, HAM_AUTOBAS_LMTO
- options after an ATOM relate to that species until the next ATOM is found
- XPOS for fractions of the lattice vectors, POS for cartesian (multiples of ALAT)
Choices:
- muffin-tin radius R – guessed manually; improve later
- basis parameters – guessed automatically using AUTOBAS
- lmax for the basis – choice
- resolution of the smooth grid, GMAX – still needs to be completed
Checking the structure with lmchk
lmchk is useful for printing geometry, neighbour distances, angles, etc. Also for determining good muffin-tin estimates:
lmchk --getwsr gan | less
recommends for R:
makrm0: initial MT radii from first estat potential maximum
site spec rmt rmt- rmt- rold lock
<spec avg> spec-min
1 1:Ga 2.0153 0.0000 0.0000 2.0000
2 1:Ga 2.0153 0.0000 0.0000 2.0000
3 2:N 1.6679 0.0000 0.0000 1.6000
4 2:N 1.6679 0.0000 0.0000 1.6000
use these to improve the input file
Setting up the basis using AUTOBAS
The basis consists of augmented smoothed-Hankel functions, up to ATOM_LMAX, described by RSMH and EH.
The default mechanism to setup the LMTO basis is based on fitting smoothed-Hankel basis functions to atomic eigenfunctions; these form one (kappa) set, giving RSMH and EH pairs for each upto LMAX. A second set is obtained from the first by adding a constant to EH. This approach is quite reliable and is fully automatic.
Automatic basis construction is activated by nonzero HAM_AUTOBAS_MTO and HAM_AUTOBAS_LMTO. If MTO=1,2 a one,two kappa set is generated. LMTO controls the size of the basis to be constructed (MTO=2 LMTO=5 ~ spd,spd).
The basis fitting is done by the free-atom code lmfa and is written to the file basp0.gan:
lmfa gan | tee lmfa_log
note: the main output for the codes is directed to standard out
Species Ga: Z=31 Qc=28 R=2.015300 Q=0
mesh: rmt=2.015300 rmax=49.443315 a=0.025 nr=389 nr(rmax)=517
Pl= 4.5 4.5 4.5 4.5
Ql= 2.0 1.0 0.0 0.0
iter qint drho vh0 rho0 vsum beta
1 31.000000 5.270E+03 155.0000 0.1542E+03 -62.2954 0.30
51 31.000000 4.475E-05 302.3465 0.3447E+05 -166.2355 0.30
sumev=-1.548744 etot=-3882.268124 eref=0.000000
Free-atom wavefunctions:
valence: eval node at max at c.t.p. rho(r>rmt)
4s -0.67398 0.807 1.882 2.876 0.584682
4p -0.20078 0.869 2.519 5.145 0.830684
4d 0.01244 1.254 33.370 49.443* 0.999995
4f 0.01890 0.000 35.790 49.443* 1.000000
core: ecore node at max at c.t.p. rho(r>rmt)
1s -750.27749 0.000 0.032 0.065 0.000000
2s -92.58142 0.065 0.183 0.285 0.000000
2p -80.67737 0.000 0.146 0.309 0.000000
3s -10.83385 0.267 0.565 0.823 0.000217
3p -7.24866 0.240 0.567 0.979 0.000870
3d -1.42658 0.000 0.531 2.017 0.022565
Optimise free-atom basis for species Ga, rmt=2.0153
l it Rsm Eh stiffR stiffE Eval Exact Pnu Ql
0 24 1.916 -0.563 266.6 191.7 -0.67366 -0.67398 4.78 2.00
1 13 1.859 -0.100 289.5 1726.4 -0.19838 -0.20078 4.63 1.00
eigenvalue sum: exact -1.54874 opt basis -1.54571 error 0.00303
Make LMTO basis parms for species Ga to lmxb=3, rmt=2.0153 vbar=0
l it Rsm Eh Eval Exact Pnu Ql Gmax
0 10 1.344* -0.251 -0.66140 -0.67398 4.78 2.00 5.5
1 11 1.344* -0.100* -0.13808 -0.20078 4.63 1.00 5.9
lmfa writes to basp0.gan in order not to overwrite basp.gan. For running lmf, we need to rename the basp file. lmfa also writes free-atom densities to disk, which are overlapped to initialise the self-consistency cycle.
cat basp0.gan
mv basp0.gan basp.gan
For lmax=2 and LMTO=5, we have a (spd,spd) basis for both Ga and N.
FREEAT: estimate HAM_GMAX from RSMH: GMAX=8.8
Once the basis has been constructed, it is possible to make an estimate of the GMAX necessary to describe the basis. This is printed at the end of the lmfa output. We need to edit the ctrl file with the recommended value (8.8) and to specify the k-mesh (say 10^3).
Running lmf
mpirun -np 24 lmf gan | tee lmf_log
This should converge in 9 iterations in about ~3s (parallelism is over k-points). Convergence progress can be monitored using the save file:
h ehf=-7983.012975 ehk=-7982.2720881
i ehf=-7982.9833865 ehk=-7982.2611507
i ehf=-7982.7731994 ehk=-7982.6602763
i ehf=-7982.7336772 ehk=-7982.7335436
i ehf=-7982.733418 ehk=-7982.7333624
i ehf=-7982.7333979 ehk=-7982.7333642
i ehf=-7982.7333797 ehk=-7982.7333467
i ehf=-7982.7333582 ehk=-7982.7333237
c ehf=-7982.7333627 ehk=-7982.7333263
^h – density started using overlapped free atoms ^i – iterating ^c – converged
ehf – Harris-Foulkes total energy (depends on input density and not the output density) ehk – Hohenburg-Kohn total energy (Ryd)
and the gap is 2.04087 eV. (grep gap lmf_log).
main output files:
- atm.gan (from lmfa) – atomic charges and potentials (ASCII)
- basp.gan (from lmfa) – basis definition (like ctrl)
- log.gan – terse output information (ASCII)
- mixm.gan – scf mixing history (binary)
- moms.gan – weights of overlap matrix (binary)
- rst.gan – restart file (binary)
- save.gan – scf progress (ASCII)
- wkp.gan – k-point weights (binary)
Adding local orbitals
HAM
AUTOBAS[MTO=2 LMTO=5 LOC=1]
Rerunning lmfa picks up the shallow Ga 3d:
Search for possible local orbitals
LO l=2: qistl=0.0226 exceeds QLOC, autogen PZ=13.91024 use PZ=13.91024
LO l=2: e=-1.43 exceeds ELOC, autogen PZ=13.91024 use PZ=13.91024
basp0.gan reads:
BASIS:
Ga RSMH= 1.344 1.344 1.344 EH= -0.1 -0.1 -0.1 RSMH2= 1.344 1.344 1.344 EH2= -0.9 -0.9 -0.9 PZ= 0 0 13.9102
N RSMH= 1.002 0.909 1.112 EH= -0.191 -0.1 -0.1 RSMH2= 1.002 0.909 1.112 EH2= -0.991 -0.9 -0.9
copy basp0 to basp and run lmfa again (change in core charge number); delete mixm, rst.gan.
grep ^c save.gan
c ehf=-7982.7333627 ehk=-7982.7333263
c ehf=-7982.7921134 ehk=-7982.7920645
the gap becomes 1.97349 eV.
Big basis – lmax
increase lmax to 3 to setup an (spdf,spd) basis set re-run lmfa mv basp0.gan basp.gan
BASIS:
Ga RSMH= 1.344 1.344 1.344 1.344 EH= -0.1 -0.1 -0.1 -0.1 RSMH2= 1.344 1.344 1.344 EH2= -0.9 -0.9 -0.9 PZ= 0 0 13.9102
N RSMH= 1.002 0.909 1.112 1.112 EH= -0.191 -0.1 -0.1 -0.1 RSMH2= 1.002 0.909 1.112 EH2= -0.991 -0.9 -0.9
lmf now gives:
grep ^c save.gan
c ehf=-7982.7333627 ehk=-7982.7333263
c ehf=-7982.7921134 ehk=-7982.7920645
c ehf=-7982.8033203 ehk=-7982.8032586
gap 1.93483 eV.
Bigger basis – combined LMTO/APW
additional ctrl flags:
HAM
PWMODE=11 PWEMIN=0.0 PWEMAX=2.0
plane-wave cut-offs in Ryd.
PWEMAX=2 Ryd
suham: q-dependent PW basis with Emin = 0 < E < 2.
Est. min,max PW dimension = 11,23. Use npwpad = 3 => ndham = 136
c ehf=-7982.8080754 ehk=-7982.8080078
gap 1.93285 eV.
PWEMAX=4 Ryd
suham: q-dependent PW basis with Emin = 0 < E < 4.
Est. min,max PW dimension = 30,51. Use npwpad = 4 => ndham = 165
c ehf=-7982.8093846 ehk=-7982.8093305
gap 1.93235 eV.
Notes:
- usually PWEMAX=1 or 2 Ryd is effective for open structures
- over-completeness becomes a problem (see OVEPS)
- for large APW component, reduce R, remove MTOs.
Ctrl preprocessor directives
The Questaal ctrl is interpreted by a preprocessor which allows variable substitution, looping, arithmetic, flow control, … Preprocessor directives are preceded by the “%” symbol, which introduces constants and variables, and braces “{}” which are then expanded.
First use of this is to shadow ctrl options and make them accessible from the command line, eg, for checking the k-convergence:
% const nk=10
BZ
NKABC={nk}
The variable nk can now be changed as an argument to lmf, eg:
lmf -vnk=12 gan
NKABC accepts either one number or three numbers to describe an M.P. grid – if one number, then all directions have the same number of k-divisions. If a negative number is given, this is interpreted as the desired number in the total BZ, and a set of divisions proportional to the length of the inverse cell vectors will be generated that is closest to this target. In this case, this reduces the number of divisions along the (longer) c-axis.
for i in -200 -400 -800 -2000 -4000 -8000
do
mpirun -np 24 lmf -vnk=$i gan |grep BZMESH
done
gives:
i nk=-200 ehf=-7982.8032021 ehk=-7982.8031549
c nk=-200 ehf=-7982.8032019 ehk=-7982.803155
i nk=-400 ehf=-7982.8032963 ehk=-7982.8032494
c nk=-400 ehf=-7982.8032962 ehk=-7982.8032495
i nk=-800 ehf=-7982.8033115 ehk=-7982.8032652
c nk=-800 ehf=-7982.8033099 ehk=-7982.8032638
i nk=-2e3 ehf=-7982.8033133 ehk=-7982.8032674
c nk=-2e3 ehf=-7982.8033133 ehk=-7982.8032676
i nk=-4e3 ehf=-7982.8033133 ehk=-7982.8032678
c nk=-4e3 ehf=-7982.8033131 ehk=-7982.8032678
i nk=-8e3 ehf=-7982.803313 ehk=-7982.8032679
c nk=-8e3 ehf=-7982.8033129 ehk=-7982.803268
Obviously metals require more attention; by default Questaal uses metallic occupations and tetrahedron integration: only in special circumstances would you need to change these. To display the as-preprocessed ctrl file, use: lmf –showp gan .
Automatic setup with blm
sample cif file:
cp /home/vol05/tmp15/gan_icsd_34476.cif .
cif2cell icsd_34476.cif
cif2cell icsd_34476.cif > out
cif2init out
cif2cell is a (very useful!) external tool (cif2cell at sourceforge); cif2init converts the output to an “init” file, which contains basic structural information than can be read by blm. blm generates R and sets various variables depending on the specific kind of calculation.
mv init init.gan
blm gan
mv actrl.gan ctrl.gan
(actrl is to avoid potentially clobbering ctrl)
vi ctrl.gan
cat site.gan
lmfa --usebasp gan
cat basp.gan
mpirun -np 24 lmf -vgmax=8.8 -vnkabc=10 gan |tee log_lmf
the result is different because the default basis is (spd,sp) and because XCFUN=2 (von Barth-Hedin) has been chosen.
mpirun -np 24 lmf -vgmax=8.8 -vnkabc=10 -vlxcf=1 gan |tee log_lmf2
Notes:
- check the preprocessor names/values generated by blm
- -v is used to change the values of preprocessor variables
- blm has lots of options – tune the setup at the point of blm invocation
- blm –help
- –express=0 gives more “old fashioned” input files