# lmgf tutorial

### Summary

This package implements the ASA local spin-density approximation using Green’s functions. The Green’s functions are contructed by approximating KKR multiple-scattering theory with an analytic potential function. The approximation to KKR is essentially similar to the linear approximation employed in band methods such as LMTO and LAPW. It can be shown that this approximation is nearly equivalent to the LMTO hamiltonian without the “combined correction” term. With this package a new program, lmgf is added to the suite of executables. lmgf plays approximately the same role as the LMTO-ASA band program lm: a potential is generated from energy moments $Q_0$, $Q_1$, and $Q_2$ of the density of states. in the same way as the lm code. You can use lmgf to make a self-consistent density as you can do with lm. lmgf is a Green’s function method: Green’s functions have less information than wave functions, so in one sense the things you can do with lmgf are more limited: you cannot make the bands directly, for example. However lmgf enables you do do things you cannot do with lm. The two most imprortant are:

• Calculate magnetic exchange interactions
• Calculate magnetic susceptibility (spin-spin, spin-orbit, orbit-orbit parts)
• Calculate properties of disordered materials, either chemically disordered or spin disorder from finite temperature, within the Coherent Potential Approximation, or CPA.
• Calculate the ASA static susceptibility at $q=0$ to help converge calculations to self-consistency.

You can find some extra information on the way lmgf works in lmgf documentation.

lmgf is a Green’s function program complementary to the ASA band code lm. For some properties, e.g. calculating moments $Q_{0..2}$ lmgf can be straightforwardly substituted for lm because both calculate the DOS. The DOS is $1/(2\pi ) {\mathrm{Im}} G$: it can be decomposed into site contributions and thus moments Q0..2 can be generated for each site and l channel, as an alternative to decomposing the eigenfunctions of the bands, as lm does. Thus it can achieve self-consistency in a manner similar to lm, but generating $Q_{0..2,{\mathbf R}l}$ by an alternate route. If the ASA hamiltionian built by lm is suitably simplified, i.e. by

• omitting the “combined correction term” (OPTIONS_ASA_CCOR)
• generating $Q_{0..2,{\mathbf R}l}$ from true power moments as the Green’s function does (HAM_QASA=0),

then lmgf and lm will produce nearly identical self-consistent solutions. When potential functions are parameterized to 2nd order in both lm and lmgf, and both methods are fully k converged, they should product nearly identical results. By default lm parameterizes the potential function to 3rd order; lmgf can do the same. The 3rd order parameterizations are similar in the two methods, but not identical. To verify this, try the following test:

gf/test/test.gf co 1 2   ← Test 1 for 2nd order parameterization; test 2 for 3rd order


lmgf is a bit messier to work with (Green’s functions are harder to stabilize than wave functions), and it a bit less accurate as the simplifications to lm amount to approximations. So, typically lm makes a better self-consistent potential.

But lmgf can do things lm doesn’t do, e.g. calculate magnetic exchange interactions through linear response as this tutorial demonstrates. Sometimes there is a need or advantage to carrying self-consistency via lmgf, e.g. when performing CPA calculations. Unless there is good reason to do otherwise, it is better use the self-consistent potential generated by lm to calculate other properties such as the magnetic exchange parameters. We follow that strategy here.

### Preliminaries

For this tutorial the blm, lmchk, lmstr, lm and lmgf executibles are required and are assumed to be in your PATH; the source code for all Questaal exectuables can be found here.

### Tutorial

##### 1. Building input file

Before starting working with this tutorial we advise you to read through the ASA-tutorial which explains building the imput file in more details (you can also look through the input described in a full-potential context). In the present tutorial we’ll focus on the part of the input specific for using with lmgf.

To get started, copy doc/demos/asa-copt/init.copt to your working directory. Inspect the init file and you will see it contains just the minimum structural information, apart from one line supplying some information about the magnetic structure:

LATTICE
ALAT=7.1866
PLAT=    1.000000   0.000000   0.000000
0.000000   1.000000   0.000000
0.000000   0.000000   1.000000
SPEC ATOM=Co MMOM=0,0,2.2
SITE
ATOM=Pt POS=  0.0  0.0   0.0
ATOM=Co POS=  0.5  0.5   0.0
ATOM=Co POS=  0.0  0.5   0.5
ATOM=Co POS=  0.5  0.0   0.5


