#!/usr/bin/perl
################################################################################
#
# state2xyz.pl: generate an XYZ formatted file from an output file of STATE
#            program package
#             with duplication of atoms and thinning of time step
#
# Usage: state2xyz.pl [ file thinning ]
#
# Written by I.Hamada @ ISIR, Osaka university (26/04/05)
# Modified by K.Inagaki @ Prec Osaka univ. (04/08/10)
#
################################################################################

if( @ARGV == 0 ){
  print STDERR "Usage: state2xyz.pl [ file ] [ thinning ]\n";
  exit;
}

if( $ARGV[0] eq "-h" || $ARGV[0] eq "--help"){
  print STDOUT "Usage: state2xyz.pl [ file ] [ thinning ]\n";
  exit;
}

if( ! -e $ARGV[0] ){
  print STDERR "$ARGV[0] does not exist.\n";
  exit;
}

 $thinning=$ARGV[1];
if( $ARGV[1] eq "" ){
  $thinning=1;
}
  

$icount=0;
$icountm=-1;
$icountxyz=0;
$katm=0;

$bohr=0.529177;

print STDERR "input file is $ARGV[0]\n";

open(DATAFILE,$ARGV[0]);

while(<DATAFILE>){
 @line=split(/\s+/,$_);
 if($icount == 1){
   $katm++;
#  print "katm=$katm\n";
 }
 if($icount == 1){
   if(@line != 11){
     $icount=0;
     $katm--;
   } 
 }
 if($line[5] eq 'TAUX'){
  $icount=1;
  $katm=0
 }

 if($icountxyz == 3){
  $x3=$line[5]*$bohr;
  $y3=$line[6]*$bohr;
  $z3=$line[7]*$bohr;
  $icountxyz=0;
    print STDERR "$x3 $y3 $z3 \n";
 }
 if($icountxyz == 2){
  $x2=$line[5]*$bohr;
  $y2=$line[6]*$bohr;
  $z2=$line[7]*$bohr;
  $icountxyz=$icountxyz+1;
    print STDERR "$x2 $y2 $z2 \n";
 }
 if($icountxyz == 1){
  $x1=$line[5]*$bohr;
  $y1=$line[6]*$bohr;
  $z1=$line[7]*$bohr;
  $icountxyz=$icountxyz+1;
    print STDERR "$x1 $y1 $z1 \n";
 }
 if($line[1] eq 'INVERSION'){
  $icountxyz=1;
 }
 if($line[0] eq 'MD:'){
   if(@line == 2){
    $icountm=$icountm+1;
    if($icountm==$thinning){
     $icountm=0;
    }
   }
   if($icountm==0){
     if(@line == 2){
#     print "$katm\n";
      printf("%s \n",$katm*4);
      print "MD STEP: $line[1]\n";
      print STDERR "MD STEP: $line[1]\n";
     }
     if(@line > 2){
      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$line[3]*$bohr            ,$line[4]*$bohr            ,$line[5]*$bohr            );
      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$line[3]*$bohr+$x1        ,$line[4]*$bohr+$y1        ,$line[5]*$bohr+$z1        );
      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$line[3]*$bohr    +$x2    ,$line[4]*$bohr    +$y2    ,$line[5]*$bohr    +$z2    );
      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$line[3]*$bohr+$x1+$x2    ,$line[4]*$bohr+$y1+$y2    ,$line[5]*$bohr+$z1+$z2    );
#      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$line[3]*$bohr        +$x3,$line[4]*$bohr        +$y3,$line[5]*$bohr        +$z3);
#      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$line[3]*$bohr+$x1    +$x3,$line[4]*$bohr+$y1    +$y3,$line[5]*$bohr+$z1    +$z3);
#      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$line[3]*$bohr    +$x2+$x3,$line[4]*$bohr    +$y2+$y3,$line[5]*$bohr    +$z2+$z3);
#      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$line[3]*$bohr+$x1+$x2+$x3,$line[4]*$bohr+$y1+$y2+$y3,$line[5]*$bohr+$z1+$z2+$z3);
     }
   }
 }
}
exit