# Spin polarized calculations

Contents

### 1 Useful info

- There is an online ASE manual on the web: http://www.fysik.dtu.dk/campos/ASE.
- There is an online Dacapo manual on the web: http://www.fysik.dtu.dk/campos/Dacapo.
- Want to learn more about the Python language? This tutorial takes you around the language: http://www.python.org/doc/current/tut/tut.html.

### 2 Kohn-Sham wavefunctions of a CO molecule

In this section, we will look at the Kohn-Sham wavefunctions of the CO molecule and compare them to results from molecular orbital theory.

- Make a script, where a CO molecule is placed in a qubic unitcell,
e.g. of 5 Å in each direction. For more accurate calculations, the
cell should definitely be bigger, but for reasons of speed, we use
this cell here. Which values for the plane-wave and density cutoffs
would you use? How many
**k**-points should you use? As this is a calculation of a molecule, one should also use a low electronic temperature. What value would you use? If`calc`is your calculator, the electronic temperature is set by`calc.SetElectronicTemperature(temp)`, where`temp`is the electronic temperature in eV units. Guess reasonable positions from the covalent radii of C and O. Then relax the CO molecule to its minimum energy position. - Plot the Kohn-Sham wavefunctions of the different wavefunctions of the CO
molecule. For this, copy the file
CO_wavefunction.py to your area.

What is the highest occupied state and the lowest unoccupied state? Which wavefunctions are bonding and which are antibonding?

### 3 Vibrational frequencies of the H_{2}O molecule

Density functional theory can be used to calculate vibrational frequencies of molecules, e.g. either in the gas phase or on a surface. These results can be compared to experimental output, e.g. from IR-spectroscopy, and they can be used to figure out, how a molecule is bound to the surface. In this example, we will calculate the vibrational frequencies for a water molecule.

- For a simple molecule, like CO, there is only one stretching mode. How would you calculate the vibrational frequency of this mode?
- For a general molecule with N atoms, how many modes are there? How many of them are vibrational modes? How would you do a calculation for the vibrational modes? Describe in detail, which steps have to be performed.
- Make a script, where a H
_{2}O molecule is placed in a cubic unitcell, e.g. of 6 Å in each direction. A good guess for the beginning is to put the O atom in (0,0,0), one H atom in (-0.77, 0, -0.6) and the other H atom in ( 0.77, 0, -0.6). Relax the H_{2}O molecule to the equilibrium positions in the same way as for the CO molecule. - Copy the file
H2O_vib.py to your area and try to understand what it does.

- Run the script and look at the output frequencies. Compare them to literature values,
which are 1595cm
^{-1}for the bending mode, 3657cm^{-1}for the symmetric stretching mode and 3756cm^{-1}for the anti-symmetric stretching mode. How good is the accuracy and what are possible error sources? - Now we want to look at the modes to see how the atoms move. For this we use the
files
`H2O_vib_mode_x.traj`where`x`is the number of the mode counted in the order they are printed out. You can look at these trajectories with the`plottrajectory`command. A nicer way to visualize them is using`xmakemol`, because you can play a movie. For this, you have to convert the trajectory files to .xyz files, which is the type`xmakemol`can read. You do this with the command`trajectory2xyz [traj filename] [xyz filename]`. Open xmakemol by typing the command`xmakemol`, open a .xyz file (Go to`File`and then`Open...`), and go to`Control`and then`Frames...`. Now you can play the movie and look at all the modes. Do they look like you expected and what would you have expected (you may have learned something about symmetry groups at one point)?

### 4 Electron spin: Fe(bcc)

As an example of spin polarized calculations, we'll study Fe(bcc) in a
two-atom unit cell, i.e. a simple cubic Bravais lattice with a basis,
where one Fe atom sits at the origin, the other in the middle of the
unit cell. We'll stick to the experimental lattice constant *a* = 2.87
Å. The atomic term of iron is [Ar]3d^{6}4s^{2}, i.e. 8
valence electrons/atom is included in the calculation.

