Basis optimization¶

Description¶

The AiiDA-siesta package offers three workchains to help users in selecting the optimal basis set for a given system:

1) The SimplexBasisOptimization that finds the minimum of a quantity (typically the basis enthalpy) varying a set of input variables (typically cutoff radii of orbitals) using the Nelder–Mead (simplex / amoeba) method. 2) The TwoStepsBasisOpt that performs a two level optimization, running simplex iterations followed by periodic restarts with new simplex hyper-tetrahedra of progressively smaller sizes. This replicates more closely the optimization util of the siesta distribution. 3) The BasisOptimizationWorkChain that performs a full optimization testing first basis cardinality and then applying the SimplexBasisOptimization to optimize all the radius of the orbitals.

The simplex optimization is performed taking advantage of the aiida-optimize package <https://aiida-optimize.readthedocs.io/en/latest/index.html>_, in particular the Nelder–Mead engine implemented in that package.

Supported Siesta versions¶

At least 4.0.1 of the 4.0 series, 4.1-b3 of the 4.1 series and the MaX-1.0 release, which can be found in the development platform (https://gitlab.com/siesta-project/siesta). For more up to date info on compatibility, please check the wiki.

SimplexBasisOptimization¶

The Nelder–Mead method (commonly known as simplex or amoeba method) is a numerical method to find the minimum of a function with N variables. A simplex is a special polytope of N+1 vertices in N dimensions (for instance a triangle in 2D, a tetrahedron in 3D and so forth). The Nelder–Mead methods in N dimensions starts from a set of N+1 test points arranged as a simplex. The value of the function is calculated at each test point and these values are then used in order to find a new test point and to replace one of the old test points with the new one in case it returns a smaller value for the function under investigation. Repeating the procedure, the technique progresses until all the N+1 test points produce values that are all within a threshold. When this happens the minimum has been reached. In the context of the basis optimization the function is usually the basis enthalpy (but also other quantities are supported) and the variable are the parameters defining the basis (typically cutoff radii of the basis orbitals).

An example of the use of SimplexBasisOptimization is in /aiida_siesta/examples/workflows/example_simplex.py

Inputs¶

Inputs are organized in two namespaces and are described in the following:

siesta_base, input namespace, Mandatory

Accepts all the inputs of a SiestaBaseWorkChain (listed here <siesta-base-wc-inputs>) with the only mandatory modification to include in the “basis” input some variables to optimize. The variable must be defined using a dollar and a string. An example:
```
basis = Dict(dict={
 '%block pao-basis': "\nSi   2\n n=3   0   2\n 4.99376      $sz2 \n n=3   1   2 P 1\n 6.2538      $pz2 \n%endblock pao-basis"
 })
```
An upper and lower value must be set for each variable and optionally one or more starting points (see next point in this list). Please note that variables are typically defined for orbitals radii in the pao-basis block, but one can also create variables for “higher level” keywords like the energy-shift or split-norm.

simplex, input namespace, Mandatory

Here all the inputs for the simplex method can be defined. They are listed in the next lines.

simplex.variables_dict class Dict, Mandatory

A dictionary containing all the info about the variables that are modified in order to find the minimum basis enthalpy. An example related to the basis block above:
```
variables_dict = Dict(dict={
    "sz2":[2.0,4.8,3.0],
    "pz2":[2.0,6.0,3.0]
    })
```
The variables names must be defined here as keys of the dictionary and must correspond to the strings defined in the basis input, but removing the dollar symbol. The list associated to each string defines in this order: 1) The lower limit for the variable, 2) The upper limit, 3) the starting value to construct the simplex hyper-tetrahedron. The up and down limit of the variables are used in such way: if the algorithm attempts the calculation of the function out of range, a huge value for the function is returned. The starting value is going to be the point from which the simplex hyper-tetrahedron is constructed. In particular, the first test point is directly formed by the specified starting points (in the example above is [3.0,3.0]). The other N test points are obtained substituing one component with num + range * simplex_inps.initial_step_fraction, where num is the defined starting point, range is the upper - lower limit and simplex_inps.initial_step_fraction is a number between 0 and 1 defined in the next point of this list. Supposing simplex_inps.initial_step_fraction = 0.2, in out example, the other two test points are [3.0,3.8] and [3.56,3.0].

When 3) is not defined, it is chosen randomly between the boundaries, but it is always suggested to set it since it will be used to construct the Alternatively to 3), N+1 values can be entered and this would correspond to define explicitly all the components of all the simplex initial points.

simplex.initial_step_fraction class Float, Optional

A fractional increment to be used in the construction of the starting simplex hyper-tetrahedron. See point above for more details. Default at Float(0.4). It is ignored if all the components af all the test points are set in the point above.

simplex.max_iters class Int, Optional

The maximum iterations for the Nelder–Mead algorithm. Please note that an iteration step usually involves more then one new test point. So the points tested at the end will be way more than the max_iters. Once the simplex.max_iters is reached, the workchain stops returning the best simplex so far, even if the threshold convergence has not been reached. Default is Int(40).

simplex.output_name class Str, Optional

The name of the output that needs to be minimized. In principle all the numerical values returned in the “output_parameters” of a SiestaBaseWorkChain are accepted, but typically the “basis_entalpy” or the “harris_energy” are of interest. Defalut is Str("basis_entalpy")

simplex.tolerance_function class Float, Optional

The tolerance accepted to define the optimization converged. If the values of the functions for all points in the simplex are all within the simplex.tolerance_function, the optimization is considered concluded. The default is Float(0.01). Please note that the choice of this parameter must be related to the variance of the output function, therefore the default might be unreasonable for your application. In the future an extension implementing a fractional tolerance will be provided.

