#!/usr/bin/perl
#
#
# This script is to sift through the output of DFT++ md dft.log file,
# and pick out ion locations
# and generate a gif movie.  You can view with with xanim or netscape.
#
#This was written by Tairan, and modified by Niko and Evan
#November 13, 1999
#


# default options
# 1 for balsac lattice, 0 for balsac cluster
$lattice = 0;

# default parameters
#supercell variables define the supercell for the cluster output
$supercell_n[0]=2.0;
$supercell_n[1]=2.0;
$supercell_n[2]=1.0;

$dt = 1.0;
$frame = 1;
$xsize = 400;
$nz1 = 4;
$nz2 = 4;
$nz3 = 40;
$ninit = 1;
$theta = 74;
$phi = 62;
$magnf = 1;
$icon = 0;
$nha = 0;

for ($i = 0; $i < 3; $i++) {
    $vector[$i] = 0.0;
}

# parse options.
if ($#ARGV < 0) {
    &help(); die;
}
for ( $i = 0; $i <= $#ARGV; $i++ ) {
    if ( $ARGV[$i] =~ /^\-h/ ) {
	&help();
	die;
    } elsif ( $ARGV[$i] =~ /^-c/ ) {
	$lattice = 0;
    } elsif ( $ARGV[$i] =~ /^-l/ ) {
	$lattice = 1;
    } elsif ( $ARGV[$i] =~ /^-param$/ ) {
	$i++;
	if ( open(PARAM, $ARGV[$i++]) ) {
	    while ( $line = <PARAM> ) {
		printf("%s", $line);
		@nums = split(' ', $line);
		# printf("nums: %s\n", $nums[0]);
		if ($nums[0] =~ /frame/) {
		    $frame = $nums[1];
		}
		elsif ($nums[0] =~ /label/) {
		    $label = $nums[1];
		    for ( $j = 2; $j <= $#nums; $j++ ) {
			$label = "$label $nums[$j]";
		    }
		}
		elsif ($nums[0] =~ /cluster_supercell/) {
		    $supercell_n[0]=$nums[1];
		    $supercell_n[1]=$nums[2];
		    $supercell_n[2]=$nums[3];
		}
		elsif ($nums[0] =~ /xsize/) {
		    $xsize = $nums[1];
		}
		elsif ($nums[0] =~ /dt/) {
		    $dt = $nums[1];
		}
		elsif ($nums[0] =~ /block/) {
		    $nz1 = $nums[1];
		    $nz2 = $nums[2];
		    $nz3 = $nums[3];
		}
		elsif ($nums[0] =~ /ninit/) {
		    $ninit = $nums[1];
		}
		elsif ($nums[0] =~ /view/) {
		    $theta = $nums[1];
		    $phi = $nums[2];
		}
		elsif ($nums[0] =~ /mag/) {
		    $magnf = $nums[1];
		}
		elsif ($nums[0] =~ /bond/) {
		    $nuca[$icon] = $nums[1];
		    $nuce[$icon] = $nums[2];
		    $rmin[$icon] = $nums[3];
		    $rmax[$icon] = $nums[4];
		    $rbond[$icon] = $nums[5];
		    $nbcol[$icon++] = $nums[6];
		}
		elsif ($nums[0] =~ /hlatom/) {
		    $hlatom[$nha] = $nums[1];
		    $hlz[$nha++] = $nums[2];
		}
		elsif ($nums[0] =~ /vector/) {
		    $vector[0] = $nums[1];
		    $vector[1] = $nums[2];
		    $vector[2] = $nums[3];
		}
	    }
	}
	close(PARAM);
    }
}

