#! /bin/sh


# Summary
# 
# - For a sequence of intermediates, each defined by a fixed angle H-C-N between
#   the linear extremes HCN and CNH, the remaining geometrical parameters are
#   optimized, giving a Linear Transit point-by-point scan of the energy curve
#   of the Hydrogen atom traveling from one end of the CN fragment to the other.
#   This is a useful way to get a reasonable first guess of the Transition
#   State.
# - At the approximate TS a Frequencies calculation is performed to obtain a
#   fairly accurate Hessian for the next calculation.
# - A TS search is carried out, using the computed Hessian. As variation, the TS
#   search is repeated, first with the automatic (internal) Hessian (based on
#   force fields) and then also with a constraint applied.
# - A full IRC scan of the full path, starting from the TS, down to the two
#   minima.


# == Linear Transit (LT) ==

# The first calculation is a Linear Transit where the Hydrogen atom moves from
# one side of CN to the other by a parametrized step-by-step change of the angle
# H-C-N. The other coordinates of the system are optimized along the path.

# In the atoms block, one coordinate value is represented by an identifier (th).
# In the geovar block this is assigned two values, implying that it is a Linear
# Transit parameter. The initial and final values for the parameter are given.

# Since the geometry block does not have OPTIM SELECTED, all other coordinates
# are optimized for each of the 10 Linear Transit points.

# The subkey iterations in the geometry block carries two arguments: the first
# is the maximum number of optimization steps (per LT point). The second is the
# number of LT points to compute in this run: 4. This implies that only a part
# of the 10-point path defined by the LT parameter(s) will be scanned. The
# remainder will be done in a follow-up run to illustrate usage of the restart
# facility.



$ADFBIN/adf <<eor
Title    HCN Linear Transit, first part
NoPrint  SFO, Frag, Functions, Computation
  
Atoms      Internal
  1 C  0 0 0       0    0    0
  2 N  1 0 0       1.3  0    0
  3 H  1 2 0       1.0  th  0
End

Basis
 Type DZP
End

Symmetry NOSYM

BeckeGrid
 Quality good
End

Geometry
  Branch Old
  LinearTransit  10
  Iterations     30  4
  Converge   Grad=3e-2  Rad=3e-2  Angle=2
END

Geovar
  th   180    0
End

eor

mv TAPE21 t21.LT
rm logfile


# The NoPRINT key turns off a lot of default output. There are several PRINT and
# NOPRINT options; see the User's Guides for details.

# Since the geometry changes from linear to planar (and finally back to linear
# again), the symmetry must be given explicitly in the input file. Otherwise the
# program would find a C(lin) symmetry for the initial geometry and assume that
# this symmetry is preserved throughout. This would of course result in an error
# abort when the first LT step is carried out, breaking the linear symmetry.

# The here specified symmetry (NOSYM: no symmetry at all) is not the true
# symmetry of the complete path C(s) but a subgroup. It is always allowed to
# specify a lower symmetry than the actually present symmetry. Such may be
# necessary (for instance when the true symmetry cannot be handled by adf) or in
# special cases required for reasons of analysis. Generally speaking, however,
# we recommend to use the highest symmetry possible (given the case at hand and
# taking into account the symmetries recognizable by ADF) to boost performance.

# Convergence thresholds in the geometry block are set less tight than the
# defaults: we need only a reasonable estimate of the path, but no highly
# converged geometries.

# At the end of the run the tape21 result file is saved and renamed t21.LT to
# serve as restart file for the follow-up calculation.

# == LT continuation ==


$ADFBIN/adf <<eor
Title    HCN Linear Transit
NoPrint  SFO,Frag,Functions,Computation

Restart 
  File t21.LT
End

Fragments 
  N   t21.N
  C   t21.C
  H   t21.H
End
  
Atoms       Internal
  1 C  0 0 0       0    0    0
  2 N  1 0 0       1.3  0    0
  3 H  1 2 0       1.0  th   0
End

symmetry NOSYM

BeckeGrid
 Quality good
End

Geometry
  Branch Old
  LinearTransit   10
  Converge        Grad=3e-2 Rad=3e-2 Angle=2
END

Geovar
  th   180    0
End

eor

rm logfile
mv TAPE21 HCN_LT.t21


# From the restart file, supplied with the key restart, the program reads off
# that the first 4 points of the LT path have been done already and the scan is
# continued with LT point #5. The same path definition is supplied again,
# including the original starting values for the coordinates. The actual
# starting coordinates (for LT point #5) are read from the restart file. The
# input values, however, serve to define and verify consistency of the defined
# LT path and must therefore be supplied correctly.

# The key noprint is used to suppress major parts of standard output: all
# information pertaining to the sfo analysis, all build-from-fragments
# information, and the lists of elementary functions in the basis sets and fit
# sets.


# == Frequencies at the estimated Transition State ==


# From the results of the Linear Transit run we can sketch the energy barrier
# that H passes over when going from one side of the molecule to the other. This
# yields a reasonable guess for the Transition State.

