3.11 Energy derivatives (forces, stress) and geometry optimization

With self-consistent Kohn-Sham orbitals, densities and total energies available, one of the primary tasks of electronic structure theory is to obtain energy derivatives with respect to the nuclear coordinates $\mathbf{R}_{\text{at}}$ . First derivatives (“forces”) allow to find the optimum structure and molecular dynamics on the Born-Oppenheimer surface. Based on the structure optimum, second derivatives then enable the calculation of the (Born-Oppenheimer) zero-point vibrational energy of the nuclei, and of vibrational excitations (phonons).

The present section deals with the options available in FHI-aims for the computation of forces in FHI-aims, and with algorithms related to geometry optimization. Specifics regarding first-principles molecular dynamics are given in Sec. 3.12, and the computation of second energy derivatives (vibrational frequencies, zero-point energies, and oscillator strengths) by a finite-difference technique is covered in Sec. 4.6.

Before coming to the full set of related keywords, we give some basic comments below. Please also consider using (and adding to) the wiki, especially if you encounter any kind of behavior that is obviously not intended. We would like to know about such cases.

One “most important” rule:

For efficient local structure optimizations, do not simply use “tight” settings from an arbitrary starting geometry and wait.

It is usually much more effective to use a two-step approach. First, use “light” settings in a pre-relaxation step to get the rough geometry right quickly. Then, use “intermediate” or “tight” settings for post-relaxation or post-processing only, e.g., based on the “geometry.in.next_step” file that is written out by the code by default (see below).

Basic handling:

For most applications, reliable structure optimizations of atomic coordinates into a local minumum of the potential energy surface can be obtained by simply setting the keyword
relax_geometry bfgs threshold
threshold denotes the minimum absolute force component in eV/Å acting on any atom. Typically, a threshold value of 5d-3, i.e., 5 $\cdot$ 10^-3 eV/Å, corresponds to a very tightly converged structure optimization. Do not use much lower values since you might spend substantial amounts of CPU time for no measurable gain.

For unit cell optimizations in periodic systems, the keyword relax_unit_cell must additionally be invoked (otherwise, the unit cell shape will remain fixed).

Finally, individual coordinates or degrees of freedom can be fixed by using the keyword constrain_relaxation.

Another method to run a constrained relaxation is the combined use of the keywords symmetry_n_params, symmetry_params, symmetry_lv and symmetry_frac. When this block is added to a geometry.in and relax_geometry (and optionally relax_unit_cell) are set in control.in, a geometry relaxation is started that is constrained to a parameter reduced space defined by the user. This can be used to enforce a symmetry-preserving relaxation by providing the exact prototype of the given structure in its analytic form. It furthermore allows for the introduction of local symmetries or local symmetry breaking like distortions. Important: When using these constraints all keywords must be specified for all lattice and atomic degrees of freedom. If you want to run with these constraints on only the lattice or the atomic positions, add the other block fixing the coordinates or with $3N_{Lattice}$ or $3N_{atomic}$ additional parameters for a free relaxation. Important: The relationship between the full coordinates and the reduced parameters has to be linear, i.e. there is a transformation matrix A and a transformation matrix B that can map the flattened 1-dimensional representation of the lattice cell/atomic positions to the parameter-reduced space. For the monoclinic and triclinic lattice systems, e.g. do not use $c\cdot\cos(\beta)$ and $c\cdot\sin(\beta)$ but replace them with independent variables e.g. $d$ and $e$ . It is on the to-do-list to make aims recognize these cases internally.
Hint: You can create geometry.in’s already containing the required input block for a symmetry-constrained relaxation using AFLOW as of version 3.1.204. It gives access to hundreds of structure prototypes in the AFLOW Prototype Library [217, 141] throughout all spacegroups. Simply add --add_equations to the command. Example:

aflow --proto=AB_cF8_225_a_b --params=5.64 --aims --add_equations

To ensure the constraints are mapped onto the positions of the atoms inside of the unit cell (AFLOW’s default behavior) use the aflow_structure_wrapper utility in utilities Example:

aflow_structure_wrapper --proto=AB_cF8_225_a_b --params=5.64 -o geometry.in.aflow

FHI-aims checks to ensure these constraints do not change the structure too much with the maximum allowed change in fractional coordinates set by symmetry_frac_change_threshold in control.in.

To obtain gradients for a given structure, but without any subsequent structure optimization, molecular dynamics run, etc., employ the keywords compute_forces, compute_analytical_stress or compute_numerical_stress instead.

relax_geometry bfgs value calls a trust radius enhanced variant of BFGS, also available as trm. This variant covers all use cases.

The bfgs algorithm can be tweaked by specifying an initial guess for the Hessian matrix, which can influence the performance. The default is the model Hessian matrix by Lindh and coworkers,[199] which leads to reliable and fast optimizations in the vast majority of scenarios.

