#!/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), 
#                              packing into primitive cell (26/07/11)
#
################################################################################

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";
# calculate 3x3 matrix inverse
    $det=($x1*$y2*$z3+$x2*$y3*$z1+$x3*$y1*$z2-$x1*$y3*$z2-$x2*$y1*$z3-$x3*$y2*$z1);
    $invx1=($y2*$z3-$y3*$z2)/$det;
    $invy1=($y3*$z1-$y1*$z3)/$det;
    $invz1=($y1*$z2-$y2*$z1)/$det;
    $invx2=($z2*$x3-$z3*$x2)/$det;
    $invy2=($z3*$x1-$z1*$x3)/$det;
    $invz2=($z1*$x2-$z2*$x1)/$det;
    $invx3=($x2*$y3-$x3*$y2)/$det;
    $invy3=($x3*$y1-$x1*$y3)/$det;
    $invz3=($x1*$y2-$x2*$y1)/$det;
    
## verify 
#    $ix1=$x1*$invx1+$y1*$invx2+$z1*$invx3;
#    $ix2=$x2*$invx1+$y2*$invx2+$z2*$invx3;
#    $ix3=$x3*$invx1+$y3*$invx2+$z3*$invx3;
#    $iy1=$x1*$invy1+$y1*$invy2+$z1*$invy3;
#    $iy2=$x2*$invy1+$y2*$invy2+$z2*$invy3;
#    $iy3=$x3*$invy1+$y3*$invy2+$z3*$invy3;
#    $iz1=$x1*$invz1+$y1*$invz2+$z1*$invz3;
#    $iz2=$x2*$invz1+$y2*$invz2+$z2*$invz3;
#    $iz3=$x3*$invz1+$y3*$invz2+$z3*$invz3;
#   
#   print STDERR "INVERSE-------------------\n";
#    print STDERR "$ix1 $ix2 $ix3\n";
#    print STDERR "$iy1 $iy2 $iy3\n";
#    print STDERR "$iz1 $iz2 $iz3\n";
#   print STDERR "INVERSE-------------------\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);
### for 3d      printf("%s \n",$katm*8);
      print "MD STEP: $line[1]\n";
      print STDERR "MD STEP: $line[1]\n";
     }
     if(@line > 2){

       $lx=$line[3]*$bohr;
       $ly=$line[4]*$bohr;
       $lz=$line[5]*$bohr;

       $lx1=$invx1*$lx+$invx2*$ly+$invx3*$lz;
       $ly1=$invy1*$lx+$invy2*$ly+$invy3*$lz;
       $lz1=$invz1*$lx+$invz2*$ly+$invz3*$lz;
       
       $lx1=$lx1+100-int($lx1+100);
       $ly1=$ly1+100-int($ly1+100);
### for 3d      $lz1=$lz1+100-int($lz1+100);
       
       $lx=$x1*$lx1+$x2*$ly1+$x3*$lz1;
       $ly=$y1*$lx1+$y2*$ly1+$y3*$lz1;
       $lz=$z1*$lx1+$z2*$ly1+$z3*$lz1;

       printf("%s %20.12f %20.12f %20.12f\n",$line[2],$lx            ,$ly            ,$lz            );
       printf("%s %20.12f %20.12f %20.12f\n",$line[2],$lx+$x1        ,$ly+$y1        ,$lz+$z1        );
       printf("%s %20.12f %20.12f %20.12f\n",$line[2],$lx    +$x2    ,$ly    +$y2    ,$lz    +$z2    );
       printf("%s %20.12f %20.12f %20.12f\n",$line[2],$lx+$x1+$x2    ,$ly+$y1+$y2    ,$lz+$z1+$z2    );
### for 3d      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$lx        +$x3,$ly        +$y3,$lz        +$z3);
### for 3d      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$lx+$x1    +$x3,$ly+$y1    +$y3,$lz+$z1    +$z3);
### for 3d      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$lx    +$x2+$x3,$ly    +$y2+$y3,$lz    +$z2+$z3);
### for 3d      printf("%s %20.12f %20.12f %20.12f\n",$line[2],$lx+$x1+$x2+$x3,$ly+$y1+$y2+$y3,$lz+$z1+$z2+$z3);

     }
   }
 }
}
exit