2.15. MedeA Transition State Search: Reaction and Diffusion Pathways, Transition States, and Activation Barriers


download:pdf

The Transition State Search module combines different computational approaches to map out the minimum energy pathway between an initial and a final atomic configuration, and to search and refine transition states on this minimum energy path, along a defined direction or in randomly chosen directions. For mapping out the minimum energy path, it uses a chain of configurations (called images) connected by elastic springs, calculates forces for each configuration and uses information from neighboring images to optimize all configurations in the chain (Nudged Elastic Band method, in brief NEB method). At each of the energy maxima along the path, transition states can be searched for, refined and optimized by elastic band methods, employing the so-called mode following techniques (Dimer method) and energy-based optimization techniques (RMM-DIIS).

The overall goal is to find one or more saddle points of the energy hypersurface, i.e. configurations that are local energy minima in all but one coordinate. That can be verified by looking at their phononic or vibrational spectrum, which has one imaginary frequency in such a case.

Based on the JobServer/TaskServer architecture, transition state search may run on any hardware configuration from single-processor desktop machines up to highly parallel computer racks.

2.15.1. Basic Tasks

In the MedeA GUI select Tools >> Transition State Search to access a new item appearing in the menu bar. Click on Transition State Search to open the menu.

../../_images/image00118.png

The menu’s first section displays three different tasks:

Map the reaction path and find transition states brings up a panel to configure a mapping and optimization of the minimum energy path between an initial configuration and a final configuration followed by a search and optimization of transition states along the path.

Find transition state in a direction brings up a panel to configure a transition state search and optimization starting from the initial configuration in a direction towards another state utilizing the Dimer method.

Find transition states in random directions brings up a panel to configure transition state searches and optimizations of created structures in several different randomly chosen directions utilizing the Dimer method.

Furthermore, Select Engine allows for selecting a specific VASP version for the calculations to search and find transition states. The three tasks above are only supported by VASP 6 and VASP 5.4. For VASP 5.2 and 4.6 a restricted set of simulation options are available from a different user interface (not shown here).

Energy Profile visualizes and analyzes the energy and structures along the path between the initial configuration and final configuration, which is created by running the task Map the reaction path and find transition states. In the below example, the blue line outlines the path between the initial configuration at the zero value of reaction coordinate the found transition state (blue diamond and energy maximum) and the final configuration at the value of one of the reaction coordinate. The transition state search started with the initial NEB path (red line) which was formed by 4 images between the initial configuration and final configuration. The initial NEB path is further refined around the energy maximum by a second NEB optimization (orange) and third NEB optimization (green line). In each of the path refinements 4 images are used. In a further step, the transition state guess (energy maximum of the green line, i.e. third NEB optimization) is further refined employing the Dimer method (not resolved in the graph). Finally, an optimization of the transition state leads to the configuration marked by the blue diamond. The basic workflow of this task and its computational options are described below, the output and analysis of results are presented in Section Output for Reaction Path Mapping. The other two tasks related to the Dimer method are discussed in sections Find Transition State in a Direction and Find Transition States in Random Directions.

../../_images/image00218.png

2.15.2. Map the Reaction Path and Find Transition States

2.15.2.1. Requirement

For this feature two structures are necessary. One structure is the initial configuration, which can reactants of a chemical reaction, the initial step of diffusion process, etc. The other structure is the final configuration, which can be the products of a chemical reaction, the final step of a diffusion process, etc. Both configurations must be a periodic models of symmetry P1 (no point group symmetry defined). Important is that both structures have identical dimensions of the simulation cell (cell lengths and cell angles) and identical compositions. If one of the two latter requirements are not fulfilled then MedeA Transition State Search cannot be used to map a path. Also it is strongly recommended that atoms in the initial configuration have a clear one-to-one correspondence to atoms in the final configuration. For example the atom with the name atom 1 should have the same name in the final configuration. With this knowledge you can construct a chemically reasonable reaction coordinate that describes the transition of the initial configuration to the final configuration via one or several energy maxima and transition states.

