1. Force relaxations and thin films
NOTE: For film calculations the inpgen k point set generation seems to be broken at the moment. Therefore replace the default k point set for the film calculations by
<kPointMesh nx="7" ny="7" nz="1" gamma="F"/>
and for the analogous bulk calculation by
<kPointMesh nx="7" ny="7" nz="3" gamma="F"/>
Note 2: The k point meshes defined here do not yield a converged calculation. In fact they are very coarse. Depending on the other numerical parameters the coarseness may lead to an error in the determination of the Fermi energy. If this is the case refine the k point mesh.
1.1. Ag film in (100) direction
In this tutorial we will set up a thin Ag film and calculate its band structure, density of states, and work function. As a preparation we first perform a calculation of the band structure and density of states for a bulk Ag fcc crystal.
Consider the conventional fcc unit cell. We cut out of this unit cell a smaller unit cell consisting of 2 atoms. In z direction the smaller, tetragonal, unit cell should feature the same length as the conventional unit cell. In the inpgen input the tetragonal unit cell is denoted as 'tP'. For the Ag fcc lattice constant we use a.
Reduce the Ag MT radius to 95% of its default value. Perform the self-consistent density calculation for this unit cell and calculate its band structure and density of states.
Now we calculate three thin Ag films with different thickness in (100) direction: 1,3, and 9 atomic layers. The setup of the respective inpgen inputs is as follows:
After the comment line insert the line
&input film=t /
to indicate that we want to set up a unit cell for a thin film. In the next line we have to provide the lattice definition. As above this is a 'tP' cell. The c parameter has to be large enough that all atoms fit in the unit cell. Besides this there are no further requirements to it. It will automatically be reduced to a reasonable value. Set it to . The a parameter should be the same one as in the bulk unit cell you set up.
The list of atoms should be arranged such that the film is centered around the (x,y,0) plane, i.e., there is a central atomic layer and every atomic layer above this central layer also has a twin below the central layer. In general, of course, there does not have to be z reflection symmetry, but for the case of the Ag crystal there is.
Generate inp.xml inputs for the 1, 3, and 9 atomic layer films. Reduce the Ag MT radius to 95% of its default value. Perform self-consistent density calculations for each of the three setups, then calculate for each setup the band structure and the density of states.
How does the band structure change when going from 1 to 3 to 9 layers? Why? Compare the three band structures to the band structure of the bulk unit cell. How does the density of states for the atom of the central layer change when going from 1 to 3 to 9 layers? Compare it to the density of states of the Ag atom in the bulk unit cell.
We next relax the atom positions in the 9 layer film such that the forces on the atoms vanish. Copy the inputs for the 9 layer film in a directory 'relax' and modify it such that for each atom group the atomGroups/atomGroup/force/@relaxXYZ XML attribute is set to "FFT". This indicates that we only want to perform atom position relaxations in z direction. On top of the self-consistent density we now set calculationSetup/geometryOptimization/@l_f to "T" to activate the calculation of all force contribution together with a structural relaxation step. Perform the relaxation step by calling fleur. Your inp.xml file will then be replaced by a new one with changed atom positions. Also a file forces.dat will be created. This file is readable but is also used as an input for further relaxation steps. It keeps the history of total energies, atom positions, and forces. With each step it gets extended with the additional data.
After the first relaxation step the l_f switch is set to "F" again because we now have to perform a new self-consistent density calculation. Just call fleur to start that calculation, then perform a relaxation step again and so on. Stop after you have the impression that the atom positions are in equilibrium. In what directions have the atoms been displaced and by what distance?
The work function of silver actually is not very sensitive to this relaxation of the atom positions. This may be different for other materials and other quantities. We obtain the work function as the difference of the vacuum potential at an infinite distance from the film (grep vzInf out.xml) and the calculated fermi energy. In our case vzInf should be set to . This may vary depending on the setup. How large is the calculated work function in eV? Compare it to the experimental value for a silver surface in (100) direction: 4.64 eV. Note that the highest occupied state of the Kohn-Sham system actually is the only state of that system that has a strict physical meaning.
2. Exercises
2.1. PbTiO3
Lead titanate (PbTiO3) is a ferroelectric material that forms a perovskite structure. Materials in perovskite structure have the chemical formular ABX3. In its basic form perovskites are related to a cubic unit cell with "A" atoms at the corners, "B" atoms in the center, and "X" atoms at the face centers. Look up a visualization of perovskite unit cells. Note that these visualizations are sometimes shifted such that the "B" atom is at a corner position of the unit cell.
As a ferroelectric material PbTiO3 features a spontaneous electric polarization. This causes a slight deviation of the actual unit cell from the basic perovskite structure. In this exercise we want to determine this unit cell under the (wrong) assumption that it is still cubic. We start with the following inpgen input:
PbTiO3
&lattice latsys='sc' a=7.50 /
5
82 0.00 0.00 0.000
22 0.50 0.50 0.475
8 0.50 0.50 0.000
8 0.50 0.00 0.500
8 0.00 0.50 0.500
&end /
As you can see the Titanium atom in this input is slightly displaced from its center position. This is done to prevent inpgen from detecting certain symmetries in the unit cell that would become a problem when we displace atoms in it in z direction.
After generating the inp.xml file place the Ti atom in it in the center of the unit cell again. To allow for a wider range of atom positions reduce the MT radii of Pb and Ti to 90% and those of O to 95%.
We now perform several self-consistent density calculations for different displacements of the Ti atom in z direction. Vary the Ti position from 1.00/2.00 to 1.05/2.00 in steps of 0.01/2.00 (or finer if you want to have a nicer result, e.g., 0.005/2.00). For each of the displacements calculate the forces on the Ti atom in z direction (1 step with l_f="T"). Plot the total energy vs. the displacements on top of a plot of the forces vs. the displacements. If everything works as expected you should see a total energy minimum for a finite but small displacement of the Ti atom.
Starting from a position of the Ti atom near the energy minimum perform a force relaxation for the Titanium and Oxygen atoms in z direction (keep Pb fixed at the predefined position). If you encounter any errors in this process follow the hints for the different errors below. Where are the atom positions for the relaxed unit cell? What total energy does the unit cell feature? Compare the displacements and total energy to the basic perovskite structure and the equilibrium position in the displacement of the Ti atom alone.
Note: In the test calculations to this exercise about 20 force relaxation steps were needed to reach the "official" message that the structure is relaxed with respect to the forces:
*****************************************
Run finished successfully
Stop message:
GEO Des woars
*****************************************
This message is reached when in one step of the calculation the inp.xml file values calculationSetup/geometryOptimization/@epsdisp for the atom displacement and calculationSetup/geometryOptimization/@epsforce for the forces on the atoms are underrun. (They were kept at the default values for the test calculation.) The values for the forces are always provided in .
3. Last weeks results
4. Notable error messages in the context of force relaxations
Note: When performing force relaxations one often encounters error messages. These can be due to an erroneous setup of the unit cell or wrong assumptions about the final atom positions or the total energy landscape. The relevant error messages are:
*****************juDFT-Error*****************
***Error message:bfgs: <p|y><0
***Error occurred in subroutine:bfgs
********************************************
This error message is due to the assumption of the underlying optimization algorithm that we are already very near to the relaxed position and the total energy rises quadratically with respect to the displacement of the atom positions from the equilibrium configuration. If this assumption is too broken, e.g., the curvature of the energy landscape is negative, one obtains this error. If one is sure that the setup is not too bad the solution to this error is to remove the forces.dat file and perform a simple linear relaxation step, like the very first relaxation step. After a few linear steps this error hopefully does not appear any more.
*****************juDFT-Error*****************
***Error message:Error checking M.T. radii
***Error occurred in subroutine:chkmt
********************************************
This error appears if due to a relaxation step atom positions get too near to each other such that the respective MT spheres overlap. In the best case this is just due to an underestimation of the required relaxation. In such a situation one may delete the forces.dat file, reduce the MT radii, and continue.
*****************juDFT-Error*****************
***Error message:Potential shift at maximum
********************************************
This error is uncommon for force relaxation (or in general) but in the test calculations to the exercises it appeared a few times. It was solved by removing the cdn.hdf file and then restarting the respective fleur run.