Лаборатория изучения механизмов органических реакций

levels2html

#!/usr/bin/perl -ws

use Data::Dump 'pp';

our($dpi,$print_levels,$index,$bg,$no_png,$jmol_dir,$G,$ZPE,$Edisp,
    $h,$help,$dir,$unit,$null,$delete_null,$no_highlight);
if ($h or $help) {
  print "Usage: $0 fig.svg

Требуется бинарник inkscape. Для работы апплета jmol нужен www-сервер.

fig.svg получаем xyz2levels *.xyz > fig.svg.

   Опции:

-dpi   Разрешение. Умолчаемое 108 (масштаб в 1.2 раза больше, чем отображается 
       в inkscape).
-dir=directory  Если *.xyz или *.xyzppm расположены в директории directory
       (относительно соответствующего index.html), то эта опция добавит путь
       в file.htm. Default -dir=.
-unit  В каких единицах энергия в *.xyz.
       Default a.u. Возможны -unit=kcal  -unit=kJ
-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)
-null=regexp  за 0 будет приниматься энергия той структуры, имя файла которой 
       (только имя, без расширения) подходит regexp. 
-delete_null  не помещать в htm-файл информацию о \"нулевых\" уровнях 
        и не создавать соответствующие html-файлы. Это нужно, если эти уровни
        на завершающем этапе удаляются с диаграммы.
-print_levels  Печатает координаты (x и y) уровней в png-файле и их надписи
-index  С этой опцией будет создан полный html-файл, готовый для запуска в 
        браузере. Без нее - только кусок тегои <MAP> в файле fig.htm (для
        вставки в свой темплат). Просто -index эквивалентно -index=index.html
        -index forces -no_highlight.
-bg[=color] do color background in png file (more size of image)
            simple -b is synonym for -b=white
-no_png     Doesn't create png file
-jmol_dir   Directory of jmol (relatively DocumentRoot)
            -jmol_dir='/jmol' by default
