Spin-Polarized Calculations for bcc Fe

by Dmitrii Nabok, Andris Gulans, Pasquale Pavone, & Benedikt Maurer for exciting oxygen

Purpose: In this tutorial you will learn how to initialize and perform spin-polarized calculations. As an example, we calculate electronic ground-state properties of bcc-Fe with different types of magnetic order.


0. Define relevant environment variables

Read the following paragraphs before starting with the rest of this tutorial!

Before starting, be sure that relevant environment variables are already defined as specified in Tutorial scripts and environment variables.

From now on the symbol $ will indicate the shell prompt.

Important note: All input parameters are in atomic units!

In order to start the tutorial create somewhere a directory spin-tutorial and move inside it.

$ mkdir spin-tutorial
$ cd spin-tutorial

1. Unpolarized ground state

The system on which we will focus is iron in the body-centered cubic phase. Therefore, inside the directory spin-tutorial create a new directory Fe-bcc and move inside it.

$ mkdir Fe-bcc
$ cd Fe-bcc

Let us first start from the spin-unpolarized case. Below is an example for the required input file (input.xml).

<input>
 
   <title>Spin-unpolarized Fe-bcc</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      <crystal scale="5.200">
         <basevect> 0.5  0.5 -0.5</basevect>
         <basevect> 0.5 -0.5  0.5</basevect>
         <basevect>-0.5  0.5  0.5</basevect>
      </crystal>
      <species speciesfile="Fe.xml" rmt="2.0">
         <atom coord="0.0 0.0 0.0"/>
      </species>
   </structure>
 
   <groundstate 
      do="fromscratch"
      xctype="GGA_PBE_SOL"
      rgkmax="8.0"    
      ngridk="8 8 8"
      stype="Gaussian"
      swidth="0.01"
      nempty="10">
   </groundstate>
 
   <properties>
      <dos
         nsmdos="2"
         ngrdos="300"
         nwdos="1000"
         winddos="-0.3 0.3"/>
      <bandstructure>
         <plot1d>
            <path steps="300">
               <point coord="0.0    0.0    0.0 " label="Gamma"/>
               <point coord="0.5   -0.5    0.5 " label="H"/>
               <point coord="0.0    0.0    0.5 " label="N"/>
               <point coord="0.0    0.0    0.0 " label="Gamma"/>
               <point coord="0.25   0.25   0.25" label="P"/>
               <point coord="0.5   -0.5    0.5 " label="H"/>
            </path>
         </plot1d>
      </bandstructure>
   </properties>
 
</input>

Make sure to set $EXCITINGROOT to the correct exciting root directory in the speciespath attribute using the command

$ SETUP-excitingroot.sh

Let us look closer at this input file. In the groundstate element some attributes are explicitly given. Be aware that for time-saving reasons, here we have chosen computational parameters which do not correspond to those devised from convergence tests, although producing reasonable description of the system under consideration. Some of the attributes are known from other tutorials. The attributes stype and swidth are specified explicitly because they are relevant for metallic systems (see Input Reference for their definition). In fact, their choice is crucial when performing convergence tests with respect to the k-point number for metals. The meaning of the attribute nempty will be clear later when we initialize spin-polarized calculations.

Furthermore, inside the input file you can find the element properties which contains the subelements bandstructure and dos. The presence of these two subelements allows for the calculation of the electronic band-structure and density of states, respectively (see Electronic band-structure and density of states).

To perform actual calculations, you may run the script EXECUTE-single.sh.

$ EXECUTE-single.sh noSPIN

Using this command all results of the calculation will be stored after the execution in the subdirectory noSPIN. The next step would be to visualize the resulting density of states. For this, you can use the script PLOT-dos.py, which is described in details in The python script "PLOT-dos.py". Execute (inside the directory Fe-bcc) the script as follows.

$ PLOT-dos.py -d noSPIN  -t 'Spin-unpolarized Fe-bcc'  -nl

This script produces a PNG (PLOT.png) output file. In this examples you will obtain the following plot.

fe-bcc-noSPIN-dos.png

In order to visualize the electronic band-structure, you can use the script PLOT-band-structure.py, described in The python script "PLOT-band-structure.py". In this case you execute the following command.

$ PLOT-band-structure.py -d noSPIN  -e -10 40  -t 'Spin-unpolarized Fe-bcc'  -nl

The output file PLOT.png will look as follows.

fe-bcc-noSPIN-bs.png

2. Spin-polarized ground state

In the previous example the spin of the electrons is not taken into account. However, it is well known that iron is a magnetic material. In this kind of materials spin up electrons behave differently from spin down ones. In this section you will learn how to include spin degrees of freedom in exciting calculations.

Inside the directory Fe-bcc, please open your favorite text editor in order to modify the input file input.xml for making it suitable for spin-polarized calculations. First of all, you may want to change the label of the calculation using

