Spin-Polarized Calculations for bcc Fe

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

(Jupyter notebook by Mara Voiculescu & Martin Kuban)

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.

Table of Contents

0. Before Starting

1. Unpolarized Ground State

2. Spin-Polarized Ground State

3. Anti-Ferromagnetic Phase

  • Comparison with the Non-magnetic (NM) Case

  • Exercises

4. Restarting Magnetic Calculations

Exercises

Literature

0. Before Starting

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

Before running any Jupyter tutorials, please refer to the 00_before_starting.md document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting's root directory:

source tools/excitingjupyter/venv/excitingvenv/bin/activate

As a first step, you may create a running directory for the notebook.

In [1]:
%%bash
mkdir -p run_spin_tutorial

1. Unpolarized Ground State

The system on which we will focus is iron in the body-centered cubic phase. Therefore, inside the running directory create a new directory noSPIN for the unpolarized calculations.

In [2]:
%%bash
cd run_spin_tutorial
mkdir -p noSPIN
cd ..

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"
         nwdos="1000"
         winddos="-0.3 0.3"
         inttype="tetra"/>
      <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>

The next step is writing the complete input as a string and saving it in your working directory as input.xml.

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

In [4]:
%%bash
cd run_spin_tutorial/noSPIN
python3 -m excitingscripts.setup.excitingroot
cd ../..

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 excitingscripts.execute.single.

In [5]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.execute.single -r noSPIN
cd ..

Plotting the DOS

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 excitingscripts.plot.dos, which is fully described in The python script "plot.dos". Execute the script as follows.

In [6]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.plot.dos -d noSPIN  -t 'Spin-unpolarized Fe-bcc'  -nl
cd ..

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

Plotting the Band Structure

In order to visualize the electronic band-structure, you can use the script excitingscripts.plot.band_structure, which is described in The python script "plot.band_structure". In this case you can execute the following command.

In [9]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.plot.band_structure -d noSPIN  -e -10 40  -t 'Spin-unpolarized Fe-bcc'  -nl
cd ..

The output file PLOT.png will look as follows.

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.

You may start by creating a subdirectory SPIN in which you will run the calculation.

In [7]:
%%bash
cd run_spin_tutorial
mkdir -p SPIN
cd ..

The next step is to modify the input file input.xml for making it suitable for spin-polarized calculations and to save it in the subdirectory SPIN. 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 by executin the script excitingscripts.execute.single.

In [1]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.execute.single -r SPIN
cd ..

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

Furthermore, you can visualize the new electronic band-structure and density of states by proceeding as for the spin-unpolarized case.

In [2]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.plot.dos -d SPIN  -t 'Spin-polarized Fe-bcc'  -e -8 8  -db -3.5 3.5  -nl
cd ..

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

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.

In [5]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.plot.band_structure -d SPIN  -e -10 40  -t 'Spin-polarized Fe-bcc'  -l 'center left'
cd ..

The result will be:

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.

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 input file input.xml in a new subdirectory AFM. The file should contain the following the following lines, corresponding to the hypothetical AFM phase of bcc Fe.

<input>

   <title>Anti-ferromagnetic Fe-bcc</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"
         nwdos="1000"
         winddos="-0.3 0.3"
         inttype="tetra"/>
      <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>
In [1]:
%%bash
cd run_spin_tutorial
mkdir -p AFM-B_4.0
cd ..
In [36]:
%%bash
cd run_spin_tutorial/AFM-B_4.0
python3 -m excitingscripts.setup.excitingroot
cd ../..

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.

In [37]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.execute.single -r AFM-B_4.0
cd ..

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

In [4]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.plot.dos -d AFM-B_4.0  -t 'Anti-ferromagnetic Fe-bcc'  -e -8 8  -db -3.5 3.5  -nl
cd ..

The electronic band-structure can be obtained as follows.

In [11]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.plot.band_structure -d AFM-B_4.0  -e -10 40  -t 'Anti-ferromagnetic Fe-bcc (SC)'
cd ..

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.

The values of the local magnetic moments can be found at the end of the INFO.OUT file in the directory AFM-B_4.0. You should obtain something like:

...
Moments :
     interstitial                           :         0.00000000
     moment in muffin-tin spheres :
                  atom     1    Fe          :        -0.70467731
                  atom     2    Fe          :         0.70467731
     total moment in muffin-tins            :         0.00000000
     total moment

As you can see the magnetic moment is localized (with opposite sign) on the two atoms. Only the non-vanishing component of the local moments along the z-axis is shown in the file.

Comparison with the Non-magnetic (NM) Case

In order to compare the energy and the other properties discussed for the antiferromagnetic system with a non-magnetic (NM) calculation, we change the input file input.xml by setting to zero the value of the initial local magnetic field:

...
         <atom coord="0.00000 0.00000 0.00000" bfcmt="0.0 0.0  0.0"/>
         <atom coord="0.50000 0.50000 0.50000" bfcmt="0.0 0.0 -0.0"/>
...
In [8]:
%%bash
cd run_spin_tutorial
mkdir -p AFM-B_0.0
cd ..

You can execute the NM calculation using the command

In [33]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.execute.single -r AFM-B_0.0
cd ..

The results contained in the directory AFM-B_0.0 can be compared with the ones obtained above for the antiferromagnetic phase. For example in the case of the density of states by typing the command

In [38]:
%%bash
cd run_spin_tutorial
python3 -m excitingscripts.plot.dos -d AFM-B_4.0 AFM-B_0.0  -t 'Anti-ferromagnetic Fe-bcc'  -e -8 8  -db -5 5
cd ..

you get the following plot

Exercises

  • Repeat the calculation by changing the length of the local magnetic field to "2.0". Analyze the results as described above. What does happen in this case?
  • Change again the length of the local magnetic field to "8.0" and "16.0". Do you see any variation?

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, ngrik) 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