# Author: Artur Michalak
# Added to database: 7/02/05 

# Purpose: Takes an ADF frequency calculation output and writes a series of
# xyz structures that correspond to a chosen vibration. The series of
# structures can be viewed as a movie as usual.
#
# To run:   awk -f freq freq=xxx.xxx output.log > vib.xyz
#
# where the xxx.xxx is the frequency (in cm-1) of the vibration that
# you are intersted in taken from the ADF output. It must be given 
# rounded to three decimal places. output.log is the name of your
# ADF output file and can of course be anything that you want.
# The program writes the structures to standard out and so in the
# example here the output is redircted to the file vib.xyz which
# again can be anything that you want.
#

BEGIN {
   infinity      = 10000000 ;
   COORDRecord   = infinity ;
   MODERecord    = infinity ;
   natoms        = 0 ;
   iatoms        = 0 ;
   freqOn        = 0 ;
   modeOn        = 0 ;
   snapshots     = 20 ;
   scale         = 3.0 ;
   }
                                                                                                        
NR==1 && (freq=="") {
   print "Missing freq= on the command line" ;
   exit 1 ;
   }
                                                                                                        
/\*  F R E Q U E N C I E S  \*/ { freqOn = 1 ; }
/^  *Coordinates \(Cartesian\)$/ && (freqOn==1){ COORDRecord = NR ; natoms = 0 ; next ; }
NR == COORDRecord + 1 && !/^ *==*$/{ COORDRecord = infinity ; next ; }
NR == COORDRecord + 3 && !/  *Atom  *bohr  *angstrom/{ COORDRecord = infinity ; next ; }
NR == COORDRecord + 4 && !/  *X  *Y  *Z  *X  *Y  *Z/{ COORDRecord = infinity ; next ; }
NR == COORDRecord + 5 && !/^ --*$/{ COORDRecord = infinity ; next ; }
NR >= COORDRecord + 6 {
   if( NF < 8 ){ COORDRecord = infinity ; next ; }
   label[natoms] = $2 ;
   x    [natoms] = $6 ;
   y    [natoms] = $7 ;
   z    [natoms] = $8 ;
   natoms++ ;
   }
                                                                                                        
/Vibrations and Normal Modes.*cartesian coordinates, NOT mass-weighted/{ modeOn = 1 }
/Normal Modes in Z-matrix coordinates/{ modeOn = 0 }
(NF<=3) && (modeOn==1) {
   for( i = 1 ; i <= NF ; i++ ){
      if( $i == freq ){
         ndisp = 0 ;
         MODERecord = NR ;
         field = i ;
         next ;
         }
      }
   }
(NR >= MODERecord + 2) && NF < 1 { MODERecord = infinity ; reportVib() ; exit 0 }
(NR >= MODERecord + 2) {
   alabel = substr($1,index($1,".")+1) ;
   if( alabel != label[ndisp] ){
      printf "Atom %d has different labels in XYZ (%s) and displacement (%s)\n", \
             ndisp+1, label[ndisp], alabel ;
      exit 1 ;
      }
   dispx[ndisp] = $(2+3*(field-1)) ;
   dispy[ndisp] = $(3+3*(field-1)) ;
   dispz[ndisp] = $(4+3*(field-1)) ;
   ndisp++ ;
   }
                                                                                                        
function dumpgeo() {
   print natoms ; print "" ;
   for( i = 0 ; i < natoms ; i++ ){
      printf "%4s %20.8f %20.8f %20.8f\n", label[i], x[i], y[i], z[i] ;
      }
   }
                                                                                                        
function reportVib() {
   if( ndisp != natoms ){
      printf "Number of atoms (%d) and displacement vectors (%d) do not agree\n", \
             natoms, ndisp ;
      exit 1 ;
      }
   for( iatom = 0 ; iatom < natoms ; iatom++ ){
      rx[iatom] = x[iatom] ;
      ry[iatom] = y[iatom] ;
      rz[iatom] = z[iatom] ;
      }
   for( isnap = -snapshots ; isnap <= snapshots ; isnap ++ ){
      s = (scale*isnap)/snapshots ;
      for( iatom = 0 ; iatom < natoms ; iatom++ ){
         x[iatom] = rx[iatom] + s*dispx[iatom] ;
         y[iatom] = ry[iatom] + s*dispy[iatom] ;
         z[iatom] = rz[iatom] + s*dispz[iatom] ;
         }
      dumpgeo() ;
      }
   }