Understanding the different force contributions

When calculating the forces that act on the atoms of a crystal's unit cell, there are several contributions of varying size and computational cost that, when summed up, yield the final result. The total force on an atom is defined as the negative derivative of the total energy by the position of said atom. In the out file of a Fleur calculation all contributions are listed for each atom type and spin independently by referring to the equation numbers in the original paper they stem from. The brunt of it is found in a paper by Rici et al. (reference 1) - with the notable exception of the enhanced force formulae introduced by Klueppelberg et al. (reference 2). The standard contributions are as follows:

  • FORCES: EQUATION A3: The most basic Hellmann-Feynman force contribution from differentiating the Coulomb energy by the atomic position. Formula A3 in (1).
  • FORCES: EQUATION A4: The corresponding contribution from the core density. Formula A4 in (1).
  • FORCES: EQUATION A8: The first so-called Pulay contribution, that stems from acknowledging that the states need to differentiated as well when looking at the matrix element of an operator. This occurs in the differention of the eigenenergies. Formula A8 in (1).
  • FORCES: EQUATION A12: The Pulay contribution from equation A12 in (1). [Only written out in non-parallel calculations, as the alorithm yields wrong values when MPI-parallelized. The total force is correct nevertheless.]
  • FORCES: EQUATION A21: The Pulay contributions from equation A17+A20 in (1). [Only written out in non-parallel calculations, as the alorithm yields wrong values when MPI-parallelized. The total force is correct nevertheless.]

Klueppelberg introduced additional so-called "force levels" that cover several effects not contained in the previous terms, namely Pulay core terms and surface integrals due to the discontinuity of quantities at the MT-boundaries:

  • Level 1: Pulay terms from the core states leaking into the interstitial region, that are added onto the A4 contribution (found in out file as FORCES: IS/MT ADDITION TO EQUATION A4).
  • Level 2: Surface integral terms of the kinetic density and eigenenergies (modifies the output of A12).
  • Level 3: Surface integral terms of the effective potential, exchange-correlation potential and exchange-correlation energy density (found in out file as FORCES: SURFACE CORRECTION).

To access these additional force contributions, the f_level switch needs to be added to the inp.xml:

<geometryOptimization l_f="T" f_level="n" forcealpha="1.00" forcemix="BFGS"
                      epsdisp=".00001" epsforce=".00001"/>

The default value hardcoded into Fleur is -1. Setting n to 0 gives the standard force calculation albeit with additional output files to be used with the phonopy code by Atsushi Togo and 1-3 obviously coincide with levels 1-3. For the impact of the enhanced force routines on the calculation time and quantities like the drift force (that should be as close to 0 as possible to enable e.g. phonon calculations in a finite displacement approach) consult (2).