Simple stuctural relaxation
Objectives
- Calculate forces in FLEUR.
- Use forces to find equilibrium positions of atoms.
- Understand important parameters for force calculations.
Introduction
The atomic stucture you provide to FLEUR might not be the ground state in several aspects. In particular, the atoms might not be in their relaxed positions, i.e. there might be forces acting on the atoms. Here we will now discuss how to calculate these forces and how to use them to find the relaxed positions.
Note
The structure you obtain here is not necessarily the atomic ground state. You might also have to change the lattice constant/lattice structure, the magnetic state or other parameters external to the discussion here.
Info
To obtain a relaxed structure a series of FLEUR calculations needs to be run. The general workflow however is simple: 1. Perform an SCF calculation for your current structure. 2. If the structure is converged, calculate the forces on the atoms. 3. Modify the structure in a way to reduce the forces. 4. These steps 1.-3. can be done in a single call to FLEUR. To eliminate the forces and to find the relaxed structure, you have to restart FLEUR several times until you consider the forces in step 3. to be small enough.
## FLEUR steps
***This tutorial uses the film structure of the previous tutorial(F5)***
```bash
cd MoS2
We start by using the inpgen to create a new setup.
inpgen -f inpMoS2.txt
Modification of inp.xml
We now modify the inp.xml file to perform the force calculations and generate the
atom displacements. The most fundamental switch, which enables the evaluation of forces and the calculation of a new structure is /fleurInput/calculationSetup/geometryOptimization/@l_f which we switch to "T".
In the geometryOptimization XML element you can also change the relaxation algorithm to something different (we keep
BFGS), change the mixing
parameter, or adjust the criteria for a successful relaxation. Note that you can limit the force calculations for each atom type to certain orientations. This is done in the /fleurInput/atomGroups/atomGroup/force/@relaxXYZ XML attribute. Each logical in this variable stands for one direction.
In this example we want to keep the Molybdenum atom fixed (having one fix point in a structural relaxation is always a good idea because otherwise and depending on the accuracy of the force calculation one may obtain a small overall drift term for the whole unit cell). Set relaxXYZ for Molybdenum to "FFF". For the Sulphur atom, we only consider a relaxation in z-direction (even though in this example the forces in the other directions are anyway). Set relaxXYZ for Sulphur to "FFT".
Multiply the MT radii for the atoms by to allow atom movements without overlapping spheres. The radii can be set in /atomSpecies/species/mtSphere/@radius. This is a general adjustment advisable for force relaxations. For this example this additional buffer between the MT spheres is actually needed to avoid that the atoms leave the interstitial volume and touch the vacuum boundary, but in other cases atoms might also touch each other. Note that Fleur understands simple mathematical expressions for many parameters in its input file. You can actually write something like (there are no spaces within such an expression allowed):
<mtSphere radius="0.95*2.57000000" gridPoints="849" logIncrement=".01500000"/>
Please check that you change the radius of both atoms!
```{attention}
Open an editor and modify the ```inp.xml``` file as sketched above.
Invoke Fleur to start calculating for this configuration.
fleur_MPI
Please note that in this step actually no new structure was calculated as the forces did not yet converge. Hence, restart the calculation
fleur_MPI
Output of a structural relaxation run
In this call of FLEUR the forces now have been calculated. However, a new structure is only generated if the forces are sufficiently converged. This implies:
- at least two iterations will be performed to check the convergence of forces.
- only if the difference between the forces is small FLEUR will create new positions and stop. The usual criteria for checking the convergence are therefore ignored. (However, the
itmaxparameter will still be honoured, so make sure it is large enough)
If new positions are calculated, you will
- get the message
Structural relaxation: new displacements generatedinstead of the usualall done. - find a file
relax.xml.
In every call of FLEUR the relax.xml file is automatically included in the inp.xml file whenever it is present. Investigate the file. You will first find a list of displacements for each atom group. This is used to modify the atom positions in the inp.xml file. The displacements are updated by every new force calculation. The other big block in the relax.xml file is the relaxation history. Here you find for every relaxation step and atom type a list of atom positions in cartesian coordinates followed by the forces on the representative atom of this atom type in . You can use this list as an output, e.g., to look up the final atom positions in it, but it is also an input for Fleur that is used by the relaxation algorithm to determine the next atom positions. An experienced user may modify this file whenever a relaxation process seems to be problematic. For example, in certain situations it may be a good idea to delete a part of the relaxation history: The BFGS algorithm constructs a quadratic model of the energy landscape around the equilibrium position. If the initial guess of the atom positions is too far away from the equilibrium positions, such a model might not be a good approximation, and then it can be advisable to delete parts of the relaxation history.
We now have to repeat the self-consistency calculation followed by the force step several times to reach the equilibrium atom positions. You can just invoke Fleur (keeping l_f="T") several times to do so. Have a look at the changes in the relax.xml file while you do so. The relaxation of the atom positions is finished once FLEUR ends with the message Structural relaxation: Done.
Repeat this step until Fleur stops with the terminal output
```text
*****************************************
Run finished successfully
Stop message:
Structural relaxation: Done
*****************************************
#Repeat invoking FLEUR until the relaxation is finished
fleur_MPI
Note that a systematic approach to an automatization of this workflow is possible by using AiiDA.
Now the forces on the atoms are smaller than the respective convergence criterion and the atoms can be considered to having reached their equilibrium positions.
The relevant convergence criteria for this step are given in the inp.xml file in the geometryOptimization tag:
- epsdisp: the minimal displacement
- epsforce: the minimal force
at which we consider a structure to be converged.
As soon as this is the case find the Sulfur z-positions of the last force step in the relax.xml file and compare them to the original positions. How much did they change? What is the energy difference between the two configurations?
cat relax.xml
Your relax.xml file should look similar to the following output from the reference calculation (modified by removing some irrelevent digits):
<!-- Attention, absolute coordinates used here -->
<relaxation>
<displacements>
<displace> 0.0000000000 -0.0000000001 0.0000000000 </displace>
<displace> 0.0000000000 0.0000000001 0.1060500828 </displace>
</displacements>
<relaxation-history>
<step energy=" -4848.0718804597">
<posforce> 0.0 -3.4367564958 0.0000000000 0.0 0.0 0.0000000000 </posforce>
<posforce> 0.0 3.4367564958 2.8500000000 0.0 0.0 0.0232689368 </posforce>
</step>
<step energy=" -4848.0727664846">
<posforce> 0.0 -3.4367564958 0.0000000000 0.0 0.0 0.0000000000 </posforce>
<posforce> 0.0 3.4367564958 2.8732689368 0.0 0.0 0.0178807525 </posforce>
</step>
<step energy=" -4848.0739479634">
<posforce> 0.0 -3.4367564959 0.0000000000 0.0 0.0 0.0000000000 </posforce>
<posforce> 0.0 3.4367564959 2.9504871752 0.0 0.0 0.0011397651 </posforce>
</step>
<step energy=" -4848.0739385584">
<posforce> 0.0 -3.4367564959 0.0000000000 0.0 0.0 0.0000000000 </posforce>
<posforce> 0.0 3.4367564959 2.9557443713 0.0 0.0 0.0000623203 </posforce>
</step>
<step energy=" -4848.0739375964">
<posforce> 0.0 -3.4367564959 0.0000000000 0.0 0.0 0.0000000000 </posforce>
<posforce> 0.0 3.4367564959 2.9560484518 0.0 0.0 0.0000003325 </posforce>
</step>
</relaxation-history>
</relaxation>
We can observe in this list of energies and forces that while the forces nearly vanish for the last relaxation step, the total energy becomes minimal for the second last step. Even though this is a very small discrepancy, one should be aware that something like this can happen due to several reasons, e.g., insufficient convergence with respect to parameters or missing force terms.
Note
The cdn.hdf file has become quite large because many densities are stored in it. You can investigate it by calling Fleur with the -info command line switch. You will obtain a list of the densities together with time stamps and the distance to the previous density. Of course, there are ways to get rid of certain densities. One way is to perform the Fleur calculations with the command line switch -last_extra. This will write another file out in every iteration that only contains the last charge density. One can then overwrite the cdn.hdf file with this additional file. Another way is to use Fleur command line options to actually delete densities from the cdn.hdf file. A list of all Fleur command line options can be obtained with the command line option -h. But in this exercise we will not manipulate the cdn.hdf file in such a way.
## Aiida-fleur
The somewhat cumbersome process of calling FLEUR over and over again is of course a clear task that can be automated using an aiida workflow. Using the aiida-fleur command line interface this can be done very easily.
We start by using the `inp.xml` but delete the `relax.xml` file created before
```bash
rm relax.xml
Now we can simply start the relaxation workflow. It will take quite some time
aiida-fleur launch relax
The workflow runs essentially the same steps as we performed manually. At the end, the following output is created:
- relax.json: The usual file with the main results of the workflow
- relax.xml: The
relax.xmlfile corresponding to theinp.xmlfile is copied from the AiiDA repo - energy/force_convergence.png: The simple plots visualize the convergence of the forces/energies during the structural relaxation
- cdn_last.hdf: Last charge density obtained.
To check the convergence of the relaxation we might want to investigate the plots and/or the relax.json file.
display < energy_convergence.png
cat relax.json
You could now try the -wf template.json option of aiida-fleur again and check the resulting wf_relax.json file for possible options of the workflow.
aiida-fleur launch relax -wf template.json
cat wf_relax.json
For example we might reduce the force_criterion to 0.0001 and restart the workflow.
Please modify the file accordingly
aiida-fleur launch relax -wf wf_relax.json
This time only a single additional relaxation step is performed until this new low convergence criterion is reached.
Learn more
To learn more about the possibilities for bandstructure calculations in FLEUR please visit: - https://aiida-fleur.readthedocs.io/en/latest/user_guide/workflows/relax_wc.htm - https://www.flapw.de/MaX-7.0/documentation/forces/