Noncollinear fcc Iron
This tutorial illustrates noncollinear magnetism in lmf, using fcc iron as an example.
Here we take a four-atom supercell of fcc Fe, with lattice vectors in a simple cubic structure. Using the ASA code and spin statics, the stable configuration was found to be one where the spin quantization axis pointed along some permutation of the (1,1,1) direction, so that the angles between each spin are arccos(1/√3), or 109.5o.
The FP code does not as yet implement spin statics, so we take the ASA spin configuration as input.
Setup: files init, eula, ctrl, site, and basp
We will use blm, the input file generator, to build the input file, ctrl.fccfe. Following this tutorial, we build it from an init file specifying structure, in this case a simple cubic cell of Fe with four atoms.
Paste the box below into init.fccfe
HEADER init file for fccfe constructed from questaal test case.
LATTICE
% const a=6.7794173
ALAT={a}
PLAT= 1.0000000 0.0000000 0.0000000
0.0000000 1.0000000 0.0000000
0.0000000 0.0000000 1.0000000
SPEC
ATOM=Fe MMOM=0,0,2
SITE
ATOM=Fe POS= 0.0000000 0.0000000 0.0000000
ATOM=Fe POS= 0.5000000 0.5000000 0.0000000
ATOM=Fe POS= 0.0000000 0.5000000 0.5000000
ATOM=Fe POS= 0.5000000 0.0000000 0.5000000
Note there is a species category; this is needed for non-structural information, in this case the initial trial Fe magnetic moment (2 μB on the d orbital).
We also need to specify the Euler angles. The following specifies Euler angles for each site. The lmchk instruction below shows that they point along permutations of [1,1,1] directions, with 109.5o between spins. Paste the box below into eula.fccfe.
pi/4 acos(1/sqrt(3)) -pi/4
-3*pi/4 acos(1/sqrt(3)) 3*pi/4
3*pi/4 pi-acos(1/sqrt(3)) -3*pi/4
-pi/4 pi-acos(1/sqrt(3)) pi/4
Make the ctrl file (input file).
blm --brief --mag~nc~spinor=3 --forces=0 --autobas~ehmx=-.4~eloc=-4.5~loc=5 '--mix~nit=30~mixp=b4,b={beta},k=7' --nk~10 --gw~nk=8 fccfe
This instruction makes the main input file, ctrl.fccfe; the site file specifying structure, site.fccfe; and finally basp.fccfe, the file defining the basis set.
Command-line switches for blm are documented on this page. Specifically:
- –brief : A moderately simple ctrl file, but one with embedded variables to enable command line control of some of the input.
- –mag~nc~spinor=3 : Calculation is to be magnetic and noncollinear.
- –forces=0 : No forces are to be calculated
- –autobas~ehmx=-.4~eloc=-4.5~loc=5 : Sets the basis set.
This switch makes up the semicore 3p local orbital, included as a valence state (not needed, but the calculation is more accurate). - ’–mix~nit=30~mixp=b6,b={beta},w=1,1,wnc=.1,k=6,n=6’ : String for the ITER_MIX tag that controls how the density is mixed.
This exercise does not converge well if blm’s default parameters are used. (This is taken up below.) - –nk~8 : Sets the number of k divisions (8×8×8 k mesh). If you want the calculation to run more quickly, you can reduce the the mesh
- –gw~nk=8 : Enables the ctrl file for a GW calculation. This tutorial doesn’t carry out GW calculations, but the switch enlarges the basis set.
The Fe moment can very sensitive to details.
To confirm that Euler angles in file eula.fccfe orient spins along the [1,1,1] directions, do the following:
lmchk --euler ctrl.fccfe
You should see that the spin quantization axis points along one permutation of [±0.577350, ±0.577350, ±0.577350].
This table should also be recorded in stdout. It shows that the angles between spins are 109.5o :
1 2 3 4
-----------------------------------
1 -- 109.47 109.47 109.47
2 109.47 -- 109.47 109.47
3 109.47 109.47 -- 109.47
4 109.47 109.47 109.47 --
Partial self-consistency
The noncollinear calculation proceeds in the same way as a collinear one. First make the atomic densities to be used as a trial starting point for the self-consistency cycle
lmfa ctrl.fccfe
This calculation is too heavy for a scalar job, so you will want to use MPI. Assign variable mn to the MPI instruction, e.g. for a machine with 16 processors,
mn='mpirun -n 16'
Make this system (partially) self-consistent
rm -r mixm.fccfe ! only needed if restarting a calculation
$mn lmf --rs=0 fccfe > llmf
The ctrl is built to iterate to a maximum 30 iterations. Confirm that the density is reasonably self-consistent:
grep DQ llmf
You should see something like this:
mixrho: sought 6 iter, 11362 elts from file mixm; read 0. RMS DQ=4.45e-2 mixrho: sought 6 iter, 11362 elts from file mixm; read 1. RMS DQ=3.66e-2 last it=4.45e-2 mixrho: sought 6 iter, 11362 elts from file mixm; read 2. RMS DQ=1.01e-2 last it=3.66e-2 mixrho: sought 6 iter, 11362 elts from file mixm; read 3. RMS DQ=5.40e-3 last it=1.01e-2 mixrho: sought 6 iter, 11362 elts from file mixm; read 4. RMS DQ=2.89e-3 last it=5.40e-3 mixrho: sought 6 iter, 11362 elts from file mixm; read 5. RMS DQ=2.09e-3 last it=2.89e-3 mixrho: sought 6 iter, 11362 elts from file mixm; read 0. RMS DQ=5.02e-3 last it=2.09e-3 mixrho: sought 6 iter, 11362 elts from file mixm; read 1. RMS DQ=1.93e-2 last it=5.02e-3 mixrho: sought 6 iter, 11362 elts from file mixm; read 2. RMS DQ=1.08e-3 last it=1.93e-2 mixrho: sought 6 iter, 11362 elts from file mixm; read 3. RMS DQ=3.36e-3 last it=1.08e-3 mixrho: sought 6 iter, 11362 elts from file mixm; read 4. RMS DQ=7.29e-4 last it=3.36e-3 mixrho: sought 6 iter, 11362 elts from file mixm; read 5. RMS DQ=7.36e-4 last it=7.29e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 0. RMS DQ=7.63e-4 last it=7.36e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 1. RMS DQ=1.73e-3 last it=7.63e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 2. RMS DQ=7.52e-4 last it=1.73e-3 mixrho: sought 6 iter, 11362 elts from file mixm; read 3. RMS DQ=8.56e-4 last it=7.52e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 4. RMS DQ=1.31e-3 last it=8.56e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 5. RMS DQ=7.38e-4 last it=1.31e-3 mixrho: sought 6 iter, 11362 elts from file mixm; read 0. RMS DQ=7.46e-4 last it=7.38e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 1. RMS DQ=1.14e-3 last it=7.46e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 2. RMS DQ=7.33e-4 last it=1.14e-3 mixrho: sought 6 iter, 11362 elts from file mixm; read 3. RMS DQ=8.23e-4 last it=7.33e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 4. RMS DQ=9.48e-4 last it=8.23e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 5. RMS DQ=7.07e-4 last it=9.48e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 0. RMS DQ=7.49e-4 last it=7.07e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 1. RMS DQ=1.02e-3 last it=7.49e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 2. RMS DQ=7.26e-4 last it=1.02e-3 mixrho: sought 6 iter, 11362 elts from file mixm; read 3. RMS DQ=7.90e-4 last it=7.26e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 4. RMS DQ=8.88e-4 last it=7.90e-4 mixrho: sought 6 iter, 11362 elts from file mixm; read 5. RMS DQ=7.00e-4 last it=8.88e-4
The density is close enough to a fully converged result to be nearly identical. (A fully self-consistent density difficult to stabilize in this system.)
It is also interesting to see what the local moments are. Do the following:
tail -40 llmf
You should see something close to the following
Sphere charges and magnetization
ib q M Mx My Mz
1 in 13.217070 1.626168 0.939366 -0.938695 0.938544
out 13.218013 1.617144 0.933801 -0.933060 0.934114
diff 0.000944 -0.009024 -0.005565 0.005635 -0.004430
mix 13.217025 1.627597 0.940229 -0.939660 0.939192
2 in 13.217073 1.626160 -0.938736 0.939612 0.938243
out 13.216254 1.618504 -0.934130 0.934653 0.934547
diff -0.000819 -0.007656 0.004606 -0.004959 -0.003696
mix 13.217026 1.627452 -0.939465 0.940567 0.938796
3 in 13.216949 1.626500 -0.938942 -0.938565 -0.939673
out 13.219039 1.617189 -0.933932 -0.933309 -0.933812
diff 0.002090 -0.009311 0.005011 0.005255 0.005861
mix 13.216876 1.627889 -0.939671 -0.939466 -0.940449
4 in 13.217268 1.621191 0.935667 0.936227 -0.936091
out 13.216993 1.613792 0.931615 0.932085 -0.931470
diff -0.000275 -0.007399 -0.004051 -0.004143 0.004621
mix 13.217211 1.622463 0.936384 0.937047 -0.936757
-----
total 52.868138 6.505401 -0.002523 -0.001512 0.000782
The local Fe moment is about 1.6μB, somewhat smaller than 2.2μB found in bcc Fe.
Energy bands with spin texture
To plot bands on the specified lines, create a symmetry lines file. You can do this by hand, or let lmchk create one for you. We have a simple cubic lattice, i.e. bcc Fe in a four-atom supercell. However, the underlying lattice is bcc, let’s make symmetry lines for it, keeping in mind that four k points in the bcc Brillouin zone will be folded into a single k point.
lmchk ctrl.fccfe --syml~n=50~lblq:G=0,0,0,L=1/2,1/2,1/2,X=1,0,0,W=1,1/2,0,K=3/4,3/4,0~lbl=LGXWGK
Create the bands using generating color weights with spin texture. The following will generate four color weights, projecting the charge and three components of spin onto the Fe d orbitals
$mn lmf fccfe --band@scolst~z=26~l=d@fn=syml
This creates an energy bands file, bnds.fccfe, in an old fashioned format. Inspect the first line of the file. You should see something close to
264 0.00107 4 spin-texture lbl=LGXWGK col1=5:9,21:25,29:33,38:42,54:58,62:66,71:75,87:91,95:99,104:108,120:124,128:132 col2= col3= col4=
This line tells you that the basis consists of 264 orbitals, with 4 color weights corresponding to spin texture. The color weights correspond to projections onto d orbitals, with charge, and the three vector components of the magnetization.
One way to view the spin texture is to use the colors to depict projections onto x, y, and z. In the following instruction spin projections are depicted as three colors (black=, red=, green=), similar to the example for Co in the plbnds documentation.
Do the following
plbnds -range=-9,5,5,10 -wscl=.6,.6 -fplot~sh~ext=z~lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0~spintexture -ef=0 -scl=13.6 bnds.fccfe
This instruction creates a file fplot.ps, which you can view with your favorite postscript viewer. You should see something similar to the following. (Note that the Γ and X points are identical because of zone folding.)
You can see how the spin quantization axis evolves with k. Particularly interesting are points where the color changes abruptly (see for example midway between L and X at −6 eV), or where two bands of different texture cross.
Fermi surface with spin texture
This section shows how to draw the Fermi surface using colors to reveal the evolution of spin texture.
lmf can create a file with selected bands plotted on a regular 2D mesh (--band option in mesh mode). Data on a regular mesh is the input needed to draw 2D contour plots of constant energy, which if chosen as the Fermi level, is the Fermi surface. This tutorial shows how to use the 2D mesh mode of lmf with fplot to make a figure with a Fermi surface.
To get started, you need to specify an area (square, rectangle, or parallelepiped) lying on a plane in the Brillouin zone), the number of divisions on each axis, and know specify which bands to store.
For the latter, 48 bands lie strictly below the Fermi surface in this system, and the next 12 cross it at some point. To see this run lmf with --minmax :
$mn lmf fccfe -vnit=1 --rs=1,0 --minmax:ef0
--minmax:ef0 causes lmf to display a table with largest and smallest eigenvalues in each band, shifted to align the Fermi level with 0. Thus any band where the maximum value is less than 0, or the minimum value is greater than 0, does not cross the Fermi level. You should find a table something like this
Lowest, highest evals by resolved band, spin 1, ef shifted to 0 Ry 1 -3.9578 -3.9532 | 89 1.0628 1.7466 | 177 4.2055 4.8588 47 -0.1440 -0.0399 | 135 2.8685 3.4373 | 223 6.1189 6.7170 48 -0.1440 -0.0397 | 136 2.8685 3.4373 | 224 6.1191 6.7188 49 -0.1327 0.0344 | 137 2.9664 3.4373 | 225 6.3931 7.3532 50 -0.1326 0.0345 | 138 2.9665 3.4373 | 226 6.3945 7.3538 51 -0.1325 0.0346 | 139 2.9749 3.4726 | 227 6.3947 7.3556 52 -0.1324 0.0347 | 140 2.9749 3.4727 | 228 6.3953 7.3560 53 -0.0890 0.0347 | 141 3.1198 3.6165 | 229 6.4532 7.4404 54 -0.0890 0.0349 | 142 3.1198 3.6166 | 230 6.4541 7.4432 55 -0.0815 0.0386 | 143 3.1202 3.6166 | 231 6.7443 7.5163 56 -0.0815 0.0389 | 144 3.1202 3.6166 | 232 6.7450 7.5182 57 -0.0135 0.0426 | 145 3.2639 3.6167 | 233 6.7551 7.7141 58 -0.0135 0.0427 | 146 3.2640 3.6168 | 234 6.7562 7.7153 59 -0.0133 0.0878 | 147 3.2654 3.6386 | 235 6.9025 7.7166 60 -0.0133 0.0878 | 148 3.2656 3.6387 | 236 6.9054 7.7182 61 0.0051 0.0951 | 149 3.3522 3.6387 | 237 7.0465 7.8693 ...
Bands 49:60 cross EF.
The --band switch in mesh mode requires an input file defining the mesh (syntax is described on this page). We will choose a square with the Γ point at the center, abscissa along x and ordinate along y, dividing each axis into 35 mesh points (more mesh points gives better definition to the figure).
Paste the box below into file fs.fccfe:
# vx range n vy range n height band
1 0 0 -.5 .5 35 0 1 0 -.5 .5 35 0.00 49:60
Generate the bands in mesh mode, specifying spin texture for color weights:
$mn lmf fccfe --band~con~scolst@atom=fe@l=2~fn=fs
Inspect the contents of bnds.fccfe. You should see a sequence of 2D arrays in the standard Questaal format for 2D arrays, stacked one after another. There are 60 such data sets: the first 12 are energy bands 49:60 on a 35×35 mesh. The next 12 are weights for the charge density; the next 12 are weights for the x component of magnetization, then y component, then z component.
You can use any facility to draw the Fermi surface. Here we use the contour mode facility in fplot. For each band, fplot will require 4 files: one for the eigenvalues, and three for x, y, z components of the magnetization.
To simplify matters, we will extract the bands and shift them by the Fermi level, to put it at zero. Extract the Fermi level with this
grep efermi= bnds.fccfe | head -1 | vextract . efermi
You could see something similar to 0.001066, which we will use for the Fermi level.
Here we will used the mcx calculator to read bnds.fccfe for each of 12 data sets, containing eigenvalues and write them into files 49, 50, …, 60. Also we will write weights into files x49, x50, …, x60, y49, y50, …, y60, z49, z50, …, z60, for x, y, z components of the the weights.
Set the following shell variables:
rf="-r:open -qr bnds.fccfe -shft=-0.001066"
rfw="-r:open -qr bnds.fccfe"
and run mcx as follows:
mcx [ i=49:60 -qr $rf -w \{i} ] [ i=49:60 -qr $rfw ] [ i=49:60 -qr $rfw -abs -s5 -w x\{i} ] [ i=49:60 -qr $rfw -abs -s5 -w y\{i} ] [ i=49:60 -qr $rfw -abs -w z\{i} ] -show
This should create 60 files, named as described above.
Finally we need an fplot script to draw the figure. Paste the box below into plot.fs :
fplot -frme 0,1,0,1 -x -.5,.5 -y -.5,.5 -tmx .1
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x49~colw2=0,1,0~wtf2=y49~colw3=0,0,1~wtf3=z49 -0.00 49
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x50~colw2=0,1,0~wtf2=y50~colw3=0,0,1~wtf3=z50 -0.00 50
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x51~colw2=0,1,0~wtf2=y51~colw3=0,0,1~wtf3=z51 -0.00 51
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x52~colw2=0,1,0~wtf2=y52~colw3=0,0,1~wtf3=z52 -0.00 52
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x53~colw2=0,1,0~wtf2=y53~colw3=0,0,1~wtf3=z53 -0.00 53
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x54~colw2=0,1,0~wtf2=y54~colw3=0,0,1~wtf3=z54 -0.00 54
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x55~colw2=0,1,0~wtf2=y55~colw3=0,0,1~wtf3=z55 -0.00 55
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x56~colw2=0,1,0~wtf2=y56~colw3=0,0,1~wtf3=z56 -0.00 56
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x57~colw2=0,1,0~wtf2=y57~colw3=0,0,1~wtf3=z57 -0.00 57
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x58~colw2=0,1,0~wtf2=y58~colw3=0,0,1~wtf3=z58 -0.00 58
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x59~colw2=0,1,0~wtf2=y59~colw3=0,0,1~wtf3=z59 -0.00 59
-lt 1,bold=4 -con~colwt=1,0,0~wtfn=x60~colw2=0,1,0~wtf2=y60~colw3=0,1,0~wtf3=z60 -0.00 60
This script draws 12 contour plots for contours at zero energy, one for each band, with three weights
Invoke fplot with
fplot -f plot.fs
The instruction should create a file fplot.ps, which you can view with your favorite postscript viewer. You should see something similar to the following. Red, green, blue correspond to spin oriented along x, y, z, and yellow is a combination of x and y.
It is interesting to see where bands nearly cross (e.g. around (0.2, 0.3)) and the spin quantization axis abruptly changes around the crossing point.
Additional Exercises
Shifting the Fermi level slightly radically alters the Fermi surface. Try modifying the energy in script plot.fs a little (say from (-0.00 to -0.01) see how the Fermi surface changes.
You can confirm that the noncollinear total energy is deeper than the collinear one. To see this, remove eula.fccfe and repeat the self-consistent calculation. You should see that the collinear case is less bound.
Getting the density to go fully self-consistent is exceedingly difficult. You can do it with some special mixing modes. One way accomplish this is to perform iterations in two steps. In the first, modify the mixing mode so that it mixes the spin density keeping the charge density fixed. Insert the following line just above the line beginning with mix= b6,b …
mix= b4,b={beta},w=0,1,k=10,n=10(Place it on the line before the existing tag; the ctrl file parser reads the first occurence of this tag.) Run lmf again
$mn lmf fccfe -vnkabc=10 fccfe > out.2This instruction increases the k mesh a little, which stabilizes convergence. You can see how lmf proceeds to convergence by invoking:
grep DQ out.2
Next, remove or comment out the new constrained mixing line just inserted and run lmf again.
$mn lmf fccfe -vnit=200 -vnkabc=10 fccfe > out.3
lmf proceeds to convergence, albeit very slowly. Monitor it with
grep DQ out.3