If, for any reason, the BFGS algorithm does not perform well for a given structure initially, a number of options are available. The standard first step should always be to simply restart the optimization from the updated geometry and Hessian matrix in geometry.in.next_step (see below). Alternatively, one may preset a simple, diagonal Hessian using the keyword init_hess diag value (see the actual keyword description for more details). Finally, the energy_tolerance and harmonic_length_scale keywords decide when the bfgs algorithm will abort a relaxation because a presumed deviation between the predicted and the real energy landscape. They can be set to more forgiving values if needed.

For the relax_geometry bfgs algorithm, the handling of restarting relaxations is formalized by writing out a file geometry.in.next_step after each relaxation step. By default, this file contains both the geometry and the current estimate of the Hessian matrix of the system that are used by the BFGS algorithm. For an organized restart, simply copy this file to geometry.in, and the stored geometry and Hessian will be used to reinitialize the BFGS algorithm if the control.in file requests a structure relaxation. If the keyword distributed_hessian is in use, the current estimate of the Hessian matrix will be stored in a separate binary file, which should not be modified by hand.

Use the ’get_relaxation_info.pl’ utility script to monitor the progress of an ongoing or finished relaxation run. This information can be immensely helpful to make sure that you are not, e.g., spending your time optimizing the last 10^-5 eV out of an already converged structure relaxation.

To stop an ongoing structure relaxation in an organized way, create the abort_opt file in the respective directory.

MP2 and RPA total energy derivatives with pz-lda and pbe potential:

For running particle-hole RPA total energy derivatives calculation following tags in control.in file are necessary: One should set (keyword rpa_force freq_formula_method or matrix_diag_method) depending on two different approaches. The first one use frequency integral formula for evaluation of correlation energy which scale $N^{4}$ , with $N$ being the system size. The later one scales $N^{6}$ , comparatively consume more time as compared to first one. Others additional tags required are as RI_method lvl, use_2d_corr .false., DFPT vibration or vibration_reduce_memory if doing calculation with pz-lda, with pbe only use vibration_reduce_memory, total_energy_method rpa, this is optional does not effect anything. Calculation with bigger Gaussian basis set, the accuracy is affected with the choice of number of auxiliary basis sets control by tag prodbas_acc, it is highly recommended to use larger values like 1.E-2. Currently, MP2 and RPA forces are evaluated for closed shell systems only. For MP2 force, same (keyword rpa_force mp2_force) invoke. Frozen-core approximation can be invoked by using additional tag frozen_core_postSCF 1.
Some useful tips to run RPA force calculations are as

frequency_points 12 # for bigger systems
rpa_force  freq_formula_method
use_2d_corr  .true. # for bigger systems

There are several version of RPA forces which allows the additional tags as

post_SCF_force rpa_parametric
force_reduce_memory  .true.

post_SCF_force rpa_parametric
least_memory  .true.

post_SCF_force rpa_parametric
least_memory_1  .true.

and two more options are least_memory_2 and least_memory_3.

The RPA+SE and RPA+rSE total energy derivatives:

The forces at the RPA+SE (RPA+single excitation) and RPA+rSE (re-normalized single excitation) level can also be calculated by using additional tag to RPA calculations by setting post_SCF_force to SE or rSE. Additionally, for RPA, RPA+SE and RPA+rSE calculations for big molecules and with bigger basis set the tag force_reduce_memory .true. will be quite helpful.

The RPA single point energy with RI-LVL can be computed by using tag rpa_ene_lvl .true.. For this purpose total_energy_method rpa will only work.

The RPA total energy derivatives with PBE potential for solids
For running RPA total energy derivative calculations for solids rpa_force freq_formula_method is needed together with DFPT phonon_reduce_memory, total_energy_method rpa, freq_grid_type minimax(optional), RI_method LVL, and use_2d_corr .false. is working now. If use_2d_corr is not set manually then, if k-points are more than number of chosen core then it set use_2d_corr .false. automatically and vice versa. If you are only interested in force calculation then compute_forces .true. or sc_accuracy_forces needed only.

Stress tensor:

For unit cell optimization (keyword relax_unit_cell), an analytically computed stress tensor will be used where available, thanks to work by Viktor Atalla, Christian Carbogno, and Franz Knuth [171]. For some density functionals, the analytical stress tensor is not available. In these cases, the unit cell itself can still be relaxed, but the stress tensor must be computed numerically from finite differences of total energies.

The numerical finite-difference stress tensor is always available as a fallback.

Finally, we note that, in some cases, just optimizing the unit cell shape blindly may not be what you want. For example, in high-symmetry structures a fit to Murnaghan’s Equation of state—see the utility provided in the utilities directory—will be more accurate and give you more information about the system (bulk moduli and, in case of more than one phase, also access to transition pressures).

Hybrid functionals:

For hybrid functionals, analytic energy gradients and the analytic stress tensor are only available together with the RI_method LVL_fast version of “resolution of identity” of the two-electron Coulomb operator.[151]