Then use the blm tool (described in more details in ASA-tutorial and full-potential tutorial. )

blm --mag --nk=8 --asa --gf copt


blm should generate file actrl.copt, which should be essentially the same as doc/demos/asa-copt/ctrl.copt (commented lines might be different though). If this is not the case, something is wrong with your configuration. At the moment if actrl looks slightly different from the one provided just move on to the next steps.

The command-line arguments are not required, but they supply quantities blm cannot determine automatically, that you will have to supply at some point. If you supply them on the command-line they are folded into the ctrl file at the outset; or, you can edit the ctrl file after it is generated. Command-line switches blm recognizes are summarized in Building FP input file.

###### The --asa switch

This switch tailors the ctrl file for the ASA. To see how it affects the ctrl file, try running blm without --asa. For more details see the ASA-tutorial.

###### The --mag switch

This switch tells blm that you plan on doing a spin polarized calculation. All it does is change the preprocessor variable nsp to 2. This turns on the spin polarization through NSPIN={nsp}.

Without any other information the spin polarized calculation will proceed with zero magnetic moment. The system needs a “push” in the initial direction to find the magnetic state. You have to supply some initial information about the magnetic structure. Since we know that the magnetization is concentrated on the Co (Pt is paramagnetic, though it has a high magnetic susceptibility), the init file supplies an initial magnetic moment on the Co site of about 2 Bohr on the Co d orbital, in the SPEC category (SPEC ATOM=Co MMOM=0,0,2.2 in the initial file). The precise value 2.2 is not important: this quantity is determined self-consistently later. Choosing it rather large (the bulk moment is 1.8 $\mu_B$) gives it a strong initial push so to encourage it not revert to a (metastable) nonmagnetic state in the course of a self-consistent calculation.

###### The --gf switch

When --gf is used, blm prepares the input file for the Green’s function program lmgf. This tutorial uses lmgf to calculate magnetic exchange interactions. Adding --gf to the blm command line argument modifies actrl.copt in two ways:

1. The GF category is created:

% const gfmode=1 c3=t
GF      MODE={gfmode} GFOPTS={?~c3~p3;~p2;}


To see the purpose of GF_MODE, do:

lmgf --input


and look for GF_MODE. You should see:

    GF_MODE           reqd   i4       1,  1          default = 0
0: do nothing
1: self-consistent cycle
10: Transverse exchange interactions J(q), MST
11: Read J(q) from disk and print derivative properties
...


So, if MODE=1, lmgf does a self-consistent calculation, generating the P and $Q_{0..2}$ for each l channel using Green’s functions rather than wave functions as lm does.

GFOPTS bundles a variety of lmgf-specific options, which you supply through a sequence of strings separated by semicolons. This tag:

GFOPTS={?~c3~p3;~p2;}


becomes GFOPTS=p3 after parsing by the preprocessor, because c3 is nonzero (see preprocessor documentation). p3 tells lmgf to use $3^{rd}$ order potential functions (somewhat more accurate than $2^{nd}$ order, but also prone to generating false poles not too far from the real axis).

2. EMESH is added to BZ:

% const nz=16 ef=0
EMESH={nz},10,-1,{ef},.5,.3  # nz-pts;contour mode;emin;emax;ecc;bunching


Green’s functions are energy resolved; thus physical properties such as the charge density or magnetic exchange interactions require an integration over the energy as well as over the BZ. For both density and static exchange interactions, the integration must be taken on the real axis from below the lowest eigenstate in the system to the Fermi level $E_F$. Im G is basically the density-of-states. It is very spikey on the real axis, and a very fine energy mesh would be required to integrate Im G close to the real axis. The integration can be accomplished with vastly greater ease by deforming the contour into an elliptical path in the complex plane. A gaussian quadrature is used; typically 15 or so energy points is sufficient for a well converged result.

This contour is specified through EMESH. Breaking down the constituents of EMESH as autogenerated by blm:

     EMESH={nz}          ←  number of energy points in the contour; {nz} evaluates 16 in this file
10           ←  elliptical contour
-1           ←  starting energy on the contour.  Must be deeper than the lowest state in the system (-0.776 Ry)
{ef}          ←  Fermi level determined by charge neutrality; see below
0.5           ←  eccentricity of the ellipse ranging from 0 (circle) to 1 (line)
0.3           ←  bunching parameter, bunching points near Ef. 0→no bunching