...
   <title>Spin-polarized Fe-bcc</title>
...

The main change is the inclusion of the element spin inside the groundstate element.

...
   <groundstate 
      do="fromscratch"
      xctype="GGA_PBE_SOL"
      rgkmax="8.0"    
      ngridk="8 8 8"
      stype="Gaussian"
      swidth="0.01"
      nempty="10">
 
      <spin 
         bfieldc="0.0 0.0 -4.0" 
         reducebf="0.5"
         spinorb="false">
      </spin>
 
   </groundstate>
...

The presence of the spin element will tell to exciting to initialize the spin-polarized calculations. The relevant attributes of the spin element are described in the following Table.

attribute description
bfieldc Specifies the cartesian coordinates of the external magnetic field required in spin-polarized calculations to break spin symmetry. In case of collinear magnetism, it is important to note that the preferable direction of the magnetic field should be selected along the z-axis (the axis of the spin quantization). Another important notice is that the amplitude of this field may play a crucial role to converge the system to the state with the right magnetic moment. Thus, as in the case of bulk iron, one has to specify rather high value of the field to obtain the magnetic state.
reducebf If one searches the ground state in absence of magnetic field, this attribute specifies which portion of the starting field is acting at each step during the self-consistent cycle (see more in Input Reference).
spinorb Should be set to "true" to include the spin-orbit coupling ("false" by default).

The way in which exciting deals with spin-polarized systems (the so-called "second-variational approach", described, e.g., in Singh-2006) requires also the explicit specification in the input of the attribute nempty inside the groundstate element.

attribute description
nempty Determines the size of the "second-variational" Hamiltonian. The optimized value should be chosen on the basis of a convergence test.

As before, you may run the calculation using the script EXECUTE-single.sh.

$ EXECUTE-single.sh SPIN

All results of the calculation are now stored in the subdirectory SPIN. Then, move inside the directory SPIN.

$ cd SPIN

The main output file INFO.OUT contains information on the magnetic moments of the spin-polarized phase. These moments are written at each iteration. The values of the magnetic moments at the last SCF cycle should look like the following.

...
Moments :
     interstitial                           :        -0.00911718
     moment in muffin-tin spheres :
                  atom     1    Fe          :         1.98008333
     total moment in muffin-tins            :         1.98008333
     total moment                           :         1.97096615

Check the convergence of the total moment for during the SCF cycle this calculation using the script PLOT-convmoment.py as follows:

$ PLOT-convmoment.py .

Visualize the PNG (PLOT.png) output file.

Furthermore, you can visualize the new electronic band-structure and density of state. To do this,
move back to the parent directory Fe-bcc using the commands

$ cd ../
$ pwd
.../Fe-bcc

You can now proceed as for the spin-unpolarized case by typing

$ PLOT-dos.py -d SPIN  -t 'Spin-polarized Fe-bcc'  -e -8 8  -db -3.5 3.5  -nl

Thus, you will get the following plot for the density of states:

fe-bcc-SPIN-dos.png

Notice that the DOS for the spin-down electrons is represented by the negative values in the plot.

In order to visualize the electronic band-structure, you execute the following command.

$ PLOT-band-structure.py -d SPIN  -e -10 40  -t 'Spin-polarized Fe-bcc'  -l 'center left'

The result will be:

fe-bcc-SPIN-bs.png

As you can see the main difference with the spin-unpolarized case is now the explicit distinction between electrons with spin up and spin down.

Before completing this part of the tutorial, it is useful to make a safe copy of the input file

$ mv input.xml input-spin-polarized.xml

3. Anti-ferromagnetic phase