2.15.2.2. Select Initial Configuration

Select the structure window with the initial configuration and make it active, so that its title bar turns green.

Note

If some atoms of the selected initial configuration are frozen, these degrees of freedom are omitted from the optimization process. As a consequence, these atoms are not moved and remain fixed for all images.

Click on Transition State Search >> Map the reaction path and find transition states to select the final configuration and to configure the protocol to map the minimum energy path, perform the transition state search and optimization structures of identified energy maxima.

../../_images/image00614.png

2.15.2.3. Select Final Configuration

All periodic structures in the MedeA GUI that match with the structure of the initial configuration (identical stoichiometry and cell contents) and P1 symmetry are available for being selected as the final configuration. Click on the drop-down menu next to the final configuration to change the proposed choice.

../../_images/image0079.png

2.15.2.4. Step 1: Map the Reaction Path Using the Elastic Band Method

../../_images/image0098.png

The first step of this simulation protocol is mandatory and cannot be deselected. Its objective is to explore the energy hypersurface and to identify possible locations of transition states and intermediate minima along the reaction pathway between the initial configuration and final configuration. In the subsequent steps, such candidates for transition states can be refined and optimized either for all identified saddle points or for the highest saddle point only. Three main parameters define the mapping of the reaction path using the elastic band methods:

Number of images: The number of images needed depends on the complexity of the reaction path (single or multiple saddle points), the reaction barrier and the distance between respective atoms in the initial configuration and final configuration. A smaller number of images can still deliver excellent energy barriers when using path refinement steps or the climbing image option (see the Advanced Settings below).

Convergence: sets the criterion for stopping the optimization process, also for the path refinement steps. The maximum force acting on any atom of any image including the spring contributions is tested.

Maximum number of steps allowed for the optimization process, also for each path refinement step

Clicking on the Show Advanced Settings button displays several further options to tweak the performance and improve the success rate of the elastic band technique:

../../_images/image0107.png

For the Type of calculation there is a choice between two variants of the Elastic Band method: Nudged Elastic Band or Elastic Band. For the conceptual difference between these approaches see the Section Elastic Band Method

With the Number of refinement steps you define how often the path around the energy maxima is defined, after a completed initial elastic band calculation. With the path refinement the resolution of the section of the energy path is improved, which is located between images that bracketing the expected saddle points. In each path refinement the initial paths is set up by the same number of linearly interpolated images between the images that bracketing the expected saddle points. With a small total number of images the transition barrier can be over or underestimated; creating a path with smaller distances between the images by increasing total number of images between the initial configuration and final configuration also increases the computational cost and affects convergence. The path refinement overcomes this situation by a closer examination of the region around the suspected transition states.

Spring constant: Atoms of adjacent images are connected by elastic springs to keep images close to the minimum energy path. The spring constant determines the chain stiffness. A soft chain allows adjacent images to differ substantially from each other and may sample a larger search space, however, it is more vulnerable to corner-cutting and down-sliding of images in steep and curved sections of the energy hypersurface. The recommended default value is 5 eV/\({\mathring{\mathrm{A}}}\)2.

Initial inverse Hessian matrix: the diagonal elements are set to this value and pre-condition the optimizer to take bold moves or rather small steps. Reduce the value of the initial inverse Hessian, if the initial gradient is too big and the chain of images gets off track after 3 or 5 elastic band steps. A general value for all systems cannot be provided, as the choice depends on the ratio of moving atoms to more or less static atoms, initial forces and masses of moving atoms. However, the inverse of this value should be larger than the force constant that is associated with the highest vibrational frequency of both the initial configuration and the final configuration. That is why the recommended default value is 0.001 \({\mathring{\mathrm{A}}}\)2/eV.

Additional options specific for the nudged elastic band method are:

