• modified-Becke-Johnson-Potential, vor allem für bessere Bandlücken nützlich:
• normale PBE-Rechnung konvergieren lassen
• dann:
  • init_mbj_lapw
  • run_lapw -i 1 -NI -p &
  • save_lapw pbe
  • init_mbj_lapw
• danach neue SCF-Läufe
• Konvergenzprobleme, eigene Erfahrungen dazu (nur wenige bisher)
  • PRATT-Mixing nehmen, Parameter 0.2 oder 0.1, später auch 0.4
  • nach 5-10 Zyklen ev. auch wieder auf MSR1 gehen

• !! automatisiert sind nur Achs- und Volumenvariationen möglich !!
• x optimize erzeugt eine Serie von struct-Dateien (Variationen im w2web auswähbar und im Namen kodiert)
  • Minibug: optimize hex a+c (5) nimmt die Paramter von (6)
• edit optimize.job : hier Konvergenzkriterien anpassen und Parallelisierung (-p) ergänzen
• run optimize startet dann eine Serie von Rechnungen
• eplot erlaubt Plots von E vs. vol etc. (Dateien: *vol*.scf oder *a+c*.scf etc.)
• Tests und Eingabedateien: s. ./Wien2k/AM2_Geometrie

• monokline Strukturen
  • angeblich nur bei monoC nach monoB transformieren (Bilbao!), ABER
  • bei Wien2k_2.23 platzt TETRA, wenn beta nicht 90 ist, also auch monoP-Strukturen umindizieren
  • Warnungen über falsche Winkel bleiben, können aber wohl ignoriert werden.
• rhomboedrische Strukturen
  • Eingabe der EZ in hexgonaler Aufstellung
  • Eingabe der Atomlagen in Rhomboeder-Aufstellung
  • hex2rhomb hilft

Auswertung der Ergebnisse

• DOS-Plots
  • zuerst wird lapw2 mit der option -qtl gestartet
  • im File NAME.int eintragen, welche DOS geplottet werden soll. Wer wer ist, erkennt man im Kopf der Datei NAME.qtl, der auch oben in der Maske angezeigt wird. Hier sind maximal sieben Zeilen zulässig. Der Rest sollte ausreichend kommentiert sein.
    1. -0.5 0.002 1.5 0.003 => EMIN (in Ryd.), DE, EMAX
    2. N => Zahl der DOS-Plots = Zeilen die folgen
    3. 0 1 => totale Zustandsdichte = 1. Zeile, voreingestellt
    4. X Y Kommentar => Atomnummer, Spalte in der NAME.qtl, beliebiges Label
       • Y=0 alle s-Zustände
       • Y=1 alle p-Zustände
       • Y=2 alle d-Zustände
       • der Rest in NAME.qtl erklärt
    5. X=X_{max. Atomzahl}+1 und Y=1 plottet die interstitielle Elektronendichte.
    6. usw. usw. insgesamt N-Zeilen
  • tetra sammelt die Dichten über alle k-Punkte zusammen, NAME.outputt enthält die Ausgabe!
  • dosplot (Achtung: hierzu muss das Verzeichnis /home/krist/WIEN2k/SRC_w2web/htdocs/tmp für den normalen User auch schreibbar sein, weil w2web dort seine Bilder vor der Darstellung im Browser temporär ablegt!)
  • TIPP 1: am besten hier nach dem Betrachten gleich die Files NAME.dos1ev, die spaltenweise E und die N DOS enthalten, wegsichern. Es empfiehlt sich, im Header vernünftig dokumentieren, welche Spalte genau was ist (Wien schreibt da praktisch nichts rein!). Diese Dateien sind mit xmgrace weiterverarbeitbar (DATA -> Import-ASCII -> load as NXY).
  • TIPP 2: Die PDOS beziehen sich auf die Ladungen innerhalb der MT-Radien. D.h. viele der eigentlich interessierenden Elektronen wandern in den interstitiellen Bereich zwischen die Kugeln. Diesen kann man ebenfalls plotten, wenn man die Atom-Nr. um 1 grösser macht als die maximale Atomzahl. Für die Beurteilung der Ladungsdichten kann man z.B. eine Bader-Analyse (AIM, s.u.) durchführen.
  • TIPP 3: Wenn vor der DOS Bandstrukturen gerechnet wurden, muss der LAPW1 !ohne -band! gelaufen sein -> ggf. unter single programs nochmal laufen lassen.
  • TIPP zur SUM Mit SUM: kann man (bis zu 6) DOS, die oben in der Datei mit gerechnet werden müssen, summieren. Die Summen wandern in Dateien *.dossum*. Ein Beispiel:

    Title
     -0.50 0.001 1.500 0.002  # EMIN, DE, EMAX, Gauss-broadening(>;de)
      15    N   0.000          # NUMBER OF DOS-CASES below, G/L/B broadening (Ry)
       0    1   total        # atom, case=column in qtl-header, label
       1    1   Fe11
       2    1   Fe12
       15    1   S11  
       16    1   S12
       17    1   S21
       18    1   S22
       19    1   S31
       20    1   S32
       21    1   S41
       22    1   S42
       23    1   S51
       24    1   S52
       25    1   S61
       26    1   S62
       SUM: 3 6      # Zahl der Summen, maximal Zahl von DOS-Zeilen von oben
       4 5  
       6 7 8 9
       10 11 12 13 14 15
    

