Nudged Elastic Band (NEB) for gasphase reaction with DFTB
*********************************************************

DFTB NEB Catalysis

**Requires:** AMS2026 or later

-[ Related examples ]-

* Breaking bonds in a ring molecule and modifying distances for NEB

* Automated Transition-State Workflow for a Diels-Alder Addition

-[ Related tutorials ]-

* NEB tutorial

* Tips and tricks for transition states tutorial

-[ Related documentation ]-

* NEB task

* DFTB engine

[image]


Initial and final (approximate) states
======================================

By default, the NEB task performs geometry optimizations of the
initial and final states before creating interpolated structures
(images).

But it is a good idea to optimize the molecules before starting the
NEB, to make sure that they are in reasonable positions.


Create the systems
==================

Check out the GUI tutorial or the related “building structures”
example to find out how to get the initial structures used here.

   from scm.plams import view
   from scm.base import ChemicalSystem

   initial_sys = ChemicalSystem("""
   System
      Atoms
         C 1.0065881512963437 -1.417955004240563 0.12885420257373528
         C -0.08873018725926754 -2.1350000996605933 -0.03448649693154524
         O -1.5081108856471217 1.002080375323799 -0.2045570136880532
         N -0.4803377947492119 1.514182931120784 -0.052917508314934095
         N 0.4992559331764258 2.001851603909277 0.09160944741412266
         H 1.3565548327663348 -1.1143852514314225 1.1071131753303691
         H 1.6148659535066119 -1.095826278535031 -0.7065218137603257
         H -0.6985763255773102 -2.454234856430313 0.8007371044735704
         H -0.4402537210929026 -2.435675244924285 -1.0129788309893883
      End
   End""")

   view(initial_sys, width=150, height=150, direction="tilt_pca3", guess_bonds=True)

[image]

   final_sys = ChemicalSystem("""
   System
      Atoms
         C 0.8176518265162523 0.006185364744877703 0.11652129166388352
         C -0.5531099061901967 -0.6363379504647356 -0.08529176454211394
         O -1.4164729952124286 0.48607450912407413 -0.19677333330677277
         N -0.6859030532374636 1.632537947036214 -0.08098563281192021
         N 0.5074709667559635 1.4283821326486141 0.08689807666792455
         H 1.2798721321383 -0.22791534308315498 1.0777663111036928
         H 1.5305100901334066 -0.2099050616627963 -0.6819652771198405
         H -0.866292408825807 -1.2485848049850283 0.7685478972017538
         H -0.6137266520780276 -1.2304367933580558 -1.0047175688566132
      End
   End""")

   view(final_sys, width=150, height=150, direction="tilt_pca3", guess_bonds=True)

[image]


NEB and DFTB settings. Run the Job.
===================================

   from scm.plams import Settings, AMSJob

   settings = Settings()

   # Input options for the AMS driver
   settings.input.ams.Task = "NEB"
   settings.input.ams.Properties.PESPointCharacter = "Yes"

   # Input options for the engine (DFTB in this case)
   settings.input.DFTB.Model = "GFN1-xTB"

   # You can pass multiple systems to an AMSJob by using a dictionary of molecules.
   # The key of the dictionary will be used as the header of the 'System' block
   molecules = {"": initial_sys, "final": final_sys}

   # Create and run the job
   job = AMSJob(settings=settings, molecule=molecules, name="NEB")
   print(job.get_input())

   Properties
     PESPointCharacter Yes
   End

   Task NEB

   System
      Atoms
         C 1.0065881512963437 -1.417955004240563 0.12885420257373528
         C -0.08873018725926754 -2.1350000996605933 -0.03448649693154524
         O -1.5081108856471217 1.002080375323799 -0.2045570136880532
         N -0.4803377947492119 1.514182931120784 -0.052917508314934095
         N 0.4992559331764258 2.001851603909277 0.09160944741412266
         H 1.3565548327663348 -1.1143852514314225 1.1071131753303691
         H 1.6148659535066119 -1.095826278535031 -0.7065218137603257
   ... output trimmed ....
         O -1.4164729952124286 0.48607450912407413 -0.19677333330677277
         N -0.6859030532374636 1.632537947036214 -0.08098563281192021
         N 0.5074709667559635 1.4283821326486141 0.08689807666792455
         H 1.2798721321383 -0.22791534308315498 1.0777663111036928
         H 1.5305100901334066 -0.2099050616627963 -0.6819652771198405
         H -0.866292408825807 -1.2485848049850283 0.7685478972017538
         H -0.6137266520780276 -1.2304367933580558 -1.0047175688566132
      End
   End

   Engine DFTB
     Model GFN1-xTB
   EndEngine

   job.run();

   [28.02|09:45:50] JOB NEB STARTED
   [28.02|09:45:50] JOB NEB RUNNING
   [28.02|09:45:51] JOB NEB FINISHED
   [28.02|09:45:51] JOB NEB SUCCESSFUL

   # job = AMSJob.load_external("/path/to/ams.rkf")   # if loading results from disk
   assert job.ok(), "Looks like NEB calculation failed?"
   print("Successful NEB calculation!")

   Successful NEB calculation!


