4.8 Finding Transition States: the aimsChain

The chain control file is the file that governs the evaluation of MEP. It follows the same configuration convention set by the aims control.in, namely a list of key value pairs. The control file should be named chain.in.

This is a sample control file.

All the available keywords for the control are listed below:

Tag: run_aims (chain.in)
        Usage: run_aims your command Purpose: Provide the command that is used to call aims in your environment. Default value: mpiexec -ppn 8 -n $NSLOTS ~/bin/aims.081912.scalapack.mpi.x

You will need to change the default value for whatever command you are using.

There is no need to quote your command, simply type it as you would in a bash shell. Everything after the keyword are concatenated into a single string.

Tag: initial_file (chain.in)
        Usage: initial_file filename Purpose: Provide the initial geometry file for the path. Default value: ini.in

The format should be the same as the aims geometry.in. Flags such as constraints are interpreted, and will affect the evaluation of MEP. All the intermediate images will have the same constraints as the initial geometry.

Atoms in the initial and final geometry must establish a one to one correspondence. The $i$th line in the initial geometry and the final geometry must correspond to the same atom. Often there are many possible combinations, it is the user’s responsibility to use the most physical one and/or check different combinations. You can perform a test interpolation with the interpchain.py tool provided.

Tag: final_file (chain.in)
        Usage: final_file filename Purpose: The final geometry for the path. Refer to initial_file . Default value: fin.in

Tag: n_images (chain.in)
        Usage: n_images value Purpose: The number of images used for MEP calculation. Default value: 5

Suitable value for this variable is highly dependent on the desired degree of resolution and the complexity of the system. A simple single barrier reaction will require 5 images to achieve a good result, while more complex potential energy surfaces require roughly 3 or more images per hill/trough.

