NEB calculation for a reaction barrier using DFTB


Performing a spin polarized calculation with DFTB	Bibliography

Setting up the NEB object

Open the Builder .

Click the Add button in the lower left-hand panel and select Add From Database....

The Database tool will open. From the menu, select Databases → Molecules. Type NH3 in the search field, and you should see the following.

Double-click the "Ammonia" entry (or the icon), and a NH3 molecule will be added to the stash.

When running NEB calculations it is important to restrict at least a few degrees freedom. What primarily happens in the ammonia reaction considered here is that the nitrogen atom moves through the hydrogen plane. It is therefore a good idea to keep the hydrogen atom plane fixed when defining the initial reaction path. This can be acheived rather easily.

Activate the Move Tool . Select one of the hydrogen atoms as an anchor atom (marked by a red halo); this is the atom you will operate on. Since no atoms were selected when the Move tool was opened, all atoms will simultaneously be selected (marked by a yellow halo); these atoms will be affected by any operation carried out in the Move tool.

In the Move widget, enter the value 0 for the z coordinate of the selected atom, and press Enter.

Close the Move tool by closing the widget or clicking the select tool.

Now ensure that the Ammonia molecule is selected on the stash by left-clicking it, then click the Copy button. This adds an additional, identical ammonia molecule to the stash.

Open the Mirror tool from the Coordinate Tools in the right tool bar. Ensure that the mirror plane has a normal parallel to the Z axis and contains the origin (this is the default), and press Apply to invert the molecule.

You have now defined the initial and final states of the reaction path, and the next step is to set up the NEB configuration by inserting a number of images in-between these states.

To this end, open the Create Nudged Elastic Band tool under Builders in the right-hand panel, and drag the two configurations into the "Initial" and "Final" slots.

Set the maximum distance between images to 0.09 Å to have an odd number of images, which is very important since the reaction path is symmetric. Choose the Linear method and press Create NEB.

The NEB object will now be built and added to the stash. Push the play button to see the different configurations that will be used to as starting configurations when optimizing the reaction pathway.

Performing the NEB simulation

To set up the NEB calculation, drag the NEB configuration from the stash area and drop it onto the Script Generator .

In the Script Generator,

Add a New Calculator block.
Add an OptimizeGeometry block.
Change the default output name to nh3_neb.nc

Now open the New Calculator block.

Select the "ATK-SE: Slater-Koster" calculator.
Go to the "Slater-Koster basis set" and check that the "DFTB [mio]" basis set is selected.
Uncheck the No SCF iteration box.

Open the OptimizeGeometry block.

You are now ready to perform the NEB optimization. Send the script to the Job Manager and execute it.

The optimization will take 20-30 minutes depending on your computer.

Note

While most other implementations of the DFTB model use a pointwise electrostatic pair potential interaction model, ATK-SE uses a Poisson solver for calculating the electrostatic interactions. This approach is significantly slower for small systems, like the current NH3 molecule, but faster for large systems where it scales as $O(N)$ instead of a typical $O(N^2)$ scaling.

The use of a Poisson solver ensures that the DFTB model is equivalent to the other calculators in ATK, and thus can be used for modeling device geometries and allows for using implicit solvent models.

Analyzing the NEB simulation

Return to the main VNL window and select the file nh3_neb.nc. Select the last NEB configuration in the file (the one with the highest ID number). Press the Show configuration and you will see the calculated energy barrier and the corresponding configurations in the reaction path. The reaction path can also be illustrated as an animation by pressing the play button.

Figure 2: Reaction path of the ammonia inversion, computed with the DFTB [mio] method. The end-points were not optimized.

Note

The obtained reaction barrier is not fully accurate since the end points were not optimized. This can not least be seen from the fact that the second image has a lower energy than the end-points. For real work with the NEB method the end-points should always be optimized first. If you do this for the ammonia molecule, you will obtain the barrier shown in the figure below.

Reaction path of the ammonia inversion, computed with the DFTB [mio] method. In this case the end-points were first optimized, and the image distance was also adjusted (0.07 Å) to have the same number of images as above.

Figure 3: Reaction path of the ammonia inversion, computed with the DFTB [mio] method. In this case the end-points were first optimized, and the image distance was also adjusted (0.07 Å) to have the same number of images as above.

	Tip
	You can inspect an individual optimized configurations in the reaction path closer by clicking the dots in the reaction paths plot (and stopping the movie).

A recipe for faster calculations

In the calculation above you used a linear interpolation between the end points to set up the initial guess for the NEB path. In many cases the linear path is however very far from the optimized one. Comparing the movie of the converged result above to the initial NEB path you set up, it can be seen that that the hydrogen atoms need to move out to make room for the nitrogen atom as it passes through the center of the molecule, and then they move back in. If this behavior somehow could already be part of the initial path, one would save some steps in the ionic dynamics, and hence time. However, since the initial and final positions of the hydorgen atoms are identical, this cannot be captured by the linear interpolation.

There are however other methods for generating the initial path, and these are is available in ATK. Go back to the Builder, and the Nudged Elastic Band plugin. Keep the image distance but change the Method to BC Method and click Create NEB. If you now look at the initial path, in the Builder, you can see that it already contains this behavior of the hydrogen atoms.

Follow the same procedure as above to run the NEB calculation via the Scripter, and you will find that it runs about 30% faster this time, but the results are identical.

	Note
	The BC (bond conservation) method is similar to the Halgren-Lipscomb method [5] which is also available in ATK (but not suitable for this particular system).