Новосибирский институт органической химии им. Н.Н. Ворожцова СО РАН Лаборатория изучения механизмов органических реакций |
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||
xyz2levels#!/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; } |
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||