Now we will learn how to initialize calculations for materials with anti-ferromagnetic (AFM) type of spin ordering. As an example, we consider a hypothetical AFM phase of bcc iron. In oder to realize the anti-ferromagnetic order, we need at least two nonequivalent atoms per cell. To do this with we simply apply the (standard) definition of the bcc lattice as a simple cubic (SC lattice with a basis consisting of two atoms.

Let us now create a new empty input file input.xml in the actual directory (it should be Fe-bcc). Copy inside input.xml the following lines, corresponding to the hypothetical AFM phase of bcc Fe.

<input>
 
   <title>Anti-ferromagnetic Fe-bcc (SC)</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      <crystal scale="5.200">
         <basevect> 1.0  0.0  0.0</basevect>
         <basevect> 0.0  1.0  0.0</basevect>
         <basevect> 0.0  0.0  1.0</basevect>
      </crystal>
      <species speciesfile="Fe.xml" rmt="2.0">
         <atom coord="0.00000 0.00000 0.00000" bfcmt="0.0 0.0  2.0"/>
         <atom coord="0.50001 0.50001 0.50001" bfcmt="0.0 0.0 -2.0"/>
      </species>
   </structure>
 
   <groundstate
      do="fromscratch"
      xctype="GGA_PBE_SOL"
      rgkmax="8.0"    
      ngridk="6 6 6"
      stype="Gaussian"
      swidth="0.01"
      nempty="20"
      epsengy="5.d-6"
      niterconvcheck="1"> 
 
      <spin/>
 
   </groundstate>
 
   <properties>
       <dos
         nsmdos="2"
         ngrdos="300"
         nwdos="1000"
         winddos="-0.3 0.3"/>
      <bandstructure>
         <plot1d>
            <path steps="300">
               <point coord="0.0  0.0  0.0" label="Gamma"/>
               <point coord="0.0  0.5  0.0" label="X"/>
               <point coord="0.5  0.5  0.0" label="M"/>
               <point coord="0.0  0.0  0.0" label="Gamma"/>
               <point coord="0.5  0.5  0.5" label="R"/>
               <point coord="0.0  0.5  0.0" label="X"/>
            </path>
         </plot1d>
      </bandstructure>
   </properties>
 
</input>

Here, the relevant difference with respect to the other calculations is the presence of the additional attribute bfcmt inside the atom element. This attribute plays a similar role as the attribute bfieldc but only inducing the local atomic (muffin-tin) magnetic moment. The presence of this magnetic moment breaks the symmetry of the crystal and make the two basis atoms indeed nonequivalent. Since in this example we do not desire any external magnetic field or spin-orbit coupling, we can skip the explicit specification of the attributes inside the spin element, keeping them to the default values.

Notice also that the path for the electronic band-structure in this example corresponds to the simple cubic supercell.

Following the steps as in the previous sections, we can calculate the total energy of the system and the local atomic magnetic moments.

As before, after updating the species path

$ SETUP-excitingroot.sh

you may run the calculation using

$ EXECUTE.single.sh AFM

Once the run is finished, proceed as in the other sections to visualize the density of states.

$ PLOT-dos.py -d AFM  -t 'Anti-ferromagnetic Fe-bcc'  -e -8 8  -db -3.5 3.5  -nl
fe-bcc-AFM-dos.png

The electronic band-structure can be obtained as follows.

$ PLOT-band-structure.py -d AFM  -e -10 40  -t 'Anti-ferromagnetic Fe-bcc (SC)'
fe-bcc-AFM-bs.png

Here, you can notice that, due to the symmetry of the anti-ferromagnetic phase, both the DOS and the electronic band-structures corresponding to the two spins (up and down) are identical.


4. Restarting magnetic calculations

It is possible to restart a magnetic calculation using a previously saved potential in STATE.OUT as the initial guess. To do so, you have to specify the following attribute in the element groundstate:

   <groundstate
      do="fromfile"
      ...
   </groundstate>

The previous spin-polarized calculations in this tutorial have relied on the initial magnetic field to make sure that the self-consistent procedure converges to a spin-polarized solution. But you do not need this trick for a restarted calculation, because STATE.OUT contains an initial potential that is already close to a self-consistent one. You should not set the initial field exactly to 0, because it will trigger a non-magnetic calculation. Instead, you can set it a small value, such as 1d-6. In the calculation with spin-polarized iron, you would have to modify the input given in the Section 2 as follows:

...
      <spin 
         bfieldc="0.0 0.0 -1d-6"
         reducebf="0.5"
         spinorb="false">
      </spin>
...

In the anti-ferromagnetic case, make the following changes:

...
      <species speciesfile="Fe.xml" rmt="2.0">
         <atom coord="0.00000 0.00000 0.00000" bfcmt="0.0 0.0  1d-6"/>
         <atom coord="0.50001 0.50001 0.50001" bfcmt="0.0 0.0 -1d-6"/>
      </species>
...

If you repeat an already converged calculation using this restart option, you will able to reach self-consistence in 4 steps. But it will take more of them if you change some other parameters (for example, rgkmax, ngridk) before the restart.

Exercises
  • At the given lattice constants, perform complete convergence test of the total energy and (if applicable) magnetic moments for the different phases with respect to the values of the attributes ngridk, swidth, rgkmax, and nempty.
  • Perform the lattice optimization (see Volume optimization for cubic systems) for the paramagnetic, ferromagnetic, and anti-ferromagnetic phases. Compare the total energies for the optimized structures of the different phases. Which phase has the minimum total energy per atom? How does the energetics of the different phases depend on the main computational parameters?
  • Calculate the magnetic moments of the optimized structure of the ferromagnetic phase. How does the obtained total magnetic moment compare with the one calculated in Section 2. and with the experimental value of about 2.12 ╬╝B?
  • In order to visualize the volume dependence of total magnetic moment (magnetization) in the ferromagnetic phase, use the script PLOT-totalmoment.py inside the directory in which the volume optimization for this phase was performed.

Literature

  • Singh-2006: David J. Singh and Lars Nordstr├Âm, "Planewaves, Pseudopotentials, and the LAPW Method, Second Edition", Springer 2006
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License