#!/usr/bin/perl

my $atatdir;
BEGIN {
    open(my $cf,"<",join('/',$ENV{"HOME"},".atat.rc")) or die "Can't open ~/.atat.rc\n";
    $atatdir = <$cf>;
    chomp $atatdir;
    $atatdir =~ s/.*atatdir\s*=\s*//g;
}

my @relaxpath;

$relaxpath{"bcc-fcc"}=q{
1 1 {1+x*(sqrt(2)-1)} 90 90 90
-0.5 0.5 0.5
0.5 -0.5 0.5
0.5 0.5 -0.5
0 0 0 {A}
};

$relaxpath{"bcc-hcp"}=q{
#const a=1; b=sqrt(2); c=sqrt(2);
#const d=x/12.; th=109.47*(1-x)+x*120;
#const deg=3.14159/180; e=sqrt(sin(th*deg/2)/cos(th*deg/2)/sqrt(2));
{a/e} {b*e} {c} 90 90 90
0.5 -0.5 0
0.5 0.5 0
0 0 1
0 {0.25+d} 0.25 {A}
0 {0.75-d} 0.75 {A}
};

if ( (@ARGV == 0) || ($ARGV[0] eq "-h") ) {
  print q{Usage: symbrklib [-s fraction] [element] [initial_structure] [final_structure]
Currently, the initial and final structures can be one of bcc,fcc,hcp.
The -s option (optional) indicates to just generate a small distortion (a given fraction along the instability path).
This command is intended to be run before robustrelax_xxxx
};
exit;
}

my $x=1;
my $x0=0;

my @myargv=@ARGV;
if ($ARGV[0] eq "-s") {
  $x=$myargv[1];
  shift @myargv;
  shift @myargv;
}

my $pathname;
if ($myargv[1] gt $myargv[2]) {
  $pathname=$myargv[2] . '-' . $myargv[1];
  $x=1-$x;
  $x0=1-$x0;
}
else {
  $pathname=$myargv[1] . '-' . $myargv[2];
}

if ( $relaxpath{$pathname} eq "" ) {
  print "Unkown instability path: $pathname\n";
  exit;
}

{
    open(my $fh,">","path.sspp");
    print $fh "#const x=$x;\n";
    print $fh qq{#const A="$myargv[0]";};
    print $fh $relaxpath{$pathname};
}
system("sspp < path.sspp > path.in; nntouch -l=path.in -r=$atatdir/data/radii.in > str_hint.out");

{
    open(my $fh,">","path.sspp");
    print $fh "#const x=$x0;\n";
    print $fh qq{#const A="$myargv[0]";};
    print $fh $relaxpath{$pathname};
}
system("sspp < path.sspp > path.in; nntouch -l=path.in -r=$atatdir/data/radii.in > str_sup.out");
system("rm -f path.in path.sspp");