-no_highlight не использовать подсветку уровней на диаграмме (для подсветки 
              нужен mapper.js (http://www.netzgesta.de/mapper/))
";
  print "-help  дополнительная информация\n" unless $help;
  if ($help) {
    print "
   Для правильной работы программы *.xyz или *.xyzppm должны находится 
в текущей директории (или в директории, задаваемой опцией -dir=), 
и их имена должны совпадать с подписями уровней. Чтобы программа правильно
определила координаты уровней, уровни не должны входить в состав групп.
   Берем подходящий index.html и меняем в нем соответствующий кусок на
содержимое file.htm и имя png-файла.  
   Программа также создаст на каждый xyz-файл свой html-файл, который 
будет запускаться в браузере по клику на соответствующий уровень
(необходим ява-апплет jmol). 
   Если в xyz-файле есть хим.сдвиги (пятой колонкой), то будет возможность 
их посмотреть (как метки атомов), активировав чек-бокс внизу всплывающего 
окна.
   Если вместе с *.xyz (*.xyzppm) обнаружатся *.mol с теми же именами, 
то из html-файлов будут вызываться *.mol. 
Делать mol-файлы, в которых будут показаны кратные связи, можно командой
babel *.xyz -omol -m

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

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

my $inkscape_x = 'inkscape';

$dpi ||= 108; # 
$dir ||= '.';

$index = 'index.html' if ($index && $index eq '1');
$jmol_dir ||= '/jmol'; # from DocumentRoot

$no_highlight = 1 if $index;

if ($delete_null && ! defined($null)) {
  die "You must define -null option\n";
}

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

my $file_svg = $ARGV[0];
(my $filename = $file_svg) =~ s/\.svg$//;

my @levels = svg2levels_inkscape($dpi,$file_svg);
warn $#levels+1, " levels\n\n";
#pp @levels; exit;

my @mols;
foreach my $level (@levels) {
  my $name = $level->[1][0];
  my $ext;
  
  if (-f "$name.xyzppm") {
    $ext = 'xyzppm';
    push @mols, read_molden("$name.xyzppm");
  }
  elsif (-f "$name.xyz") {
    $ext = 'xyz';
    push @mols, read_molden("$name.xyz");
  }
  else {
    print STDERR "No $name.xyzppm or $name.xyz\n";
    next;
  }
  
  @{$mols[-1][0]}{qw(id rx ry rw rh)} = @{$level->[0]};
  @{$mols[-1][0]}{qw(name tx ty tw th)} = @{$level->[1]};
  $mols[-1][0]{'ext'} = $ext;
}
#pp @mols; exit;

my %isomers = isomers(@mols);
#pp %isomers; exit;

my $factor = 627.5;
$factor = 1 if $unit && $unit eq 'kcal';
$factor = 1/4.184 if $unit && $unit eq 'kJ';

my $mapname = ($dir eq '.') ? $filename : $dir;
$mapname =~ s/_/UNDERSCORE/g;
my $htm = $index ? $index : "$filename.htm";
open HTM, '>', $htm or die "Can't write $htm: $!\n";
warn "Create $htm\n\n";
if ($index) {
  print HTM <<HTM;
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">
<html>

<head>
  <meta http-equiv="content-type" content="text/html; charset=utf-8">
  <title>$filename</title>
  <SCRIPT TYPE="text/javascript">
    function popup(mylink, windowname)
    {
      if (! window.focus)return true;
      var x = window.open(mylink.href, windowname, 'width=100 height=100');
      x.resizeTo(380, 410);
      return false;
    }
  </SCRIPT>

</head>

<body>
<h1>Energy levels</h1>
<small>
  Click on a level title to download xyz file<br>
  Click on energy level to view 3D structure (run jmol java-applet)<br>
</small>
<br>
HTM
}
print HTM <<HTM;
<MAP NAME="$mapname">
HTM


foreach my $formula (keys %isomers) {
	my @mols = @{$isomers{$formula}};
  my $min_en = $mols[0][0]{Energy};
  if (defined $null) {
    foreach my $mol (@mols) {
      if ($mol->[0]{name} =~ /$null/) {
        $min_en = $mol->[0]{Energy};
        last;
      }
    }
  }
  foreach my $mol (@mols) {
    my ($name,$ext,$en) = ($mol->[0]{name},$mol->[0]{ext},$mol->[0]{Energy});
	  next if $delete_null && $name =~ /$null/;
    my $coord  = ($mol->[0]{rx}) . ',' . 
                 ($mol->[0]{ry}-3) . ',' . 
                 ($mol->[0]{rx}+$mol->[0]{rw}) . ',' . 
                 ($mol->[0]{ry}+$mol->[0]{rh}+3);
    my $coord2 = ($mol->[0]{tx}-2) . ',' . 
                 ($mol->[0]{ty}-2) . ',' . 
                 ($mol->[0]{tx}+$mol->[0]{tw}+2) . ',' . 
                 ($mol->[0]{ty}+$mol->[0]{th}+2);
 	  $en = sprintf '%.2f',($en-$min_en)*$factor;
    #print "$formula\t $name\t $en\n";
    
    print HTM <<HTM;
<AREA HREF="$dir/$name.html" TITLE="$en kcal/mole" ALT="$en" COORDS="$coord" SHAPE=RECT onClick="return popup(this, '$name')">
<AREA HREF="$dir/$name.$ext" TITLE="Download $name.$ext" ALT="$en" COORDS="$coord2" SHAPE=RECT>
HTM

	  $ext = 'mol' if -f "$name.mol";
    open F, '>', "$name.html" or die "Can't open : $!\n";
	  warn "Create $name.html\n";
    print F <<FFF;
<head>
  <META HTTP-EQUIV="CACHE-CONTROL" CONTENT="NO-CACHE">
  <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
  <title>$name</title>
  <script src="$jmol_dir/Jmol.js" type="text/javascript"></script>
</head>
<body>
  <pre>$name    $en kcal/mole</pre>
  <form>
    <script type="text/javascript">
      jmolInitialize("$jmol_dir");
      jmolApplet(["100%","90%"], "load $name.$ext");
FFF
	  my @ppm = map {$mol->[$_][4]} 1..$#{$mol};
    #pp @ppm;
    if (grep {$_} @ppm) {
      print F "      jmolBr();\n      jmolCheckbox('\\\n";
      for (my $i=1; $i<=@ppm; $i++) {
        print F "      select (*)[$i];set labeloffset 5 0; label \"$ppm[$i-1]\"; \\\n" if $ppm[$i-1];
      }
      print F "      ', 'select all; label off', \"show chemical shifts\");\n";
    }
    print F <<FFF;
    </script>
  </form>
</body>
FFF
    close F;
  }
}

my $highlight = $no_highlight ? '' : 'class="mapper"';
print HTM <<HTM;
<AREA SHAPE=DEFAULT NOHREF>
</MAP>
<DIV><IMG $highlight SRC="$dir/$filename.png" border=0 USEMAP="#$mapname" hspace="5" vspace="5"></DIV>
<br>
HTM
if ($index) {
  print HTM <<HTM;
</body>

</html>
HTM
}
close HTM;

# Make png
if (! $no_png) {
  my $file_svg_for_png = $file_svg;
  if ($delete_null) {
    warn "\nGenerate svg without null level.\n";
    $file_svg_for_png = "PNG_$file_svg";
    open IN, '<', $file_svg or die "Can't open $file_svg: $!\n";
    local $/;
    my $svg = <IN>;
    close IN;
    open OUT, '>', $file_svg_for_png or die "Can't write $file_svg_for_png: $!\n";
    foreach my $l (@levels) {
      if ($l->[1][0] =~ /$null/) {
        my $text = $l->[1][0];
        my ($id1,$id2,$id3) = $l->[0][0]=~/(rect\S+?)(rect\S+?)(rect\S+)/;
        #warn "$text,$id1,$id2,$id3\n";
        $svg =~ s/<text[^>]*>\s*(<tspan[^>]*>)?\s*\Q$text\E\s*(<\/tspan>)?\s*<\/text>//s;
#        $svg =~ s/<text[^>]*>\s*\s*\Q$text\E\s*\s*<\/text>//s;
        foreach my $id ($id1,$id2,$id3) {
          $svg =~ s/<rect(.(?!id=))*(.(?=id=))id="$id"[^\/]*\/>//s;
        }
      }
    }
    print OUT $svg;
    close OUT;
  }
  
  $bg = 'white' if $bg && $bg eq '1';
    
  my @args = ($inkscape_x, 
    "--export-dpi=$dpi", "--export-area-drawing", "--export-png=$filename.png",
    "--without-gui", $file_svg_for_png);
    splice(@args, -2, 0, "--export-background=$bg") if $bg;
  warn "\nRun `@args`\n";
  system(@args) == 0 or warn "System\n@args\nfailed: $?\n";
  unlink $file_svg_for_png if $delete_null;
}

if ($index) {
  warn "\nTest if jmol directory exists (DocumentRoot$jmol_dir) and run $index in browser\n";
} else {
  warn "\nInsert $htm in your html-template\n";
}

######################################################################
# 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 = 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 {
  open L, '<', shift or do {warn "Can't open $_[0]: $!\n"; return};
	return undef if eof(L);
  my @mol;
  # 1-st string
	if (<L>=~/^\s*(\d+)\s*$/) {
		$N = $1;
	} else {
	  return undef;
	}
	# 2-nd string
	my $line = <L>;
  ($mol[0]{Energy}) = $line=~/(-?\d+\.\d+)/;
  if ($G) {
    $line=~/\sG\(($ddd)\)\s+($ddd)/ && ($mol[0]{Energy} += $2/627.51);
  }
  if ($ZPE) {
    $line=~/\sZPE\s+($ddd)/ && ($mol[0]{Energy} += $1);
  }
  if ($Edisp) {
    $line=~/\sEdisp\s+($ddd)/ && ($mol[0]{Energy} += $1/627.51);
  }
# 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;
  return \@mol;
}

# Печатает xyz. Параметры -- список молекул, в конце м.б. имя файла.
# Если последний элемент списка - имя файла (не ссылка на массив),
# то печать в этот файл, иначе - на stdout
sub write_molden {
	my $fh = \*STDOUT;
  if (ref($_[-1]) ne 'ARRAY') {
    my $file = pop @_;
    open $fh, '>', $file or die "Can't write to $file: $!\n";
  }
  
  foreach my $mol (@_) {
		my $N = $#{$mol};
    print $fh " $N\n";
		print $fh " Energy $mol->[0]{Energy} " if $mol->[0]{Energy};
		print $fh " Symmetry $mol->[0]{Symmetry} " if $mol->[0]{Symmetry};
		print $fh "\n";
		for (my $i=1; $i<=$N; $i++) {
      my ($atom,$x,$y,$z,$ppm) = @{$mol->[$i]};
      printf $fh " %-2s %12.8f %12.8f %12.8f", $atom, $x, $y, $z;
      printf $fh uc($atom) eq 'H' ? " %10.3f" : " %9.2f"  , $ppm if $ppm;
      print $fh "\n";
		}
	}
  #close $fh;
}

sub svg2levels_inkscape {
 #Получает dpi, имя svg-файла.
 #Через inkscape извлекает координаты и размеры всех элементов.
 #Парсит svg-файл и находит ближайшие к уровням надписи (имена файлов).
 #Возвращает массив ([ [IDIDID,X,Y,width,height], [text,X,Y,width,height] ], ...)
 #где IDIDID - ID прямоугольников, составляющих уровень, text - ближайшая надпись.
 #Расстояние определяется по серединам левых сторон.

  my ($dpi,$svg_file) = @_; # input
  my @levels; # output
  
  # Через inkscape получает координаты и размеры всех элементов в массив @elements
  # и в отдельные массивы - прямоугольники и текст.
  # @elements = ([ID,x,y,width,height], ...)
  my @arg = ($inkscape_x, '--query-all', "--without-gui", $file_svg);
  warn "Run `@arg`\n";
  open L, '-|', @arg or die "Can't run $inkscape_x: $!\n";
  my @elements;
  while (<L>) {
    chomp;
    my @ar = split ',';
    push @elements, [@ar];
    push @rect_elements, [@ar] if $ar[0] =~ /^rect/;
    push @text_elements, [@ar] if $ar[0] =~ /^text/;
  }
  close L;
  
  my $X = $elements[0][1];
  my $Y = $elements[0][2];
  my $R = $dpi/90;
  
  # Из массива прямоугольников выделяем в подмассивы расположенные на одинаковом уровне
  my %hhh;
  foreach (@rect_elements) {
    push @{$hhh{$_->[2]}}, $_;
  }
  #pp %hhh; exit;
  
  # сортируем по уровню (по Y)
  @rect_elements = sort {$a->[0][2] <=> $b->[0][2]} values %hhh;
  #pp @rect_elements; exit;
  
  # Выделяем тройки примыкающих друг к другу прямоугольников
  my @res; # в нем накапливаем подходящие тройки
  foreach (@rect_elements) {
    next if @$_< 3;
    @$_ = sort {$a->[1] <=> $b->[1]} @$_;
    my @t = @$_; # массив прямоугольников на одном уровне
    #pp @t;
    LOOP:
    while (1) {
      my $triples; # счетчик троек
      for (my $i=0; $i<@t-2; $i++) {
        for (my $j=$i+1; $j<@t-1; $j++) {
          #print "abs($t[$j][1]-$t[$i][1]-$t[$i][3])>1\n";
          next if abs($t[$j][1]-$t[$i][1]-$t[$i][3])>1;
          for (my $k=$j+1; $k<@t; $k++) {
            #print "abs($t[$k][1]-$t[$j][1]-$t[$j][3])<1\n";
            if (abs($t[$k][1]-$t[$j][1]-$t[$j][3])<1) {
              if (@res &&
                  $res[-1][0][2] == $t[$i][2] &&
                  $res[-1][2][1]+$res[-1][2][3] > $t[$i][1])
              {
                #pp @res;
                #print "$res[-1][2][1]+$res[-1][2][3] > $t[$i][1]\n";
                warn "Levels $res[-1][0][0]-$res[-1][1][0]-$res[-1][2][0] and $t[$i][0]-$t[$j][0]-$t[$k][0] are interlaced\n";
              }
              push @res, [ $t[$i],$t[$j],$t[$k] ];
              $triples++;
              splice @t, $k, 1;
              splice @t, $j, 1;
              splice @t, $i, 1;
              next LOOP;
            }
          }
        }
      }
      last LOOP unless $triples;
    }
  }
  #pp @res; exit;
  @rect_elements = @res;
  undef @res;

  open L, '<', $svg_file or die "Can't open $svg_file: $!\n";
  my $svg = join "", <L>;
  close L;

  my ($H) = $svg =~ /<svg\s(?:.(?!height="))*\sheight="([^"]+)"/s;
  
  my %htext;
  while ($svg =~ /<text[^>]*id="([^"]+)"[^>]*>(.*?)<\/text>/gs) {
    my ($id,$text) = ($1,$2);
    $text =~ s/<tspan[^>]*>//; $text =~ s/<\/tspan>//; # if text in tspan tag
    $text =~ s/^\s+//s; $text =~ s/\s+$//s; # spaces
    $htext{$id} = $text;
  }
  #pp %htext; exit;
  foreach (@text_elements) {
    if (! exists $htext{$_->[0]}) {
      warn "Can't find text elements $_->[0]\n";
      next;
    }
    push @$_, $htext{$_->[0]};
  }
  #pp @text_elements; exit;

  foreach my $r (@rect_elements) {
    my ($rx,$ry) = ($r->[0][1], $r->[0][2]+$r->[0][4]/2);
    my $dist=1e9;
    my @nearest;
    foreach my $t (@text_elements) {
      my ($tx,$ty) = ($t->[1], $t->[2]+$t->[4]/2);
      my $d = sqrt(($rx-$tx)**2 + ($ry-$ty)**2);
      if ($d<$dist) {
        @nearest = (@$t,$d);
        $dist = $d;
      }
    }
    
    if ($nearest[-1] > 20) {
      warn "Level $r->[0][0]$r->[1][0]$r->[2][0] ", 
            sprintf("(X,Y %d,%d: ",$r->[0][1],$H-$r->[0][2]),
            "nearest text is ",
            sprintf("$nearest[5] (%d px away)", $nearest[6]),
            "\n";
    }
    push @levels, [ ["$r->[0][0]$r->[1][0]$r->[2][0]",
                      sprintf("%d", ($r->[0][1]-$X)*$R),
                      sprintf("%d", ($r->[0][2]-$Y)*$R),
                      sprintf("%d", ($r->[0][3]+$r->[1][3]+$r->[2][3])*$R),
                      sprintf("%d", ($r->[0][4])*$R) ],
                    [$nearest[5],
                     sprintf("%d", ($nearest[1]-$X)*$R),
                     sprintf("%d", ($nearest[2]-$Y)*$R),
                     sprintf("%d", ($nearest[3])*$R),
                     sprintf("%d", ($nearest[4])*$R)] ];
  }
  #pp @levels; exit;
  
  if ($print_levels) {
    print STDERR "Levels:\n pngX  pngY |  svgX  svgY | Title\n"; 
    foreach (@levels) {
      printf STDERR "%5d %5d | %5d %5d | %s\n",
        $_->[0][1], $_->[0][2], $_->[0][1]/$R+$X, $H-($_->[0][2]/$R+$Y), $_->[1][0];
    }
    #print STDERR "@$_\n" for @levels;
    #exit;
  }

  return @levels;
}