For perturbative methods (MP2 perturbation theory, RPA and beyond), analytical gradients are not yet available.

Tags for geometry.in:

Tag: constrain_relaxation(geometry.in)

Usage: constrain_relaxation constraint
Purpose: In geometry.in, fixes the position of the last specified atom/lattice_vector in a structure optimization.
constraint is a string, indicating what exactly will be constrained. Default: .false.

Allows to relax only parts of a structure, while keeping the rest at fixed positions. Currently, the following simple options for constraint are possible:

•

.true.: All coordinates for this atom will be constrained.
•

.false.: The relaxation of this atom will not be constrained.
•

x: The $x$ coordinate of this atom is not allowed to move.
•

y: The $y$ coordinate of this atom is not allowed to move.
•

z: The $z$ coordinate of this atom is not allowed to move.

Attention: If you wish to constrain more than one coordinate, the required constraints must be specified as separate lines, like this:

  atom 0. 0. 0. Fe
    constrain_relaxation x
    constrain_relaxation y

In contrast, specifying two constraints in one line will not work. The second constraint would simply be ignored!

Tag: hessian_block(geometry.in)

Usage: hessian_block i_atom j_atom block
Purpose: In geometry.in, allows to specify a Hessian matrix explicitly, with one line for each 3 $\times 3$ block. The option block consists of 9 numbers in column-first (Fortran) order. The 3 $\times$ 3 block corresponding to j_atom, i_atom is initialized by the transposed of block. The Hessian matrix is input in units of eV/Å².

If at least one hessian_block line is found in the file, the Hessian is constructed using this mechanism. So far there is no safe-guard from overriding Hessian blocks with subsequent lines with equal i_atom, j_atom.

There are two scripts in the utilities directory to automatically construct such Hessian matrix approximations. First, conversions/thess2aims.py converts a Tinker generated Hessian matrix. Second, Lindh.py constructs a general purpose model matrix [199]. Please note that the Lindh model matrix is now also directly available with init_hess Lindh.

Tag: hessian_block_lv(geometry.in)

Usage / purpose: Like hessian_block, but for Hessian matrix elements between lattice vector degrees of freedom.

Tag: hessian_block_lv_atom(geometry.in)

Usage / purpose: Like hessian_block, but for Hessian matrix elements between lattice vector degrees of freedom.

Tag: hessian_file(geometry.in)

Usage: hessian_file
Purpose: In geometry.in, this keyword indicates that there exists a hessian.aims file to be used to construct the Hessian.

If hessian_file is found in geometry.in, the Hessian is constructed using the data in hessian.aims, which is a binary file generated by geometry optimization calculations with FHI-aims (please see distributed_hessian). The user should not try to create such a file by hand.

Tag: trust_radius(geometry.in)

Usage: trust_radius value
Purpose: In geometry.in, allows to specify the initial trust radius value for the trm relaxation algorithm.
value is the initial value set for trust radius, in Å. Default: 0.2 Å.

This keyword is a significant exception – it is the only algorithmic keyword found in the geometry.in file. Conceptually, it would belong into control.in, but since it is written out as part of the geometry.in.next_step file which can be used to restart a structure relaxation, we keep it here.

The keyword specifies the initial value of the “trust radius” used to limit structure relaxation steps determined by the bfgs (synonymous with trm) relaxation algorithm.

The default is set to 0.2 Å. It will be overridden by the default value of max_atomic_move, which can be set in control.in.

Specifically, the trust radius is compared to combined step length $D X$ as proposed by the bfgs / trm algorithm across all active coordinates $X_{l}$ , i.e., all coordinates that are not explicitly constrained:

DX=\sqrt{\sum_{l=1}^{N_{\mathrm{coords}}}(X_{l}^{\mathrm{new}}-X_{l}^{\mathrm{% prev}})^{2}}.

(3.40)

$X_{l}^{\mathrm{new}}$ are the new coordinates that were proposed by the bfgs / trm algorithm, starting from the last accepted coordinates $X_{l}^{\mathrm{prev}}$ for which the total energies and energy gradients were calculated. The active coordinates $X_{l}$ include all $N_{\mathrm{coords}}$ unconstrained atomic positions ( $x_{I},y_{I},z_{I}$ , where $I$ =1, …, $N_{\mathrm{atoms}}$ ) and, if applicable, all unconstrained latticed vector coordinates ( $a_{i,x},a_{i,y},a_{i,z}$ , where $i$ =1, 2, 3).

Since the trust_radius will be dynamically adjusted during a structure optimization run, an allowed minimum value can be set using the min_trust_radius keyword in control.in. If the trust_radius attempts to fall below this minimum value, the computation will stop since the forces likely deviate from the underlying energy surface.

Tag: symmetry_n_params(geometry.in)