Climbing image close to the saddle point: The almost equidistant images of a nudged elastic band may be inconveniently aligned off the saddle point positions. As an alternative to path refinement steps, a selected image close to the saddle point can be made climbing up into the transition state, once the convergence reaches a value specified by Convergence before climbing. This may yield transition states directly by a single optimization without an intermediate second step for further refinement.

Note

It is recommended to apply the latter two options in the step 2 of the transition state search protocol and do not add a mark to the check box of the option Climbing image close to the saddle point.

Add double nudging: This uses a refined scheme for calculating the forces perpendicular to the chain of spring connected images (see Section Elastic Band Method). This option is especially very useful in case of long straight energy paths that have small energy barriers.

VASP Settings: This brings up the MedeA VASP interface to customize computational parameters of the VASP calculations. If the property Energy of formation is selected from the Calculation tab of the MedeA VASP interface, the formation energy for initial configuration and final configuration, as well as all transition states is reported in Job.out.

2.15.2.5. Image Initialization

As a default, initial images are obtained from linear interpolation, i.e. without user interaction, a set of linearly interpolated images will be used to create the initial chain between the initial configuration and final configuration filling the full range.

../../_images/imageTSS08.png

Translation criterion deals with some ambiguity in setting up the initial chain of images due to the periodicity of the cell: When the distance between matching atoms in the initial configuration and final configuration is larger than half of the respective cell coordinate, the movement can be either within the cell or between the atom and its periodic copy in a neighboring cell.

Table 34 Examples for the Translation criterion
image12 The value of 0.5 always uses the shortest distance between the initial configuration and final configuration, whereas a larger value allows for longer paths within the cell. In the illustration on the left, a value of 0.75 is needed to allow for the path indicated by the dashed red arrow between the initial configuration (0) and configuration (1) image. For the default value of 0.5, the shorter path along the solid blue arrow is pursued.

Suppress net atom movements takes care that the sum of all movements of atoms between the initial configuration and final configuration vanish. This is accomplished by an overall shift of the final image, and improves convergence due to vanishing spring forces at the boundary images, in particular for smaller systems or systems with strong relative drift.

Setting the Interpolation range to in the middle, parts of the reaction path can be excluded from optimization by specifying the Starting boundary image and Ending boundary image path parameters. For example, when a reaction occurs close to a surface, but the initial configuration or final configuration are in quite some distance from the surface, these options allow the omission of the physisorption part of the reactants and products. It is noted, however, that the new boundary images are affected by linear interpolation and do not correspond to local minima. It is recommended to avoid these options by finding initial configuration and final configuration close to the surface before transition state search.

../../_images/imageTSS10.png

Setting the Interpolation range to close to initial and final, the interpolated images are accumulated close to the initial configuration and final configuration. Specifically, half of the images are created to be close to initial configuration, i.e. between the Starting boundary image and Upper initial area image path parameters. The other half if the created images is interpolated between the Lower final area image and Ending boundary image path parameters. This setting is recommended if unphysical images would occur in the range between Upper initial area image and Lower final area image upon linearly interpolating the full range. Since the nudged elastic band technique enforces equidistant images during optimization, the images will move into the middle range avoiding close contacts between atoms, which would occur interpolating linearly in the entire range.

../../_images/imageTSS11.png

To speed up the redistribution of images, an order of magnitude larger spring constants and low precision VASP parameters are recommended. Such a simulation could be run before the actual Transition State Search simulation, with the sole objective to obtain physically reasonable initial images, which can be introduced into the subsequent simulation in the following way:

You can click on View initial images to inspect the images generated by any of the interpolation modes.

../../_images/image0162.png

You can select Initial images from specified systems to define a specific pathway with structures that you created by manual modifications of the structure of the initial configuration, the final final configuration, or both. Also you can define a specific pathway by specifying images with structures of a previously performed NEB calculation to, e.g. continue a non-converged path mapping.

The amount of structures that should be specified depends on the number of images that you want to use in the path mapping step. In this example the number of images is set to 4, hence, for 4 structures can be specified.

../../_images/image0137.png

