Spin polarized calculations

Contents

1 Useful info
2 Kohn-Sham wavefunctions of a CO molecule
3 Vibrational frequencies of the H₂O molecule
4 Electron spin: Fe(bcc)
- 4.1 Density of states (DOS)
- 4.2 Projected Density of states (PDOS)

1 Useful info

There is an online ASE manual on the web.
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.

Save the relaxation in a trajectory and look at it:
```
plottrajectory -p "[[Distance(0,1)]]" CO-traj.nc
```
Plot the Kohn-Sham wavefunctions of the different wavefunctions of the CO molecule. For this, copy the file plotwavefunction.py to your area. (Alternatively, you can plot the wavefunctions with VTKplotwavefunction.py;

usage: python -i VTKplotwavefunction.py [bandnumber].) 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₂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₂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₂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^-1for the bending mode, 3657cm^-1for the symmetric stretching mode and 3756cm^-1for 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 VMD, because you can play a movie. For this, you have to convert the trajectory files to .xyz files, which is the type VMD can read. You do this with the command trajectory2xyz [traj filename] [xyz filename]. Open VMD by typing the command vmd filename.xyz. 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⁶4s², 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:
1. A non-magnetic calculation
2. A ferro-magnetic calculation (aligned atomic moments)
3. 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.
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.IO.Cube import WriteCube
>>> atoms = Dacapo.ReadAtoms('anti.nc')
>>> calc = atoms.GetCalculator()
>>> spin0 = calc.GetDensityArray(spin=0)
>>> spin1 = calc.GetDensityArray(spin=1)
>>> zeta = (spin0 - spin1) / (spin0 + spin1)
>>> WriteCube(atoms,zeta,'magnetization.cube')

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.

4.2 Projected Density of states (PDOS)

The projected density of states is usefull for for analyzing chemical bonding. There exist several studies where the density projected onto the d states of a given surface atom is used. This short excercise demonstrates how to construct the PDOS of Fe. Take a look at the manual on local density of states. Take notice of the speciel numbering of atoms which start from 1.

We will get a feel for the local density of states by plotting the PDOS for the non magnetic Fe crystal. Look at LDOS.py. Use it to plot the s,p,d states on one of the Fe atoms.

PDOS can also be used to investigate the magnetic moment on the individual atoms. Look at localmagmoment.py script and use it to determine the magnetic moment on each Fe atom in the ferro and anti ferro magnetic case. Comment on the dependence of the magnetic moments on the cutoff radius.