Example: PES scan and transition state search for H2 on graphene
****************************************************************

"Download PESScan_and_TS_H2_on_Graphene.run"

   #! /bin/sh

   # First we do a 2D PES scan varying the z-coordinate of the two hydrogen atoms
   # In this example we will keep the graphene slab fixed. From a physical/chemical
   # standpoint this is not a good approximation. The graphene slab is
   # intentionally not perfectly symmetric.

   AMS_JOBNAME=PESScan $AMSBIN/ams << EOF

      Task PESScan

      System
         Atoms
            H   0.0      1.53633037    1.1
            H   0.0     -0.11341359    1.1
            C   0.001    1.42028166    0.0
            C   1.230    2.13042249    0.0
            C   1.230   -0.71014083    0.0
            C   2.460    0.00000000    0.0
            C   2.460    1.42028167    0.0
            C   0.000    0.00000000    0.0
         End
         Lattice
            3.69   -2.13042249   0.0
            0.00    4.26084499   0.0
         End
      End

      PESScan
         ScanCoordinate
            nPoints 10
            Coordinate 1 Z 1.1 2.0
         End
         ScanCoordinate
            nPoints 10
            Coordinate 2 Z 1.1 2.0
         End
      End

      GeometryOptimization
         Convergence Step=1.0e-3
      End

      Constraints
         # Fix the entire graphene slab.
         Atom 3
         Atom 4
         Atom 5
         Atom 6
         Atom 7
         Atom 8
      End

      Engine DFTB
         Model DFTB
         ResourcesDir DFTB.org/3ob-3-1
         DispersionCorrection D3-BJ
         KSpace
           Type Symmetric
           Symmetric KInteg=3
         End
      EndEngine

   EOF


   # A human looks at the PES scan and picks a reasonable starting point for the
   # TS search. (Normally you would do that in AMSMovie by looking at the PES and
   # then exporting the geometry into an xyz file.)

   #           _           ____
   #        ___))         [  | \
   #        ) //o          | | ]
   #     _ (_    >         | | ]
   #    (O)  \__<          | | ]
   #    [/] /   \)        [__|/_
   #    [\]|  ( \         __/___\_____
   #    [/]|   \ \__  ___|            |
   #    [\]|    \___E/%%/|____________|_
   #    [/]|=====__   (_________________)

   cat << EOF > initial_geometry_for_TS.xyz
   8

   H    0.4145668856457391      1.72927656037925      1.100000023839768     region=H2
   H   -0.05533871972549955    -0.06805093626643093   1.500000013242627     region=H2
   C    0.001                   1.42028166            0.0
   C    1.230                   2.13042249            0.0
   C    1.230                  -0.71014083            0.0
   C    2.460                   0.00000000            0.0
   C    2.460                   1.42028167            0.0
   C    0.000                   0.00000000            0.0
   VEC1 3.69 -2.13042249 0.0
   VEC2 0.0   4.26084499 0.0
   EOF


   # Compute the partial initial Hessian to be used in the transition state
   # search. (The Hessian will be computed only for the hydrogen atoms.)

   AMS_JOBNAME=Hessian $AMSBIN/ams << EOF

      Task SinglePoint

      System
         # Load the geometry we just saved.
         GeometryFile initial_geometry_for_TS.xyz
      End

      Properties
         # Calculate the Hessian (implied when calculating normal modes) ...
         NormalModes True
         # ... but only the part related to the hydrogen atoms.
         SelectedRegionForHessian H2
      End

      Engine DFTB
         Model DFTB
         ResourcesDir DFTB.org/3ob-3-1
         DispersionCorrection D3-BJ
         KSpace
           Type Symmetric
           Symmetric KInteg=3
         End
      EndEngine

   EOF

   echo "Extract the frequencies from the kf file using amsreport:"
   $AMSBIN/amsreport Hessian.results/dftb.rkf -r "Vibrations%Frequencies[cm-1]##1"


   # Do a transition state search using the initial Hessian just computed (the
   # Graphene slab is constrained). Also compute the final Hessian for the
   # hydrogen atoms to validate the TS.

   AMS_JOBNAME=TS $AMSBIN/ams << EOF

      Task TransitionStateSearch

      System
         # Load the geometry we just saved.
         GeometryFile initial_geometry_for_TS.xyz
      End

      GeometryOptimization
         Quasi-Newton
            Step TrustRadius=0.05
         End
         Convergence Gradients=1.0e-4
         InitialHessian
             # Load previously calculated Hessian as initial Hessian for a
             # transition state search with the Quasi-Newton optimizer.
             Type FromFile
             File Hessian.results/dftb.rkf
         End
      End

      TransitionStateSearch
          # Follow the mode with the smallest frequency.
          ModeToFollow 1
          # (This is also the default, we wouldn't need to specify this.)
      End

      Constraints
         # Fix the entire graphene slab.
         Atom 3
         Atom 4
         Atom 5
         Atom 6
         Atom 7
         Atom 8
      End

      Properties
         NormalModes Yes
         SelectedRegionForHessian H2
      End

      Engine DFTB
          Model DFTB
          ResourcesDir DFTB.org/3ob-3-1
          DispersionCorrection D3-BJ
         KSpace
           Type Symmetric
           Symmetric KInteg=3
         End
      EndEngine

   EOF

   echo "Extract energy from the rkf file using amsreport:"
   $AMSBIN/amsreport TS.results/dftb.rkf -r "AMSResults%Energy"