We don’t know what $E_F$ is a priori. In the ASA, a general reasonable guess is 0. If we perform the band calculation, see ASA-tutorial, we get $E_F$ generated by lm: it is −0.12927 Ry. $E_F$ is fixed by charge neutrality. If lmgf generated exactly the same spectrum as lm, and the k-integration were fully converged (or at least identical in the two cases) $E_F$ would be the same for lm as for lmgf. However we can expect that the charge neutrality points will slightly different in the two methods. Further we’ll find $E_F$ using lmgf.

If you want to practice finding $E_F$ using lm use the following commands (for the details see ASA-tutorial) (We advise you to do those steps since you’ll need some of them further anyway but the following lm-part contains some useful comments):

Invoking blm with the switches given above is sufficient to make a working input file. Normally you can copy actrl.copt to ctrl.copt as it is.

All the ASA electronic structure codes (lm, lmgf, and lmpg) use a tight-binding form of the LMTO basis, where the envelope functions are screened to make them short ranged. This information is carried through screened structure constants, which in this package are precomputed and stored using lmstr. Run this setup to make the structure constants:

 lmstr ctrl.copt                             ← Make and store structure constants


It should store str.copt and sdot.copt on disk. (If not, something is wrong and you should not proceed.)

Invoking blm with the switches given above is sufficient to make a working input file. Normally you can copy actrl.copt to ctrl.copt as it is.

As of yet we have no starting density or potential. You can see this immediately by trying to run the band code straight off:

lm ctrl.copt


The program will stop with this message:

LM:  Q=ATOM encountered or missing input


In usual LDA calculations, a trial density is obtained by generating densities for free atoms, and superposing them (Mattheis construction). While the ASA code could have been written to do just this, it does something different. This code takes advantage of the simplification the ASA offers, namely that the sphere density is completely determined by a small number of parameters, namely the log derivative parameters P and energy moments of the charge density $Q_{0..2}$ for each l channel. We can supply reasonable guesses through the ctrl file, or let the program pick some defaults as a first guess. Defaults are typically assigned so that $Q_0$ is the charge in the l channel of the atom and $Q_{1..2}$ are taken to be zero. While this is a pretty crude guess (cruder than the Mattheis construction) usually it is good enough that the program can find its way to the proper self-consistent solution.

The ASA code can either start from “potential parameters”, which gives it enough information to generate energy bands and calculate moments ($P_l$, $Q_{0..2,l}$) , or from the moments ($P_l$, $Q_{0..2,l}$)  which is sufficient for the sphere routine to fix the potential and calculated potential parameters. The band and sphere blocks of the code are thus complementary: one takes the input of the other and generates output required by the other. The cycle is described in the LMTO-ASA documentation.

The ctrl file is built with the following START category:

START CNTROL={nit==0} BEGMOM={nit==0}


If BEGMOM is nonzero, lm will start from potential parameters (which don’t exist yet, in the present case). If BEGMOM=0 lm will start from the ($P_l$, $Q_{0..2,l}$). These haven’t been given either, but lm can pick defaults for them. We get an initial potential by doing:

lm ctrl.copt -vnit=0         ← Because -vnit=0, BEGMOM={nit} is preprocessed into BEGMOM=0


lm will start from (default) moments and generate a trial density for each sphere, together with potential parameters corresponding to potential generated.

The output should generate a table of potential parameters like this:

PPAR:  Pt        nl=4  nsp=2  ves=  0.00000000
l     e_nu          C        +/-del     1/sqrt(p)      gam         alp
...
1 -0.33739981  0.66438341  0.17542338   6.2239786  0.13462478  0.13462478
2 -0.21536750 -0.17914256  0.02841817   1.1299420  0.01358564  0.01358564
...
1 -0.33739981  0.66438341  0.17542338   6.2239786  0.13462478  0.13462478
2 -0.21536750 -0.17914256  0.02841817   1.1299420  0.01358564  0.01358564
...


and a similar table for Co. Particularly important are C, the band center of gravity C and the bandwidth del. You can see that $C_p$ sits far above zero while $C_d$ is a little below. It tells you that the Pt d orbital is important for bonding while the Pt p orbital is pretty far above from the Fermi level and of much less importance. del is the bandwidth parameter; a little more detail is given in the LMTO-ASA documentation. A disk file is created for each class. It contains the ($P_l$, $Q_{0..2,l}$), the potential parameters, and possibly other things. Take a look at files co.copt and pt.copt. You can see what defaults were chosen for ($P_l$, $Q_{0..2,l}$).

