xdat2xyz.pl

eval '(exit $?0)' && eval 'exec perl -S $0 ${1+"$@"}' && eval 'exec perl -S $0 $argv:q' if 0;
#;-*- Perl -*-

# Open the XDATCAR. If no argument is supplied XDATCAR is open
# by default. If argument is supplied that file will be opened.

  $zip = $ENV{'VTST_ZIP'} ;
  if($zip eq ''){ $zip = 'gzip' ; }

  if(@ARGV == undef){
    if(-e "XDATCAR"){
      $xdatzipped = 0 ; }
    elsif(-e "XDATCAR.gz"){
      print " XDATCAR IS ZIPPED, UNZIPPING \n" ;
      system "gunzip XDATCAR.gz" ;
      $xdatzipped = 1 ; }
    elsif(-e "XDATCAR.bz2"){
      print " XDATCAR IS ZIPPED, UNZIPPING \n" ;
      system "bunzip2 XDATCAR.bz2" ;
      $xdatzipped = 1 ; }
    else{ die " NO XDATCAR IN THIS DIRECTORY \n" ; }
    print " OPEN XDATCAR ... \n" ;
    open XDAT, "XDATCAR"
      or die " ...  XDATCAR MISSING FROM THIS DIRECTORY \n" ;
  }else{
    print " OPEN ",@ARGV[0]," ... \n" ;
    open XDAT, @ARGV[0]
      or die " ... ",@ARGV[0]," MISSING FROM THIS DIRECTORY \n" ;
  }

# Get the forces and energies to print in the xyz-file
#
#  if(-e "OUTCAR"){
#    $zipped = 0 ; }
#  elsif(-e "OUTCAR.gz"){
#    print " OUTCAR IS ZIPPED, UNZIPPING \n" ;
#    system "gunzip OUTCAR.gz" ;
#    $zipped = 1 ; }
#  elsif(-e "OUTCAR.bz2"){
#    print " OUTCAR IS ZIPPED, UNZIPPING \n" ;
#    system "bunzip2 OUTCAR.bz2" ;
#    $zipped = 1 ; }
#  else{ die " NO OUTCAR IN THIS DIRECTORY \n" ; }
#
#  $forces = `grep 'FORCES: max atom, RMS' OUTCAR` ;
#  $energy = `grep 'energy  without entropy' OUTCAR` ;
#  @forces = split /\n/ , $forces ;
#  @energy = split /\n/ , $energy ;
#  if($zipped){print " ZIPPING OUTCAR AGAIN \n" ; system "$zip OUTCAR &" ;}
#
# Get information about, number and types of atoms, box lengths ets.
# from the POSCAR file.

  open POS, "POSCAR" or die " NO POSCAR IN THIS DIRECTORY \n";
  for($i=0; $i<8; $i++){
    $line = <POS> ;
    chomp($line) ; $line=~s/^\s+//g ; @line=split /\s+/,$line ;
    if($i == 0){
      @elements = split /\s+/ , $line ;
      $nel = @elements ;
    }
    if($i == 1){$latt = $line[0] ;}
    if($i == 2){@sidex[0..2] = @line[0..2] ;}
    if($i == 3){@sidey[0..2] = @line[0..2] ;}
    if($i == 4){@sidez[0..2] = @line[0..2] ;}
    if($i == 5){
      @not[0..$nel-1] = @line[0..$nel-1] ;
      while($not[$k] != undef){$natoms+=$not[$k++] ;}   # Calculate the number of atoms
    }
  }
  close POS ;

# Assign a type (element) to each atom.

  $n = 0 ;
  $j = 0 ;
  for($i=1; $i<=$natoms; $i++){ 
    $j++ ;
    if($j <= $not[$n]){$type[$i-1] = $elements[$n] ;}
    else{$n++ ; $type[$i-1] = $elements[$n] ; $j = 1 ;}  
#    print $i,"  ",$j,"  ",$type[$i-1];
#    $IN = <STDIN> ;
  }
 
# Read the XDAT file and make .xyz files

  for($i=0; $i<5; $i++){$line = <XDAT> ; }      # Jump over the first few lines
  $n=0 ;  
  open MOV, ">TRAJEC.xyz" ;
  while($line = <XDAT>){
    chomp($line) ;
    print MOV $natoms,"\n" ;
    $f = $forces[$n] ; chomp($f) ; $f=~s/^\s+//g ; @f=split /\s+/,$f ;
    $e = $energy[$n] ; chomp($e) ; $e=~s/^\s+//g ; @e=split /\s+/,$e ;
#    print MOV "FORCE:  $f[4]  ...  ENERGY:  $e[6]","\n" ;
    print MOV $n, "\n";
    for($i=0; $i<$natoms; $i++){
      $line = <XDAT> ;
      chomp($line) ; $line=~s/^\s+//g ; @line=split /\s+/,$line ;
#  Transform from direct coordinates to cart. coordinates.
      $x = $line[0] ; $y = $line[1] ; $z = $line[2] ;
#  Periodic boundaries
#      while($x<-0.5){$x+=1.0;} while($x>0.5){$x-=1.0;}
#      while($y<-0.5){$y+=1.0;} while($y>0.5){$y-=1.0;}
#      while($z<-0.5){$z+=1.0;} while($z>0.5){$z-=1.0;}
      $xt=$x*$sidex[0]+$y*$sidey[0]+$z*$sidez[0];
      $yt=$x*$sidex[1]+$y*$sidey[1]+$z*$sidez[1];
      $zt=$x*$sidex[2]+$y*$sidey[2]+$z*$sidez[2];
      $x=$latt*$xt; $y=$latt*$yt; $z=$latt*$zt;
      printf MOV "%2s %18.13f %18.13f %18.13f \n",$type[$i],$x,$y,$z ;  
    }
    $n++ ;
  }
  close XDAT;
  if($xdatzipped){print " ZIPPING XDATCAR AGAIN \n" ; system "$zip -9 XDATCAR &" ;}