#!/usr/bin/perl -ws

#use Data::Dump 'pp';

our($h,$help,$width,$fig_h,$fig_w,$level_w,$level_h,$zazor,
    $unit,$no_isomers,$null,$sort,$group,$G,$ZPE,$Edisp,$delete_null);
if ($h || $help) {
  print "Usage: $0 *.xyz > file.svg

Из xyz-файлов (2-я строка содержит энергию)
создается картинка для inkscape с уровнями и осью энергии.
Каждый уровень состоит из 3 прямоугольников, 2 маленьких по краям --
для связи уровней (инструмент \"Создавать линии соединения в диаграммах\").
Подписи уровней -- названия файлов *.xyz (берется последняя геометрия).

-unit  в каких единицах энергия в xyz-файлах. 
       По умолчанию a.u. (хартри). Другие значения -unit=kcal -unit=kJ
-kJ  с этой опцией ось энергии будет в кДж, а не в ккал.
-no_isomers  без этой опции (т.е. по умолчанию) в каждой группе изомеров 
             за 0 будет приниматься энергия самого стабильного в этой группе.
             С опцией -no_isomers будет единая шкала.
-ZPE    К энергиям прибавляются ZPE-поправки
        Во 2-й строке xyz-файлов д.б. подстрока типа ZPE 0.095000 (a.e.)
-G      Аналогично, термические поправки 
        Во 2-й строке xyz-файлов д.б. G(298.25) 33.98 (kcal/mol)
        Д.б. только одна из -G и -ZPE (т.к. в терм. поправки ZPE уже входит)
-Edisp  Дисперсионные поправки Grimme. См. скрипт Edisp.
        Во 2-й строке xyz-файлов д.б. подстрока типа Edisp -9.79 (kcal/mol)
-sort  Обычно на диаграмме структуры группируются по изомерам. С этой опцией
       будет сквозная сортировка по относительной энергии всех изомеров.
-null=regexp  за 0 будет приниматься энергия той структуры, имя файла которой 
              (только имя, без расширения) подходит regexp
-delete_null  не помещать на диаграмму \"нулевые\" уровни
-width=744  Ширина страницы A4 в пикселах (в 744 помещается 
            12 уровней шириной 40 и 15 уровней шириной 30)
-fig_w=0.8  Рисунок будет примерно на ширину страницы A4
-fig_h=0.6  Рисунок будет занимать ~0.6 страницы A4 по высоте
-level_w=40 Ширина уровня в пикселах
-level_h=2  Толщина уровня в пикселах
-zazor=2    Зазор м/д уровнем и надписью (толщин уровней)

-group  Уровни и их подписи будут сгруппированы. Не применяйте эту опцию, если 
        хотите потом соединять уровни.
";
  if ($help) {
    print "

Разработка этой программы поддержана РФФИ (грант 13-03-00427a).

";
  }
  exit;
}

$kJ = $kJ ? 4.18605 : 1;

if ($G && $ZPE) {
  die "-G and -ZPE are nonconsistent options\n";
}

$unit ||= 1;
$unit = 627.5 if $unit eq 'kcal';
$unit = 627.5*4.184 if $unit eq 'kJ';

my $ddd = qr/-?\d+(?:\.\d+)?/;

# Получаем названия файлов и энергии в массив @arr
my @arr;
my @mols = read_molden(@ARGV);
my %isomers = isomers(@mols);
my $min_E = 1e12;
my $max_E = -1e12;
foreach my $formula (keys %isomers) {
  my @mols = @{$isomers{$formula}};
  my $min_en = $mols[0][0]{Energy};
  if ($null) {
    foreach (@mols) {
      $_->[0]{File} =~ m/$null/ and $min_en = $_->[0]{Energy} and last;
    }
  }
  foreach my $mol (@mols) {
    (my $name = $mol->[0]{File}) =~ s/\.xyz(ppm)?$//;
    push @arr, [$name, 
                ($mol->[0]{Energy}-$min_en)/$unit*627.5*$kJ,
                $mol->[0]{Symmetry}];
    $min_E = $arr[-1][1] if $min_E > $arr[-1][1];
    $max_E = $arr[-1][1] if $max_E < $arr[-1][1];
  }
}
#pp @arr;
my $delta_E = $max_E - $min_E;
#warn "$delta_E = $max_E - $min_E\n";

# Сортируем массив по возрастанию энергий
@arr = sort {$a->[1] <=> $b->[1]} @arr if $sort;
#@arr = sort {$a->[0] cmp $b->[0]} @arr;

# Всякие параметры
$width ||= 744; # Ширина страницы в пикселах
$fig_h ||= 0.6; # Рисунок будет занимать ~0.6 страницы A4 по высоте
$fig_w ||= 0.8; # Рисунок будет примерно на ширину страницы A4
$level_w ||= 40; # Ширина уровня в пикселах
$level_h ||= 2; # Толщина уровня в пикселах
$zazor ||= 2; # Зазор м/д уровнем и надписью (толщин уровней)
$level_w = $level_h if $level_w < 3*$level_h;

my $height = sprintf '%d', $width*sqrt(2); # page A4
#my $x0 = int($width/10); # 10% слева
my $x0 = 50;
my $y0 = int($height*(1+$fig_h*$min_E/$delta_E))-50; # 10% снизу
#warn( ($height-$y0), "\n");
my $px_kJ = $height*$fig_h/$delta_E;

# Вычисляем положения уровней в пикселах и добавляем в массив
my $dx = 0;
foreach (@arr) {
  push @{$_}, $x0 + $dx, int($y0 - $_->[1]*$px_kJ);
  $dx += $level_w+5*$level_h;
  $dx = 0 if $dx > $width*$fig_w;
}

#pp @arr;

# Шапка svg
print <<PPP;
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<svg
   xmlns:svg="http://www.w3.org/2000/svg"
   xmlns="http://www.w3.org/2000/svg"
   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
   width="$width"
   height="$height"
   id="svg2"
   sodipodi:version="0.32"
   inkscape:version="0.46"
   inkscape:output_extension="org.inkscape.output.svg.inkscape"
   version="1.0">
  <defs />
  <sodipodi:namedview
     inkscape:document-units="px"
     inkscape:window-width="1030"
     inkscape:window-height="894"
     inkscape:zoom="1.2"
     inkscape:cx="372"
     inkscape:cy="266"
     inkscape:window-x="68"
     inkscape:window-y="54"
     inkscape:current-layer="layer1" />
  <g
     inkscape:label="Layer 1"
     inkscape:groupmode="layer"
     id="layer1">
PPP

# Уровни и подписи
my $text_ID = 1;
my $rect_ID = 1;
my $path_ID = 1;
foreach (reverse @arr) {
  next if $delete_null && $null && $_->[0] =~ /$null/;
  printf("    <g id=\"g%04d\">\n", $group++) if $group;
  print '
    <rect width="',$level_h,           '" height="',$level_h,'" x="',$_->[3],                  '" y="',$_->[4],'" id="rect', $rect_ID++, '" />
    <rect width="',$level_w-2*$level_h,'" height="',$level_h,'" x="',$_->[3]+$level_h,         '" y="',$_->[4],'" id="rect', $rect_ID++, '" />
    <rect width="',$level_h,           '" height="',$level_h,'" x="',$_->[3]+$level_w-$level_h,'" y="',$_->[4],'" id="rect', $rect_ID++, '" />
    <text
       style="font-size:10px;font-style:normal;font-family:DejaVu Sans"
       x="',$_->[3],'" y="',$_->[4]-$level_h*$zazor,'" id="text', $text_ID++,'">
       ',$_->[0],'
    </text>
  ';
  if ($_->[2] && $_->[2] ne 'C1') {
    print '
    <text
       style="font-size:12px;font-style:oblique;fill:#0000ff;font-family:DejaVu Sans"
       x="',$_->[3],'" y="',$_->[4]-$level_h*$zazor-12,'">
       ',$_->[2],'
    </text>
  ';
  }
  print "    </g>\n" if $group;
}

#pp @arr;
# Ось энергии и ее подпись
my ($lx0,$ly0,$lx1,$ly1);
$lx0 = int($x0 - $level_w/2);
$ly0 = $arr[0][4]+10;
$lx1 = $lx0;
$ly1 = $arr[-1][4]-10;
print "    <g id=\"axis\">\n";
print '
    <path
       style="stroke:#000000;stroke-width:1px"
       d="',"M $lx0,$ly0 C $lx1,$ly1 $lx1,$ly1 $lx1,$ly1",'" id="path', $path_ID++,'" />
      <text
         style="font-size:12px;font-style:normal;font-family:DejaVu Sans"
         x="',$lx1+5,'" y="',$ly1+3,'" id="text', $text_ID++,'">
         ',$kJ==1 ? 'kcal' : 'kJ','
      </text>
'; 

# Шаг тиков
my $step = $delta_E/10;
sprintf("%e",$step) =~ /^(\d).*e(.*)/;
if    ($1 >= 5) { $step = 10*10**$2 }
elsif ($1 >= 2) { $step =  5*10**$2 }
else            { $step =  2*10**$2 }

# Тики
$lx0 -= 6; # Ширина тика
# 0 и положительные энергии
my $en_tic = 0;
while ($en_tic < $max_E) {
  $ly0 = int($y0 - $px_kJ*$en_tic);
  $ly1 = $ly0;
  print '
      <path
         style="stroke:#000000;stroke-width:1px"
         d="',"M $lx0,$ly0 C $lx1,$ly1 $lx1,$ly1 $lx1,$ly1",'" id="path', $path_ID++,'" />
      <text
         style="font-size:12px;font-style:normal;font-family:DejaVu Sans;text-align:end;text-anchor:end"
         x="',$lx0-2,'" y="',$ly0+3,'" id="text', $text_ID++,'">
         ',$en_tic,'
      </text>
  '; 
  $en_tic += $step;
}
# Отрицательные энергии
$en_tic = -$step;
while ($en_tic > $min_E) {
  $ly0 = int($y0 - $px_kJ*$en_tic);
  $ly1 = $ly0;
  print '
      <path
         style="stroke:#000000;stroke-width:1px"
         d="',"M $lx0,$ly0 C $lx1,$ly1 $lx1,$ly1 $lx1,$ly1",'" id="path', $path_ID++,'" />
      <text
         style="font-size:12px;font-style:normal;font-family:DejaVu Sans;text-align:end;text-anchor:end"
         x="',$lx0-2,'" y="',$ly0+3,'" id="text', $text_ID++,'">
         ',$en_tic,'
      </text>
  '; 
  $en_tic -= $step;
}
print "    </g>\n";

# Конец svg
print "
  </g>
</svg>
";

sub read_molden0 {
  my $inp = shift;
  $inp = 'STDIN' unless $inp;
  return () if eof($inp);
  my ($N,$energy,$symmetry,@atom,@x,@y,@z,@charge);
  # 1-st string
  <$inp>=~/^\s*(\d+)\s*$/ ? $N = $1 : return ();
  # 2-nd string
  my $str2 = <$inp>;
  ($energy) = $str2=~/(-?\d+\.\d+)/;
  ($symmetry) = $str2=~/Symmetry\s+(\S+)/;
  $symmetry ||= '';
  # Geometry
  for ($i=1;$i<=$N;$i++) {
    my $str = <$inp>;
    #($atom[$i],$x[$i],$y[$i],$z[$i],$charge[$i]) = split ' ', $str;
    #return undef unless $atom[$i]=~/^\w\w?/ &&
    #                    $x[$i]=~/^-?\d+\.\d+/ &&
    #                    $y[$i]=~/^-?\d+\.\d+/ &&
    #                    $z[$i]=~/^-?\d+\.\d+/;
    next;
  }
  return [$energy,$symmetry],\@atom,\@x,\@y,\@z,\@charge;
}

######################################################################
# Get array of molecules
# Return hash of isomers in which keys is formulas 
# and isomer is ref to array of molecules with this formula. 
sub isomers {
  my @mols = @_;
  my %isomers;
  foreach my $mol (@mols) {
    my $formula = $no_isomers ? 'no_isomers' : get_formula($mol);
    push @{$isomers{$formula}}, $mol;
  }
  foreach (values %isomers) {
    $_ = [sort {($a->[0]{Energy}||1e6) <=> ($b->[0]{Energy}||1e6)} @$_];
  }
  return %isomers;
}

sub get_formula {
	my @mol = @{$_[0]};
	my %f;
	my $formula = '';
	$f{ucfirst lc $_->[0]}++ foreach @mol[1..$#mol];
	foreach (sort by_Hill keys %f) {
	  $formula .= $_ . ($f{$_}==1 ? '' : $f{$_});
	}
	return $formula;
}

sub by_Hill {
	return -1 if uc($a) eq 'C' and uc($b) ne 'C';
	return  1 if uc($a) ne 'C' and uc($b) eq 'C';
	return -1 if uc($a) eq 'H' and uc($b) ne 'H';
	return  1 if uc($a) ne 'H' and uc($b) eq 'H';
	uc($a) cmp uc($b);
}

sub by_formula {
  my $aa = $a;
	my $bb = $b;
	$aa =~ s/([A-Z](?:[a-z]*))(?!\d)/${1}1/g;
  $bb =~ s/([A-Z](?:[a-z]*))(?!\d)/${1}1/g;
	$aa =~ s/(\d+)/sprintf("%03d",$1)/eg;
	$bb =~ s/(\d+)/sprintf("%03d",$1)/eg;
	$aa cmp $bb;
}
######################################################################

# Читает xyz. Параметр - имя xyz-файла.
# Возвращает массив найденных молекул.
sub read_molden {
  my @mols;
  foreach my $file (@_) {
    my $mol;
    $mol->[0]{File} = $file;
    open L, '<', $file or do {warn "Can't open $file: $!\n"; return};
    return undef if eof(L);
    # 1-st string
    if (<L>=~/^\s*(\d+)\s*$/) {
      $N = $1;
    } else {
      return undef;
    }
    # 2-nd string
    $_ =<L>;
    ($mol->[0]{Energy}) = m/(-?\d+\.\d+)/;
    if ($G) {
      /\sG\(($ddd)\)\s+($ddd)/ && ($mol->[0]{Energy} += $2/627.51);
    }
    if ($ZPE) {
      /\sZPE\s+($ddd)/ && ($mol->[0]{Energy} += $1);
    }
    if ($Edisp) {
      /\sEdisp\s+($ddd)/ && ($mol->[0]{Energy} += $1/627.51);
    }
    ($mol->[0]{Symmetry}) = m/Symmetry\s+(\w+)/;
    # Geometry
    for ($i=1;$i<=$N;$i++) {
      my ($atom,$x,$y,$z,$ppm) = split ' ', <L>;
      return undef unless $atom=~/^\w\w?/ &&
                          $x=~/^-?\d+\.\d+/ &&
                          $y=~/^-?\d+\.\d+/ &&
                          $z=~/^-?\d+\.\d+/;
      $mol->[$i] = [$atom,$x,$y,$z,$ppm];
    }
    close L;
    push @mols, $mol;
  } 
  return @mols;
}