We are now ready for a self-consistent calculation. Doing:

lm ctrl.copt -vnit=30 --pr31,20                ← NIT={nit} is preprocessed into NIT=30.
--pr31,20 sets verbosity fairly low


will perform up to 30 self-consistent cycles, that is

lm will continue until the RMS change in ($P_l$, $Q_{0..2,l}$) falls below tolerance ITER_CONVC, or until 30 iterations is reached.

If it’s converged you’ll get the following phrase at the end of your output:

Jolly good show! You converged to ...


In this demo convergence should be reached in 21 iterations.

Interpreting the output:

The output can provide some very useful information. For example, the self-consistent Co moment is 1.8 $\mu_B$; Pt has a moment (induced by the neighboring Co) of 0.356 $\mu_B$. You can see it in the line of the following form

ATOM=PT Z=78 Qc=68 R=2.928343 Qv=-0.009029 mom=0.356 a=0.025 nr=491
ATOM=Co ...


Scrolling up you can find the density-of-states at the Fermi level $E_F$ is D(Ef)=79.729 (units of $Ry^{-1}$ per unit cell), or about 1.3 $eV^{-1}/atom$. Had the calculation been done without spin polarization, D(Ef) would be ~187, more than twice larger. This is a very large number and suggests there is a likely instability. Indeed, the system can lower its energy by spontaneously magnetizing. Consider the Stoner criterion for spontaneous magnetization, I D(EF) > 1. In 3d transition metals I is about 1 eV. Thus the Stoner criterion is easily satisfied and the system should spontaneously magnetize. In magnetizes so strongly that the Co moment (1.8 $\mu_B$) is larger than that for elemental Co (1.6 $\mu_B$).

The same line also provides you with the Fermi energy:

BZINTS: Fermi energy:   -0.129264; ...

##### 2. Making ctrl file and structure constants

If you’ve done the lm-part above go straight to ‘The Green’s function program lmgf’

Invoking blm with the switches given above is sufficient to make a working input file. Normally you can copy actrl.copt to ctrl.copt as it is.

For a fuller description of the ctrl file, see the the ASA-tutorial, the FP tutorial, and also Building FP input file.

All the ASA electronic structure codes (lm, lmgf, and lmpg) use a tight-binding form of the LMTO basis, where the envelope functions are screened to make them short ranged. This information is carried through screened structure constants, which in this package are precomputed and stored using lmstr. Run this setup to make the structure constants (for more detailes see see the the ASA-tutorial):

 lmstr ctrl.copt                             ← Make and store structure constants


It should store str.copt and sdot.copt on disk. (If not, something is wrong and you should not proceed.)

As of yet we have no starting density or potential.We get an initial potential by doing:

lm ctrl.copt -vnit=0


We are now ready for a self-consistent calculation. Do:

lm ctrl.copt -vnit=30 --pr31,20

##### 3. The Green’s function program lmgf

###### a) Finding $EF$

If GF_MODE=1, lmgf will generate the $Q_{0..2,{\mathbf R}l}$ for whatever $E_F$ you give it. However there is only one physically meaningful $E_F$ – the one that satifies charge neutrality. The input file is constructed so you can supply $E_F$ through command-line argument -vef=expr: the preprocessor evaluates  ef  from  expr, substitutes it for  {ef}  in the input file (see preprocessor documentation). We’re going to use the one obtained by running lm (see above).

The simplest way to find the charge neutrality point is to run lmgf interactively in the self-consistent mode (GF_MODE=1). By running lmgf interactively you can monitor convergence. Do:

lmgf ctrl.copt -vnit=30 --pr31,20 --iactiv -vef=-.1293


Since we’re using --iactiv switch the code is going to stop and ask us to make some choices. At first you’ll see

QUERY: max it (def=30)?


Just hit Enter (return) to confirm. Output contains two tables the first of which looks like

GFASA:  integrated quantities to efermi = -0.1293
PL      D(Ef)      N(Ef)       E_band      2nd mom      Q-Z
spin 1   11.450485   21.527413   -7.139582    2.654697    3.027413
spin 2   79.610384   15.411861   -4.696173    1.699842   -3.088139
total   91.060868   36.940675  -11.835755    4.354540   -0.060726
deviation from charge neutrality: -0.060726