• spaghetti
  • Generate k-mesh legt die Richtungen im k-Raum fest, entlang derer die Bandstruktur geplottet werden soll. Diese k-Liste kann mit dem Programm XCrySDen festgelegt werden. In diesem Programm die Richtungen anklicken und als xcrysden.klist speichern. XCrySDen macht häufig merkwürdige k-Listen, bei als Muster verfügbaren Bravaisgittertypen lieber diese verwenden. Falls keine verfügbar sind, ein Bild beim Bilbao-Server betrachten und die Liste selber verzeugen.

    ACHTUNG: Bei zentrierten Gittern ist die BZ des konventionellen Gitters kleiner (reziproker Raum!) als bei primitive setting. D.h. man sieht zwar bei conventionell die Richtungen besser, aber es fehlen Teile des k-Raums! Am besten erstmal genau das anschauen. Dann in primitive setting die Richtungen so auswählen, dass sie von den Richtungen zu den conventionell passen (z.B. wir bei einem I-Gitter aus 0,1/2,0 (d.h. der konventionellen b-Achse) 1/4,-1/4,1/4. Geplottet werden muss dann aber bis 1/2,-1/2,1/2!!. Unbedingt support-k-path dazu anschauen!!

  • Create NAME.klist_band kopiert entweder die mit XCrySDen erzeugte k-Listen-Datei oder eine entsprechende Datei aus der Sammlung bekannter k-Listen nach NAME.klist_band. Die Datei NAME.klist vom SCF-Zyklus, die den IBZ erfaßt wird nicht überschrieben!
  • lapw1 -band
  • Die nächsten beiden Punkte der Liste (Programme irrep und lapw2 -band -qtl werden bei "normalen" Bandstrukturplots ausgelassen (nur für Fat-Band Plots bzw. wenn die Charaktere der Bänder gekennzeichnet werden sollen).
  • NAME.insp enthält diverse Parameter für den Bandstrukturplot. Diese sind in der Datei gut dokumentiert. Von Hand muß die Fermikante EF (wird oben aus der Datei NAME.scf eingeblendet) eingesetzt werden. Außerdem ggf. den Energiebereich (in eV voreingestellt) für den Plot auswählen.
  • spaghetti
  • plot bandstructure
  • TIPP: Der Plot landet in einer Datei NAME.spaghetti_ps. Diese läßt sich z.B. mit pstoedit in einer verarbeitbare xfig-Datei konvertieren. Der nicht ausgewählte Bereich wird in der PS-Datei durch die Bounding-Box abgeschnitten, d.h. die ganz tiefliegenenden Zustände sind dann in der xfig-Datei mit drin!
  • fat-band Plots sind extrem nützlich, um den Charakter der Bänder zu erkennen. Die Auswahl der jeweiligen Zustände erfolgt wie bei der DOS (s.o.) durch Angaben der Atome und der Nummer der Zustände aus der qtl-Datei (werden oben angezeigt).
  • WICHTIG: Vor weiteren Rechnungen (DOS, Elektronendichten usw.) muß nach spaghetti mit lapw1 wieder für den ganzen k-Raum gerechnet werden.
  • TIPP 2: Nachdem irrep gelaufen ist, sehen (manchmal?immer???) normale Bandstrukturen chaotisch aus. Dann einfach alles was im Verzeichnis irrep im Namen trägt, weglöschen und spaghetti neu laufen lassen.
• LOCAL ROTATION MATRIX oder Who is who?
  Sowohl für die Bandstrukturen als auch für die DOS-Plots ist wichtig, wie die lokale Symmetrie an den einzelnen Atomen (beschrieben durch die Punktgruppe der Lage) mit den Kristallachsen zusammenhängt, wie also z.B. die einzelnen p-Orbitale im Kristall orientiert sind. Die Transformationsmatrix steht für alle Atomlagen einzeln (und ANDERS!) in der Datei NAME.outputs (und lt. Blaha NICHT in NAME.outputsgroup).
  
  • Die Symmetrieelemente werden analysiert bzgl. ihrer Richtungen, z.B.
    • m n z = Spiegelebene senkrecht (normal) zur z-Achse
    • 2 || z = 2-zählige Achse parallel zur z-Achse
  • und dann auf Standards der Punktgruppe umdefiniert
  • mit der als LOCAL ROTATION MATRIX angegebenen Matrix nach:

    		M11	M12	M13
    [ x y z ]	*	M21	M22	M23	=	[x-neu y-neu z-neu]
    		M31	M32	M33

  • Wichtig: Die Richtungen/Orbitale müssen dabei von vorne als Zeile mit der Matrix multipliziert werden.
  • Am besten zur Sicherheit prüfen, ob die Symmetrieelemente nach der Transformation der Reihenfolge der gewünschten Punktgruppe entsprechen.
• Elektronendichten
  Mit Wien2k lassen sich verschiedene Analysen der Elektronendichten im Realraum durchführen. Da die Dichte natürlich an den Atomkernen derart groß ist, dass man nicht mehr viel von den Bindungselektronen erkennen könnte, werden am besten die gewünschten Zustände über das Energiefenster ausgeschnitten (sog. PED, partielle Elektronendichten).
  • Vor allen Elektronendichteanalysen in der DOS (bzw. in *.outputt) schauen, welche Elektronen man erfassen will.
  • Untermenü electron density auswählen
  • Die gewählten minimalen und maximalen Energien (Default-Maximum ist E_F) in die Datei NAME.in2 (bzw. in die Maske) eintragen.
    !!! Diese Werte sind in Ry und nicht relativ zu E_F=0 !!!
    Achtung: Vor irgendwelchen anderen Rechnungen nie vergessen, diesen Wert wieder zurückzusetzen auf die voreingestellten -9.0 Ry!
  • lapw2 -p laufen lassen
  • Zum der Plotten der Elektronendichtekarten gibt es verschiedene Möglichkeiten:
    1. XcrysDen
       
       • Calculate Electron density with xcrysden ist am einfachsten, weil man in diesem Programm sowohl Schnitte (2D) als auch eine 3D-Dichte (3D) sowie totale bzw. Valenz-Dichte rechnen kann. Das Programm bietet gleich zu Anfang die beiden Möglichkeiten zur Auswahl:
         • 2D Plot:
           • Durch Auswahl dreier Atome und ggf. Erweiterung der rechteckigen Ebene in alle Richungen wird bei eine bestimmte Auflösung (0.01 A reicht i.a. aus) vorgegeben.
           • Im folgenden Menü kann zwischen VALenz (Valenzdichte) und RHO (Gesamtdichte?) (!? nur VAL sieht korrekt aus, bei RHO ergibt sich eine Unstetigkeit in der Elektronendichte am Kugelrand?!) und der Einheit e^-/A³ bzw. e^-/Bohr³ ausgewählt werden.
           • Dann wird lapw5 gestartet.
           • Nach dessen Ende erhält man das Menü zurück mit der Möglichkeit, die Ebene zu colorieren mit der Dichte. Also: Ebene einblenden und Bereich für die Elektronendichte auswählen (z.B. für typische Bindungselektronendichten: 0.01 bis 0.3 e/A³).
         • 3D Plot
       • TIPP: Die Verarbeitung der Elektronendichteplots für Publikationen (keine bunten Blubberbildchen möglich), am besten mit 2D-Plot das Optimale rausholen (Ausschnitt, Höhenlinienzahl/Abstand usw.). Dann die noch bunte PS-Datei mittels pstoedit -f fig name.ps name.fig in eine xfig-Datei konvertieren. Rest im xfig weiterverarbeiten. Vorteil: Es ist und bleibt ein Vektorformat!
    2. wien2venus und DrawXTL/VESTA
       • Das Phython-Skript wien2venus.py wirft eine Serie von lapw5-Jobs an.
       • Aufruf direkt gefolgt von den Voxeln, also z.B.
         wien2venus.py 200 200 200.
       • Die Ausgabedatei *.rho3d kann dann im DrawXtl unter Map Parameters eingelesen werden.
       • Plots der Dichten für bestimmte Energiebereiche erfordern die zusätzliche Option -V, trigonale Strukturen zusätzlich noch -t.
       • Betrachten der Maps geht auch mit VESTA ganz gut.
    3. 3ddens und DrawXTL/VESTA
       • Ab WIEN2k_19 ist die Berechnung der Dichte (Ausgabe in *.xsf-Datei) möglich mit x 3ddens (-val)
       • Dazu braucht man eine Eingabedatei *.in3d:

         300 300 300 # number of 3D-gridpoints in a-, b- and c-direction
         0.0 0.0     # extend (pos.)/reduce (neg.) ploting grid by some fraction of the lattice vectors
         0.0 0.0     # 2nd line: a-direction, 3rd line: b-direction, 4th line: c-direction
         0.0 0.0     # all values set to 0.0: use full unit cell
         add         # no: use file 9; add: sum of files 9 and 11; sub: difference of 9 and 11
         conv        # calculate density in conventional (conv) or primitive (prim) cell
         

       • Die Ausgabe *.xsf [mit den Dichten in e^-/Bohr³] kann ...
         • ... in XcrysDen eingelesen werden: xcrysden --xsf *.xsf Im XcrysDen dann Tools -> Datagrid verwenden.
           !! Die aus XcrysDen exportieren *.xsf-Dateien sehen in VESTA nicht o.k. aus!??
         • ... mit VESTA dargestellt werden.
           • Properties -> Isosurface -> dort Level anpassen
         • ... mit DrawXTL geplottet werden (als MAP zusätzlich zum *.struct-file einlesen).
         • Dazu ggf. in *.xsf noch die Syntax ändern:

           BEGIN_BLOCK_DATAGRID_3D
            9:clmval+++11:clmvaldn
           BEGIN_DATAGRID_3D_
           

           und am Ende analog:

           END_DATAGRID_3D_
           END_BLOCK_DATAGRID_3D 
           

         • Achtung! die MAP-INFO ist aktuell noch falsch, kann aber einfach ignoriert werden.
• Für die Darstellung von Fermiflächen startet man xcrysden von der Kommandozeile und importiert dann WIEN2k-Files. Zum Plot der Fermifläche wählt man die Bänder an der Fermienergie aus. Die Darstellung erfolgt für jedes Band extra oder zusammen (MERGE), per Reiter auswählbar).
• Mit AIM (Bader Analysis) (unter Single Program) kann man die Ladungsverteilung nach Bader berechnen. Den Input muss man von Hand zusammenstricken, z.B.

  SURF
  1                                (atom in center of surface !! including MULT)
  200 0.0 3.1415926536             (theta, 200 points, from zero to pi)
  200 0.0 6.283185                 (phi)
  0.04 1.0 4                       (step along gradient line, rmin, check)
  1.45 0.04                        (initial R for search, step (a.u))
  3 3 3                            (nshell)
  IRHO                             ("INTEGRATE" rho)
  WEIT                             (WEIT (surface weights from case.surf), NOWEIT)
  30                               (30 radial points outside min(RMIN,RMT))
  END
  

  Da zur Berechnung von AIM die bindungskritischen Punkte berechnet werden, kann man diese unter :PC in der Datei *.outputaim finden (in Bohr-basierten Einheiten) bzw. mit dem Skript extractaim_lapw in die Datei critical_points_ang (in Angstrom-basierten Einheiten) rausschreiben lassen.
  Alle Informationen stehen auch in der Zeile, die mit :PC beginnt, z.B.:
  

  :PC  0.5852 0.1704 0.2538 -0.1709E-01 -0.1543E-01 0.5221E-01 -1 0.1969E-01 0.3520E-01 -2 2.521 -3 2.530

  Die letzten vier Zahlen sind die Atome mit Abständen zum BCP, direkt davor steht die Dichte (in Bohr³), d.h. *6.755 ergibt die Dichte in A³ (hier: 0.03520 * 6.755 = 0.238).

  Eine nette Alternative zur Bestimmung von bindungskritischen Punkten, Ladungsverteilungen und sonstigen Analysen der Elektronendichte (oder auch anderem) ist das Programm critic2, das auch mit Elk-Output klarkommt.

• Aufteilungsraster der Brillouin-Zone: Zahl der k-Punkte, einmal angegeben als Gesamtraster über die gesamte Brilloiun-Zone (sog. Monkhorst-Pack (MP)-Gitter [X*Y*z]) und - als eigentlich wichtige Größe - die Zahl der k-Punkte im irreduziblen Teil der BZ (k/IBZ). Motto: Je mehr desto gut (s.o.)
• Als Kriterien für die Güte der Behandlung innerhalb der MT-Kugeln gibt es verschiedene sogenannte lm-Parameter (für jede kristallographische Atomlage separat!)
  • l_max: größte Drehimpulsquantenzahl für die Entwicklung der Wellenfunktion.
  • l_ns,max größte Drehimpulsquantenzahl für Nichtdiagonalelemente der Matrix des Hamiltonoperators.
  • LM_max größte Dreimpulsquantenzahl für die Entwicklung des Potentials.
• Als Kriterien für das Interstitium zwischen den MT-Kugeln, wo ebene Wellen (PW) verwendet werden, sind zu beachten:
  • Abschneideenergie E_max^pot. für das Potential (als G²_max im Programm enthalten; G_max ist voreingestellt auf 14, das sollte eigentlich immer reichen! Der Einfluß von Änderungen dieses Parameters auf die Ergebnisse der Rechnungen ist außerdem i.a. sehr klein.
  • Dagegen ist die Abschneideenergie E_max^wf (gelegentlich auch als cutoff bezeichnet) für die Entwicklung der Wellenfunktion im Interstitium genauer zu beachten. Sie steckt in Wien in der Größe RMT*KMAX. Voreingestellt ist 7.0 in der 2. Zeile der Datei NAME.in1. RMT*KMAX wird ev. vom Programm wegen Feldeinschränkungen runtergesetzt, daher sollte der real verwendete Wert in der Datei NAME.??? überprüft werden. Aus diesem Wert kann man bei bekannten RMT die Abschneideenergie berechnen. Es gilt mit guter Genauigkeit: K_max² = Abschneideenergie in Ry. Der Wert von RMT*KMAX durchläuft ein Minium, d.h. man kann nicht einfach sagen, dass beliebig groß immer gut ist. Das liegt daran, dass bei zu großem RMT*KMAX entsprechend mehr Parameter auftauchen und damit die Gefahr falscher Minima und von Korrelationen groß ist. Wenn alle Atome ähnliche RMT haben und bei metallischen Systemen reicht meistens 7.0 aus. Am besten den Einfluß von RMT*KMAX durch Ändern des Wertes ausprobieren! Unter 7.0 sollte man nicht gehen. Als kritische Obergrenze gilt 9.