# To check that the so-obtained estimate is adequate we compute the frequencies
# in that geometry: one of them should be imaginary.

# Apart from serving as a check that the TS estimate is not too bad, the
# computed Hessian will also serve in the follow-up calculation to obtain the
# true TS.



$ADFBIN/adf <<eor
Title    HCN Frequencies in LT max (approx), moderate precision
NoPrint  SFO,Frag,Functions,Computation

BeckeGrid
 Quality good
End
  
Fragments 
  N   t21.N
  C   t21.C
  H   t21.H
End
  
Atoms      Internal
  1 C  0 0 0      0       0   0
  2 N  1 0 0      1.186   0   0
  3 H  1 2 0      1.223  70   0
End

Geometry
  Frequencies
  End
End

eor

rm logfile
mv TAPE21 HCN_Freq.t21

# Inspection of the output file shows that one of the frequencies is imaginary,
# as expected (printed as negative), signaling the proximity of the Transition
# State.

# The TAPE21 result file of the calculation is renamed and saved. Later we will
# use it as a 'restart' file for a TS search, namely to supply the computed
# Hessian as the initial 'guess' of the Hessian in the (TS) optimization run.


# == TS search ==

# Now carry out the Transition State search, starting from the lt-derived guess.

# In this first attempt to find the TS, no use is made of the tape21 result file
# from the Frequencies run. That will be done in the next calculation.



$ADFBIN/adf <<eor
Title    HCN Transition State, automatic initial Hessian
NoPrint  SFO,Frag,Functions,Computation

BeckeGrid
 Quality good
End

Atoms      Internal
  1 C  0 0 0      0       0   0
  2 N  1 0 0      1.186   0   0
  3 H  1 2 0      1.223  70   0
End

Fragments 
  N   t21.N
  C   t21.C
  H   t21.H
End
  
Geometry
  TransitionState
  End
  Converge Grad=0.0001
End

eor

rm TAPE21 logfile

# The TS-search run type is specified in the geometryblock.

# No symmetry is specified; the program determines the symmetry to be C(s) and
# consequently carries out the ts search in that symmetry.


#  == TS search, using the Hessian from the Frequencies run ==


$ADFBIN/adf <<eor
Title    HCN Transition State,  initial Hessian from Freq run
NoPrint  SFO,Frag,Functions,Computation

Restart  
  File HCN_Freq.t21
End
Save     TAPE13

BeckeGrid
 Quality good
End

Atoms          Internal
  1 C  0 0 0      0       0   0
  2 N  1 0 0      1.186   0   0
  3 H  1 2 0      1.223  70   0
End

Fragments 
  N   t21.N
  C   t21.C
  H   t21.H
End
  
Geometry
  TransitionState
  End
  Converge Grad=0.0001
End

eor

mv TAPE13  t13.TS
rm logfile
mv TAPE21 HCN_TS.t21


# The CheckPoint file TAPE13, at normal termination automatically deleted by the
# program, is here saved, using the SAVE key. TAPE13 is as good a restart file
# as TAPE21 is, but it is a lot smaller. TAPE21 contains a large amount of
# information for analysis purposes, while TAPE13 contains essentially only
# restart-type data.

# The input is identical to the previous one, except for the restart file. This
# is used here to provide the Hessian computed in the Frequencies run as the
# start-up Hessian for the ts optimization. At the same time the atomic
# coordinates are read off from the restart file and override the values in the
# input file. This latter aspect could have been suppressed; see the User's
# Guide for using the restart key.



#  == Constrained TS search ==

# Finally the ts search where one coordinate is kept frozen, to illustrate a
# constrained optimization.



$ADFBIN/adf <<eor
Title    HCN  constrained TS search
NoPrint  SFO,Frag,Functions,Computation

Restart 
  File HCN_Freq.t21
End

BeckeGrid
 Quality good
End

Atoms     Internal
  1 C  0 0 0      0       0   0
  2 N  1 0 0      rNC     0   0
  3 H  1 2 0      1.223  70   0
End

GeoVar
  rNC=1.186 F
End

Fragments 
  N   t21.N
  C   t21.C
  H   t21.H
End
  
Geometry
  TransitionState
  End
  Converge Grad=0.0001
End
FULLSCF

eor


echo ""
echo "Final Zmatrix and energy at TS"
$ADFBIN/adfreport -i TAPE21 'Geometry%zmatrix#12.6f##3' bonding
rm TAPE21 logfile

# The geovar key specifies that the nc distance, rNC has the initial value 1.15
# and remains frozen ('F').

# The fact that the optimization is now carried out in a different subspace of
# atomic coordinates does not prevent us from using the t21.Freq restart file to
# supply the initial Hessian.

$ADFBIN/adf <<eor
Title    HCN Frequencies at TS
NoPrint  SFO,Frag,Functions,Computation

BeckeGrid
 Quality good