Note

All periodic structures in the MedeA GUI that match with the structure of the initial configuration (identical stoichiometry and cell contents) and P1 symmetry can be selected.

2.15.2.6. Configurations Tab

Click on the Configurations tab to inspect how individual atoms of the initial configuration and the final configuration are connected and apply corrections if needed. Each atom of the initial configuration and final configuration needs to connected and correspond to each other, i.e. appearing in the same line of the table. If this is not the case, the sequence of atoms of the final configuration can be modified by choosing atoms to be matched and applying the change. These modifications directly update the final configuration model beneath.

../../_images/image0185.png

For example, if you want to connect atom #36 O_CO2_2 of the initial configuration with atom #38 O_CO2_1 of the final configuration then select these names in the drop-down menus for Link the initial configuration atom and with the final configuration atom and confirm with Apply.

The value of Total Distance is the sum of the values in the Distance column in the below table.

To determine which atoms move most between the initial configuration and final configuration right-click on the header of the Distance column >> Sort Descending.

2.15.2.7. Step 2: Refine the Transition State(s) Using Elastic Band Methods

A second simulation step can be selected to refine the transition state candidate(s) obtained in the first step by two different methods. With the Elastic Band Methods the Transition State Search module performs a Climbing Image Nudged Elastic Band (CI-NEB) calculation between the images that are found adjacent to the transition state position in the first step with Number of images and Convergence criterion is chosen independently from those in the first simulation step. All other parameters for the elastic band method are applied as set for step 1, the path mapping.

../../_images/image0203.png

2.15.2.8. Step 2: Refine the Transition State(s) Using Dimer Method

Alternatively, the transition state candidate can be further refined with the Dimer Method, which is a minimization technique following the lowest eigenmode of the Hessian. The latter is the matrix with the force constants.

../../_images/image0214.png

The following parameters trigger the performance of the Dimer method:

The Convergence criterion and the Maximum number of steps are set independently of those set in step 1.

Dimer separation: the step size for numerical differentiation

Minimum rotation: the dimer is rotated only if the predicted rotation angle is greater than this value

Trial step size: the trial step size for the optimization step

Maximum step size: the trust radius (upper limit) for the optimization step

Differentiation scheme: whether to use a forward finite difference or a central finite difference formula for the numerical differentiation to compute the curvature along the dimer direction

2.15.2.9. Step 3: Refine the Transition State(s) Using the RMM-DIIS Minimizer

Check the box next to Refine the transition state(s) using RMM-DIIS minimizer to apply RMM-DIIS optimization to the approximate saddle point structure obtained by the search process specified in step 2. The convergence radius of this method is rather small and requires an initial structure very close to the optimized transition state.

Ensure that the transition state is well approximated by the search process, e.g. based on one or more refinement steps or the CI-NEB procedure and sufficiently strict convergence criteria in the first step, and/or by including a second refinement step. The attempted optimization may fail, if a structure of the transition state guess converges towards a local minimum energy configuration. If the optimization succeeds the optimized transition state can be used for subsequent phonon calculations to derive temperature-dependent thermodynamic functions and calculate free energies.

2.15.3. Output for Reaction Path Mapping

2.15.3.1. Job.out

The detailed output file describes the method and parameters chosen for the transition state search as well as the computational parameters of the underlying VASP calculations, followed by one or more (e.g. for path refinements) tables of energies with additional notes of accepted steps and climbing of certain images summarizing the optimization process.

The results of transition state refinement with the Dimer Method in step 2 and RMM-DIIS optimization in step 3 is indicated by headers as shown below, followed by the output as obtained after a standard VASP geometry optimization. The results for transition state refinement with the Elastic Band Methods are presented like a refinement step, however, with a different number of images.

2.15.3.1.1. Settings for the Transition State Search Method

Nudged Elastic Band for mapping the minimum energy path between
the initial system min image10 36772
and the final system neb0_image18
with 4 intermediate images and a spring constant of 5 eV/Ang^2 and 2 refinement steps.