NEB Results
===========

   print(job.results.engine_names())

   ['NEB_TS_final']

   from pprint import pprint

   pes_point_character = job.results.readrkf("AMSResults", "PESPointCharacter", file="NEB_TS_final")
   print(f"{pes_point_character=}")

   pes_point_character='transition state'

   neb_res = job.results.get_neb_results()
   pprint(neb_res)  # energies in hartree

   {'Climbing': True,
    'Energies': [-17.333016762960952,
                 -17.33294947958635,
                 -17.33256809516803,
                 -17.331527201901654,
                 -17.328829238015103,
                 -17.319753903625962,
                 -17.307297193368917,
                 -17.322575012605878,
                 -17.364900430989824,
                 -17.393606196413522],
    'HighestIndex': 6,
    'HistoryIndices': [309, 301, 302, 303, 304, 305, 306, 307, 308, 310],
    'LeftBarrier': 0.025719569592034475,
    'Molecules': [Molecule('C2H4N2O' at 0x7e3cea270370),
                  Molecule('C2H4N2O' at 0x7e3cea270670),
                  Molecule('C2H4N2O' at 0x7e3cea270790),
                  Molecule('C2H4N2O' at 0x7e3cea270280),
                  Molecule('C2H4N2O' at 0x7e3cea270b20),
                  Molecule('C2H4N2O' at 0x7e3cea270dc0),
                  Molecule('C2H4N2O' at 0x7e3cea6c9160),
                  Molecule('C2H4N2O' at 0x7e3cf20742b0),
                  Molecule('C2H4N2O' at 0x7e3ce9ffb0a0),
                  Molecule('C2H4N2O' at 0x7e3ce9ffb310)],
    'ReactionEnergy': -0.06058943345257006,
    'RightBarrier': 0.08630900304460454,
    'nImages': 8,
    'nIterations': 216}

   import matplotlib.pyplot as plt

   neb_res = job.results.get_neb_results()
   fig, ax = plt.subplots()
   ax.plot(neb_res["Energies"])
   ax.set_ylabel("Energy (hartree)")
   ax.set_xlabel("Image (zero-based indexing)");

[image]

   import matplotlib.pyplot as plt
   from scm.plams import view

   neb_res = job.results.get_neb_results()
   mols = neb_res["Molecules"]  # these are PLAMS Molecule, not ChemicalSystem
   indices = [0, 5, 6, 7, 9]  # show only these
   fig, axes = plt.subplots(1, len(indices), figsize=(16, 16))
   for ax, i in zip(axes, indices):
       img = view(mols[i], width=160, height=160, direction="tilt_pca3", guess_bonds=True)
       ax.imshow(img)
       ax.axis("off")
   axes;

[image]


See also
========


Related examples
----------------

* Breaking bonds in a ring molecule and modifying distances for NEB

* Automated Transition-State Workflow for a Diels-Alder Addition


Related tutorials
-----------------

* NEB tutorial

* Tips and tricks for transition states tutorial


Related documentation
---------------------

* NEB task

* DFTB engine