Usage:symmetry_n_params n_total n_lv n_coords
Purpose: In geometry.in, specifies the number of parameters to be optimized in a symmetry-constrained relaxation. n_total is the total number of parameters, n_lv is the number of parameters that define the lattice cell, n_coords is the number of parameters defining the fractional atomic positions.
n_lv + n_coords = n_total
Example for an orthorhombic cell and no parameters for the atomic positions:
symmetry_n_params 3 3 0

This keyword also serves as a flag, setting use_symm_const_geo_relaxation internally to True. If n_lv is 0, then no unit cell relaxation is done. For unit cell relaxations make sure to set relax_unit_cell to full.

Tag: symmetry_params(geometry.in)

Usage:symmetry_params [list of variables for parameters]
Purpose: In geometry.in, list all variables that are used as parameters separated by blanks. Always list the lattice parameters first. The number of variables used for parameters has to be equal to n_total specifies in symmetry_n_params.
Example for an orthorhombic cell and no parameters for the atomic positions:
symmetry_params a b c

Tag: symmetry_lv(geometry.in)

Usage:symmetry_lv x , y , z
Purpose: Specifies the analytic form of the lattice vectors. Use exactly the same cell as given in lattice_vector (same order of lattice vectors, same orientation), and simply replace entries that are free to relax by their analytic expression in terms of parameters specified by symmetry_params. Example for an orthorhombic cell:
symmetry_lv a , 0 , 0
symmetry_lv 0 , b , 0
symmetry_lv 0 , 0 , c

Note the comma between the components of the lattice vectors.

Tag: symmetry_frac(geometry.in)

Usage:symmetry_frac $n_{1}$ , $n_{2}$ , $n_{3}$
Purpose: Specifies the analytic form of the fractional atomic positions. Use exactly the same form as given in atom_frac (same order atoms), and simply replace entries that are free to relax by their analytic expression in terms of parameters specified by symmetry_params.
Example:
symmetry_frac 0.0 , 0.0 , x1
symmetry_frac 1./2 , 1./2 , x1 + 1./2

Note the comma between the fractional coordinates of the atoms and that no species entry is needed.

Tag: symmetry_frac_change_threshold(control.in)

Usage:symmetry_frac_change_threshold max_change
Purpose: Specifies the maximum allowed change in the initial structure to any fractional coordinate after applying the parametric constraints. If set to 1.0 then these tests will be ignored. max_change is the maximum allowed change in the fractional coordinates of an atom Default: 0.01

Tags for general section of control.in:

Tag: aggregated_energy_tolerance(control.in)

Usage: aggregated_energy_tolerance tolerance
Purpose: Sets the energy amount by which the energy across an entire relaxation trajectory may ever go uphill, based on the lowest known energy so far.
tolerance is a positive real number, in eV. Default: 3 $\cdot$ 10^-3 eV.

Small uphill steps of a relaxation trajectory are allowed up to the keyword energy_tolerance, but a relaxation trajectory should never go uphill for an extended number of steps in small uphill increments. The keyword aggregated_energy_tolerance sets an overall cap for any accepted uphill steps across an entire relaxation trajectory. The default is much larger than the allowed energy_tolerance in a single step and should, in principle, never be breached. If the aggregated_energy_tolerance criterion triggers, please contact us.

Tag: calc_analytical_stress_symmetrized(control.in)

Usage: calc_analytical_stress_symmetrized flag
Purpose: If .false., calculates all 9 components of the analytical stress tensor. If .true. calculates only the upper triangle (6 components) of the tensor and copies the result to the lower triangle.
flag is a logical string, either .true. or .false. Default: .true.

Generally, it is sufficient to calculate the upper triangle of the tensor. This flag is mainly for debugging purposes.

Tag: clean_forces(control.in)

Usage: clean_forces type
Purpose: Can remove small unitary force components (rotation and translation of the whole structure due to residual numerical noise) in relaxations.
type is a string, specifying whether and how any overall rotations / translations are removed.

The default for type depends on the exact circumstances (see below). The following choices exist:

•

none : No removal of residual rotations and translations. This is the default if any external embedding fields or charges are specified in geometry.in.
•

sayvetz : Non-periodic structures: Removal of rotations and translations by a formal projection [84, 270]. In periodic systems, only translations are removed.
•

fixed_plane : experimental Simple alternative algorithm by constraining three atoms into a plane (implicitly constraining all others).

Tag: compute_analytical_stress(control.in)

Usage: compute_analytical_stress flag
Purpose: If .true., switches on the computation of the analytical stress tensor.
flag is a logical string, either .true. or .false.

Default: .true. if a unit cell relaxation was requested and computation is possible. Otherwise, .false.

This flag allows to request an explicit analytical stress tensor computation for an otherwise explicit single-point calculation.

The calculation of the analytical stress is limited to LDA, GGA, Meta-GGA and hybrid functionals and is not possible with load_balancing. The vdW correction based on Hirshfeld partitioning (vdw_correction_hirshfeld) is included into the analytical stress tensor.

Tag: compute_forces(control.in)

