Personal tools
You are here: Home Documentation Help files for v5.7 tutorial ABINIT. Tutorial spin


ABINIT, spin lesson of the tutorial:

Properties related to spin: spin polarized calculations and spin-orbit coupling.


This lesson aims at showing how to get the following physical properties:

  • the total magnetisation of a ferromagnetic material
  • the magnetisation of an antiferromagnetic material
  • analyse the total density of states per spin direction 
  • analyse the density of states per atom and per spin direction
  • look at the effect of spin-orbit coupling for a non magnetic system
  • non-collinear magnetism (not yet)
  •  spin-orbit coupling and magnetocristalline anisotropy (not yet)
You will learn to use of  features of ABINIT which deal with spin.

This lesson should take about (to be provided) hours to be done.

Copyright (C) 2005-2009 ABINIT group (GZ)
This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

Goto : ABINIT home Page | Suggested acknowledgments | List of input variables | Tutorial home page | Bibliography
Help files : New user's guide | Abinis (main) | Abinis (respfn) | Mrgddb | Anaddb | AIM (Bader) | Cut3D | Optic | Mrgscr

 

Content of spin lesson

  • 1 A ferromagnetic material: the  bcc Fe
  • 2 An antiferomagnetic example: fcc Fe
  • 3 Another look at fcc Fe
  • 4 A ferrimagnetic example (Mn4N)(Not yet)
  • 5 The spin-orbit coupling
  • 6 Rotation of the magnetization and Spin-orbit (not yet)


1. A ferromagnetic material: the  bcc Fe

 

Before beginning, you might consider to work in a different subdirectory as for the other lessons. Why not "Work_spin" ?

The file ~abinit/tests/tutorial/Input/tspin_x.files lists the file names and root names. You can copy it in the Work_spin directory (and change it, when needed, as usual).
You can also copy the file ~abinit/tests/tutorial/Input/tspin_1.in in Work_spin. This is your input file. You can run the calculation, then you should edit the input file, and read it carefully.
Because we are going to perform magnetic calculations, there a two new types of variables:

Please, read their description in the help file.
You will work at fixed ecut (=18Ha) and k-point grid, defined by kptrlatt (the 4x4x4 Monkhorst-Pack grid). It is implicit that in "real life", you should do a convergence test with respect to both convergence parameters.
This takes about 30 seconds on a 2.8 Ghz PC, but one needs a sufficiently large cut-off to exhibit magnetic effects.
We will compare the output with and without magnetisation. (Hence, there are two datasets in the run)
We now look at the output file.
In the magnetic case, the electronic density is split in two parts, the "Spin-up" and the "Spin-down" parts to which correspond different Kohn-Sham potentials and different sets of eigenvalues whose occupation are given by the Fermi-Dirac function (without the ubiquitous factor 2)

For the first k-point, we get for instance:
(no magnetization)
       occ     2.00000   1.99989   1.99989   1.22915   1.22915   0.28676   0.00000   0.00000
(magnetic case)
       occ   1.00000   0.99999   0.99999   0.98396   0.98396   0.69467   0.00000   0.00000 (spin-up)
             1.00000   0.99730   0.99730   0.00898   0.00898   0.00224   0.00000   0.00000 (spin-down)

