# NEB via File-IO¶

[Input: recipes/interfaces/ase/neb/fileio/]

In order for ASE to find the DFTB+ executable and Slater-Koster files, environment variables must be set. Assuming the use of the BASH shell, this is done as follows (in general an adaptation to your specific computing environment will be necessary):

$DFTB_PREFIX=~/slakos/mio-0-1/$ DFTB_COMMAND=~/dftbplus/bin/dftb+


## Geometry¶

To calculate the energy barrier of the umbrella inversion process, it is necessary to provide geometries for the initial and final state.

A possible choice for the .gen files defining the initial and final state of the calculation could be the following:

initial state:

  4  C
N  H
1 1    0.3282523436E+00   -0.6075947736E+00    0.1062852063E+00
2 2    0.1193013364E+01   -0.3265926020E+00    0.5715962197E+00
3 2    0.2770647466E+00   -0.1801902091E+00   -0.8199941035E+00
4 2   -0.4767004541E+00   -0.3215324152E+00    0.6662026774E+00


final state:

  4  C
N  H
1 1    0.3206429517E+00   -0.6319775238E+00    0.1527307910E+00
2 2    0.1179077740E+01   -0.1010702015E+01    0.5563961495E+00
3 2    0.2628891413E+00   -0.8794981362E+00   -0.8365549890E+00
4 2   -0.4905898328E+00   -0.1001302325E+01    0.6515180485E+00


Figure 56 visualizes the defined geometries. The NH3 molecule passes through a planar transition state, thus inverting the pyramidal shape:

Figure 56 Initial a) and final b) state of the NH3 (ammonia) umbrella inversion.

## Main script¶

The Python script shown below is a suitable and simple implementation to run a calculation with NIMAGES as the number of intermediate images (in this case set to 13).

The main method then first reads in the geometry and immediately writes it out for the reasons mentioned in Socket-Comm.: Geometry Optimization by ASE. The trajectories (obtained from previous calculations) of the initial and final state are read in and a list with the individual images is created, whereby the first and last image is determined by the states that have just been read in. True, independent copies of the first image are created for the 13 intermediate images. A NEB object is instantiated and a linear interpolation of the path from the intial and final state is calculated.

Subsequently, the ASE driver BFGS() for geometry optimization is specified and where to write the resulting trajectory (i2f.traj). A list of instantiated DFTB calculators is created and its elements attached to the single intermediate images. Finally the calculation is started by calling the run command, while passing the maximum force component fmax (mind the ASE Units!):

from ase.io import read, write
from ase.neb import NEB
from ase.optimize import BFGS
from ase.calculators.dftb import Dftb

NIMAGES = 13

def main():
'''Main driver routine.'''

images = [initial]
images += [initial.copy() for ii in range(NIMAGES)]
images += [final]

neb = NEB(images)
neb.interpolate()

opt = BFGS(neb, trajectory='i2f.traj')

calcs = [Dftb(label='NH3_inversion',
Hamiltonian_SCC='Yes',
Hamiltonian_SCCTolerance='1.00E-06',
Hamiltonian_MaxAngularMomentum_N='"p"',
Hamiltonian_MaxAngularMomentum_H='"s"')
for ii in range(NIMAGES)]

for ii, calc in enumerate(calcs):
images[ii + 1].set_calculator(calc)

opt.run(fmax=1.00E-02)

if __name__ == "__main__":
main()


## Analysis¶

The results of the NEB calculation can be extracted out of the trajectory file i2f.traj using corresponding ASE tools. If only the last 15 images (intermediate images plus initial and final state) are of interest, a construction like the following one would be suitable:

ase gui i2f.traj@-15:


Since ASE is capable of quickly visualizing the reaction path via the Tools ‣ NEB tab, you should obtain something similar to Figure 57:

Figure 57 Energy barrier of the umbrella inversion of a single NH3 (ammonia) molecule.

Since the initial and final state are identical after a suitable unitary transformation, the energetic difference $$\Delta E$$ between the states obviously vanishes.