Nudged Elastic Band (NEB) Examples

Here are a few examples showing how the NEB method can be used to obtain the the path and transition state of a reaction.

See also

Nudged Elastic Band (NEB) Documentation

HCN isomerization reaction with NEB

Download NEB_HCN.run

#!/bin/sh

# This example demonstrates the use of the Nudged Elastic Band method in AMS for
# finding a transition state of the HCN isomerisation reaction. A shell script
# used to run the example calculation is shown below.


$AMSBIN/ams <<eor

Task NEB

System
    Atoms
        C         0.000000    0.000000    0.000000
        N         1.180000    0.000000    0.000000
        H         2.196000    0.000000    0.000000
    End
End
System final
    Atoms
        C         0.000000    0.000000    0.000000
        N         1.163000    0.000000    0.000000
        H        -1.078000    0.000000    0.000000
    End
End

NEB
    Images 9
    Iterations 100
End

Engine DFTB
    Model DFTB3
    ResourcesDir DFTB.org/3ob-3-1
    DispersionCorrection D3-BJ
EndEngine
eor

H2 dissociation on graphene

Download NEB_H2_on_graphene.run

#!/bin/sh

for map1 in no yes
do

for map2 in no yes
do

AMS_JOBNAME=test.map1=$map1.map2=$map2 $AMSBIN/ams << eor

Task NEB
System
    Atoms
        C -1.23079526  1.45426666 -0.60065817 
        C  1.25628125  1.47378921 -0.00438520 
        C -2.44885253 -2.15934252 -0.51047194 
        C -0.01262403 -2.15933012 -0.51045137 
        C -1.23076689 -0.05455308 -0.60068696 
        C  1.25630128 -0.07404160 -0.00441877 
        C -2.44887495 -0.74178904 -0.51048970 
        C -0.01264534 -0.74177652 -0.51049375 
        H  1.25629472 -0.34738650  1.07475852
        H  1.25630533  1.74708724  1.07481086
    End
    Lattice
        4.97414302 0.0 0.0
        0.0 4.30083513 0.0
    End
    MapAtomsToUnitCell $map1
End

System final
    Atoms
        C -1.24330284  1.42711444 -0.25240444 
        C  1.23494804  1.42739453 -0.25493930 
        C -2.48233069 -2.12769513 -0.25358631 
        C -0.00426996 -2.12770161 -0.25349945 
        C -1.24330266  0.00893844 -0.25240978 
        C  1.23494831  0.00864800 -0.25496555 
        C -2.48233386 -0.70957077 -0.25358531 
        C -0.00426696 -0.70956458 -0.25349840 
        H  1.27176634  0.72925445  2.26688554 
        H  1.28137254  0.73053214  3.01215519 
    End
    Lattice
        4.95651766 0.0		  0.0
        0.0        4.27331845 0.0
    End
    MapAtomsToUnitCell $map2
End

Properties
    NormalModes Yes
End

GeometryOptimization
    Convergence Gradients=3.0e-3
    HessianFree
        Step TrustRadius=0.5
    End
End

NEB
    ClimbingThreshold 0.01
    OptimizeLattice Yes
End

Engine DFTB
    Model DFTB3
    ResourcesDir DFTB.org/3ob-3-1
    DispersionCorrection D3-BJ
    KSpace Quality=Basic
EndEngine

eor
done
done

Running multiple NEB calculations using PLAMS

Download NEB.plms

This example should be executed using PLAMS.

See also

PLAMS documentation and tutorial

import os

# This PLAMS script locates TS for several reactions using the Nudged Elastic Band (NEB) method.
# For each reaction two xyz files are required: one corresponding to the reactant state, called
# '{}_initial.xyz', and one corresponding to the product state, called '{}_final.xyz' (e.g.
# 'MyMolecule_initial.xyz' and 'MyMolecule_final.xyz')

# The folder containing the xyz files:
xyz_folder = os.path.join(os.environ["AMSHOME"], "examples", "AMS", "NEB_PLAMS", "xyz")

names = [name.rsplit("_initial.xyz")[0] for name in os.listdir(xyz_folder) if name.endswith("_initial.xyz")]

# Settings for the AMS driver
neb_sett = Settings()
neb_sett.input.ams.Task = "NEB"
neb_sett.input.ams.Properties.NormalModes = "Yes"
neb_sett.input.ams.GeometryOptimization.Convergence.Step = 1.0e-3

# Settings for the engine (here we use the DFTB engine with semiempirical GFN1-xTB method)
engine_sett = Settings()
engine_sett.input.DFTB.Model = "GFN1-xTB"

for name in sorted(names):

    # For NEB we need two system blocks in the AMS input (the initial and final molecule).
    # In PLAMS you can have multiple system blocks by passing the the AMSJob a dictionary of molecules.
    # The 'keys' of the dictionary will be used as the headers of the System block
    mols = {}
    # The initial molecule should be in the main 'System' block (the main system block has no header).
    # The key of the mols dictionary should therefore be an empty string
    mols[""] = Molecule(os.path.join(xyz_folder, name + "_initial.xyz"))
    # The final molecule should be specified in a system block with the header 'final'
    mols["final"] = Molecule(os.path.join(xyz_folder, name + "_final.xyz"))

    # Create and run the job:
    job = AMSJob(molecule=mols, name=name, settings=neb_sett + engine_sett)
    job.run()

    print("")
    print("System name: {}".format(name))
    if job.ok():
        pes_point_character = job.results.readrkf("AMSResults", "PESPointCharacter", file="NEB_TS_final")
        print("NEB calculation converged. PES Point character: {}".format(pes_point_character))
        print("Left  TS barrier: {0:.6f} [Hartree]".format(job.results.readrkf("NEB", "LeftBarrier")))
        print("Right TS barrier: {0:.6f} [Hartree]".format(job.results.readrkf("NEB", "RightBarrier")))
    else:
        print("Unsuccesful NEB calculation")
    print("")