Performing structural relaxations
The geometry of crystal unit cells features many degrees of freedom. The most basic of these are the Bravais lattice vectors, which can be described by the lattice type together with the lattice parameters. Next the atom positions in the unit cell also represent geometrical degrees of freedom and it may even be that the groundstate structure of the crystal is only describable in a supercell because the atom positions in neighboring unit cells interact.
The Fleur code provides a feature to semi-automatically optimize the atom positions for a given unit cell shape by minimizing the forces acting on the atoms. The other degrees of freedom have to be investigated by the user explicitely if they are of interest. Especially the lattice parameters have to be optimized manually by performing a series of calculations with test lattice parameters and extracting from these the optimal parameters minimizing the total energy. Note that all geometrical degrees of freedom may interact with each other. Also the magnetic ground state configuration may affect the optimal geometrical parameters.
The total energy for each iteration of an SCF calculation can be extracted from the out.xml file by
invoking grep "<totalEnergy" out.xml. A listing of the individual energy contributions is also
provided within the totalEnergy section.
The workflow to optimize the atom positions essentially works such that Fleur is invoked several times and each time an SCF calculation with a subsequent calculation of the forces and the generation of suggested atom displacements is performed. In a preprocessing step to this workflow and other structural relaxations the Fleur input file should be adapted such that the MT spheres are slightly reduced in size to allow a movement of the atoms without obtaining overlapping spheres.
For each atom species the MT radius is specified in atomSpecies/species/mtSphere/@radius. Default radii are automatically
determined by the Fleur input generator and may subsequently by modified by the user. Since the default parameters lead to nearly
overlapping spheres it is a good strategy to reduce them to about before performing structural relaxations. The input file
understands mathematical expressions so that the user does not have to perform the involved calculations by manually. An example
for the reduction of a MT radius to is provided below.
<mtSphere radius="0.95*2.2200" gridPoints="733" logIncrement=".0160"/>
The reduction of the MT radii should be performed for the optimization of the atom positions as well as for manual optimizations of the lalltice parameters. Note that total energies from calculations with different MT radii are typically not comparable since the quantity of interest are energy differences between the calculations and not total energies of individual calculations. Total energy differences are typically converged with smaller basis set sizes and thus the total energies themselves are typically underconverged. If it turns out that the reduction of the MT radii was not significant enough and one still obtains overlapping MT spheres in a certain calculation probably all calculations have to be performed again with smaller MT radii. It is therefore advisable to start with a geometry of the unit cell that is already near the expected optimum.
The storage of the charge density also depends on the MT radii. If the radii or other geometrical parameters are changed it is therefore important to also delete the charge density file if there already is one in the directory.
The next step for the optimization of the atom positions is to select the atoms for which the forces are to be calculated and the directions in which the force relaxations may be performed.
In the Fleur input file the selection of atoms and directions for force relaxations is performed for each group of symmetry
equivalent atoms separately. The switch to calculate forces for the atoms of a certain group is found in atomGroups/atomGroup/force/@calculate.
Note that even if this switch is set to T all force contributions are only calculated if a structural optimization of
the atom positions is activated. The directions for which the relaxations are to be performed are selected
in atomGroups/atomGroup/force/@relaxXYZ by three logical values for the x, y, and z direction. An example force tag
used for the relaxation of the z coordinate only is presented below. Such a relaxation is often used when calculations
for thin film systems have to be performed.
<force calculate="T" relaxXYZ="FFT"/>
Different iterative optimization schemes can be used to find the atom positions featuring vanishing forces. In each force iteration they generate a new displacement of the atom positions from the original input file based on the position-dependent forces calculated up to that point. A very efficient optimization scheme for this purpose is the BFGS algorithm. But also a straight linear mixing may be useful in certain situations. These optimization schemes are parametrized by a mixing parameter and convergence criteria relating to the maximal change of the positions between the last two force iterations and the maximal forces.
The force optimization algorithm is selected by the calculationSetup/geometryOptimization/@forcemix parameter. The
default BFGS entry selects the BFGS algorithm, straight selects a straight linear mixing. The mixing parameter
is set in calculationSetup/geometryOptimization/@forcealpha, the convergence criterion for the change of the atom
displacements in calculationSetup/geometryOptimization/@epsdisp, and the convergence criterion for the
forces in calculationSetup/geometryOptimization/@epsforce.
After selecting and parametrizing the force optimization scheme. The only missing step is the activation of the force optimization procedure.
The force optimization procedure is activated by the calculationSetup/geometryOptimization/@l_f switch. An example for
the geometryOptimization tag with activated force optimization is provided below.
<geometryOptimization l_f="T" forcealpha="1.00" forcemix="BFGS"
epsdisp=".00001" epsforce=".00001"/>
After activating the force optimization starting Fleur yields an SCF calculation. When this calculation reaches the convergence criterion for the minimal distance between the charge densities the forces are calculated and a new atom displacement is calculated.
Fleur stores the atom displacements relative to the atom positions provided in the Fleur input file in the file relax.xml.
By default this file is included into the Fleur input file in terms of an XML include tag. If no such file is present in the
working directory or the inlude tag has been removed no displacements of the atom positions are taken into account. The file
includes a list of diplacements for the representative atoms of each atom group in relaxation/displacements/displace. The
order of displacements has to be consistent with the order of the respective atom groups. For every force relaxation step the
file also contains a relaxation/relaxation-history/step section. This section provides the total energy related to
the SCF run for the respective atom positions in relaxation/relaxation-history/step/@energy and also for each atom group
a relaxation/relaxation-history/step/posforce element containing the x, y, z coordinates of the representative atom followed
by the forces on this atom in x, y, and z direction. The relaxation/relaxation-history is considered for the mixing part
of the chosen force optimization algorithm. After performing the force optimization the final atom positions can also easily be
extracted from this section.
Performing a force optimization of the atom positions for a given unit cell shape therefore mainly reduces to invoking Fleur several times with activated geometry optimization. In each force calculation step the atom displacements as well as the relaxation history are updated. They are automatically included in the following Fleur run.
If Fleur is started with an old charge density file in the working directory the starting density is taken from this file, even if the atom positions have changed. This may be a very good starting point if the change of the atom positions is small. In such a situation the number of required iterations to reach a self-consistent density is strongly reduced. If the shifts in the atom positions is large this approach can be a catastrophy. It may lead to a slower convergence or even to an error in the calculation because the potential obtained from the stored density is too inconsistent to the actual atom positions. If a problem like this is observed it may be a good idea to delete the charge density file and automatically construct a new starting density for the current atom positions.
It can also be that the relaxation-history gives rise to large atom displacements. This is because the model assumption in optimization algorithms typically is that the parameters of the function to be optimized are near an optimum and that the function can be approximated by a parabolic model near this optimum. If this is not fulfilled the hints extracted from the relaxation-history may be contradictory ending in an extreme shift of the next atom position. In such a situation the advanced user can delete parts of the relaxation-history by hand or switch to a straight linear optimization scheme until the parabolic model assumptions are better fulfilled.
Besides these challenges during the force optimization the choice of the initial atom positions is also crucial. As already mentioned it is of benefit to start from atom positions that are already near their optimum. Beyond this the user should pay attention that the starting atom positions do not preserve symmetries that are not featured by the possible optimized positions. Otherwise one may end up with vanishing forces because the initial positions mark a local maximum in the energy landscape, not because they mark a minimum.
Note that forces are defined as the negative variation of the total energy with respect to the atom positions. This implies that the demands for high-quality force calculations are stricter than for the simple calculation of total energy differences. Especially due to the atom position dependence of the LAPW basis and the partitioning of the unit cell care has to be taken to ensure that the forces are carefully converged with respect to the calculation setup. This implies that the need to describe semicore states by LOs as valence electrons is increased. Also the demands for the cutoff parameters of the LAPW basis may be increased.