If the geometry is not computationally expensive(tens of seconds per calculation with a reasonable # of CPUs), the simple principle of “the more the merrier” can always be applied.

Tag: external_geometry (chain.in)
        Usage: external_geometry file_name Purpose: If this variable is set, the initial set of images will be obtained from a external source instead of a linear interpolation. Default value: None

This flag should be set to the filename of a text file listing all of the intermediate geometries in order, e.g.

./geo/1.in
./geo/2.in
./geo/3.in

It is advised that if you have anything better than a direct linear interpolation, you should use it as the initial guess. Linear interpolation can be highly inefficient, and would result in cases which never converges. When in doubt, you can always use the growing string method.

Tag: periodic_interpolation (chain.in)
        Usage: periodic_interpolation true/false Purpose: Whether or not interpolation is done while considering periodic boundary. Default value: False

This setting affects periodic systems only. If set to true, each atom’s final geometry will be checked against all periodic counterparts to find the shortest path between the initial and final geometry. This can resolve some problems caused by the periodic boundary condition. Consult sample 4 for detailed usage. Always use the provided interpchain.py tool to check the resulting interpolation and confirm that it’s physical before running.

Tag: resample (chain.in)
        Usage: resample true/false Purpose: Resample the input path into arbitrary number of images. Default value: False

This tag is effective only when external_geometry is set. When resample is set to true, the package will resample the path you provided and interpolate it to have # of images equal to n_images . This can be very useful if you find your calculation to have an inadequate number of images. You can extract the geometries from the final iteration, and start a new job by resampling them.

Tag: aims_restart (chain.in)
        Usage: aims_restart value Purpose: Reuse the wave function restart file to speed up the calculation. This should set to be the same as restart in control.in. WARNING! This function is not compatible with aims as of July 2013, and only works occasionally. Default value: None

Value should be the same as restart in your control.in. The wave function restart file will be copied from iteration to the next, so that the scf cycles can start from an existing configuration. When working correctly, this can cut the computation time required by half or more.

However, at this stage restart has very limited check on the input restart file, resulting in an error when the wave functions don’t match. Do not use this keyword unless you are sure that your version of aims will start a new scf initialization instead of reporting errors when the wrong restart file is provided. If certain of which, however, the keyword should always be used to speed up the calculation.

Setting restart will consume a decent amount of disk space. You may want to remove the restart files if you are planning to store the calculation permanently.

Tag: restart (chain.in)
        Usage: restart true/false Purpose: Whether the package will allow restarts. When set to true, you can simply submit the job again to continue from the previous calculations. Default value: true

When set to true, a restart file will be written as the path iterates. When the job is submitted/resubmitted, it will first check for the restart file to see if it’s possible to start from a previous calculation. This is useful, when, for example, your job is killed due to time limitations on the cluster. Please refrain from changing settings between restarts, which is always error prone. If you are planning to change some settings, it’s always safer to extract the most recent geometries and use them as the initial path for a new calculation.

There is still the risk that the job is killed while the restart file is been written, in which case anything could happen when the job is resumed. However, considering the rarity of such event, it should not be a problem in practice.

Tag: method (chain.in)
        Usage: method neb/string Purpose: Pick the method to use for the MEP calculation. Default value: string

Currently only string method and NEB are supported. In our testing, string method have shown to be slightly more stable than NEB, and therefore is our recommended method.

Tag: neb_spring_constant (chain.in)
        Usage: neb_spring_constant value Purpose: This sets the spring constant used in NEB calculation, and has no effect if method is string. Default value: 20.0

It should be set to have the same magnitude as the force felt by the geometry, which is hard to determine before hand. However a generic value of 20 is good in general. Please try method string first if you suspect that the spring constant have to be changed for the system to converge.

Tag: force_thres (chain.in)
        Usage: force_thres value Purpose: The threshold for convergence. Default value: 0.2

Optimization will stop when the residual forces in the system is smaller than this preset value. If use_climb is true, climbing-image calculation will start after this threshold has reached. $0.2eV$ is always a good starting point, and for more rigorous calculations $0.1eV$ or $0.05eV$ can be used. Do not set it too small, since some numerical noises in the forces will always exist even when the path is well converged.

Also note that adequate value is dependent on the system. A reaction with a $3eV$ barrier can be much more well converged at $0.05eV$ threshold than a reaction with a $0.3eV$ barrier at the same threshold, because the former, very often, has a steeper barrier and hence larger force by nature.

Tag: optimizer (chain.in)
        Usage: optimizer dampedBFGS/BFGS/LBFGS/trm/CG/FIRE Purpose: This picks the optimizer used for optimization. Default value: dampedBFGS

BFGS is a textbook implementation of BFGS optimizer, which was observed to have wild behaviour is some situations. The fact that BFGS series of optimizer can be used for this type of calculations is in fact purely coincidental. The series of force projections involved in the chain of states methods can severely hamper the effectiveness of quasi-Newton optimizers. Please resort to FIRE optimizer whenever you observe wild behaviour in the optimization process.

dampedBFGS damps the original BFGS optimizer using several techniques, and is slightly more stable in those situations. (But they still get stuck from time to time!) This is the default setting that should always be tried first.

LBFGS is the limited memory version of BFGS, which uses less memory for extremely large systems. However, given that in general FHI-aims is applied to systems with hundreds of atoms at most, this optimizer has little advantage in this front. However, due the fact that LBFGS uses only recent iterations for approximation, the results can be quite different from BFGS in a some PES. Therefore it’s worth trying if BFGS fails.

trm is the same trust-radius method ported from FHI-aims. In our tests it has proven to be reasonably stable and efficient. We would recommend this as another alternative along with FIRE when dampedBFGS fails.

CG is the textbook implementation of conjugate-gradient algorithm using finite difference scheme. In our tests it has not shown any advantage over other algorithms, ans is provided for the sake of completeness.

FIRE is the Fast Inertial Relaxation Engine, one of the better non-Newton type optimizers. It is slower than BFGS series of optimizers in general, but is immune to the aforementioned instability because it does not approximate the Hessian. If dampedBFGS fails, try FIRE with global_optimizer off as an alternative.

Tag: global_optimizer (chain.in)
        Usage: global_optimizer true/false Purpose: Whether all images are optimized as a single object, or if individual images are optimized as separated object. Default value: true

The global version, when coupled with BFGS-series of optimizer, is often faster than the non-global optimizer. However, this combination is less stable, and can lead to wild results. The non-global optimizer coupled with FIRE seems to be the most stable one in our tests, but is much slower than the former. Non-global optimizer combined with BFGS, on the other hand, slows the calculation significantly due to overestimation, and should be avoided.

Tag: xyz_lattice (chain.in)
        Usage: xyz_lattice a b c Purpose: This key only affects the .xyz file written in paths directory. It governs how the lattice is repeated in each lattice vector. The default is a standard 2 2 1 setting, valid for most surfaces. It has no effect on clusters. Default value: 2 2 1

Tag: map_unit_cell (chain.in)
        Usage: map_unit_cell true/false Purpose: This key only affects the .in file written in paths directory. It determines whether atoms in these geometries are mapped back to the central unit cell. It has no effect on clusters Default value: false

Tag: use_climb (chain.in)
        Usage: use_climb true/false Purpose: This will turn on the climbing image[134] feature. Default value: false

The points with highest energies will move toward a higher energy location along the path until the saddle point is reached. Please consult climb_mode for a detailed explanation of this process. If you are only interested in finding the energy barrier, it’s possible to set force_thres to some larger value, such as 0.2, and set climb_thres smaller(e.g. 0.05) for a reasonably tight convergence. In this case decreasing force_thres will not affect the accuracy of climbing image, so long as the climbing image converges. Transition state finding can be sped up significantly this way by reducing the number of single point calculations required.

Tag: climb_thres (chain.in)
        Usage: climb_thres value Purpose: Set the convergence criterion for the climbing image. Default value: same as force_thres

The climbing image will stop when the residual forces in the system is smaller than this value. It is normally, although not necessarily, set to a value smaller or equal to force_thres . In simple systems, the climbing image is capable of pushing the system done to $meV/{\AA }^2$ range, but that’s not as plausible in larger systems.

Tag: interpolated_climb (chain.in)
        Usage: interpolated_climb true/false Purpose: Determine whether the climbing image is chosen from one of the existing images or interpolated from known energies and geometries. Default value: true

When set to true, the current energies and geometries will be fitted with a cubic spline to identify the point with highest energy. If this geometry is sufficiently far from existing nodes, then the interpolated geometry is used. If the geometry is close to one of the existing images, than that image is set to be the climbing image.

When set to false, the image with highest energy is set to be the climbing image.

Tag: climb_mode (chain.in)
        Usage: climb_mode 1/2/3 Purpose: This determines the “mode” of the climbing image. Default value: 2

Increasing mode corresponds to increasing stability and decreasing efficiency. The main difference lies in the number of images that are allowed to move during the climbing process.

climb_mode 1: Only the image with highest energy is allowed to move. This setting can manage many circumstances, provided that force_thres is small enough so that the converged path provides a stable basin.

climb_mode 2: The image with highest energy and its two neighbouring images are allowed to move. The two additional image provides an evolving tangent estimate as the central image climbs. This mode can tackle nearly all cases where climb_mode 1 fails under the same force_thres setting. It is safe, in general, to set force_thres to 0.3 if you are only interested in energy barrier.

climb_mode 3: All of the images excluding the end points are allowed to move. This is the most stable, but rarely necessary setting. The efficiency drop is highly dependent on the system and the force_thres setting. In general this should be used as a last resort when force_thres can be converged to a tight setting but climbing image on other modes fails.

Tag: climb_global_optimizer (chain.in)
        Usage: climb_global_optimizer true/false Purpose: The same as global_optimizer , except this keyword is dedicated for climbing image. Default value: true

The tag has no effect when climb_mode is 1, in which case the global and non-global optimizer are equivalent. Unlike global_optimizer , BFGS+Non-global setting can be a good combination for climb_mode 2 if the default setting fails, mainly due to the well behaving local basin for a roughly converged string.

Tag: climb_optimizer (chain.in)
        Usage: climb_optimizer dampedBFGS/BFGS/LBFGS/CG/TRM/FIRE Purpose: The same as optimizer , but dedicated for climbing image. Default value: fire

Tag: climb_control (chain.in)
        Usage: climb_control file_name Purpose: Specify a different control.in for the climbing image. Default value: control.in

It’s possible, for example, to converge the entire path with a light setting, and then converge the climbing image with tight setting, which consumes far less computational power.

However, there are a few things to note if this feature is utilized. First, force_thres must be set smaller than normal, perhaps even the same value as climb_thres . Since a different control file can produce a radically different PES, you would like to be as accurate as possible when converging the string, so that the error is not amplified with a tighter setting. Secondly, it’s better to use climb_mode 2 or 3 for these kind of computations, since the single image climbing mode can rarely recover from a bad starting point.

This feature is best suitable for production works where a tighter setting which will dramatically increase the time required for single point calculations.

As an alternative, you can also kill your running process and change the control.in file manually before restarting the process.

Tag: use_gs_method (chain.in)
        Usage: use_gs_method true/false Purpose: This controls whether the growing string method is used or not. Default value: False

The growing string method(GSM) is explained in sample 6. It is best suitable for cases where it is certain that linear interpolations between initial and final images will not lead to a correct path. This can be caused by movements such as rotations. It is also useful when doing large-scale automated scanning, where the user does not have the time to look at geometries on a case by case basis.

When set to true, the path will start from the two end point and slowly grow inward to generate a physical path. When the path is completely grown, it will be passed onto NEB/SM for further calculations.

Beware that when direct interpolation can lead to correct paths, e.g. surface dispersion, using growing string method may reduce efficiency.

Tag: gs_thres (chain.in)
        Usage: gs_thres true/false Purpose: The GSM counterpart for force_thres . Default value: force_thres *1.5

When the forces in the system goes below this preset value, a new node will be added to the path. It should not be set too small, since the purpose of GSM is to generate a physical path, not a well-converged path. The default value gives a good guideline for this key.

Tag: gs_n_images (chain.in)
        Usage: gs_n_images value Purpose: The GSM counterpart for n_images . Default value: n_images

In some circumstances, you may want to specify a different number of images for the growing stage of the calculation. This key serves this purpose. After the growing process, the path will be re-sampled to match the value of n_images

Tag: gs_optimizer (chain.in)
        Usage: gs_optimizer dampedBFGS/BFGS/LBFGS/trm/CG/FIRE Purpose: The GSM counterpart for optimizer . Default value: trm

For GSM, dampedBFGS, trm, and FIRE are optimizers worth trying.

Tag: gs_global_optimizer (chain.in)
        Usage: gs_global_optimizer true/false Purpose: The GSM counterpart for global_optimizer . Default value: false

We have not done enough testings to determine conclusively the better set of optimizers. The default provided here have performed well in our benchmarks, but feel free to try other combinations.

Tag: lbfgs_alpha (chain.in)
        Usage: lbfgs_alpha value Purpose: Set the curvature used to initialize the Hessian(in $eV/{\AA }^2$) in LBFGS. In general any value that is not magnitudes off are acceptable. Default value: 120.0

Tag: lbfgs_memory (chain.in)
        Usage: lbfgs_memory value Purpose: Set the number of past iteration that LBFGS is going to remember. A larger value will increase memory consumption, and a small value will decrease accuracy. Default value: 30

Tag: lbfgs_maxstep (chain.in)
        Usage: lbfgs_maxstep value Purpose: The maximum step in $\AA$ that an atom can take in a single iteration. This is similar to max_atomic_move, but defaulted to a much smaller value due to the increasing complexity. Default value: 0.04

Tag: bfgs_alpha (chain.in)
        Usage: bfgs_alpha value Purpose: Same as lbfgs_alpha , but used by BFGS, trm, and dampedBFGS optimizers. Default value: 120

Tag: bfgs_maxstep (chain.in)
        Usage: bfgs_maxstep value Purpose: Same as lbfgs_maxstep , but used for BFGS and dampedBFGS optimizers. Default value: 0.04

Tag: fire_dt (chain.in)
        Usage: fire_dt value Purpose: The initial time step used by the FIRE optimizer. Default value: 0.02

It is set to a very conservative value because FIRE is intended to be used as a fall back when BFGS fails. Values up to 0.1 can be used to speed up the calculation, provided that the PES is smooth enough. Internally, the time step is dynamically adjusted, and this key only serves to initialize the value.

Tag: fire_maxstep (chain.in)
        Usage: fire_maxstep value Purpose: Same as lbfgs_maxstep , but used for FIRE optimizer. Default value: 0.04

You should create a directory for each and every aimsChain calculation you are planning to run. It should contain the following files.

 -Project directory
 |
 |--chain.in
 |
 |--control.in
 |
 |--ini.in*
 |
 |--fin.in*
 |
 |--extgeo.lst*+
 |
 |--images*+
      |
      |-image1.in*+
      |
      |-image2.in*+
      |
      |-...
*The filename can be set by the user
+The files are only necessary for starting from external geometry.

The initial and final geometry should be in the standard aims input format. Keywords such as constraints will be interpreted. The constrain_relaxation tag should be used with caution–it acts as a double-edged sword in terms of MEP evaluation. You can try to remove the constraint tag from sample 2 and observe the efficiency boost for example. In other cases, setting it may help convergence by reducing degree of freedoms.

In addition, you should pre-align your initial and final geometries to remove any rotational/translational component in the geometry, this will reduce the computational effort required.

The ordering of atoms in the geometry file is very crucial. The atoms in the initial and final geometry must establish a one to one correspondence, so that the $i$th line in the initial and final geometry represents the same atom. Changing the order can change the result from the evaluation dramatically.

For example, consider the diffusion of benzene on a surface. By changing the ordering of atoms you may force the diffusion to be done via rotational or translational motion, which would result in very different MEP. If you can visualize more than one possible combination, try all of them to find the lowest barrier. It is a matter of fact that there often exists more than one MEP between two geometries, but the overall barrier is only determined by the MEP with lowest barrier. Chain of states methods will evolve toward the MEP that’s closest to the initial path, but not necessarily the lowest in energy.

The aims control file should contain two lines.

compute_forces .true.
final_forces_cleaned .true.

This will turn on the force evaluation for aims, which is required for any chain of states method. relativistic, when required, must be set to atomic_zora scalar, since we require force evaluations.

Sometimes additional configurations can require more than the two lines listed above. For example, when using b3lyp you need to set

RI_method lvl_fast

for a correct force evaluation. Generally you can find these information in the aims output. A rule of thumb is to relax the geometry using your configurations and skim through the aims output for any warnings. This will force aims to provide all warnings related to the calcualtion of forces. When you have confirmed that the configuration if warning-free, just comment out the relax_geometry for aimsChain calculation.

You can also set sc_iter_limit to a lower value, which should still be much higher than the normal number of cycles that your system consumes. As the path evolves, it can be the case that during a few iterations the geometry become non-physical, and its scf cycle will never converge. Such cases often recover itself within a few iteration, and will not affect the final result. However, a default limit of 1000 will consume lots of computational power in these cases, which is a complete waste especially for larger systems.

The control file should not have relax_geometry, any molecular dynamic keywords, and etc. The geometry should not be altered by aims. It is advised that you always use a light setting for the first run, so that tunning settings and confirming convergence can be done efficiently. If a more accurate result is desired, you can submit a tight run using the resulting geometries from the light run, which is always faster than a direct tight run from linear interpolation.

To ensure efficiency, you should configure your control file so that each single point calculation takes at most few minutes to finish. If a tighter setting is desired, consider using climb_control for a different control file during climbing image.

It is unlikely that you will need a long list of keywords in your control. Tags should be added gradually if you find that the default settings is not sufficient for your purpose. Keys that should be included on the first run are run_aims , initial_file , final_file , n_image , force_thres , restart , and external_geometry if so desired. Simply grabbing the control from the samples will also work most of the times.

If you believe that your system follows a rather straightforward path geometrically(i.e. diffusion, small rotation, etc.), you can set the climbing image on your first run and see if it works out correctly. For any complex paths, such as non-trivial isomerization, fine tunning of the control file can be required for the system to converge.

The job script should contain a call to runchain.py, as well as loading the necessary modules if that is not done in .bashrc. No post processing should be included unless you are certain that the system will converge and terminate within the time limit.

If you have any information about the path you are trying to find, please include them here. This can be a guess for transition state, geometry reported by papers, or your own interpolation by tweaking atoms in a visualization software. Coupled with resample , any number of external geometries can help speed up the calculation process, making the process very flexible.

It is a common practice to take the intermediate path from a previous calculation (perhaps before it has gone wild) and use them as the initialization path for a new job. This is perfectly fine, and aimsChain outputs intermediate paths just for this purpose. However, when doing so please remember to remove the first and last image from the path, which are the same geometry as the initial and final states. The list of external geometries should only contain intermediate images, entering initial/final states in there will most likely lead to a dead end.

One of the most important step you can take to ensure a successful MEP evaluation is to ensure that the initial path is physical. This is especially true for periodic systems where the periodic boundary condition plays an unwanted part in this.

We have provided a tool in aimsChain/tools, the interpchain.py. This gadget performs the resampling and interpolation process exactly the same way that the actual script is doing it, and therefore is a good way to check if the initial path is reasonable.

When ran, the program will create the interpolation directory in your project directory (and will clear that folder if it already exists!).

 -Project directory
 |
 |--interpolation
      |
      |-image001.in
      |
      |-image002.in
      |
      |-...
      |
      |-path.xyz

There will be n_images +2 aims geometry created, including initial and final geometry. A multi-frame .xyz file is also created in the directory for use with visualization softwares that doesn’t support multi-file animation. If the geometry is periodic, the .xyz file is repeated in each lattice vector according to xyz_lattice to make visualization easier. (because the standard xyz does not encode periodic information) Always check the interpolation if you are working with a new project, which can save you lots and lots of time if you happen to spot a bad interpolation(as shown in sample 4).

Running aimsChain on periodic systems requires extra precaution, because initial files from a periodic system can often provide misleading initial path.

The default interpolation algorithm when periodic_interpolation works as follows. For each atom in the final geometry, the coordinates are offset by each lattice vector in both the positive and negative direction. The results are compared with the same atom in the initial geometry, and the coordinate with the shortest distance is used as the actual coordinate for the final atom. This method is good in general.

However, it’s still possible to imagine cases where a longer path is the actual path. For example, when trying to simulate the effect of an edge on one end. The shortest path may corresponds to climbing over the edge to reach the other end, while the true path, moving in the other direction, is to simply walk away from the edge.

In these cases, the coordinate of the final geometry must be adjusted manually, so that its coordinate corresponds to its true coordinate after the move, and does not involve any periodic boundary condition. periodic_interpolation should be turned off, and the interpolated path should be looked in details to ensure that the atoms in the base are also correctly interpolated.