More advanced options in the XML based input file

Setup of non-collinear magnetism in the XML input file

To configure a Fleur calculation incorporating non-collinear magnetism, some parameters have to be set in the calculation setup section and further parameters have to be set for each atomGroup in the atom groups section. In each of these sections nocoParams elements have to be added. Templates with default parameters are generated by using the input generator with the -explicit command line option.

An example for the nocoParams element in the calculation setup section is:

      <nocoParams l_ss="F" l_mperp="F" l_constr="F" l_disp="F" sso_opt="FFF" mix_b=".00000000" thetaJ=".00000000" nsh="0">
         <qss>.0000000000 .0000000000 .0000000000</qss>
      </nocoParams>

The following attributes have to be set here:

Attribute Description
l_ss This boolean switch is used to activate spin-spiral calculations.
l_mperp Here the output of the magnetization perpendicular to the chosen axis can be activated.
l_constr Switch this on to constrain the magnetic moments.
l_disp This is the dispersion switch. If l_disp is set to true, the Force theorem is used to calculate the sum of eigenvalues for each predefined qss.
sso_opt This attribute sets three logical switches used in the context of [[SsoDetails
mix_b This is a mixing factor. If l_constr is set to true then the constraint field is mixed. In this case mix_b="0.5" should work fine. In the case of an atom with l_relax being set to true the input/output-directions of the moments are mixed. Here you can choose mix_b > 1 (e.g. 4).
thetaJ Used for the [[HeisenbergInteractionParametersJij
nsh Used for the calculation of Heisenberg J_ij interaction coefficients. This is the number of shells of neighbors to be considered.

The enclosed XML element is used to define the spin spiral vector in reciprocal lattice vectors.

An example for the nocoParams element in each atomGroup element of the atom groups section is:

         <nocoParams l_relax="F" l_magn="F" M=".0000000000" alpha=".0000000000" beta=".0000000000" b_cons_x=".0000000000" b_cons_y=".0000000000"/>

The following attributes have to be set here:

Attribute Description
l_relax This has to be set to true to relax the directions of the magnetic moments at the associated atoms.
l_magn Used for the calculation of Heisenberg J_ij interaction coefficients. Set this to true iff the atoms in the group are magnetic.
M Used for the calculation of Heisenberg J_ij interaction coefficients. The value of the magnetic moment(including sign) has to be set here.
alpha The 1st angle that determines the magnetic structure. It is equal to "phi" in spherical coordinates.
beta The 2nd angle that determines the magnetic structure. It is equal to "theta" (measured from the z axis) in spherical coordinates.
b_cons_x,b_cons_y These are the constraint fields in x and y direction. They are determined self-consistently if l_constr is set to true.

Setup of LDA+U calculations

To amend the description of electron correlations in local and semilocal XC functionals, up to 4 Hubbard U parameters can be defined for each species in the atom species section. For this optional ldaU XML elements have to be inserted into the respective section below the energyParameters, electronConfig, and nocoParams entries and above the lo entries. The following example demonstrates how an ldaU element looks like:

         <ldaU l="2" U="8.0" J="0.9" l_amf="F"/>
Attribute Description
l The angular momentum quantum number of the orbital for which the U parameter is set.
U The U parameter in eV.
J The J parameter in eV.
l_amf This logical switch determines whether the "around mean field" limit (if true) or the atomic limit (if false) is used.

Mixing of the density matrix

Whenever a Hubbard U parameter is added to an atom not only the density has to be part of the mixing from iteration to iteration but the density matrix, too. For this additional parameters can be set in an optional ldaU XML element (different from the one above) in the calculation setup section. Such an element looks like:

<ldaU l_linMix="F" mixParam="0.05" spinf="1.00"/>
Attribute Description
l_linMix This switch determines whether a straight mixing algorithm is applied to the density matrix (if true) or the mixing of the density matrix will be performed like the mixing of the density (if false). The switch is optional and set to false by default.
mixParam This is the optional mixing parameter that is used for the straight mixing of the density matrix. By default this parameter is 0.05.
spinf Optional, default ist 1.0.

If the ldaU XML element in the calculation setup section is not present all parameters that can be set in it have their default values.

Further reading

Setup of core spectrum calculations for EELS in the XML input file

In the output section it can be specified that a core spectrum has to be calculated. For this the attribute coreSpec has to be set to "T" and in the output section the optional xml element coreSpectrum has to be defined analogously to the following example:

      <coreSpectrum eKin="300.0" atomType="1" lmax="2" edgeType="L" eMin="-1.0" eMax="15.0" numPoints="17" nqphi="10" nqr="10" alpha_Ex="0.024" beta_Ex="0.050" I_initial="155" verbose="T">
         <edgeIndices> 3 </edgeIndices>
      </coreSpectrum>
Attribute Description
eKin The kinetic energy of the incoming electron in units of keV. The relativistic correction term, which was shown to be important (see PhD thesis of Kevin Jorissen ), is proportional to eKin.
atomType This is the index of the atom group for which the EELS is to be calculated.
lmax Maximum value of the angular momentum to be considered in the EELS calculation. By setting it to a reasonably low value (e.g. 2 or 3), the calculation is significantly faster in comparison with no lmax restriction.
edgeType The edge selection: K, L, M, ...
eMin The bottom edge of the energy interval with respect to the Fermi energy for which the EELS is calculated.
eMax The top edge of the energy interval with respect to the Fermi energy for which the EELS is calculated.
numPoints Number of points in the [eMin,eMax] interval for which the EELS is calculated..
nqphi Number of angular divisions of the q-mesh, optional, standard value: 10
nqr Number of radial divisions of the q-mesh, optional, standard value: 10
alpha_Ex Experimental convergence semi-angle in rad, optional, standard value: 0.024
beta_Ex Experimental collection semi-angle in rad, optional, standard value: 0.05
I_initial Initial intensity of the incoming electron beam, optional, standard value: 155
verbose Verbose output is produced, optional, standard value: F

In the edgeIndices element a list of space separated indices has to be provided.

Using the force-theorem in FLEUR calculations

Fleur has the option to calculate some properties using the force-theorem. In such calculations the last self-consistency iteration (or the only iteration if itmax=1) is replaced by a loop in which the eigenvalue sum for different configurations at a fixed potential are calculated.

This is controlled by a special section in the inp.xml file to be inserted after the output-section:

   <forceTheorem>
      <MAE theta="0.0 0.1*Pi" phi="0.0 0.2*Pi" />
   </forceTheorem>

Here several different modes are possible. Exactly one should be present:

Magnetic Anisotropy Energy MAE

      <MAE theta="0.0 0.1*Pi" phi="0.0 0.2*Pi" />

This is a simple mode to estimate the magnetocyrstalline anisotropy energy. A loop over different spin-quantization directions is performed. A list of angles should be given. Please note that the number of "theta" and "phi" angles must of course be the same.

Spin-Spiral Dispersion

      <spinSpiralDispersion>
          <q> 0.0 0.0 0.0 </q>
          <q> 0.1 0.0  0.0 </q>
     </spinSpiralDispersion>

This is a simple mode to calculate a spin-spiral dispersion for several q-values. Please note that the first q-Vector given will overwrite the q-vector specified in the '''qss'''-tag given in the '''nocoParams'''.

Dzyaloshinskii Moriya Interaction

      <DMI theta="0.0 0.1*Pi" phi="0.0 0.2*Pi" >
         <qVectors>
             <q> 0.0 0.0 0.0 </q>
             <q> 0.1 0.0  0.0 </q>
         </qVectors>
     </DMI>

This mode is actually slightly different from the modes above as it actually does not only calculate the eigenvalue sum for different setups, but also employs first order perturbation theory to estimate the effect of spin-orbit coupling on a spin-spiral calculation. Hence you can specify here a list of q-vectors as well as different angles for the magnetisation.

Structure relaxations with FLEUR

ATTENTION: this documentation describes work in progress. Your version might not support it yet

General notes:

  1. Structural relaxations should be performed using the HDF5-version of the code to enable better handling of charge densities after a relaxation step
  2. You should be aware that good forces in LAPW are only obtained if you perform very accurate calculations. In particular you should:
    • Use high cutoffs for your LAPW-basis set (high Kmax and high lmax).
    • Use the core-tail correction ctail=T possibly checking the influence of coretail_lmax.
    • Use LOs for semi-core states as tails or core states can give contributions to forces otherwise not covered.

Switching on the calculation of forces

To calculate forces on an Atom use the ''force'' tag in the ''atomgroup'' tag.

<force calculate="T" relaxXYZ="TTT"/>

For each atom you specify if forces are calculated and which of the directions should be used. Please be aware that setting this tag alone only enables the calculation of the contributions to the force simple to obtain. The full force including the Pulay terms are calculated by setting l_f="T".

Relaxing the structure

To calculate all forces (including Pulay terms) and to perform a structure relaxation you have to specify the ''geometryOptimization'' tag in the input:

<geometryOptimization l_f="F" epsdisp="0.001" epsforce="0.001" forcemix="2"  forcealpha="1.0" qfix="0" force_converged="0.00001" />

All attributes except the l_f are optional with defaults as specified above.

Attribute Description
l_f This switches on the calculation of the full force including the more expensive to obtain Pulay terms
epsdisp This is the minimal displacement at which the relaxation is considered converged
espforce The minimal force at which the relaxation is considered converged
forcemix The scheme to use for calculating displacements
forcemix=0 Simple relaxation by moving the atoms in the direction of the force (forcealpha gives the corresponding scaling factor)
forcemix=1 A CG scheme for relaxations
forcemix>1 A BFGS scheme to determine new displacements.
forcealpha This is the scaling factor used to shift the atoms in the ''forcemix=0'' case or in the case of no history.
qfix The qfix parameter determines how to deal with non-charge neutral densities. These occur if you want to reuse the charge after new displacements have been calculated. In addition, this parameter determines the number of times the code checks of charge neutrality in a self-consistency run. The logic is as follow:
even number leads to less frequent checks. This should be OK in general.
odd number leads to frequent checks as in older FLEUR versions
numbers 0/1 will just fix the charge after a displacement
numbers 2/3 will fix the charge by adjusting the interstitial charge only. This is a better approximation if your charge was generated by a displaced configuration.
numbers 4/5 will try to use an additional decomposition of the charge to get a better estimate for displaced positions (experimental)
force_converged the value given here determines a criterion if an actual relaxation step is performed. Only if the maximal change of force between two SCF iterations is below this value a structural relaxation step is done.

Output

If a relaxation step is performed the code stops with the message of either Structual relaxation: new displacements generated or Structual relaxation: Done if the forces or displacements are smaller than epsforce or epsdisp, respectively.

In addition a file relax.xml is generated. This file contains the current displacements and the history of the positions and forces from previous relaxation steps.

Shifting positions

During structural relaxation the atomic positions will update. These updates should be preformed by specifying additional tags in the inp.xml file. FLEUR creates such a tag-structure in the relax.xml file which looks similar to the one given here:

<relaxation>
   <displacements>
    <displace>   0.0000000000    0.0000000000   -0.0122134660 </displace>
    <displace>   0.0000000000    0.0000000000    0.0031814755 </displace>
    <displace>   0.0000000000    0.0000000000   -0.0044667051 </displace>
    <displace>   0.0000000000    0.0000000000    0.0179966687 </displace>
  </displacements>
   <relaxation-history>
    <step energy="   -22005.9436915333">
      <posforce>   0.0000000000    0.0000000000    0.5000000000    0.0000000000    0.0000000000   -0.0111067330 </posforce>
      <posforce>   0.5000000000    0.5000000000    0.0350000000    0.0000000000    0.0000000000    0.0030907378 </posforce>
      <posforce>   0.5000000000    0.0000000000    0.1130000000    0.0000000000    0.0000000000   -0.0042333525 </posforce>
      <posforce>   0.5000000000    0.5000000000   -0.3960000000    0.0000000000    0.0000000000    0.0189983343 </posforce>
    </step>
   </relaxation-history>
 </relaxation>

You should include this file into inp.xml by incorporating a line like

<xi:include xmlns:xi="http://www.w3.org/2001/XInclude" href="relax.xml"> <xi:fallback/> </xi:include>

just before the closing </fleurInput> tag. This is done automatically by the input-generator now. The empty fallback enables you to use this line even if no relax.xml is present. The code will check that no overlapping MT-spheres arise from your displacements and it will decrease them if necessary. Hence you have to make sure that no such warning for overlapping MT-spheres still remain if your structure is considered converged or you might be stuck with a configuration limited by the size of your MT-spheres.

Magnetic fields

There are two options to apply B-Fields:

  1. You can specify an overall B-field in the fields-tag of the calculation setup:
<calculationSetup>
      ....
      <fields b_field="0.1"/>
</calculationSetup>
  1. If you want to have a field only applied within the MT-sphere of a single atom (old mfee-file) you should use the tag in the species definition:
<species ....>
      ....
      <special b_field_mt="0.1"/>
     ....
</species>

Electric Fields in FLEUR

Here, we describe how to put a metallic system in an electric field. For the actual application you might want to jump to see the [[#file_settings | input file settings ]]. In case of nonmetallic systems, you might be interested to find the description [[#asymmetricfield | how to place an asymmetric field]] or in an electric field with [[#Dirichlet|Dirichlet boundary conditions]].

The screening charge

The ability to include applied electric fields in our electronic structure calculations is a very powerful tool. We can then begin to study magnetic tunnel junctions as these systems are a junction under the influence of an applied E-field or we can study field induced reconstructions at surfaces. It opens the door for us to a whole new class of systems and effects.

OK so how do we do it then?

In the FLEUR code we have a film, so the simplest method of introducing a field is to place a sheet of charge on either side of the film in the vacuum, this then provides an electric field perpendicular to the surface of the film. We need to ensure charge neutrality of the system so we must subtract the same amount of charge from the slab as we have placed on the charged sheets; that's where our screening charge comes from! The figure below shows schematically our film, the charge sheet and the induced screening charge. To answer the question that everyone asks, yes we could just put equal and opposite charges on either side of the film and then not have to subtract any charge from the film, we just have not implemented it yet! [Note: That has changed meanwhile, see below.]

Can we test that it works OK?

Yes, we can test it. We can show that there is a relationship between the change in electronic occupation and the the chemical potential: Remembering that applying an electric field corresponds to changing the occupancy we now have a test. We have two ways of calculating the change in total energy as a function of the change in occupancy (which corresponds to the change in field). We can calculate it directly from the total energy of two calculations (one with an applied field, one without) and we can calculate it by integrating the curve that represents the chemical potential as a function of applied field. The results of the test are shown below.

We see from the above figures that the open squares that represent the results obtained from integrating the chemical potential curve sit almost perfectly centered upon the results calculated directly from the total energy calculations, so we can be confident that the method works well.

Example results

The first results show the screening charge induced in the Ag(001)c(2x2)-Xe structure when a field is applied. We can see that the screening still occurs at the metal surface and not where the adsorbate sits and that the adsorbate atom becomes polarised.

The second example shows the spin-resolved screening charge in Fe(001). Here we see that even though the screening charge in the separate spin channels penetrates right to the center of the slab, the total screening charge only penetrates as deep as the first layer - the field is screened out.

Input file settings

There are two basic ways of specifying electric fields:

  1. The values of the sheets of charge for external electric fields is set by requiring charge neutrality. Thus, most easily you can add or subtracting a fractional charge by changing the number of valence electrons. The resulting field is shown in the external electric field section of the out file.
<bzIntegration valenceElectrons="8.00100000" ...
  1. More complex settings are possible using the <fields> tag:
<calculationSetup>
  ....
  <fields zsigma="0.0" sig_b_1="0.0" sig_b_2="0.0" autocomp="T" dirichlet="F" eV="F">
      <shape> shapestring </shape>
  </fields>
</calculationSetup>

The following attributes are provided: * zsigma: the position of the sheets of charge relative to the vacuum boundary (set by default to 10 a.u. (= 5.291772 Å)). * sig_b_1/2 for charges on the upper and one for the lower plate (default 0.0). Setting these to different values enables to place an asymmetric field. * the autocomp switch makes sure that overall charge neutrality is automatically calculated. * the dirichlet switch enables use of Dirichlet boundary conditions. * the eV switch is used to modify units in the dirichlet="T" case. * in addition an (unlimited) number of <shape> tags can be given to specify inhomogenous fields.

Schematic representation of a two-dimensional film in Fleur

---------------------------------           -.
                                              \
================================= -.           \ delz*nmz  = 25 a.u.
                                    \          /
.................................    > zsigma /       -.
                                    /        /          \
********************************* \          \           \
*********************************  > z1 = D/2 \           \
********************************* /            > D = Dvac  > Dtilde
*********************************             /           /
*********************************            /           /
                                                        /
................................                      -'


================================

--------------------------------

Placing asymmetric fields

Sometimes, it might be useful to place a thin film in an antisymmetrc or even asymmetric field. Since this requires to place charges of opposite sign on both sides of the film, it is necessary to provide additional input in the tag as mentioned above. It allows to place the film into a condenser that creates the electric field like that:

In the image above, we omitted the film for clarity and just show the effect of the plates put at +/-10 a.u. in the vacuum. For polar films this is the correct setup, since it allows to have a flat potential in the vacuum and a vanishing electric field inside the film. Notice, however, that for large applied fields the potential in the vacuum can drop below the Fermi-level and you can get field-emission (i.e. the program stops with a message like vacuz).

Placing fancy inhomogeneous fields

Since version 0.26b inhomogeneous fields can be generated: * 'sig_b_1/2' contain the additional (homogeneous) charge for the top and bottom sheet. By default, excess (positive/negative) charge of the film is compensated by equally charging the charge sheets; if 'autocomp' is false, the user has to do this manually. Note: Fleur requires an overall (film plus top plus bottom sheet) charge neutrality. * the inhomogeneous charge can be places using the in which stings are supplied. These strings specify the inhomogenous charges using the key-words rect, circ, rectSinX, polygon, and datafile. Their detailed syntax is:

  rect <sheet> <x>,<y> <width>,<height> <charge> [options]

  circ <sheet> <x>,<y> <radius> <charge> [options]

  rectSinX <sheet> <x>,<y> <width>,<height>  <amplitude> <n> <delta> [options]

  polygon <sheet> <n_points> <x1>,<x2> ... <x_n>,<y_n> <charge> [options]

  datafile <filename> [zero_avg] [options]

Note that all positions and lengths are currently relative values (i.e. between 0 and 1). The sheet to be used can be set using , which can be either top, bot or topbot/bottop. Options are: add (default) to add the charge to the charge of previous tags or replace to use the new charge instead; zero to place the charge only to areas which were before zero, nonzero to place it at areas which where before nonzero or all (default) to place it in the whole area covered by the tag.

Note: The regions can exceed the unit cell plane and then cut off, e.g. circ top 0,0 0.25 0.5 places half an electron in a quarter circle with origin (0,0). Note, however, that circ creates a perfect circle only on the grid; this generates a circle and not an ellipse only if the ''k1d''/''k2d'' ratio matches the crystal's ''a''/''b'' ratio.

rectSinX creates a sinodal potential in ''x'' direction (constant in ''y'' direction for any ''x'' value), i.e. {$V(x,y) = A \sin(2\pi nx + 2\pi\delta)$}, where ''A'' is the amplitude; however, the argument in apwefl is not ''A'' directly but {$A' = A L_z$}, where {$L_z$} is the number of points in ''z'' direction. Contrary to circ and rect, charges are mask out without being redistributed to non-masked positions. It is {$\int_0^{2\pi} A|\sin(x)| {\rm d}x = 4A'$}, n is the order and delta ({$\delta$}) the offset.

polygon creates a polygon-shaped charge distribution; note that the currently used algorithm does not always give the perfectly shaped polygon - and the edge points are not always included in the polygon.

The datafile reads the data from a file; if a zero_avg has been given, the charge is averaged to zero, i.e. only the inhomogeneous contributions are taken into account. The option replace/add is supported, but zero/not_zero is not. The syntax of the data file itself is as follows. First line: number of ''x'' and ''y'' points; second line: charge for point (x=1,y=2), third: (x=1,y=2) etc. The number of points must be 3*k1d and 3*k2d.

Example 1: To have two top plates (segments):

  rect top 0,  0  0.5,1.0  0.2
  rect top 0.5,0  0.5,1.0 -0.2

Example 2: To have a charged ring with 0.2e and -0.2e of charge evenly distributed outside this ring:

  circ topBot 0.5,0.5 0.2  1           ! Create temporary an inner ring
  circ topBot 0.5,0.5 0.3  0.2 zero    ! Create outer ring
  circ topBot 0.5,0.5 0.2  0   replace ! Delete inner ring
  rect topBot 0,0     1,1 -0.2 zero    ! Place smeared opposite charge

Using Dirichlet boundary conditions

In the scheme above, a (locally) constant surface charge density has been assumed, which matches a Neumann boundary condition. Another common boundary condition is a constant potential (Dirichlet boundary condition) -- such as provides by a metal plate, which can be either grounded or kept at a certain potential. In order to use Dirichlet boundary condition, add 'dirichlet="T" ' to the tag, which has been introduced above. The value given for 'sig_b_1/2' and in the different shapes are regarded as potential in Hartree or, if 'eV="T"' is set, as potential in (electron) volts.