Spin polarized calculations

1   Useful info

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 H2O 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 H2O 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 H2O 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.

System Message: WARNING/2 (<string>, line 86)

Bullet list ends without a blank line; unexpected unindent.

System Message: WARNING/2 (<string>, line 87)

Explicit markup ends without a blank line; unexpected unindent.
  • 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 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]3d64s2, 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.

System Message: WARNING/2 (<string>, line 132)

Bullet list ends without a blank line; unexpected unindent.
  • 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.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) ~ E1/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.