The initial images are created from linear interpolation.

In a second step, transition states are searched for the highest saddle point only by the Elastic Band methods.
In a last step, optimization of transition states by the RMM-DIIS algorithm is attempted.

Optimization parameters for the first step:
 Convergence: 0.1 eV/Ang
 Number of steps: 50
 Diagonal elements of the inverse Hessian are initially set to 0.001 Ang^2/eV.
------------------------------------------------------------------------

2.15.3.1.2. VASP Parameters for Energy and Force Evaluation

VASP parameters
===============
This is a calculation based on density functional theory and the GGA-PBEsol exchange-correlation functional for describing the interactions.

Since no magnetic moments are in the model, this is a non-magnetic calculation using 'normal' precision
and a default planewave cutoff energy of 282.853 eV.

The electronic iterations convergence is 1.00E-07 eV using the Normal (blocked Davidson) algorithm
and real space projection operators.

The requested k-spacing is 1.5 per Angstrom, which leads to a 1x1x1 mesh.
This corresponds to actual k-spacings of 0.823 x 0.582 x 0.419 per Angstrom.
The k-mesh is forced to be centered on the gamma point.

Using Gaussian smearing with a width of 0.05 eV.

2.15.3.1.3. Progress Report on Initial Elastic Band Calculation

VASP energy of initial and final boundary images in kJ/mol per cell:

      Image         Energy (kJ/mol)
------------------- -------------
     neb0_image00     -32595.347
     neb0_image05     -32504.154

Total and image energies below are given with respect to the energy of the initial boundary image in kJ/mol per cell

Iter Energy_total max grad     image01      image02      image03      image04   Iter_accepted
---- ------------ --------- ------------ ------------ ------------ ------------ -------------

1      2082.97  41.24831      171.106      642.226      982.047      287.594
2      1509.00  28.35392      146.905      453.671      661.775      246.645            1
3       849.84  12.25940       98.416      254.135      331.455      165.831            2

...

49       400.88   0.58698       60.776      106.787      129.541      103.772          48
50       405.70   0.66861       60.954      112.379      128.754      103.615          49

The columns of the table are (from left to right)

  • the NEB step (iteration) index
  • the sum of all relative energies of the images between the two end configurations,
  • the maximum force (energy gradient) including the contributions due to the spring forces that connect atoms between images
  • the relative energies of each image, and whether an NEB iteration is accepted.

In this case example, mapping the reaction path was performed in 50 NEB steps (iterations) but without reaching the convergence criterion of the forces. Nevertheless, the mapped path shows an energy maximum for image03 which is about 128 kJ/mol energetically less stable than the initial configuration. Hence, Transition State Search expects that the saddle point (transition state) is between image 2 and image 4 and continues with a first path refinement between the two latter images. For the first path refinement four new images are created between image 2 and image 4 by linear interpolation.

2.15.3.1.4. Report on 1st Path Refinement by the Elastic Band Method:

Refinement step 1 between image 2 and image 4:
------------------------------------------------

Iter Energy_total max grad     image01      image02      image03      image04   Iter_accepted
---- ------------ --------- ------------ ------------ ------------ ------------ -------------

1       711.40  11.02872      167.304      155.900      232.549      155.645
2       726.87   9.02329      159.339      212.192      207.551      147.783            1
3       555.77   2.44558      135.778      150.963      145.161      123.872            2

...

49       461.00   0.44334      117.406      121.924      116.397      105.271          48
50       457.29   0.32663      117.308      122.086      116.403      101.489          49

In this case example, for the first path refinement also 50 NEB steps (iterations) are performed but without reaching the convergence criterion of the forces. Nevertheless, the path refinement yielded an energy maximum for image02 which is about 122 kJ/mol energetically less stable than the initial configuration. Hence, Transition State Search expects that the saddle point (transition state) is between image 1 and image 3 of the refined path and continues with a second path refinement between the two latter images. For the second path refinement four new images are created between image 1 and image 3 by linear interpolation.