Usage: compute_forces flag
Purpose: If .true., switches on the computation of forces.
flag is a logical string, either .true. or .false.

Default: .true. if a geometry optimization or molecular dynamics run was requested, or if the sc_accuracy_forces convergence criterion was set. Otherwise, .false.

This flag allows to request an explicit force computation for an otherwise explicit single-point calculation. This is necessary for use with external tools that require forces, such as a finite-difference calculation of vibrational frequencies (see Sec. 4.6) or a transition state search (see Sec. 4.8). In these cases, keyword final_forces_cleaned should also be set.

Tag: compute_numerical_stress(control.in)

Usage: compute_numerical_stress flag
Purpose: If .true., switches on the computation of the numerical stress tensor based on central finite differences.
flag is a logical string, either .true. or .false.

Default: .true. if a unit cell relaxation was requested and the computation of the analytical stress is not possible. Otherwise, .false.

If not further specified (by delta_numerical_stress) the default value for the scaling factor delta is set to $10^{-4}$ .

Tag: delta_numerical_stress(control.in)

Usage: delta_numerical_stress value
Purpose: Specifies the scaling factor delta in the computation of the numerical stress tensor (compute_numerical_stress).
value is a dimensionless real number $>0$ . Default: $10^{-4}$ .

Tag: distributed_hessian(control.in)

Usage: distributed_hessian flag
Purpose: If .false., each MPI task holds a complete copy of the Hessian matrix. If .true., the Hessian matrix is distributed across tasks.
flag is a logical string, either .true. or .false.. Default: .true. if both MPI and ScaLAPACK are available, .false. otherwise. If flag is .true., the Hessian matrix is written to the file hessian.aims (cf.hessian_file). If flag is .false., the Hessian matrix is written into geometry.in.next_step (cf.hessian_block).

This keyword is particularly useful when relaxing a large structure. Please note that it only works if FHI-aims is built with both MPI and ScaLAPACK. If init_hess reciprocal_lattice is found in control.in, distributed storage of the Hessian will be automatically turned off.

Tag: energy_tolerance(control.in)

Usage: energy_tolerance tolerance
Purpose: Sets the energy amount by which a relaxation step can move upwards and is still accepted.
tolerance is a small positive real number, in eV. Default: 1.5 $\cdot$ 10^-3 eV.

Small upward steps during relaxation may occur as a result of a slightly mis-guessed bfgs Hessian matrix somewhere along the path, or as a result of some residual numerical noise that leads to a discrepancy between energies and forces. In the present code version, such noise is always safely below the default energy_tolerance for reasonable settings. However, be sure to check that the total energy does not go up across several successive steps in a relaxation run. For the trm optimizer, also see harmonic_length_scale.

Tag: external_force(geometry.in)

Usage: external_force x y z
Purpose: Experimental – Applies an external force to the atom previous to this keyword.
x,y,z are the force components in eV/Å applied to the atom in $x$ , $y$ , $z$ direction.

When an external force is applied it is necessary to contstrain the relaxation of at least on other atom to avoid a constant shift of the geometry. Also the value has to be reasonably chosen. Tear apart geometries can result in very flat energy landscapes, which take a large amount of time to optimize. Typically this feature should be used in cases where a small external force, e.g. a STM-tip is applied on an atomic layer and the geomtry response of this external force is of interest.

This feature is experimental since no extensive testing was done for it.

Tag: external_pressure(control.in)

Usage: external_pressure value
Purpose: Applies external pressure to the unit cell during a unit cell relaxation. Setting pressure_unit is mandatory when using this keyword.
value is the pressure in units of pressure_unit applied to the unit cell.

In the periodic case, it is possible to apply hydrostatic pressure to the unit cell. To actually see the effect of the external pressure, a unit cell relaxation is required (see relax_geometry and relax_unit_cell). The crystal is then relaxed with the external pressure added to the stress tensor.

Tag: pressure_unit(control.in)

Usage: pressure_unit unit
Purpose: Sets the unit for the value of external_pressure.
unit: Presently supported options are GPa or eV_per_A^3 for units of GPa or eV/Å³, respectively.
Default: No default.

Tag: final_forces_cleaned(control.in)

Usage: final_forces_cleaned flag
Purpose: Decides whether spurious unitary transformations of the complete system (translations and rotations) are removed before the final output.
flag is a logical string, either .true. or .false. Default: .true.

This option affects directly the long-format (15 decimal) output of total energies and forces at the end of the s.c.f. cycle in the standard output file. If flag is .true., the final output forces are “cleaned” using the sayvetz [84, 270] mechanism of keyword clean_forces (removal of translations and rotations for cluster geometries; only translations removed for periodic systems).

final_forces_cleaned .true. should be set for use with external tools that require forces, such as a finite-difference calculation of vibrational frequencies (see Sec. 4.6) or a transition state search (see Sec. 4.8).

Tag: force_correction(control.in)