# input data file is the last parameter
$input = $ARGV[$#ARGV];
@nums = split('.log', $input);
$output = "$nums[0].gif";
printf("%s\n", $output);

open(INF, $input) || die "Can't open $ARGV[$#ARGV]: $!\n";
#open(STDOUT, ">/dev/null");
#qopen(STDOUT, ">>OUT");

# Skip to lattice vector section.
#while ( $line = <INF> ) { last if $line =~ /Ion kT =/; }
#@nums = split(' ', $line);
#$kT = $nums[3];
#while ( $line = <INF> ) { last if $line =~ /dt =/; }
#@nums = split(' ', $line);
#$dt = $nums[2];
while ( $line = <INF> ) { last if $line =~ /setup_basis()/; }
while ( $line = <INF> ) { last if $line =~ /latvec =/; }
$line = <INF>; @nums = split(' ', $line); 
$latt{0.0} = $nums[1]; $latt{0.1} = $nums[2]; $latt{0.2} = $nums[3];
$line = <INF>; @nums = split(' ', $line); 
$latt{1.0} = $nums[1]; $latt{1.1} = $nums[2]; $latt{1.2} = $nums[3];
$line = <INF>; @nums = split(' ', $line); 
$latt{2.0} = $nums[1]; $latt{2.1} = $nums[2]; $latt{2.2} = $nums[3];



$len[0] = sqrt($latt{0.0}*$latt{0.0} + $latt{0.1}*$latt{0.1} 
	       + $latt{0.2}*$latt{0.2});
$len[1] = sqrt($latt{1.0}*$latt{1.0} + $latt{1.1}*$latt{1.1} 
	       + $latt{1.2}*$latt{1.2});
$len[2] = sqrt($latt{2.0}*$latt{2.0} + $latt{2.1}*$latt{2.1} 
	       + $latt{2.2}*$latt{2.2});

printf("lattice vectors:\n");
printf("%10.6f %10.6f %10.6f\n", $latt{0.0},  $latt{0.1}, $latt{0.2});
printf("%10.6f %10.6f %10.6f\n", $latt{1.0},  $latt{1.1}, $latt{1.2});
printf("%10.6f %10.6f %10.6f\n", $latt{2.0},  $latt{2.1}, $latt{2.2});

# Skip to ion position section.
#while ( $line = <INF> ) { last if $line =~ /setup_ioninfo()/; }
#while ( $line = <INF> ) { last if $line =~ /nspecies/; }
#@nums = split(' ',$line); $nspecies = $nums[2];
$nspecies = 4;
printf("Found %d species \n", $nspecies);

while ( $line = <INF> ) {last if $line =~ /Setting Vpseudo/; }
for ( $sp = 0; $sp < $nspecies; $sp++ ) {
    while ( $line = <INF> ) { last if $line =~ /species/; }
    @nums = split(' ', $line); $line = <INF>;
    if ( $z_in[$sp] > 0 ) {
	$z[$sp] = $z_in[$sp];
    } else {
	$z[$sp] = $nums[5];
    }
    $natoms[$sp] = $nums[8];
    printf("Detect %d charge %d ions\n", $natoms[$sp], $z[$sp]);
}

$t = 0.0;
#$t_hyp = 0.0;
for($count = 0; $line = <INF>; $count++) {
    $t = $t + $dt;
    #while ( $line = <INF> ) { last if $line =~ /dv:/;}
    # printf("[%s]\n", $line);
    #@nums = split(' ',$line);
    # printf("[%s]\n", $nums[1]);
    #$t_hyp = $t_hyp + exp($nums[1]/$kT)*$dt;
    # printf("t_hyp: %f\n", $t_hyp);
    while ( $line = <INF> ) { last if $line =~ /iondyn()/;} 
# This is the frame filter
   if (($line =~ /iondyn()/) && !($count%$frame)) {
#    if (($line =~ /iondyn()/) ) {
	$ion_count = 0;
	for ( $sp = 0; $sp < $nspecies; $sp++ ) {
	    while ( $line = <INF> ) { last if $line =~ /species/; }
	    $line = <INF>;
	    for ($nat = 0; $nat < $natoms[$sp]; $nat++) {
		$line = <INF>;
#		printf("%s", $line);
		@nums = split(' ', $line);
		$nums[2] -= $vector[0];
		$nums[3] -= $vector[1];
		$nums[4] -= $vector[2];
#		printf("%15.2f%15.2f%15.2f\n", $nums[2], $nums[3], $nums[4]);
		for ($i = 2; $i < 5; $i++) {
		    $nums[$i] = $nums[$i] - int($nums[$i]);
		    if ($nums[$i] < 0.0) {
			$nums[$i]++;
		    }
		}
		$atomx{$sp.$nat} = $nums[2];
		$atomy{$sp.$nat} = $nums[3];
		$atomz{$sp.$nat} = $nums[4];
#radius of balsac atom balls
#		$r_out{$sp.$nat} = $z[$sp]/8.0;
                $r_out{$sp.$nat} = .25;
		$z_out{$sp.$nat} = $z[$sp];
		for ($i = 0; $i < $nha; $i++) {
		    if($hlatom[$i] == $ion_count) {
			$z_out{$sp.$nat} = $hlz[$i];
		    }
		}
		$ion_count++;
	    }
	}

#	if ($t_hyp == $t) {
	    $label_out = sprintf("%s  t = %f fs", 
				 $label, $t);
#	} else {
#	    $label_out = sprintf("%s  t = %8.3e s  t_hyp = %8.3e s %s", 
#				 $label, $t*10e-15, $t_hyp*10e-15);
#	}
	if ($lattice) {
# Produce lattice balsac input file.
	    $bfname = "lattice.balsac";
	    printf("Writing %s\n", $bfname);
	    open(BF, ">$bfname") || die "Can't open $bfname: $!\n";
# Card 1
	    printf(BF "%6d%2d%-20s%2d%s\n", 0, 6, "balsac.eps", 3, "YO,");
# Card 2
	    printf(BF "%5d%15.9f%-40s\n", 10, 1.0, $label_out);
# Card 3, lattice basis vectors
	    printf(BF "%15.9f%15.9f%15.9f\n", 
		   $latt{0.0}, $latt{0.1}, $latt{0.2});
	    printf(BF "%15.9f%15.9f%15.9f\n", 
		   $latt{1.0}, $latt{1.1}, $latt{1.2});
	    printf(BF "%15.9f%15.9f%15.9f\n", 
		   $latt{2.0}, $latt{2.1}, $latt{2.2});
# Card 4, number of basis-vectors
	    printf(BF "%5d%5d\n", $ion_count,0);
# Card 5, ionic positions in absolute coordinates, (NBG<0)
	    for ($sp = 0; $sp < $nspecies; $sp++) {
		for ($nat = 0; $nat < $natoms[$sp]; $nat++) {
		    printf(BF "%15.9f%15.9f%15.9f%15.9f%10d\n",
			   $atomx{$sp.$nat}*$latt{0.0} 
			   + $atomy{$sp.$nat}*$latt{0.1} 
			   +  $atomz{$sp.$nat}*$latt{0.2}, 
			   $atomx{$sp.$nat}*$latt{1.0} 
			   + $atomy{$sp.$nat}*$latt{1.1} 
			   +  $atomz{$sp.$nat}*$latt{1.2},
			   $atomx{$sp.$nat}*$latt{2.0} 
			   + $atomy{$sp.$nat}*$latt{2.1} 
			   +  $atomz{$sp.$nat}*$latt{2.2}, 
			   $r_out{$sp.$nat}, $z_out{$sp.$nat});
		}
	    }
# Card 6,  miller indices of stacking plane
	    printf(BF "%5d%5d%5d%5d\n",0,0,-1,0);
# Card 7
	    printf(BF "%5d%5d%5d%5d%10.3f%10.3f%10.3f%10.3f%10.3f\n",
		   $nz1, $nz2, $nz3, $ninit, 0.0, 0.0, 0.0, 0.0, 0.0);
# Card 8
	    printf(BF "%5d\n", 0);
# Card 10
	    printf(BF "%12.5f%12.5f%12.5f%12.5f%12.5f%12.5f\n",
		   $theta, $phi, 0.0, $magnf, 0.0, 0.0);
# Card 11
	    printf(BF "%5d%5d%5d%5d%5d%5d%5d%5d\n",
		   19, 20901, 1, 5, 0, 0, $icon, 0 );
# Card 11a
	    for ($i = 0; $i < $icon; $i++) { 
		printf(BF "%5d%5d%15.7f%15.7f%15.7f%5d\n", $nuca[$i], 
		       $nuce[$i], $rmin[$i], $rmax[$i], $rbond[$i], 
		       $nbcol[$i]);
	    }
	    close(BF);
	} else {
# Beginning of cluster balsac output file




# These lines fold the unit cell to make the supercell
if ($count == 0) {
	for ($i = 0; $i < 3; $i++) {
		$len[$i] = $len[$i]*$supercell_n[$i];
		}
	}

for ($i = 1; $i < $supercell_n[0]; $i++) {
	for ($sp = 0; $sp < $nspecies; $sp++) {
	$natt = $i*$natoms[$sp];
	for ($nat = 0; $nat < $natoms[$sp]; $nat++) {
		$atomx{$sp.$natt} = $atomx{$sp.$nat} + $i;
		$atomy{$sp.$natt} = $atomy{$sp.$nat};
		$atomz{$sp.$natt} = $atomz{$sp.$nat};
		$r_out{$sp.$natt} = $r_out{$sp.$nat};
		$z_out{$sp.$natt} = $z[$sp];
		$natt++;
		$ion_count++;
		}
		}

	}
for ($sp = 0; $sp < $nspecies; $sp++) {
	$natnew[$sp]=$natoms[$sp]*$supercell_n[0];
	}


for ($i = 1; $i < $supercell_n[1]; $i++) {
	for ($sp = 0; $sp < $nspecies; $sp++) {
	$natt = $i*$natnew[$sp];
	for ($nat = 0; $nat < $natnew[$sp]; $nat++) {
		$atomx{$sp.$natt} = $atomx{$sp.$nat};
		$atomy{$sp.$natt} = $atomy{$sp.$nat} + $i;
		$atomz{$sp.$natt} = $atomz{$sp.$nat};
		$r_out{$sp.$natt} = $r_out{$sp.$nat};
		$z_out{$sp.$natt++} = $z[$sp];
		$ion_count++;
		}
		}

	}
for ($sp = 0; $sp < $nspecies; $sp++) {
	$natnew[$sp]=$natnew[$sp]*$supercell_n[1];
	}


for ($i = 1; $i < $supercell_n[2]; $i++) {
	for ($sp = 0; $sp < $nspecies; $sp++) {
	$natt = $i*$natnew[$sp];
	for ($nat = 0; $nat < $natnew[$sp]; $nat++) {
		$atomx{$sp.$natt} = $atomx{$sp.$nat};
		$atomy{$sp.$natt} = $atomy{$sp.$nat};
		$atomz{$sp.$natt} = $atomz{$sp.$nat} + $i;
		$r_out{$sp.$natt} = $r_out{$sp.$nat};
		$z_out{$sp.$natt++} = $z[$sp];
		$ion_count++;
		}
		}

	}
for ($sp = 0; $sp < $nspecies; $sp++) {
	$natnew[$sp]=$natnew[$sp]*$supercell_n[2];
	}

		



# Add extra atoms
# These lines add atoms that lie on cell edges  
	    for ($sp = 0; $sp < $nspecies; $sp++) {
#		$natnew[$sp] = $natoms[$sp];
#		for ($nat = 0; $nat < $natoms[$sp]; $nat++) {
#		    if ($atomx{$sp.$nat} == 0.0) {
#			$atomx{$sp.$natnew[$sp]} = $atomx{$sp.$nat} + 1.0;
#			$atomy{$sp.$natnew[$sp]} = $atomy{$sp.$nat};
#			$atomz{$sp.$natnew[$sp]} = $atomz{$sp.$nat};
#			$r_out{$sp.$natnew[$sp]} = $r_out{$sp.$nat};
#			$z_out{$sp.$natnew[$sp]++} = $z[$sp];
#			$ion_count++;
#		    }
#		    if ($atomy{$sp.$nat} == 0.0) {
#			$atomx{$sp.$natnew[$sp]} = $atomx{$sp.$nat};
#			$atomy{$sp.$natnew[$sp]} = $atomy{$sp.$nat} + 1.0;
#			$atomz{$sp.$natnew[$sp]} = $atomz{$sp.$nat};
#			$r_out{$sp.$natnew[$sp]} = $r_out{$sp.$nat};
#			$z_out{$sp.$natnew[$sp]++} = $z[$sp];
#			$ion_count++;
#		    }
#		    if ($atomz{$sp.$nat} == 0.0) {
#			$atomx{$sp.$natnew[$sp]} = $atomx{$sp.$nat};
#			$atomy{$sp.$natnew[$sp]} = $atomy{$sp.$nat};
#			$atomz{$sp.$natnew[$sp]} = $atomz{$sp.$nat} + 1.0;
#			$r_out{$sp.$natnew[$sp]} = $r_out{$sp.$nat};
#			$z_out{$sp.$natnew[$sp]++} = $z[$sp];
#			$ion_count++;
#		    }
#		}
	    }


# Add atoms at corners of cell to keep cell fixed when things move 	    
	    $sp = $nspecies;
	    $natnew[$sp] = 0;
	    for ($i = 0; $i < 2; $i++) {
		for ($j = 0; $j < 2; $j++) {
		    for ($k = 0; $k < 2; $k++) {
			$atomx{$sp.$natnew[$sp]} = $supercell_n[0]*$i;
			$atomy{$sp.$natnew[$sp]} = $supercell_n[1]*$j;
			$atomz{$sp.$natnew[$sp]} = $supercell_n[2]*$k;
			$r_out{$sp.$natnew[$sp]} = 0.0;
			$z_out{$sp.$natnew[$sp]++} = 0.0;
			$ion_count++;
		    }
		}
	    }
	    $nspeciesnew = $nspecies + 1;	



	    
# Produce cluster balsac input file.
	    $bfname = "cluster.balsac";
	    printf("Writing %s\n", $bfname);
	    open(BF, ">$bfname") || die "Can't open $bfname: $!\n";
# Card  1.      (A60,I5)  TITLE,NBATCH
	    printf(BF "%-60s%5d\n", $label_out, 2);
# Card  1a.     (48A1)  (BOPT(I),I=1,NBATCH)
	    printf(BF "%-48s\n", "O,");
# Card  2.      (I5,3F15.9)  NKTOT,XTOT,YTOT,ZTOT	    
	    printf(BF "%5d%15.9f%15.9f%15.9f\n", $ion_count, $len[0], 
		   $len[1], $len[2]);
# Card(s)  3.   (4F15.9,I4)  X(I),Y(I),Z(I),RAD(I),NUC(I),I=1,NKTOT
	    for ($sp = 0; $sp < $nspeciesnew; $sp++) {
		for ($nat = 0; $nat < $natnew[$sp]; $nat++) {
		    printf(BF "%15.9f%15.9f%15.9f%15.9f%4d\n",
			   $atomx{$sp.$nat}*$latt{0.0} 
			   + $atomy{$sp.$nat}*$latt{0.1} 
			   +  $atomz{$sp.$nat}*$latt{0.2}, 
			   $atomx{$sp.$nat}*$latt{1.0} 
			   + $atomy{$sp.$nat}*$latt{1.1} 
			   +  $atomz{$sp.$nat}*$latt{1.2},
			   $atomx{$sp.$nat}*$latt{2.0} 
			   + $atomy{$sp.$nat}*$latt{2.1} 
			   +  $atomz{$sp.$nat}*$latt{2.2},
			   $r_out{$sp.$nat}, $z_out{$sp.$nat});
		}
	    }
# Card  4.      (5F12.5,12X,I5)  THETA,PHI,ROT,MAGNF,PERSP,NFOC
	    printf(BF "%12.5f%12.5f%12.5f%12.5f%12.5f%12s%5d\n",
		   $theta, $phi, 0.0, $magnf, 0.0, " ", 0);
# Card  5.      (8I5)  M,N,ICLC,IBFT,NSHX,NSHY,ICON,INITSV
	    printf(BF "%5d%5d%5d%5d%5d%5d%5d%5d\n",
		   19, 20901, 1, 5, 0, 0, $icon, 0 );
# Card(s)  5a.  (2I5,3D15.7,I5)  NUCA(I),NUCE(I),RMIN(I),RMAX(I),
#                                               RBOND(I),NBCOL(I),I=1,ICON
	    for ($i = 0; $i < $icon; $i++) { 
		printf(BF "%5d%5d%15.7f%15.7f%15.7f%5d\n", $nuca[$i], 
		       $nuce[$i], $rmin[$i], $rmax[$i], $rbond[$i], 
		       $nbcol[$i]);
	    }
# Card 1'       (I2,A72)  NFRM,FIL2O
	    printf(BF "%2d%s\n", 6, "balsac.eps"); 
	    close(BF);
	}	

	$name = sprintf("balsac%04d.gif", $count);
	printf("name: %s\n", $name);

	if ($lattice) {
	    system("./balsac -l lattice.balsac");
	} else {
	    system("./balsac -c cluster.balsac");
	}
	#exit();
	system("pstopnm -xsize $xsize  balsac.eps");
	open(STDOUT, ">balsac.ppm");
	system("pnmrotate -90 balsac.eps001.ppm");
	open(STDOUT, ">$name");
	system("ppmtogif balsac.ppm");
	open(STDOUT, ">>OUT");
	#open(STDOUT, ">/dev/null");
    }
}
system("whirlgif -loop -o $output balsac*gif");

system("rm *.balsac balsac.eps balsac.eps001.ppm balsac.ppm");
system("rm balsac????.gif");

close(INF);

#
# Subroutine help
#  print some help information.
#
sub help {
    $name = $0;
    $name =~ s/.*\///g ;
    printf(" Usage: $name [options] datafile.\n\n");
}