End
  
Fragments 
  N   t21.N
  C   t21.C
  H   t21.H
End
  
Atoms      Internal
   1 C         0    0    0    0.000000    0.000000    0.000000
   2 N         1    0    0    1.185055    0.000000    0.000000
   3 H         1    2    0    1.208984   69.196764    0.000000
End

AnalyticalFreq
End

eor

mv TAPE21 HCN_Freq_TS.t21


# == IRC scan of the reaction path ==


# The IRC calculation is split in three steps, to illustrate the Restart
# facility applied to the IRC functionality.

# In the first only a few points are computed, along one of the two paths
# leading from the TS to the adjacent minima. Since no explicit directives are
# given in the input to specify the direction of the first path, the so-called
# 'forward' path is taken. The definition of which is 'forward' and which is
# 'backward' is in fact quite arbitrary and is determined by the program. See
# the User's Guide for details.

# The saved TAPE13 file from one of the TS calculations is used as restart file.
# This provides (a) the optimized coordinates of the TS as starting point, (b)
# the initial Hessian to guide the point-by-point optimizations along the IRC
# path, and (c) the eigenvector of the lowest Hessian eigenvalue to define the
# initial direction of the IRC path.

# The TAPE13 file from this partial IRC scan is saved to serve as start-up file
# for the next calculations, which will continue the IRC scan.

# In the Geometry key block, the run type is set to IRC and the 'Points' option
# is used to limit the number of IRC points to compute.



$ADFBIN/adf  <<eor
Title   HCN  IRC partial path (forward)
NoPrint SFO,Frag,Functions, Computation

BeckeGrid
 Quality good
End

Restart  
  File HCN_Freq_TS.t21
End
Save     TAPE13

Atoms          Internal
  1 C  0 0 0      0       0   0
  2 N  1 0 0      1.186   0   0
  3 H  1 2 0      1.223  70   0
End

Fragments
  N   t21.N
  C   t21.C
  H   t21.H
End

Geometry
  IRC  Points=5
  Converge Grad=0.001
End
SCF
  Converge 0.00000001
End
FULLSCF

eor

mv TAPE13  t13.IRC_1
rm TAPE21 logfile
echo ""
echo "Final Zmatrix after IRC"
$ADFBIN/adfreport -i t13.IRC_1 'Geometry%zmatrix#12.6f##3'



# The IRC is continued in the next calculation, using the TAPE13 file from the
# previous one as restart file. From this file, the program reads the IRC path
# information computed sofar. By default, it would continue on the 'forward'
# path, since that was not yet finished. However, in the Geometry key block, we
# now specify not only that a limited number of points is to be computed in this
# run (5 again), but we instruct the program also to compute only points on the
# 'backward' path.



$ADFBIN/adf  <<eor
Title   HCN  IRC partial part (backward)
NoPrint SFO,Frag,Functions, Computation

Restart  
  File t13.IRC_1
End
Save     TAPE13

BeckeGrid
 Quality good
End

Atoms          Internal
  1 C  0 0 0      0       0   0
  2 N  1 0 0      1.186   0   0
  3 H  1 2 0      1.223  70   0
END

Fragments
  N   t21.N
  C   t21.C
  H   t21.H
End

Geometry
  IRC  Points=5  Forward=False
  Converge Grad=0.001
End
SCF
  Converge 0.00000001
End
FULLSCF

eor

mv TAPE13  t13.IRC_2
rm TAPE21 logfile
echo ""
echo "Final Zmatrix after IRC"
$ADFBIN/adfreport -i t13.IRC_2 'Geometry%zmatrix#12.6f##3'


# In the third IRC run, the IRC scan is finished. We start with the TAPE13 file
# from the previous run and set a maximum of 70 IRC points to compute (which
# turns out to be sufficient for the complete IRC scan). The program starts on
# the forward path, continuing where the first (not the previous) had stopped
# after 5 points, completes the forward path, and then continues on the backward
# path, starting where the second IRC run had stopped. Both paths are finished
# and a summary of the path characteristics is printed in the final part of the
# output.


$ADFBIN/adf  <<eor
Title    HCN  IRC completion
NoPrint  SFO,Frag,Functions, Computation

Restart 
  File t13.IRC_2
End

BeckeGrid
 Quality good
End

Atoms          Internal
  1 C  0 0 0      0       0   0
  2 N  1 0 0      1.186   0   0
  3 H  1 2 0      1.223  70   0
End

Fragments
  N   t21.N
  C   t21.C
  H   t21.H
End

Geometry
  IRC  Points=100 StepMax=0.2 StepMin=0.15
  Converge Grad=0.0001
End
SCF
  Converge 0.00000001
End
FULLSCF

eor

mv TAPE21 HCN_irc.t21

echo ""
echo "Final Zmatrix and energy after IRC"
$ADFBIN/adfreport -i HCN_irc.t21 'Geometry%zmatrix#12.6f##3' bonding