2.15.3.1.5. Report on 2nd Path Refinement by the Elastic Band Method:

Refinement step 2 between image 1 and image 3:
------------------------------------------------

Iter Energy_total max grad     image01      image02      image03      image04   Iter_accepted
---- ------------ --------- ------------ ------------ ------------ ------------ -------------

 1       499.35   2.26269      122.751      127.484      127.220      121.891
 2       496.54   2.02346      122.329      126.604      126.134      121.477            1

 ...

 49        93.47   2.66443      115.954      117.716     -254.812      114.610          45
 50       502.23   3.70621      124.838      130.187      130.768      116.439

In this case example, for the second path refinement also 50 NEB steps (iterations) are performed but without reaching the convergence criterion of the forces. Nevertheless, the path refinement yielded for image02 and image03 almost identical energies, which are about 130 kJ/mol energetically less stable than the initial configuration. Hence, Transition State Search expects that the saddle point (transition state) is somewhere between image 2 and image 3 of the refined path and continues with a CI-NEB calculation to refine the transition state. The initial guess of the actual structure that is further refined in the subsequent CI-NEB calculation is created by linearly interpolating the structures of image 2 and image 3.

2.15.3.1.6. Report on Transition State Refinement by the Climbing-Image NEB Method:

Search the transition state between image 2 and image 3 by the climbing image NEB method:
----------------------------------------------------------------------------------------

Iter Energy_total max grad     image01   Climbing Iter_accepted
---- ------------ --------- ------------ -------- -------------

1       129.66   2.86203      129.664  ---
2       127.57   2.45484      127.575  ---      1
3       119.91   0.48743      119.907  ---      2

...

26        93.86   0.47533       93.860  1
27       117.32   0.08600      117.324  1      25

Iterations: 25 using 27 calls to the function
Energy: -32478.02325647

The CI-NEB method stopped after 27 steps (iterations) because the maximum gradient was smaller than the convergence criterion of 0.1 eV/\({\mathring{\mathrm{A}}}\). The final structure of image01, i.e. that of iteration 27 of the CI-NEB method is about 117 kJ/mol less stable than the initial configuration.

In the refinement of the transition state the CI-NEB approach (variant of the elastic band method), image02 and image03 of the second path refinement step are set as new boundaries states.

2.15.3.1.7. Report on Transition State Optimization by the RMM-DIIS Minimizer:

The final section of Job.out summarizes the VASP parameters and results of the optimization with VASP using the RMM-DIIS structure optimize (step 3). The initial structure for step 3 is the final structure of step 2.

Results for the attempted transition state optimization:
------------------------------------------------------------

There are 1 symmetry-unique k-points
The plane wave cutoff is 282.85 eV

VASP energy:            -336.687935 eV for Ce12CO25H cell

Initial VASP energy:    -336.477950 eV for Ce12CO25H cell
Relaxation energy:        -0.209985 eV gained after 79 optimization steps.

Electronic contributions:
                       Empirical Formula          Cell
                          Ce12CO25H       Ce12CO25H
                       ----------------- -----------------
           VASP Energy        -336.687935     -336.687935 eV
                     =      -32485.439      -32485.439 kJ/mol

2.15.3.2. Energy Profile

Select Transition State Search >> Energy Profile to visualize results from completed transition state searches running the task Map the reaction path and find transition states.

The Energy Profile panel graphically summarizes the results of the transition state search in terms of energy profiles for each refinement step and provides access to all structures of the images as well as trajectories of the optimization process. Also, images can be selected and removed from the graph, and thereby custom energy profiles can be created for publication and presentation purposes. For the custom profiles, the reaction pathway can be animated, visualizing the structural evolution during the reaction. The capabilities and functionality of the Energy Profile panel are discussed in more detail below.