Outputs¶

The following outputs are returned:

last_simplex class List

The output containing the values of the last simplex. Always returned, even if the optimization does not reached the required tolerance. It is a list of lists. The first element of the list is always the best choice of the parameters obtained by the optimization so far.

optimal_process_input class List

This output contains the optimal set of parameters obtained after optimization. This corresponds to the first entry of the list return by the last_simplex, however it is returned only if the optimization succeed.

optimal_process_output class Float

The value of the function for the optimal set of parameters obtained with the optimization. Returned only if the optimization succeed.

optimal_process_uuid class List

The uuid of the SiestaBaseWorkChain that has the optimal_process_input as variables and that returned the optimal_process_output. Returned only if the optimization succeed.

It is important to note that the optimization is entirely an AiiDA process, therefore the provenance of all calculation called is preserved. We can have a look at the attempted variables values and the obtained basis entalpy in this simple way. In the verdi shell:

node=load_node(<PK>)  #PK of your SimplexBasisOptimization
for wc in node.called[0].called:
     print(wc.inputs.the_values.get_list(),wc.outputs.ene.value)

And many more info can be extracted from the inputs and outputs of each run wc. These wc are SiestaBaseWorkChain wrapped into a thin layer that attach to each calculation the information needed by the optimizer.

TwoStepsBasisOpt¶

This workchain uses the SimplexBasisOptimization, but it adds a step in the optimization, which consists in restarting the simplex with a subsequently smaller simplex.initial_step_fraction. This is implemented in the original simplex optimization code that can be found in the Util of the SIESTA package. There the fractional step is called “lambda” and we will follow the same notation here.

Inputs¶

All the inputs of SimplexBasisOptimization are inputs of this workchain except the simplex.initial_step_fraction. This include the way to specify the optimization variables in the siesta_base.basis input. This workchain adds a further called macrostep. This allows:

macrostep.initial_lambda class Float

The value of lambda to be used as simplex.initial_step_fraction in the first iteration. Default Float(0.4),

macrostep.lambda_scaling_factor class Float

The rate at which lambda decreases between from a macrostep to the other. Default Float(0.5)

macrostep.minimum_lambda class Float

When this value for lambda is reached, the macrostep iteration stops. Default Float(0.01).

Outputs¶

Same outputs of SimplexBasisOptimization.

BasisOptimizationWorkChain¶

This workchain manages entirely the optimization of the basis sets for a SIESTA calculation. It first run calculations with different basis sizes (using the “PAO-BasisSize” option of SIESTA) and gets the size that gives minimum of the monitored quantity (e.g. basis enthalpy).

NOTE: This does not include yet the possibility to test different basis sizes for different species.

It then allow to add extra orbitals to the calculation manually and see if this leads to a further decrease in the monitored quantity.

Then automatically sets up a SimplexBasisOptimization according to an optimization schema defined by the user.

Inputs¶

All the inputs of SimplexBasisOptimization are inputs of this workchain except the simplex.variables_dict. Please note that whathever is specified in siesta_base.basis will be copied in every calculation. So we prevwnt in this keyword to set the basis bloch or the basis sizes since the alghoritm will take care of it. In siesta_base.basis can put keywords like the “pao-non-perturbative-polarization-schema” or choices on the pseudopotential grid.

Few more inputs are allowed:

basis_sizes class List Optional

The list of basis sizes to try out. Default List(list=["DZ", "DZP", "TZ"]).

add_orbital class List Optional

A dict of lists, the key of the dict must be the name of an element of the periodic table, the list must list the orbitals to add at that atom, for instance:
```
add_orbital = Dict(dict={
  "Ca":["3d1","4f1"],
  "O" :["4f2"]
  })
```
This would add a f orbital with two zetas for O and a d and f orbital to Ca (one zeta each). As already specified, the presence of this input implies an extra step between the check of basis cardinality and the actual SimplexBasisOptimization.

sizes_monitored_quantity List Optional

The quantity to monitor in the check of the cardinality. If not specified is going to be the same specified in simplex.output_name.

optimization_schema.global_energy_shift List

If set to true, the energy shift and the pao-split-norm are used as optimization variables, not the explicit radius of the basis block. Default is False

optimization_schema.global_split_norm List

If set to true, the pao-split-norm is optimized as a global variable. Please note that this can be used in conjunction with global_energy_shift in order to optimize only global variables and not the pao block, but it can be also used alone to set that the first zeta radii of the orbitals are optimized, but the second zetas no! If optimization_schema.global_split_norm is True and optimization_schema.global_energy_shift is False the basis block is created putting all the second and further zetas to zero and the globas pao-split-norm is a variable for optimization. Default False.

optimization_schema.charge_confinement List

If set to true, the empty orbitals will receive a charge confinement and the charge of the confinement is a variable for optimization. Default False

To conclude, the inputs allow to do various type of optimizations. As default all the radia are optimized, but this can be modified using the optimization_schema keywords

Outputs¶

Only one output is produced:

optimal_basis_block class Dict

Returning the optimal pao block, meaning the one that gives the minimum of the monitored quantity.

Protocol system¶

The protocol system is not directly available for this WorkChain. However inputs of the SiestaBaseWorkChain can be obtained in a dictionary in this way:

inp_gen = SiestaBaseWorkChain.inputs_generator()
inputs = inp_gen.get_inputs_dict(structure, calc_engines, protocols)

The inputs of get_inputs_dict are explained in the protocols documentation. Then the user can place these inputs in the siesta_base namespace.