Usage: force_correction flag
Purpose: When computing the forces, determines whether or not to include the Hartree potential force correction term. Consideration of this keyword can be helpful to speed up processes such as geometry relaxations or a molecular dynamics run.
flag is a logical string, either .true. or .false. Default: .true.

As described elsewere[32, 60], omission of this term is one of the main reasons why the Hellmann-Feynman forces require a high level of self-consistency before their values can be trustworthy. In fact, this term is only equal to zero at full self-consistency.

Because this correction is only meaningful at a low level of self-consistency, for an appropriate use of this keyword, sc_accuracy_rho must also be set within a reasonable value, i.e., a too tight threshold for the density can lead to an insignificant correction to the forces. Therefore, unless specified by the user, FHI-aims lowers the default value of sc_accuracy_rho accordingly if a calculation involves force_correction.

Although force_correction and sc_accuracy_forces can be used together inside a particular calculation, its joint use is strongly discouraged due to the considerably high cost that each forces computation imply.

It should be noted that currently, in case of relax_unit_cell or any kind of analytical stress computation, as well as usage of ouput, force_correction becomes ineffective in determining the default value of sc_accuracy_rho.

Tag: harmonic_length_scale(control.in)

Usage: harmonic_length_scale length
Purpose: The trm/bfgs algorithm of relax_geometry judges a step by its energy gain. Usually, one simply uses the energy difference. For very short steps and rather light grids, however, it turns out that the qualitiy of the energy is inferior to the quality of the forces. For steps shorter than length, do not look at the energy but use the harmonic approximation $-\Delta\tilde{E}=(\bm{X}_{2}-\bm{X}_{1})\cdot(\bm{F}_{2}+\bm{F}_{1})/2$ as an estimate for the energy gain. If this procedure willfully accepts a step which increases the energy by more than energy_tolerance, the code stops to warn the user about the inconsistency between energy functional and forces.
length is a length scale in Å. Default: 0.025

Effectively, this flag switches from a real energy minimizer to a search for a stable zero of the force field for short step lengths.

Tag: hessian_to_restart_geometry(control.in)

Usage: hessian_to_restart_geometry flag
Purpose: Exports the current approximation to the Hessian matrix to geometry.in.next_step during a relaxation restart using hessian_block or hessian_file.
flag is a logical string, either .true. or .false. Default: .true.

Note: The geometry.in.next_step file is written out by default when the relax_geometry bfgs algorithm is used.

Tag: init_hess_lv_diag(control.in)

Usage: init_hess_lv_diag value
Purpose: In a geometry relaxation with a unit cell optimization, allows to specify the initial Hessian matrix elements used to estimate relaxation steps that are associated with the lattice vector degrees of freedom.
length is a length scale in eV/Å². Default: 25 ev/Å²

See the init_hess keyword for more information on the initial Hessian matrix used during a structure optimization.

Tag: init_hess(control.in)

Usage: init_hess type [value]
Purpose: Defines the initial Hessian matrix that is used by the bfgs structure relaxation algorithms (synonymous with trm) of the relax_geometry keyword.
type: Presently supported options are diag [value], Lindh [value].
Default: init_hess Lindh 2. unless a specific Hessian is given in the geometry.in file.

If no explicit initial Hessian matrix is given in the geometry.in file, the keyword init_hess can be used to specify the initial Hessian matrix used in a structure optimization.

With the option diag, a diagonal initial Hessian matrix is assumed. Then, the number given by the value option sets the diagonal elements between all atomic coordinates directly. The default is 25 eV/Å². Larger values lead to a more conservative start, smaller values lead to more aggressive initial relaxation steps.

With Lindh the Lindh model matrix [199] is used to initialize the Hessian between all atomic coordinates (usually a very efficient guess). For stability reasons, add-value (in eV/Å², defaults to 2.0 eV/Å²) is added to all matrix elements on the diagonal. The parameter thres (default: 15.0) can be used to specify the accuracy of the Lindh matrix; only terms estimated to be larger than $e^{-\text{{thres}}}$ are taken into account. If you are planning a large number of complex unit cell optimizations of a similar type, we do recommend to check whether this default value is any good.

If a unit cell relaxation is additionally requested, the Hessian for the lattice degrees of freedom is set to be the square of the reciprocal lattice matrix and normalized to 25 eV/Å². This mimics a diagonal initial Hessian if strain coordinates for the representation of the lattice would be used. [251]

If the initial Hessian is specified explicitly by hessian_block or hessian_file in geometry.in, this explicit hessian overrides any information requested by the init_hess keyword.

Tag: max_atomic_move(control.in)

Usage: max_atomic_move value
Purpose: Maximum allowed step length taken during relaxation.
value is a real positive upper bound for the maximum allowed change in single atomic coordinate, in Å. Default: 0.2 Å.

If the bfgs-predicted change in an atomic coordinate exceeds max_atomic_move, the length of the entire step (all coordinates) will be scaled down to not exceed the maximum allowed displacement.

