#!/usr/bin/perl
# read normal modes from a Spartan output file
# determine bond length change between atoms $iatm and $jatm
# for each normal mode
$narg = $#ARGV ;
#$file = shift(@ARGV);    # Extract root file name into pattern
$iatm = -1 ; $jatm = -1 ; $imode = -1 ; # no stretch and none to print
$iatm = shift(@ARGV) ;    # Extract atom number
$jatm = shift(@ARGV) if $narg>0;    # Extract atom number
$imode = shift(@ARGV) if $narg>1;    # Extract mode to print
$file = "output" ;
open(ARC,$file) || die "Can't find file $file\n";
while (<ARC>) {
	if (/coordinates/) { goto data };
	}
data:
   $_ = <ARC> ;
   $_ = <ARC> ;
  $natm = -1 ;
while ( <ARC> ) {
 $natm++ ;
($atom[$natm],$atmno[$natm],$ax[$natm],$ay[$natm],$az[$natm]) = split ;
# print "$atom[$natm];$atmno[$natm];$ax[$natm];$ay[$natm];$az[$natm]\n" ;
   if ( $atom[$natm] eq '' ) { goto modes } ;
}
modes:
while ( <ARC> ) {
   if (/coordinates/) { goto data } ; # get optimized set if present
   if (/amu/) { goto freq } 
}
freq:
 $i = 0 ;
while ( <ARC> ) {
   $_ = <ARC> ; # blank line
   if (/total time/) { goto endspec } 
($f[$i],$f[$i+1],$f[$i+2],$f[$i+3],$f[$i+4],$f[$i+5]) = split ;
#print "f: $f[$i],$f[$i+1],$f[$i+2],$f[$i+3],$f[$i+4],$f[$i+5] \n" ;
  $_ = <ARC> ; # reduced mass
  $_ = <ARC> ; # symmetry
  $_ = <ARC> ; # blank line
for ($j=0 ; $j<$natm ; $j++ ) {
$_ = <ARC> ; # first data line
($a,$n,$x[$j][$i],$x[$j][$i+1],$x[$j][$i+2],$x[$j][$i+3],$x[$j][$i+4],$x[$j][$i+5]) = split;
$_ = <ARC> ;
($y[$j][$i],$y[$j][$i+1],$y[$j][$i+2],$y[$j][$i+3],$y[$j][$i+4],$y[$j][$i+5]) = split;
$_ = <ARC> ;
($z[$j][$i],$z[$j][$i+1],$z[$j][$i+2],$z[$j][$i+3],$z[$j][$i+4],$z[$j][$i+5]) = split;
}
$i+=6 ;
  }
endspec:
# get number of modes
for ($j=0 ; $j<=$i ; $j++ ) {
if ( $f[$j]<1 ) { last }
}
$nmode = $j ; # number of modes
# print "found $nmode modes____________________\n" ;
# list mode $iatm
format MODE_TOP =
atom   #      x        y        z 
.
format MODE =
@<<< @>>> @#.####  @#.####  @#.####
$atom[$j],$j,$x[$j][$imode],$y[$j][$imode],$z[$j][$imode]
.
  $^ = "MODE_TOP" ; $~ = "MODE" ; $- = 0 ; $= = 2000 ;
if ( $imode > 0 ) {
print "mode $imode:  $f[$imode] cm-1\n" ;
#print "atom   #        x            y            z \n" ;
for ($j=0 ; $j<$natm ; $j++ ) {
#print "$atom[$j] $j $x[$j][$imode]  $y[$j][$imode]  $z[$j][$imode] \n" ;
write ;
  }
print "______________________________________________________________\n" ;
# now make animation XYZ file
open(XYZ,">/web/data/chemistry/NMR/scripts/vib/xyz.xyz");
  for ($fr=0 ; $fr<20 ; $fr++) {
    $inc = $fr*0.6/19.0 ;
   print XYZ "$natm\n" ;
   print XYZ "b\n" ;
    for ($j=0 ; $j<$natm ; $j++ ) {
    $fx = $ax[$j]+$inc*$x[$j][$imode] ;
    $fy = $ay[$j]+$inc*$y[$j][$imode] ;
    $fz = $az[$j]+$inc*$z[$j][$imode] ;
	print XYZ "$atom[$j] $fx $fy $fz\n";
    } ; # end for $j
  } ; # end for $fr
 close(XYZ) ;
}
# determine bond length change
  for ($k=0 ; $k<5 ; $k++ ) { $idb[$k] = 0 ; $db[$k] = 0 }
if ( ($iatm > 0)&&($jatm > 0) ) {
print ">>> Total Relative Distance and Bond Length Changes (ang)<<<\n" ;
# unit bond direction vector
$d = ($ax[$jatm]-$ax[$iatm])**2 ;
$d += ($ay[$jatm]-$ay[$iatm])**2 ;
$d += ($az[$jatm]-$az[$iatm])**2 ;
$d = sqrt($d) ;
# get unit vector in bond direction
$ux = ($ax[$jatm]-$ax[$iatm])/$d ;
$uy = ($ay[$jatm]-$ay[$iatm])/$d ;
$uz = ($az[$jatm]-$az[$iatm])/$d ;
print "$atom[$iatm] ($iatm)- $atom[$jatm] ($jatm) " ;
  $d = substr($d,0,6) ;
print "bond length=$d \n" ;

format BOND_TOP =
                 total      bond
mode  freq cm-1 relative    <-> 
.
format BOND =
@>>> @####.##   @#.####    @#.####
$j,$f[$j],$d,$b 
.
  $^ = "BOND_TOP" ; $~ = "BOND" ; $- = 0 ;
for ($j=0 ; $j<$nmode ; $j++ ) {
$dx = $x[$jatm][$j]-$x[$iatm][$j] ;
$dy = $y[$jatm][$j]-$y[$iatm][$j] ;
$dz = $z[$jatm][$j]-$z[$iatm][$j] ;
$d = $dx**2 + $dy**2 + $dz**2 ;
$d = sqrt($d) ;
$b = $dx*$ux + $dy*$uy + $dz*$dz ;
#print "mode $j:  $f[$j] cm-1 rel: $d <-> $b \n" ;
write ;
# keep 5 biggest changes
  for ($k=0 ; $k<5 ; $k++ ) {
   if ( $b >= $db[$k]) {
      for ($m=4 ; $m>$k ; $m--) { $db[$m]=$db[$m-1] ; $idb[$m]=$idb[$m-1] }
   $db[$k]=$b ; $idb[$k]=$j ;
   last ;
   } ; # end if $d>=
  } ; # end for $k
} ; # end for $j=0
print "==============================================================\n" ;
# print top 5
print " Top five changes:\n" ;
format BIGGEST_TOP =
mode  freq cm-1   <-> 
_________________________
.
format BIGGEST =
@>>> @####.##   @#.####
$idb[$k],$f[$idb[$k]],$db[$k]
.
  $^ = "BIGGEST_TOP" ; $~ = "BIGGEST" ; $- = 0 ;
  for ($k=0 ; $k<5 ; $k++ ) {
#     print "mode $idb[$k]:  $f[$idb[$k]] cm-1 <-> $db[$k]\n" ;
   write ;
  }
print "==============================================================\n" ;
} ; # end if $iatm>0 && $jatm>0