We note that the occupation are very different for up and down spins, which means that the eigenvalues are shifted which is in turn due to a shift of the exchange-correlation potential, and therefore of the effective potential. You can indeed have a look at the output file to compare spin-up and down  eigenvalues:
  -0.48411  -0.38615  -0.38615  -0.30587  -0.30587  -0.27293   0.33747   0.33747 (up, kpt#1)
 -0.46638  -0.32383  -0.32383  -0.21767  -0.21767  -0.20371   0.36261   0.36261 (dn, kpt#1)


The magnetization density (in unit of mu_B - Bohr's magneton) is the difference between the up and down densities. The magnetization density, divided by the total density is noted "zeta". This quantity "zeta" can vary between -1 and 1. It is everywhere zero in the non-magnetic case. In the magnetic case, we can read for instance its minimal and maximal values in the output file:

Min spin pol zeta= -4.8326E-02 at reduced coord.  0.7222  0.5000  0.2222
     next min= -4.8326E-02 at reduced coord.  0.5000  0.7222  0.2222
Max spin pol zeta=  5.7306E-01 at reduced coord.  0.0000  0.8889  0.8889
     next max=  5.7306E-01 at reduced coord.  0.8889  0.0000  0.8889

The total magnetization, its integral on the unit cell, is now in the magnetic case:

  Magnetisation (Bohr magneton)=  1.96743487E+00
 Total spin up =  4.98371743E+00   Total spin down =  3.01628257E+00

We observe that the total density (up+down) yields 8.000 as expected.

The magnetization density is not the only changed quantity. The energy is changed too, and we get:
 
    etotal1  -2.4661707269E+01  (no magnetisation)
    etotal2  -2.4670792869E+01  (with magnetisation)

The energy of the magnetized system is the lowest and therefore energetically favoured, as expected since bcc iron is a ferromagnet.
Finally, one notes also that the stress tensor is affected by the magnetization. This would also be true for  the forces for a less symmetric material.

It is interesting to consider in more details the distribution of eigenvalues for each direction of magnetization, which is best done by looking  at the respective densities of state.
To this end we have set prtdos= 1 in the input file, in order to obtain the density of states corresponding to spin-up and spin-down electrons (as soon as nsppol=2).
The values of the DOS are in the files tspin_1o_DS1_DOS and tspin_1o_DS2_DOS for the magnetic and non magnetic case respectively. We can extract the values for use in a plotting software.
Traditionally, in order to enhance visibility, one affects  with a negative sign the DOS of  minority  electrons. If we compare the DOS of the magnetized system and the non magnetized system, we observe that the up and down DOS have been  "shifted" with respect each other.
The integrated density of states yields the number of electrons for each spin direction and we see the magnetization which arises from the fact that there are more up than down electrons at the Fermi level.

That the magnetization points upwards is fortuitous, and we can get it pointing downwards by changing the sign of the  initial spinat.
Indeed, in the absence of spin-orbit coupling, there is no relation between the direction of magnetization and the cristal axes.
If we start with a spinat of 0, the magnetisation remains 0. spinat gives indeed a way to initially break the spin symmetry.



2. An antiferromagnetic example: fcc Fe

 

Well sort of....
Actually, Fcc Fe, displays many complicated structures, in particular spin spirals. A spiral is characterized by a direction along an axis, an angle of the magnetization with respect to this axis and a step after which a rotation is complete.
A very simple particular case is when the angle is 90°, the axis is <100> and the step is the unit cell side, spin directions are alternating in every two planes perpendicular to the <100> axis ("spiral stairway")
For instance, if the atom at [x,y,0] possesses an "up" magnetization, the atom at [x,y,1/2] would possess a down magnetization etc...
To describe such a structure, we merely need two atoms, [0,0,0] and [1/2,0,1/2] and a unit cell comprising those two atoms.
Also, each atom will be given opposite magnetization with the help of the variable spinat.

You can copy the file ~abinit/tests/tutorial/Input/tspin_2.in in Work_spin. This is your input file. Modify accordingly the tspin_x.files file. You can run the calculation, then you should edit the tspin_2.in file, and briefly look at the two changes with respect to the file tspin_1.in : the change of unit cell basis vectors rprim, the new spinat.
Note also we use now nsppol=1 and nspden=2 : this combination of values is only valid when performing an antiferromagnetic calculation. It uses the so-called Shubnikov symmetries, to perform calculations twice faster than with nsppol=2 and nspden=2
. The symmetry of the crystal is not the full fcc symmetry anymore since the symmetry must now preserve the magnetization of each atom.  Abinit is nevertheless able to detect such symmetry belonging to the Shubnikov groups and correctly detects that the cell is primitive, which would not be the case if we had the same value of spinat on each atom.
If we now run again the calculation, this calculation will take approximately 2 minutes on a 2.8 Ghz cpu.
If we look at the eigenvalues and occupations, they are again affected by a factor 2, which results from the symmetry considerations alluded to above, and not from the "usual" spin degeneracy : the potential for spin-up is equal to the potential for spin-down, shifted by the antiferromagnetic translation vector. Eigenergies are identical for spin-up and spin-down, but wavefunctions are shifted one with respect to the other.

  kpt#   1, nband= 16, wtk=  0.05556, kpt=  0.0833  0.0833  0.1250 (reduced coord)
  -0.60539  -0.47491  -0.42613  -0.39022  -0.35973  -0.34377  -0.28893  -0.28827
  -0.25318  -0.24042  -0.22944  -0.14218   0.20264   0.26203   0.26641   0.62158
      occupation numbers for kpt#   1
   2.00000   2.00000   2.00000   1.99997   1.99945   1.99728   1.50517   1.48016
   0.15706   0.04649   0.01576   0.00000   0.00000   0.00000   0.00000   0.00000



How do we know we have a magnetic order?
The  density of state used for bcc Fe  will not be useful since the net magnetization is zero and we have as many up and down electrons.
Though, the magnetization is reflected in the existence of an up and down electronic density whose sum is the total density and whose difference yields the net magnetization density at each point in real space.
In particular, the integral of the magnetization around each atom will give an indication of the magnetic moment carried by this particular atom.
To perform this integration, we will use the utility cut3d which yields an interpolation of the magnetization at any point in space. Reminder: you should build cut3d by issuing in ~abinit : make cut3d
As of now cut3d is interactive, we will use it through a very primitive script (written in Python) to perform a rough estimate of the magnetization on each atom.
You can have a look at the ~abinit/doc/tutorial/lesson_spin/gz2.py program, and note (or believe) that it does perform an integration of magnetization on a  cube of side acell/2 around each atom.
Copy it in your Work_spin directory. If you run the program, by typing "python gz2.py" , you will see the result:
magnetization of atom 1=  0.38211
magnetization of atom 2=-0.38202
which shows that the magnetization of each atom are actually opposite.
In a next version of the tutorial, we might use the cut3d library.
With the next input file tspin_3.in, we will consider this same problem, but in a different way, and note, for future record that the total energy is: Etotal=-4.92489586027187E+01





3. Another look at fcc Fe

Instead of treating fcc Fe directly as an antiferromagnetic material, we will not make any hypotheses
on its magnetic structure, and run the calculation like the one for fcc Fe, anticipating only that the two
spin directions are going to be different. We will not even assume that the initial spins are of the same magnitude.
You can copy the file ~abinit/tests/tutorial/Input/tspin_3.in in Work_spin. This is your input file. You can modify the file tspin_x.files, and immediately start running the calculation. Then, you should edit it. Note the values of spinat. In this job, we wish again to characterize the magnetic structure.
We are not going to use zeta as in the preceding calculation, but we will here use another feature of abinit: atom and angular momentum projected densities of state.
These are densities of state  weighted by the projection of the wave functions on angular momentum channels (that is spherical harmonics) centered on each atom of the system.
Note that theses DOS are computed with the thetrahedron method, which is rather time consuming and
produces less smooth DOS than the smearing method. The time is strongly dependent on the number of kpoints and we use here only a reduced set.
(This will take about 4 minutes on a 2.8 Ghz computer)
To specify this calculation we need new variables, in addition to prtdos set now to 3 :

This will specify the atoms around which the calculation will be performed, and the radius of the sphere.
We specifically select a new dataset for each atom, a non self-consistent calculation being run to generate
the projected density of states.
First, we note that the value of the energy is:
Etotal=-4.92489557316370E+01 which shows that we have attained presumably the same state as above.

The density of states will be in the files tspin_3o_DS2_DOS_AT0001 for the first atom, and tspin_3o_DS3_DOS_AT0002 for the second atom.
We can extract the density of d states, which carries most of the magnetic moment and whose integral up to the Fermi level will yield an estimate of the magnetization on each atom.
We note the Fermi level (recalled in the file tspin_3o_DS1_DOS)
 Fermi energy :      -0.27799778
If we have a look at the integrated site projected density of states, we can compute the total moment on each atom. To this end, one can edit the file tspin_3o_DS3_DOS_AT0002, which contains information pertaining to atom 2. This file is self documented, and describes on each line, for spin up and spin down:

 index  energy(Ha)  l=0      l=1      l=2      l=3      l=4   (integral=>) l=0    l=1    l=2    l=3    l=4

If we look for the lines containing  an energy of "-0.27800",  we find
6440   -0.27800   0.2070   1.3021  47.5017   0.3928   0.1723   0.30    0.34    3.22    0.04    0.01 (up)
6440   -0.27800   0.4656   1.8227  26.0414   0.5495   0.1748   0.30    0.35    3.54    0.04    0.01 (dn)


There are apparently changes in the densities of states for all the channels, but besides the d-channels, these are indeed fluctuations. This is confirmed by looking at the integrated density of states which is different only for the d-channel. The difference between up and down is .32, in good agreement (regarding our very rough methods of integration) with the preceeding calculation. Using a calculation with the same number of k points for the projected DOS, we can plot the up-down integrated dos difference for the d-channel. Note that there is some scatter in this graph, due to the finite number of digits (2 decimal places) of the integrated dos given in the file tspin_3o_DS3_DOS_AT0002.

If we now look at the up and down DOS for each atom, we can see that the corner atom and the face atom possess opposite magnetizations, which exactly cancel each other. The density of states computed with the tetrahedron method is not as smooth as by the smearing method, and a running average allows for a better view.








4. A ferrimagnetic example (Not yet)

Some materials can display a particular form of ferromagnetism, which also can be viewed as non compensated antiferromagnetism, called ferrimagnetism.
Some atoms possess up spin and other possess down spin, but the total spin magnetization is non zero.
This happens generally for system with different type of atoms, and sometimes in rather complicated structures.
Here, we will perform  a calculation on a "simple" structure .




5. The spin-orbit coupling

 For heavy atoms a relativistic description of the electronic structure becomes necessary, and this can be accomplished through the relativistic LDA approximation. For atoms, the Dirac equation is solved and the 2(2l+1) l-channel degeneracy is lifted according to the eigenvalues of the L+S operator (l+1/2 and l-1/2 of degeneracy 2l+2 and 2l). After pseudization, the  associated wave functions can be recovered by adding to usual pseudo-potential projectors a spin-orbit term of the generic form v(r)|l,s>L.S<l,s| . Not all potentials include this additional term, but it turns out that it is the case for all the HGH type pseudopotentials. In a plane wave calculation, the wavefunctions will be two components spinors, that is will have a spin-up and a spin-down component but these components will be coupled. This mean the size of the Hamiltonian matrix is quadrupled. We will consider here a heavier atom than Iron: Tantalum. You will have to change the "files" file accordingly, as we here have used the potential: 73ta.hghsc . It is a HGH pseudopotential, with semicore states. Replace the last line of the tspin_x.files by
../../Psps_for_tests/73ta.hghsc
You can copy the file ~abinit/tests/tutorial/Input/tspin_5.in in Work_spin. Change accordingly the file names in tspin_x.files, then run the calculation. It takes about 90 secs on a 2.8 GHz PC machine.

The input file contains new variables:
Have a look at these. In this run, we check that we recover the splitting of the atomic level by performing a calculation in a big box. Two calculations are launched with and without spin-orbit.
We can easily follow the symmetry of the different levels of the non spin orbit calculation:
  kpt#   1, nband= 26, wtk=  1.00000, kpt=  0.0000  0.0000  0.0000 (reduced coord)
-2.44772 
-1.46448  -1.46448  -1.46448 
-0.17051 
-0.10859  -0.10859  -0.10859  -0.10746  -0.10746

That is, the symmetry: s,p,s,d
After application of the spin-orbit coupling, we now have to consider twice as many levels:
 kpt#   1, nband= 26, wtk=  1.00000, kpt=  0.0000  0.0000  0.0000 (reduced coord)
-2.43309  -2.43309 
-1.67343  -1.67343  -1.35515  -1.35515  -1.35515  -1.35515

-0.16816  -0.16816 
-0.11662  -0.11662  -0.11662  -0.11662  -0.09254  -0.09254
-0.09152  -0.09152  -0.09152  -0.09152

The levels are not perfectly degenerate, due to the finite size of the simulation box. We can nevetheless compute the splitting of the levels, and we obtain, for e.g. the p-channel: 1.67343-1.35515=0.31828 Ha
If we now consider the NIST table of atomic data, we obtain :

5p splitting, table: 1.681344-1.359740=0.321604 Ha
5d splitting, table:   .153395-.131684=0.021711 Ha

We obtain a reasonable agreement.
A more converged (and more expensive calculation) would yield:

5p splitting, abinit: 1.64582-1.32141=0.32441 Ha
5d splitting, abinit:   .09084-.11180=0.02096 Ha








6. Rotation of the magnetization and Spin-orbit (not yet)

The most spectacular manifestation of the spin-orbit coupling is the energy associated with a rotation of the magntisation with respect with the cristal axis. It is at the origin of the magneto cristalline anisotropy of paramount technological importance.


This ABINIT tutorial is now finished...
GZ would like to thank B. Siberchicot for useful comments.


Goto : ABINIT home Page | Suggested acknowledgments | List of input variables | Tutorial home page | Bibliography
Help files : New user's guide | Abinis (main) | Abinis (respfn) | Mrgddb | Anaddb | AIM (Bader) | Cut3D | Optic | Mrgscr