The non-zetro deviation from charge neutrality means that ef=-.1293 results in a slight electron deficiency. lmgf will estimate a constant shift to crystal potential to make the system neutral, and interpolate G to contour adjusted by this shift using a Pade approximant.

Note lmgf shifts the average crystal potential: ef  is kept fixed.

Then lmgf prints out some results of the Pade correction in a subsequent table.

Corrections to integrated quantities estimated by Pade interpolation
PL      D(Ef)      N(Ef)       E_band      2nd mom      Q-Z
spin 1    8.860480   21.530197   -7.151787    2.654754    3.030197
spin 2  110.802293   15.469803   -4.712034    1.700808   -3.030197
total  119.662773   37.000000  -11.863960    4.355561   -0.000000
deviation from charge neutrality: 0


At the prompt you should see

QUERY: redo gf pass (def=F)?


It is asking you whether you want to accept the Pade approximant, or redo the GF calculation with the potential shift added. Let’s do the latter, to see how good the estimate was. At the prompt type ‘st’ and hit return twice

QUERY: redo gf pass (def=F)?  st <RET> <RET>


After the cycle you should see

deviation from charge neutrality: 0.015781


The Pade correction reduces the deviation from neutrality but overestimates the shift. A new estimate is made for the potential shift and the prompt reappears. You can repeat the GF cycle as many times as you like. (If you see QUERY: modify vbar (def=…)? just hit return) If you iterate enough you should see something like:

gfasa:  potential shift this iter = 0.000001.  Cumulative shift = -0.000433


The shift of this last iteration is negligible.  Cumulative shift’  is the net shift accumulated over all the iterations.

You can now proceed to self-consistency but we will instead use the potential generated by lm in order to make the exchange parameters. (If you do proceed to self-consistency using lmgf, note that it writes the potential shift and Fermi level to file vshft.copt This shift is automatically read when lmgf is restarted. After self-consistency is reached you can either keep vshft.copt, or remove the file and modify  ef  so charge neutrality is satisfied without the shift.)

At the prompt enter

QUERY: redo gf pass (def=F)?  a <RET>


to prevent lmgf from continuing its self-consistency cycle. The constant potential shift is just the negative the the requisite Fermi level shift to achieve charge neutrality:  ef needs to be adjsted to -0.1293−(-0.000433) = −0.128867 Ry.

To confirm that this is the correct  ef, repeat the interactive lmgf calculation with  -vef=-0.1289.

###### b) Magnetic Exchange Interactions

As we mentioned before lmgf requires a GF-specific category (look into the ctrl.copt).

GF  MODE=1 GFOPTS=options


Token MODE= controls what lmgf calculates. Options are MODE=1, MODE=10, MODE=11, MODE=26.

Look into ctrl.copt. Two lines are important here:

% const gfmode=1 c3=t
GF      MODE={gfmode} GFOPTS={?~c3~p3;~p2;}


MODE={gfmode} means that you can define MODE in the command line by adding -vgfmode=1/10/11/26; if you don’t it will be set to 1 (from const gfmode=1). In the previous example we used MODE=1 now we’ll need MODE=10 that invokes a special branch computing magnetic exchange interactions using a linear response technique.

The Heisenberg model is an empirical model that postulates a set of interacting rigid local spins. The Hamiltonian is

