#!/usr/bin/perl
# the above is the standard linux/unix path for perl: change if your computer is nonstandard

use POSIX;

sub print_help{
    print STDERR <<EOF;

Syntax:  xyz_lattice_translate [options] input_xyzfile

generate an xyz file with a set of Bravais translations of the atoms in input

Options:
-n number_of_translations_along_each_basis_vectors (default: $nmax)
-n1 number_of_translations_along_a1 (default: $nmax)
-n2 number_of_translations_along_a2 (default: $nmax)
-n3 number_of_translations_along_a3 (default: $nmax)
-a1 a1x a1y a1z
-a2 a2x a2y a2z
-a3 a3x a3y a3z
-sc avalue
-fcc avalue
-bcc avalue
-hex avalue cvalue
    
							by Nick 13 Mar 2023
EOF
    return;
}

$BIG=99.e99;
$nmax=2;
$nmax1=$nmax;
$nmax2=$nmax;
$nmax3=$nmax;

# default lattice translations:

# this is face-centered orthorombic (fco):
$a=15.205;
$b=15.179;
$c=13.503;
@a1=($a/2,0,$c/2);
@a2=($a/2,$b/2,0);
@a3=(0,$b/2,$c/2);


while($p=shift(@ARGV)){
    if($p=~/^[^\-]/){ # all command line options start with "-", otherwise it a filename
        if(-e $p){
            push(@inputfiles,$p);
        }else{
            print STDERR "ERROR: file $p does not exist\n";
        }
        next;
    }
    if($p eq "-n"){ # number of translations
	$nmax=shift(@ARGV)+0;
	$nmax>=1 or die "Invalid number of translations $nmax\n";
	$nmax1=$nmax;
	$nmax2=$nmax;
	$nmax3=$nmax;
	next;
    }
    if($p eq "-n1"){ # number of translations
	$nmax1=shift(@ARGV)+0;
	$nmax1>=1 or die "Invalid number of translations $nmax1\n";
	next;
    }
    if($p eq "-n2"){ # number of translations
	$nmax2=shift(@ARGV)+0;
	$nmax2>=1 or die "Invalid number of translations $nmax2\n";
	next;
    }
    if($p eq "-n3"){ # number of translations
	$nmax3=shift(@ARGV)+0;
	$nmax3>=1 or die "Invalid number of translations $nmax3\n";
	next;
    }
    if($p eq "-a1"){
	$x=shift(@ARGV)+0.; $y=shift(@ARGV)+0.; $z=shift(@ARGV)+0.;
	@a1=($x,$y,$z);
	next;
    }
    if($p eq "-a2"){
	$x=shift(@ARGV)+0.; $y=shift(@ARGV)+0.; $z=shift(@ARGV)+0.;
	@a2=($x,$y,$z);
	next;
    }
    if($p eq "-a3"){
	$x=shift(@ARGV)+0.; $y=shift(@ARGV)+0.; $z=shift(@ARGV)+0.;
	@a3=($x,$y,$z);
	next;
    }
    if($p eq "-sc"){
	$avalue=shift(@ARGV)+0.;
	@a1=($avalue,0.,0.);
	@a2=(0.,$avalue,0.);
	@a3=(0.,0.,$avalue);
	next;
    }
    if($p eq "-fcc"){
	$avalue=shift(@ARGV)+0.;
	@a1=(0.,0.5*$avalue,0.5*$avalue);
	@a2=(0.5*$avalue,0.,0.5*$avalue);
	@a3=(0.5*$avalue,0.5*$avalue,0.);
	next;
    }
    if($p eq "-bcc"){
	$avalue=shift(@ARGV)+0.;
	@a1=(-0.5*$avalue,0.5*$avalue,0.5*$avalue);
	@a2=(0.5*$avalue,-0.5*$avalue,0.5*$avalue);
	@a3=(0.5*$avalue,0.5*$avalue,-0.5*$avalue);
	next;
    }
    if($p eq "-hex"){
	$avalue=shift(@ARGV)+0.;
	$cvalue=shift(@ARGV)+0.;
	@a1=($avalue,0.,0.);
	@a2=(-0.5*$avalue,sqrt(3.)*0.5*$avalue,0.);
	@a3=(0.,0.,$cvalue);
	next;
    }
    if($p eq "-h"){ # help
	print_help;
	die "\n";
    }
    print_help;
    die "invalid option $p\n";
}

$nm1=-POSIX::floor(($nmax1-1)/2);
$nM1=POSIX::floor($nmax1/2);
$nm2=-POSIX::floor(($nmax2-1)/2);
$nM2=POSIX::floor($nmax2/2);
$nm3=-POSIX::floor(($nmax3-1)/2);
$nM3=POSIX::floor($nmax3/2);


if($#inputfiles<0){
    push(@inputfiles,"-");      #       adds stdio for lack of arguments
}



while ($filename=shift(@inputfiles)){
    $natom=0;
    $count=0;
    $totp=0;
    $frames=0;
    open(INPUT, $filename) || die "Cannot open $filename\n";
    while(<INPUT>){
	s/^\s*//;	# remove leading spaces
	@li = split('\s+');
	if($#li==0 && $count!=1){
	    $count=0;
	    $frames++;
	    @ato=();
	    @atpo=();
	    $natom=$li[0]+0;
	}
	if($count==1){
	    $comment=$_;
	    $comment=~s/\n//;;
	}
	if($count>1){
	    $totp++;
	    push(@ato,$li[0]);
	    for($i=0;$i<3;$i++){
		$ri=$li[$i+1]+0.;
		push(@atpo,$ri);
	    }
	}
	$count++;
	if($natom+2==$count){
	    print $natom*$nmax1*$nmax2*$nmax3,"\n";
	    print "#".$comment." - $natom atoms translated ",$nmax1,"*",$nmax2,"*",$nmax3," times\n";
#	    print "QUII",@ato,"\n";
#	    print "QUII",@atpo,"\n";
	    for($n1=$nm1;$n1<=$nM1;$n1++){
		for($n2=$nm2;$n2<=$nM2;$n2++){
		    for($n3=$nm3;$n3<=$nM3;$n3++){			
			for($j=0;$j<$natom;$j++){
			    print $ato[$j];
			    for($i=0;$i<3;$i++){
				print "\t",$atpo[$i+3*$j]+$n1*$a1[$i]+$n2*$a2[$i]+$n3*$a3[$i];
			    }
			    print "\n";
			}
		    }
		}
	    }
	}
    }
    close(INPUT);
}