Use Hunds rule (maximum polarization rule) to calculate the magnetic moment of an isolated Fe atom. Draw schematically the one-electron eigenvalues for spin up and down on an energy axis, along with electron populations.

We'll make three calculations for bulk iron:

A non-magnetic calculation

A ferro-magnetic calculation (aligned atomic moments)

A anti ferro-magnetic calculation (anti parallel atomic moments).

Calculation 3 is only possible, if we have at least two atoms in the unit cell.

Copy the following script

Fe-ferro.py to your area.

- Now fill in the missing slots (marked with ?); the
following is useful to notice:
- How many bands are needed? Assuming the atoms polarize
maximally (as the isolated atoms). For metals, always have at least
5 extra bands to allow for uneven filling of states for different
**k**-points. - One may
*help*a magnetic calculation by providing an initial magnetic moment on an atom like`Atom('Fe', ..., magmom=?)`This option is necessary to find magnetic states. Choose the magnetic moment close to the expected/desired magnetic state of your system (the expeimental value is 2.22 per atom). The initial magnetic moment is relaxed during the self consistency cycles.

- How many bands are needed? Assuming the atoms polarize
maximally (as the isolated atoms). For metals, always have at least
5 extra bands to allow for uneven filling of states for different
- For each of the three magnetic phases ferro, antiferro and nonmagnetic, write down sensible guesses for initial magnetic moment parameters: magnetic moment for each of the two atoms in the unit cell.
- Run the
`Fe-ferro.py`calculation, and write two scripts for the other spin configuartions. If you run the calculation in parallel (`parrun -n N file.py`), where`N`is the number of processors, it will be most efficient to use a value for`N`that is a divisor in the the product*KS*, where*K*is the irreducible number of**k**-points in your calculation (i.e. number of**k**-points after symmetrization has been applied) and*S*the number of electron spins (1 or 2). Why may this be a good idea?

- Tip::
- To find the irreducible number of
**k**-points*K*, start your script as a serial calculation and look for the lines starting with`KPT`in the text output file. Then stop the serial calculation, select`N`and start it in parallel.

Compare the energies of the three magnetic phases.

Experimentally, the ferromagnetic phase is most stable. Is this reproduced for LDA and GGA?

Compare the calculated magnetic moment for the ferromagnetic phase with the experimental value. In the text output file, look for the last line starting with

`MOM`:MOM total_moment (spin_up_moment spin_down_moment )

(unit = Bohr magnetons per unit cell)

#### 4.1 Density of states (DOS)

Take a look at the

dos.py program (written in Python) and try to get
a rough idea of what it can do for you. Use it to plot DOS for the
three Fe configurations (on the *x*-axis you have the energy relative to
the Fermilevel).

Do the DOS plots integrate to the correct numbers? (i.e. number of valence electrons).

The DOS for the anti-ferromagnetic phase looks a bit like that for the non-magnetic phase - is it magnetic at all?! Try to visualize the magnetization like this:

>>> from Dacapo import Dacapo >>> from ASE.Visualization.VTK import VTKPlotArray >>> atoms = Dacapo.ReadAtoms('anti.nc') >>> calc = atoms.GetCalculator() >>> spin0 = calc.GetDensityArray(spin=0) >>> spin1 = calc.GetDensityArray(spin=1) >>> zeta = (spin0 - spin1) / (spin0 + spin1) >>> plot = VTKPlotArray(zeta, atoms.GetUnitCell())

Calculate the DOS for bulk Aluminum and compare it (qualitatively) to the DOS for the non-magnetic calculation. The DOS for a simple metal has this shape:

*g*(*E*) ~*E*^{1/2}. Explain the qualitative difference.Plot also the DOS for bulk Si and the CO molecule. Identify the bandgap between valence and conduction bands for Si and the HOMO-LUMO gap for CO.