The $J_{RR'}$ are called “Heisenberg exchange parameters”. The Heisenberg applies to a system of rigid spins undergoing small excursions about equilibrium. R and R’ are any pair sites and $J_{RR'}$ is a kind of magnetic analog to the dynamical matrix describing small oscillations of nuclei around their equilibrium point. In a crystal with periodic boundary conditions $J_{RR'}$ can be Bloch transformed to read:

$J_{RR'}(q) = \sum_T \exp(iqT) J_{R+T,R'}$, where R and R’ are now confined to sites within a unit cell.

lmgf calculates $J_{RR'}(q)$ from the “Lichtenstein formula.” This famous expression (J. Magn. Magn. Mater. 67, 65 (1987)), closely related the static transverse magnetic susceptibility $\chi^{+-}$, is derived from density functional perturbation theory. It establishes a first-principles basis for the Heisenberg model. One elegant (though approximate) feature of the ASA is that the magnetization is everywhere associated with an atomic sphere. For local moment systems, the magnetization is well confined inside a sphere; thus associated with every site R there is a well defined local moment. If sufficiently localized it rotates rigidly under the influence of an external perturbation.

When you set GF_MODE=10, lmgf will generate $J_{RR'}(q)$, and then perform an inverse Bloch transform (by Fast Fourier Transform) to make $J_{R+T,R'}$ for as many lattice translation vectors T as there are k-points.

Do

lmgf -vgfmode=10 ctrl.copt -vef=-.1289


Results are saved in file jr.copt (see below).

Most of the analysis is done in the next step, but already the output from gfmode=10 contains some useful information. In the first of this pair of tables you see J_0 and 2/3 J_0. J_0 is the net Weiss magnetic field from the surrounding neighbors; 2/3 J_0 would be the (classical) mean-field estimate for the critical temperature $T_c$ if there were one atom/cell. Since the Pt moment is very small it is weakly magnetic and has little effect on $T_c$. In the second table (J_0 resolved by L) J_0 is decomposed into lm contributions. As expected, the contributions to J_0 originates almost entirely from the d states.

Now if you run lmgf with GF_MODE=11, it reads jr.copt (which means you have to run MODE=10 first) and does some analysis with the parameters. Invoke lmgf with

lmgf -vgfmode=11 ctrl.copt -vef=-.1289


A unit cell of N sites has $N^2$ pairs. Thus jr.copt holds a succession of $N^2$ tables of J, one array for each RR’ pair in the unit cell. Each array has $n_1n_2n_3$ exchange parameters, corresponding to the lattice translation vectors that follow from the Fast Fourier Transform of a k mesh of $n_1n_2n_3$ points. You can find the headers for each array (headers follow a standard format this package uses) by doing, e.g.

grep rows jr.copt


to see:

% rows 64 cols 8 real  rs  ib=1  jb=1
% rows 64 cols 8 real  rs  ib=1  jb=2
...


ach array has 64×8 entries, for T vectors derived from 8×8×8 k-points (the 3D array is stored in a 2D format). lmgf unpacks these (GFMODE=11) and prints them out in a sequence of tables, e.g. this one coupling all pairs of atoms belong to sites 2 and 3 in the unit cell. Pairs are ordered by separation distance d. Interactions fall off rapidly with d, and oscillate around 0, as might be expected from RKKY theory. Then follow estimates for the critical temperature $T_c$. $T_c$ is estimated in Weiss mean-field theory, and also according to a spin-waves theory by Tyablikov (sometimes called the “RPA”). Mean-field tends to overestimate $T_c$; RPA tends to be a little more accurate but tends to underestimate it. From these two estimates $T_c$ should be around 1000K (see the GFMODE=11 output).

Next follows an estimate for the spin wave stiffness. We need a symmetry lines file, let’s copy it from /lm/startup/ (make sure you have the correct path there. The folder has symmetry files for different structures, e.g. syml.fcc, syml.hcp, etc):

cp startup/syml.sc syml.copt


You can also create syml.copt

41  0.5 0 0    0 0 0                  X to Gamma
41   0 0 0    .5 .5 0                 Gamma to M
41  .5 .5 0   .5 .5 .5                M to R
41  .5 .5 .5   0 0 0                  R to Gamma
0 0 0 0 0 0 0
... For fermi surface plotting
vx     range  n     vy     range  n   height  band
1 0 0   -1 1   51   0 1 0   -1 1   51   0.00    4


lmgf reads this file and calculates the spin wave spectrum from the Heisenberg model, along the lines specified. Results are saved in bnds.copt (the energy scale is now mRy). So let’s run GFMODE=11 again to get it (now we have the symmetry lines file):

lmgf -vgfmode=11 ctrl.copt -vef=-.1289


You can plot magnon spectra using the same technology you use for plotting energy bands, see ASA-tutorial. If you have the plbnds and fplot packages installed do:

echo 0 350 5 10 | plbnds -scl=13.6 -fplot -lbl=X,G,M,R,G bnds.copt
fplot -f plot.plbnds
`

Now you have the fplot.ps file with the spectra (you can rename this file to some-name.ps) and use your favorite postscript reader to view it. You should see something close to what is shown in the figure:

.

Magnon energies are in meV.

Note: The 8×8×8 mesh is a bit coarse. Use a finer k mesh for a smoother and more accurate magnon spectrum.