In the graph pane, the energy of the individual images and transition states are plotted against the reaction coordinate x, which is a measure for the progress of the reaction from the initial configuration (x=0) to the final configuration (x=1). All energies are given relative to the energy of the initial configuration. Two definitions for the reaction coordinate can be applied and are toggled by the choice Type in the Display option frame below the graph:

1. initial coordinate / energy refers to a nominal coordinate solely based on the number of intermediate images not involving a distance metric. As an example, if 3 images are specified they appear at coordinates 0.25, 0.5, and 0.75, if 4 images are specified the formal coordinates are 0.2, 0.4, 0.6, and 0.8. The coordinates for images of refinements are defined in the same manner as subdivisions. This type of reaction coordinate indicates how images are set up initially and does not depend on the structural result of the optimization.

(1)\[d_{i} = \mid R_{i+1} - R_{i} \mid x_{k} = \dfrac{\sum_{i=0}^{k} d_{i}}{\sum_{i=0}^{n} d_{i}}\]

2. optimized coordinate / energy involves a distance metric and thereby takes into account the structural evolution during the optimization. The distance between two adjacent images is calculated as the norm of the generalized distance vectors involving differences in coordinates for all pairs of matched atoms in both images. The reaction coordinate of image k is defined as the fraction of the sum of distances between initial configuration and image k to the total sum of distances between all n images of the chain.

For refinements, the reaction coordinates of boundary images are kept at their values and the intermediate images of the refinement are calculated by the above definition between these values.

(2)\[x_{TS} = \dfrac{R_{TS} - R_{0}}{R_{n} - R_{0}}\]

For transition states, the initial coordinate is defined as shown above whereas the optimized coordinate is defined as the relative distance between transition the state and initial configuration to the overall distance between initial configuration and final configuration

Differences in energy can be measured by dragging the black measurement lines vertically, for instance, to measure the activation energy or the reaction energy. The vertical black bar indicates the energy difference and can be dragged horizontally to a convenient position in the graph. The display of the measurement lines can be switched off by clicking the Measuring lines checkbox in the Display options frame.

../../_images/image0281.png

In the barrier shown above the red markers indicate image energies of the initial nudged elastic band run. Results from the further path refinement steps are shown in orange (neb1)and green (neb2). In case of the first path refinement images are placed between previous images with the highest energy; that is between images 2-3-4 of the initial chain to capture the path through the transition state in more detail. In the frame titled Select images for custom profiles, display of the curves connecting images is toggled by the checkboxes named neb0 neb1, and neb2 for the initial nudged elastic band run and the two path refinement steps. The images to be displayed as markers are selected from the table of checkboxes to the right.

The selection of images from the table of checkboxes does not only toggle the display of markers but allows for the creation of custom energy profiles. Images pulled off the minimum energy path, e.g. by spring forces acting across large distances, can be removed from the graph and for a subset of images a new energy profile is created with the option Create profile button above the center Custom profiles frame. This adds a custom profile colored in blue connecting only the selected images and a button Profile 1 for toggling its display is created in the Custom profiles frame. Multiple custom profiles can be created bringing up further buttons Profile 2 etc. All colored markers can be displayed by clicking on Select all, or can be removed with the Unselect all button (showing up only if all images are selected). To recover the information, which images are used to create a custom profile, check on the Show images of selected profile which displays these images as blue markers.

../../_images/image0292.png

The structural evolution along the reaction pathway can be visualized as trajectory animation. Over-barrier animation is only accessible for selected (i.e. displayed) blue custom profiles and can be launched by clicking on Animate selected custom profiles button at the bottom of the Display options frame.

Transition states are indicated in the energy profile panel by blue diamonds after attempting optimization by RMM-DIIS in step 3. The results of The buttons named for instance TS final allow toggling display as well as to select them for being included in custom profiles.

Structure information and collected VASP output for each image and refinement step is stored on the JobServer. The optimized results are contained in the main job directory and the structural optimization history can be deduced from trajectory files. VASP output for all iterations for each image is stored in the iterations subfolder. You can use optimized images as a starting point for further refinement of energy pathways and transition states, respectively.