Python Script
=============

   #!/usr/bin/env python
   # coding: utf-8

   # ## Initial and final (approximate) states
   #
   # By default, the NEB task performs geometry optimizations of the initial and final states before creating interpolated structures (images).
   #
   # But it is a good idea to optimize the molecules before starting the NEB, to make sure that they are in reasonable positions.
   #
   # ## Create the systems
   #
   # Check out the GUI tutorial or the related "building structures" example to find out how to get the initial structures used here.

   from scm.plams import view
   from scm.base import ChemicalSystem

   initial_sys = ChemicalSystem(
       """
   System
      Atoms
         C 1.0065881512963437 -1.417955004240563 0.12885420257373528
         C -0.08873018725926754 -2.1350000996605933 -0.03448649693154524
         O -1.5081108856471217 1.002080375323799 -0.2045570136880532
         N -0.4803377947492119 1.514182931120784 -0.052917508314934095
         N 0.4992559331764258 2.001851603909277 0.09160944741412266
         H 1.3565548327663348 -1.1143852514314225 1.1071131753303691
         H 1.6148659535066119 -1.095826278535031 -0.7065218137603257
         H -0.6985763255773102 -2.454234856430313 0.8007371044735704
         H -0.4402537210929026 -2.435675244924285 -1.0129788309893883
      End
   End"""
   )

   view(initial_sys, width=150, height=150, direction="tilt_pca3", guess_bonds=True, picture_path="picture1.png")


   final_sys = ChemicalSystem(
       """
   System
      Atoms
         C 0.8176518265162523 0.006185364744877703 0.11652129166388352
         C -0.5531099061901967 -0.6363379504647356 -0.08529176454211394
         O -1.4164729952124286 0.48607450912407413 -0.19677333330677277
         N -0.6859030532374636 1.632537947036214 -0.08098563281192021
         N 0.5074709667559635 1.4283821326486141 0.08689807666792455
         H 1.2798721321383 -0.22791534308315498 1.0777663111036928
         H 1.5305100901334066 -0.2099050616627963 -0.6819652771198405
         H -0.866292408825807 -1.2485848049850283 0.7685478972017538
         H -0.6137266520780276 -1.2304367933580558 -1.0047175688566132
      End
   End"""
   )

   view(final_sys, width=150, height=150, direction="tilt_pca3", guess_bonds=True, picture_path="picture2.png")


   # ## NEB and DFTB settings. Run the Job.

   from scm.plams import Settings, AMSJob

   settings = Settings()

   # Input options for the AMS driver
   settings.input.ams.Task = "NEB"
   settings.input.ams.Properties.PESPointCharacter = "Yes"

   # Input options for the engine (DFTB in this case)
   settings.input.DFTB.Model = "GFN1-xTB"

   # You can pass multiple systems to an AMSJob by using a dictionary of molecules.
   # The key of the dictionary will be used as the header of the 'System' block
   molecules = {"": initial_sys, "final": final_sys}

   # Create and run the job
   job = AMSJob(settings=settings, molecule=molecules, name="NEB")
   print(job.get_input())


   job.run()
   # job = AMSJob.load_external("/path/to/ams.rkf")   # if loading results from disk
   assert job.ok(), "Looks like NEB calculation failed?"
   print("Successful NEB calculation!")


   # ## NEB Results

   print(job.results.engine_names())


   from pprint import pprint

   pes_point_character = job.results.readrkf("AMSResults", "PESPointCharacter", file="NEB_TS_final")
   print(f"{pes_point_character=}")


   neb_res = job.results.get_neb_results()
   pprint(neb_res)  # energies in hartree


   import matplotlib.pyplot as plt

   neb_res = job.results.get_neb_results()
   fig, ax = plt.subplots()
   ax.plot(neb_res["Energies"])
   ax.set_ylabel("Energy (hartree)")
   ax.set_xlabel("Image (zero-based indexing)")
   import matplotlib.pyplot as plt
   from scm.plams import view

   neb_res = job.results.get_neb_results()
   mols = neb_res["Molecules"]  # these are PLAMS Molecule, not ChemicalSystem
   indices = [0, 5, 6, 7, 9]  # show only these
   fig, axes = plt.subplots(1, len(indices), figsize=(16, 16))
   for ax, i in zip(axes, indices):
       img = view(mols[i], width=160, height=160, direction="tilt_pca3", guess_bonds=True)
       ax.imshow(img)
       ax.axis("off")
   axes