Tag: max_relaxation_steps(control.in)

Usage: max_relaxation_steps number
Purpose: A structure optimization will be aborted after exceeding a prescribed maximum number of steps.
number is the prescribed maximum step number. Default: 1000 .

Tag: min_trust_radius(control.in)

Usage: min_trust_radius value
Purpose: The minimum trust radius value that the relaxation algorithm is allowed to reach before aborting the structure optimization since no path forward can be found.
value is the allowed minimum trust radius in Å. Default: 10^-4 Å.

In FHI-aims’ standard relaxation algorithm, relax_geometry trm or bfgs, the trust_radius is the limit up to which a geometry optimization step is trusted and therefore allowed. The trust_radius fluctuates during execution. It can be reduced or it can increase, depending on the success or failure of consecutive geometry optimization steps. Notably, if the relaxation algorithm encounters a scenario in which it cannot find an admissible step, the trust_radius could, in principle decrease to arbitrarily low values, never making progress while spending significant computational resources.

The min_trust_radius value sets the lower bound that is considered acceptable for the trust_radius. If the trust_radius decreases below this bound, the relaxation algorithm concludes that there does not appear to be a sensible pathway forward and stops. This can happen, for example, if the calculated energy gradients (“forces”) are inconsistent with the actual energy surface, i.e., if the forces point in a direction in which the actual energy goes uphill. Several possible strategies exist to escape this scenario – typically, denser integration grids or, perhaps, a tighter scf convergence criterion sc_accuracy_rho. However, it is also possible that a particular total energy method (e.g., a new development) does not yet have entirely consistent forces. Obviously, we’d like to know about such cases. In that case, contacting the developers would be helpful and appreciated.

Tag: numerical_stress_save_scf(control.in)

Usage: numerical_stress_save_scf flag
Purpose: Controls if constrain_relaxation directives are used to determine implicitly if a component of the numerical stress has to be calculated. This greatly accelerates unit cells with high symmetries (e.g. orthorhombic).
flag is a logical string, either .true. or .false. Default: .true.

Tag: orthonormalize_eigenvectors(control.in)

Usage: orthonormalize_eigenvectors flag
Purpose: Specifies whether or not the wave function coefficients from the previous geometry will be re-orthonormalized before initializing a new relaxation step.
flag is a logical string, either .true. or .false. Default: .true.

The orthonormalize_eigenvectors keyword allows to reorthnormalize the converged self-consistent Kohn-Sham orbitals $c_{jl}$ after a relaxation step. These are then used to reinitialize the electron density for the next relaxation step.

Due to the change in atomic positions, the wave function coefficients $c_{jl}$ for the earlier geometry are no longer orthonormal after the relaxation step. The consequence is an initial electron density which no longer satisfies the correct electron count (i.e., the system may appear to be charged immediately after a relaxation step, although a neutral system was requested). In principle, this does not matter for the outcome of a calculation, since the self-consistent solution will be independent of the starting point. In many cases, there is no clear benefit in terms of the s.c.f. convergence duration from orthonormalizing the $c_{jl}$ prior to the reinitialization; however, some cases with unstable s.c.f. convergence may benefit significantly.

Tag: relax_geometry(control.in)

Usage: relax_geometry type tolerance
Purpose: Specifies if a structure optimization (geometry relaxation) is requested, and which.
type specifies the type of requested structure optimization. Default: none.
tolerance: Specifies the maximum residual force component per atom (in eV/Å) below which the geometry relaxation is considered converged.

Finds the nearest minimum of the Born-Oppenheimer potential energy surface for the nuclei.

For periodic calculations: If you are looking to relax not just atomic coordinates but also the unit cell shape (lattice vectors), you do need to specify an additional keyword: relax_unit_cell.

The presently supported options for type are none and trm (synonymous with bfgs), as well as trm_2012 for reference with older FHI-aims versions.

•

bfgs or trm is the recommended default. It uses a trust radius method enhanced version of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization algorithm (see Ref. [231], which was the basis for Jürgen Wieferink’s code effort in this area). In our tests, this version appears to give the fastest convergence reliably.

As of December 2018, it additionally implements preconditioning of the lattice-lattice Hessian as explained in init_hess and preserves fractional coordinates of the atomic positions when predicting a new lattice. This mimics the optimization in strain coordinates. [251]
•

trm_2012 Implements the former trm method without effective strain coordinates.

A reliable force convergence criterion tolerance for most structures is 10^-2 eV/Å or 5 $\cdot$ 10^-3 eV/Å. Do not use significantly smaller values unless you have a specific reason. Smaller values may cost much computer time for essentially no further measurable total energy minimization.

Going to a much smaller tolerance value may only be useful for some very specific purposes, for example, high-accuracy finite difference calculations for vibrational properties. In other scenarios, if tolerance is set to a too small value by default, 80% or more of your CPU time may be spent groping around in the last meV of the structure optimization.