Move the cursor over an image marker in the energy profile close enough such that it becomes highlighted by a green rectangle. A right-click then provides additional items for the context menu such

  • View Image (initial): displays initial structure of the selected image as a MedeA structure in the framework viewer
  • View Image (final): displays final structure of the selected image as a MedeA structure in the framework viewer

If the selected Image is a transition state, there are two options:

  • View Transition State (initial): display the structure of the transition state before RMM-DIIS optimization
  • View Transition State (final): display the structure of the transition state after RMM-DIIS optimization

In case of every image the context menu has the item Animate Trajectory to open a window and display the trajectory of the selected image showing the structural evolution during the optimization process. The animation can be stopped at any frame and the current structure can be extracted from the trajectory by Edit >> Duplicate. This provides access to any structure of any image during optimization.

2.15.4. Find Transition State in a Direction

The menu item Find transition state in a direction raises the below panel to configure a transition state search and optimization starting from the initial configuration in a direction towards another state - which doesn’t need to be a final configuration of a reaction - through the Dimer method. The parameters triggering the performance of the Dimer method are identical to those described in Section Step 2: Refine the Transition State(s) Using Dimer Method and also their values are identical. In the compute protocol step 1 is mandatory and cannot be deselected.

../../_images/image030.png

The subpanel for Dimer initialization provides two options. If the initial configuration is a local minimum, the Dimer method would not be able to move towards a transition state and the optimization would finish after one step. To start the optimization process, mark the check box of the option Shift the initial dimer such that the dimer is shifted towards the final configuration and enables the optimization process to start. The default shift of 0.1 is reasonable and works for most of the test cases. However, in cases with expected energy barriers of larger than 0.5 eV it is recommended to increase the shift to 0.2.

../../_images/image0314.png

Furthermore, it might be beneficial to exclude certain degrees of freedom from the dimer initialization, which can be achieved by freezing atomic coordinates of the initial configuration. If this is the case, the Dimer initialization subpanel allows you to either take into account these constraints for the initial dimer orientation only or to keep these atoms fixed for the entire dimer optimization.

Mark the check box next to Step 2: Refine the transition state using RMM-DIIS minimizer to apply RMM-DIIS optimization to the approximate saddle point structure obtained by the Dimer method search.

VASP Settings: This brings up the MedeA VASP interface to customize computational parameters of the VASP calculations. If the property Energy of formation is selected from the Calculation tab of the MedeA VASP interface, the formation energy for the initial configuration and the transition state is reported in Job.out.

The Configurations tab is identical to that of the other computational tasks and is described in Section Configurations Tab.

2.15.5. Find Transition States in Random Directions

The menu item Find transition states in random directions brings up a panel to configure several transition state searches and optimizations in different randomly chosen directions employing the Dimer method. The number of directions can be specified by Number of random search directions.

For random search directions, it is particularly important to restrict possible search directions only to atoms that are directly involved in chemical reactions or diffusion processes. The coordinates of all other atoms should be frozen for the entire dimer orientation.

../../_images/image032.png

If the initial configuration is a local minimum, the Dimer method would not be able to move towards random directions and the optimization would finish after one step. To start the optimization process, mark the check box of the option Shift the initial dimer such that the dimer is shifted towards the final configurations and enables the optimization processes to start. The default shift of 0.1 is reasonable and worked for most of the test cases. However, in cases with expected energy barriers of larger than 0.5 eV it is recommended to increase the shift to 0.2.

Check the box next to Step 2: Refine the transition state using RMM-DIIS minimizer to apply RMM-DIIS optimization to each of the approximate saddle point structures obtained by the dimer searches of the first step.

VASP Settings: This brings up the MedeA VASP interface to customize computational parameters of the VASP calculations. If the property Energy of formation is selected from the Calculation tab of the MedeA VASP interface, the formation energy for the initial configuration and all transition states is reported in Job.out.

download:pdf