Новосибирский институт органической химии им. Н.Н. Ворожцова СО РАН
Лаборатория изучения механизмов органических реакций
 |
 |
 |
 |
 |
 |
xyz2jmod
#!/usr/bin/perl -ws use List::Util qw(max min sum reduce); #use Data::Dump qw(pp); our ($h, $help, $terminal, $range, $linewidth, $BB, $nuc, $label, $join, $fontsize, $jcamp); if ($h || $help) { (my $program = $0) =~ s/^.*[\/\\]//; print " Usage: $program [options] files.xyzppm Зависимости: Perl, gnuplot Рисует спектры ЯМР 13С по расчитанным хим. сдвигам с помощью gnuplot. Исходный файл д.б. в формате xyz с 5-й колонкой (в ней хим.сдвиги). Может быть и 6-я колонка, в ней обозначения атомов, если они есть, то будут использоваться вместо порядковых номеров атомом. Опции: -h help -terminal=x11 Терминал для gnuplot. По умолчания - x11. Если терминал - файловый формат (f.e. postscript), то перенаправить в файл. Список терминалов: echo 'help set terminal' | gnuplot -nuc=C Спектр на углероде, -nuc=H - на протонах (без расщеплений) -BB Обычный спектр (все сигналы вверх). По умолчанию - jmod (для С) -range=min,max Диапазон спектра (в м.д.). Если нет, то помещаются все сигналы -linewidth=0.05 Ширина линии (в м.д.). Можно давать в виде Hz/MHz, например 5/100. -label=0 Не надписывать пики. Это по умолчанию. -label Надписывать пики обозначениями атомов. -label=1 То же самое -label=2 Надписывать пики их хим. сдвигами. -label=3 Надписывать пики их хим. сдвигами и обозначениями атомов. -join=N Объединять надписи пиков, между которыми менее N ширин линий. По умолчанию N=1. -fontsize=10 Размер шрифта для надписей пиков (используется шрифт cour). -jcamp На каждый спектр делать JCAMP-DX (файл с расширением .dx). -terminal=null, если не нужен рисунок от gnuplot. "; exit; } ############################################################################# # Терминал для gnuplot если не задан, то x11 $terminal ||= 'x11'; # Кол-во точек в спектре будет определяться автоматически # $SI = 1000; # можно задать фиксированное значение $nuc ||= 'C'; $BB = 1 if $nuc ne 'C'; unless ($linewidth) { $linewidth = 5/100; # Ширина линии (Гц/МГц) $linewidth = 1/200 if $nuc eq 'H'; } if ($linewidth =~ /^(\d+(?:\.\d+)?)(?:\/(\d+(?:\.\d+)?))?$/) { die "Invalid linewidth\n" if $2 && $2==0; $linewidth = $1/$2 if $2; } else { die "Invalid linewidth\n"; } # Для линий, отстоящих менее чем на $join полуширин, объединять метки $join ||= 1; ## Умолчаемое значение для $label #$label = 1 unless defined $label; #$label = 0 if $nolabel; # размер шрифта для надписей пиков (используется шрифт cour) $fontsize ||= 10; # Относительная интенсивность синглета, дублета, триплета и квартета в JMOD @intensivity = (-0.3, 0.7, -0.7, 0.7); #@intensivity = (-0.2, 0.5, -0.5, 0.5); # @intensivity = (0.7, 0.7, 0.7, 0.7); # все сигналы одинаковой интенсивности # gnuplot (программа для рисования графиков) $gnuplot = `which gnuplot`; chomp $gnuplot; ############################################################################## $minppm = 1e5; $maxppm = -1e5; my $label_ppm_format = '%.1f'; $label_ppm_format = '%.2f' if $nuc eq 'H'; # Берем из опции диапазон, в котором рисовать спектр if (defined $range) { ($minppm, $maxppm) = split /[^\d.-]+/, $range; m/^-?\d+(\.\d+)?$/ || die "Invalid range\n" for ($minppm, $maxppm); die "Invalid range\n" if $minppm == $maxppm; ($minppm, $maxppm) = ($maxppm, $minppm) if $minppm > $maxppm; } # Зачитываем xyz-файлы с хим. сдвигами, каждый в свой датасет $dataset = 0; while (&read_molden) { $title[$dataset] = $ARGV; # Название спектра - название xyzppm-файла #$title[$dataset] =~ s/\.xyz(ppm)?$//; # Атомы углерода, у которых есть хим.сдвиги и их мультиплетность for ($i=1; $i<=$N; $i++) { if (uc($atom[$i]) eq $nuc && $ppm[$i]=~/^-?\d+\.\d+/) { my $mult = 1; if ($nuc eq 'C') { $mult = 0; for ($j=1; $j<=$N; $j++) { if ( uc($atom[$j]) eq 'H' && ($x[$i]-$x[$j])**2+($y[$i]-$y[$j])**2+($z[$i]-$z[$j])**2 < 1.44 ) { $mult++; } } } @{$PPM[$dataset]{$atom[$i].$i}} = (0+$ppm[$i], $mult, $user_label[$i]); unless (defined $range) { $minppm = $ppm[$i] if $ppm[$i] < $minppm; $maxppm = $ppm[$i] if $ppm[$i] > $maxppm; } } } $dataset++; } exit unless @PPM; # Если не задан диапазон хим.сдвигов unless (defined $range) { my $r = $maxppm-$minppm; $r = $linewidth*50 if $r == 0; $minppm = $minppm-$r/10; $maxppm = $maxppm+$r/10; } # Если не задано количество точек в спектре unless ($SI) { my $p = int(($maxppm-$minppm)/$linewidth/2); $SI = (1,2,4,4,5,8,8,8,10,10)[substr $p, 0, 1] . '0'x(length($p)-1); $SI = 500 if $SI < 500; $SI = 10000 if $SI > 10000; } # ppm.out - графики спектров в табличной форме для gnuplot unlink 'ppm.out' if -f 'ppm.out'; open OUT, ">>", "ppm.out" or die "Can't write to ppm.out: $!\n"; #pp @PPM; my @spc; foreach $dataset (0..$#PPM) { next unless $PPM[$dataset]; my ($spc,$peak) = specrtum( $PPM[$dataset], maxppm=>$maxppm, minppm=>$minppm, linewidth=>$linewidth, SI=>$SI, baseoffset=>$dataset, BB=>$BB, intensivity=>\@intensivity, ); # $spc - ссылка на массив точек (x,y) спектра # $spc->[x_точки, y_точки] # $peak - ссылка на хэш вершин пиков # $peak->{atom_label}[номер_точки, x_вершины, y_вершины, user_label] #pp $peak; foreach (@{$spc}) { print OUT ${$_}[0], " \t", ${$_}[1], "\n"; } print OUT "\n"; #push @peak, $peak if $label != 0; push @peak, $peak; push @spc, $spc; } close OUT; unless (-s 'ppm.out') { unlink 'ppm.out'; exit; } if ($jcamp) { my @files = unique_array(@title); for (my $i=0; $i<@spc; $i++) { write_jcamp($spc[$i],$files[$i]); # my $file = $files[$i]; # print "$file\n"; } } # Если есть gnuplot, то будем рисовать спектры, а если нет, # то рисовать не будем, но и удалять ppm.out тоже не будем if ($gnuplot && $terminal ne 'null') { open GP, "|-", "$gnuplot -persist" or die "Can't pipe to gnuplot: $!\n"; print GP " set xrange [$maxppm:$minppm] reverse set xlabel 'ppm' set xtics mirror set tics out set noytics set border 1 set terminal $terminal "; my $maxpeak = ['','',-1000]; my $minpeak = ['','',1000]; my $str_max = ''; my $str_min = ''; foreach $dataset (0..$#PPM) { $ytitle = $dataset+0.1; (my $title = $title[$dataset]) =~ s/\.xyz(ppm)?$//; print GP "set label '$title, $nuc' at $maxppm,$ytitle rotate\n"; # if ($label) { my @labels; # Массив меток [text,x,y] # Формируем массив меток foreach (keys %{$peak[$dataset]}) { #my $xlabel = $peak[$dataset]{$_}[1]-$linewidth; my $xlabel = $peak[$dataset]{$_}[1]; my $ylabel = $peak[$dataset]{$_}[2]; push @labels, [($peak[$dataset]{$_}[3] || $_),$xlabel,$ylabel]; } next unless @labels; # Объединяем метки по хим.сдвигам (у полученного массива другой формат) my @l = combine(\@labels, 1, $join*$linewidth); my @join_labels; foreach (@l) { my @a; foreach my $i (0..$#{$_->[0]}) { push @a, [ $_->[0][$i],$_->[1][$i],$_->[2][$i] ]; } push @join_labels, combine(\@a, 2, 0.2); } # Формируем новые названия и расположения меток #pp @join_labels; foreach (@join_labels) { #pp $_; #next unless defined $_->[0][0]; $_->[0] = join ',', sort {($A=$a)=~s/(\d+)/sprintf '%05d',$1/eg; ($B=$b)=~s/(\d+)/sprintf '%05d',$1/eg; $A cmp $B} @{$_->[0]}; $_->[1] = sum(@{$_->[1]})/@{$_->[1]}; $_->[2] = max(@{$_->[2]}); } # C1,C2,C3 => C1,2,3 foreach (@join_labels) { 1 while $_->[0]=~s/((?:^|,)$nuc(?:\d+,)+)$nuc(\d+(?:,|$))/$1$2/g; } my $max_peak = reduce {$a->[2]>$b->[2] ? $a : $b} @join_labels; $maxpeak = $max_peak if $maxpeak->[2] < $max_peak->[2]; my $min_peak = reduce {$a->[2]<$b->[2] ? $a : $b} @join_labels; $minpeak = $min_peak if $minpeak->[2] > $min_peak->[2]; #pp $max_peak; if ($label) { if ($label == 3) { $str_max = sprintf "%.2f (%s)", $max_peak->[1],$max_peak->[0]; $str_min = sprintf "%.2f (%s)", $min_peak->[1],$min_peak->[0]; } elsif ($label == 1) { $str_max = sprintf "%s", $max_peak->[0]; $str_min = sprintf "%s", $min_peak->[0]; } elsif ($label == 2) { $str_max = sprintf "%.2f", $max_peak->[1]; $str_min = sprintf "%.2f", $min_peak->[1]; } # Печатаем метки в gnuplot foreach (@join_labels) { my $direct = $_->[2]>$dataset ? 'left' : 'right'; if ($label == 3) { printf GP "set label '$label_ppm_format ($_->[0])' at $_->[1],$_->[2] $direct rotate font 'cour,$fontsize'\n",$_->[1]; } elsif ($label == 1) { print GP "set label '$_->[0]' at $_->[1],$_->[2] $direct rotate font 'cour,$fontsize'\n"; } elsif ($label == 2) { printf GP "set label '$label_ppm_format' at $_->[1],$_->[2] $direct rotate font 'cour,$fontsize'\n",$_->[1]; } } } } # } my $symbolwigth = $fontsize/10*1.37; # mm my $scale = 100; # Высота y-шкалы в mm my $max_y = $#PPM; my $min_y = 0; $max_y = $intensivity[1] if $max_y < $intensivity[1] && $BB; $min_y = $intensivity[0] if $min_y > $intensivity[0] && !$BB; $max_y = $maxpeak->[2] if $maxpeak->[2] > $max_y; $min_y = $minpeak->[2] if $minpeak->[2] < $min_y; #print "$min_y:$max_y\n"; if ($label) { #my $range = ($max_y - $min_y); my ($length_max,$length_min) = (0,0); $length_max += length($str_max) if $max_y > $#PPM; $length_min += length($str_min) if $min_y < 0; #print "length $length_max + $length_min, ", $symbolwigth*($length_max + $length_min)," mm\n"; my $k = ($scale - $symbolwigth*($length_max + $length_min))/($max_y - $min_y); $k = 10 if $k < 10; #print "k = $k\n"; $max_y += $symbolwigth*$length_max/$k if $max_y > $#PPM; $min_y -= $symbolwigth*$length_min/$k if $min_y < 0; } #print "$min_y:$max_y\n"; my $offset = ($max_y - $min_y)/20; print GP "set yrange [$min_y:$max_y]\n"; # print GP "set ytics\n"; print GP "set offsets 0, 0, $offset, $offset\n"; print GP "plot 'ppm.out' notitle with line\n"; close GP; # unlink 'ppm.out' if -f 'ppm.out'; } sub read_molden { return undef if eof(); # 1-st string if (<>=~/^\s*(\d+)\s*$/) { ($N,$energy,@atom,@x,@y,@z,@ppm,@user_label) = undef; $N = $1; } else { return undef; } # 2-nd string ($energy) = <>=~/(-?\d+\.\d+)/; # Geometry for ($i=1;$i<=$N;$i++) { ($atom[$i],$x[$i],$y[$i],$z[$i],$ppm[$i],$user_label[$i]) = split ' ', <>; return undef unless $atom[$i]=~/^\w\w?/ && $x[$i]=~/^-?\d+\.\d+/ && $y[$i]=~/^-?\d+\.\d+/ && $z[$i]=~/^-?\d+\.\d+/; $user_label[$i] =~ s/(['";,])/\\$1/g if $user_label[$i]; } return $N; } sub specrtum { my $ppm = shift; my %param = ( linewidth => 0.02, maxppm => 250, minppm => 0, SI => '8000', height => 0.7, baseoffset => 0, level => 0.001, BB => 0, intensivity => [-0.3, 0.7, -0.7, 0.7], @_ ); my $debug = 0; if ($debug) { print "%param:"; print "\t$_ => $param{$_}\n" foreach (sort keys %param); } my $SW = $param{maxppm} - $param{minppm}; $SW = $param{linewidth} if $SW==0; my $d = $SW / $param{SI}; print "SW = $SW ppm, ppm per point = $d\n" if $debug; my @SPC; foreach (0..$param{SI}) { $SPC[$_][0] = $param{minppm} + $_*$d; $SPC[$_][1] = $param{baseoffset}; } my @halfline; my $b = $param{linewidth}/2; my $y = 1; my $dd = $d; while ($y > $param{level}) { push @halfline, $y; my $t = $dd/$b; $y = 1 / (1 + $t*$t); $dd += $d; } print "halfline: @halfline\n" if $debug; my %peak; foreach (keys %$ppm) { next if $ppm->{$_}[0] < $param{minppm} || $ppm->{$_}[0] > $param{maxppm}; my $x = $ppm->{$_}[0]; my $mult = $param{intensivity}->[$ppm->{$_}[1]]; $mult = abs($mult) if $param{BB}; # Помещаем вершину линии на ближайшую точку my $I = sprintf "%d", ($x - $param{minppm}) / $d; foreach (-$#halfline..$#halfline) { my $j = $I + $_; next if $j<0 or $j>$SI; $SPC[$j][1] += $halfline[abs($_)]*$mult; } $peak{$_}[0] = $I if $I>=0 && $I<=$SI; $peak{$_}[3] = $ppm->{$_}[2]; } foreach (keys %peak) { $peak{$_}[1] = $SPC[$peak{$_}[0]][0]; $peak{$_}[2] = $SPC[$peak{$_}[0]][1]; print "$_: point $peak{$_}[0], x=$peak{$_}[1], y=$peak{$_}[2]\n" if $debug; }; return (\@SPC, \%peak); } sub combine { my @a = @{$_[0]}; my $field = $_[1]; my $criter = $_[2]; @a = sort {$a->[$field]<=>$b->[$field]} @a; my @b = [ [$a[0][0]],[$a[0][1]],[$a[0][2]] ]; for (my $i=1; $i<@a; $i++) { if ($a[$i][$field]-$a[$i-1][$field] <= $criter) { push @{$b[-1][0]}, $a[$i][0]; push @{$b[-1][1]}, $a[$i][1]; push @{$b[-1][2]}, $a[$i][2]; } else { push @b, [ [$a[$i][0]],[$a[$i][1]],[$a[$i][2]] ]; } } return @b; } sub write_jcamp { my ($spc, $file) = @_; my ($sec,$min,$hour,$mday,$mon,$year) = localtime(-f $file ? (stat $file)[9] : time()); my $DATE = sprintf "%02d/%02d/%02d", $year % 100, $mon+1, $mday; my $TIME = sprintf "%02d:%02d:%02d", $hour,$min,$sec; open F, '>', "$file.dx" or die "Can't write $file.dx: $!\n"; print F "##TITLE= Calculated $nuc NMR from $file ##JCAMP-DX= 4.24 ##DATA TYPE= NMR SPECTRUM ##ORIGIN= $0 ##OWNER= $ENV{USER} ##DATE= $DATE ##TIME= $TIME ##DELTAX= @{[$spc->[1][0]-$spc->[0][0]]} ##XUNITS= PPM ##YUNITS= ABSORBANCE ##XFACTOR= 1.00 ##YFACTOR= 1.00 ##FIRSTX= $spc->[0][0] ##LASTX= $spc->[-1][0] ##NPOINTS= @{[$#{$spc}+1]} ##FIRSTY= @{[$spc->[0][1]*1e8]} ##XYDATA= (X++(Y..Y)) "; ##.OBSERVE FREQUENCY= 75.476050452 ##.OBSERVE NUCLEUS= ^13C ##SPECTROMETER/DATA SYSTEM= Prirda, Dimitri Laikov foreach (@$spc) { printf F " %11.6f %12d\n", $_->[0], $_->[1]*1e8; } print F "##END= \n"; close F; } sub unique_array { my @a = @_; my %h; $h{$_}++ foreach @a; for(my $i = $#a; $i>=0; $i--) { next if $h{$a[$i]}==1; my $j = --$h{$a[$i]}; $j++ while exists $h{"$a[$i]-$j"}; $a[$i] .= "-$j"; $h{$a[$i]}++; } return @a; }