If tighter settings of the tolerance parameter are used, do not forget that tighter s.c.f. convergence accuracy settings may also be required to get accurate gradients in the first place. Ideally, use the sc_accuracy_rho keyword for this purpose, not sc_accuracy_forces or sc_accuracy_stress (see below).

In other words, use the tolerance criterion for a structure relaxation run wisely – decide what is the physical quantity you are actually interested in, and then check which value of the tolerance criterion is safe but still efficient.

The relaxation algorithm can be greatly sped up by using a somewhat intelligent guess for the Hessian matrix used in the initial step. By default, FHI-aims now sets the general purpose model matrix due to Lindh and coworkers [199] with a slight modification. If, for some reason, a particular initial geometry does not appear to play well with the Lindh Hessian, a simpler, slower, but more resilient diagonal approximation to the initial Hessian matrix can also be used. For more information see the init_hess keyword above.

The energy_tolerance and harmonic_length_scale keywords can be set to more forgiving values if the bfgs algorithm decides to abort relaxations because of a presumed deviation between the predicted and the real energy landscape.

Another important warning: Evaluating the forces and the stress tensor is much more expensive than a normal iteration during s.c.f. convergence. The current default behavior of FHI-aims avoids any double computations of forces and stress tensors, relying instead on a sufficiently tight convergence criterion sc_accuracy_rho to determine s.c.f. convergence and only then calculating forces and stresses once per geometry step.

While one can in principle check the s.c.f convergence of forces / stresses explicitly, the cost of multiple evaluations of forces / stresses for a single geometry can be very high. Therefore, we recommend to never use the keywords sc_accuracy_forces or sc_accuracy_stress in a control.in file unless there is a specific need. Do not set these keywords routinely.

Tag: relax_unit_cell(control.in)

Usage: relax_unit_cell type
Purpose: Relaxes unit cell (lattice vectors) using the structure optimization as specified in relax_geometry.
type specifies the type of requested unit cell optimization. Presently supported options: none, full, fixed_angles, fixed_volume, slab. Default: none

Allows to optimize the lattice vectors of a periodic calculation, in addition to the normal atomic coordinates inside the unit cell. This keyword is not on by default, as automatically optimizing the unit cell of (say) a surface calculation could do a lot of unintended harm. Possible settings:

•

none : Unit cell will be kept fixed, no optimization.
•

full : All lattice_vector degrees of freedom will be relaxed, except those affectd by explicit constraints.
•

fixed_angles : All angles between lattice vectors will be constrained (kept fixed), only the lengths of each lattice vector are varied. (This option used to be called shape, but that is a misunderstandable term. The shape term will be removed in the future to avoid confusion.)
•

fixed_volume : All angles and lengths are varied but the cell volume is kept fixed.
•
slab : This is for slab structures. Specifically, thelattice_vector aligned with the vacuum direction will be fixed usingconstrain_relaxation. To illustrate, in the case of a z-oriented slab structure, the relaxation constraint will be set in the following manner:
```
    lattice_vector 3.  0.  0.
      constrain_relaxation z
    lattice_vector 3.  6.  0.
      constrain_relaxation z
    lattice_vector 0.  0.  23.
      constrain_relaxation .true.
    
```

This keyword should be used only together with relax_geometry. Individual lattice vectors or its components can be constrained by using constrain_relaxation.

If the computation of the analytical stress is possible regarding the chosen computational settings, the analytical stress is used for the unit cell relaxation. Otherwise, the numerical stress is used. With stress_for_relaxation, one can explicitly choose either numerical or analytical stress for the unit cell relaxation.

If a unit cell relaxation produces strange results with the analytical stress, here are some potential remedies:

1.

A very possible reason may be because the integration grids are not dense enough. This could especially well be the case for “light” settings. One remedy is to just use the integration grids from “tight” and the basis functions from “light”.
2.

Another possible remedy is to switch the way the integration weights are calculated to a slightly slower, non-default version. E.g., change the partition_type to a spherical one like rho_r2.
3.

Finally, you may wish to set the convergence of the analytical stress with sc_accuracy_stress to an explicit, final value. Only ever set this keyword for test purposes, though, not routinely. The cost for too many analytical stress evaluations can be disproportionately large.

Tag: stress_for_relaxation(control.in)

Usage: stress_for_relaxation type
Purpose: Use either numerical or analytical stress for unit cell relaxations.
type can be either numerical or analytical. Default: Chosen automatically based on computational settings.

To perform an actual unit cell relaxation, one has to set relax_unit_cell. If one chooses analytical but the computation of the analytical stress is not possible, FHI-aims will abort.

Tag: write_restart_geometry(control.in)

Usage: write_restart_geometry flag
Purpose: During a structure optimization, exports the current geometrys and approximation to the Hessian matrix to a file geometry.in.next_step.
flag is a logical string, either .true. or .false. Default: .true.

Note: The geometry.in.next_step file is written out by default when the relax_geometry bfgs algorithm is used.