Новосибирский институт органической химии им. Н.Н. Ворожцова СО РАН Лаборатория изучения механизмов органических реакций |
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||
xyz2ppm#!/usr/bin/perl -ws #use Data::Dump qw(pp); use vars qw($h); (my $program = $0) =~ s/^.*[\/\\]//; if ($h) { print <<HELP; Usage: $program files.xyzppm cat files.xyzppm | $program Зависимости: perl, symmetry (не обязательно) Печать таблицы хим. сдвигов, подходящей для импорта в для Excel, из *.xyzppm Автоматически будут усреднены и объединены протоны метилов. HELP exit; } #$symmetry = '/usr/local/bin/symmetry'; #unless (-x $symmetry) { # warn "Program symmetry not found\n"; # undef $symmetry; #} # Вгоняем таблицу Менделеева @ATOM = qw (XX H He Li Be B C N O F Ne Na Mg Al Si P S Cl Ar K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe Cs Ba La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr ); for ($i=0; $i<@ATOM; $i++) { $ATOM{$ATOM[$i]} = $i }; my $prev_head = ''; while (&read_molden) { my @new_atom = map {$atom[$_].$_} 1..$#atom; unshift @new_atom, undef; #pp @new_atom; # Усредняем протоны метильных групп for ($i=1; $i<=$N; $i++) { if ( $atom[$i] eq 'C' ) { $CH[$i] = []; for ($j=1; $j<=$N; $j++) { if ( $atom[$j] eq 'H' && ($x[$i]-$x[$j])**2+($y[$i]-$y[$j])**2+($z[$i]-$z[$j])**2 < 1.44 ) { push @{$CH[$i]}, $j; } } } } for ($i=1; $i<=$N; $i++) { if ($atom[$i] eq 'C' && @{$CH[$i]}==3) { my $aver = ($ppm[$CH[$i][0]]+$ppm[$CH[$i][1]]+$ppm[$CH[$i][2]])/3; $ppm[$CH[$i][0]] = $ppm[$CH[$i][1]] = $ppm[$CH[$i][2]] = $aver; $new_atom[$CH[$i][0]] .= ",$new_atom[$CH[$i][1]],$new_atom[$CH[$i][2]]"; $new_atom[$CH[$i][1]] = $new_atom[$CH[$i][2]] = ''; } } #pp @CH; my $head = "file\tEnergy"; my $str = "$ARGV\t$energy"; for ($i=1; $i<@new_atom; $i++) { next unless $new_atom[$i]; next unless $ppm[$i]; $head .= "\t$new_atom[$i]"; $str .= sprintf "\t%.3f", $ppm[$i]; } #pp @new_atom; print "$head\n" if $head ne $prev_head; print "$str\n"; $prev_head = $head; # if ($symmetry) { # open TINP, ">symm.TINP" || warn "Can't write temporary file: $!\n" , next; # print TINP " $N\n"; # for ($i=1; $i<=$N; $i++) { # print TINP " $ATOM{$atom[$i]} $x[$i] $y[$i] $z[$i]\n"; # } # close TINP; # # $generator = undef; # open OUT, "$symmetry -final 0.0009999 symm.TINP |"; # while (<OUT>) { # if (/^Molecule has the following abelian symmetry generators: (\S+)/) { # $generator = $1; # last; # } # } # while (<OUT>) { # if (m/^Abelian symmetrized coordinates of all atoms$/) { # <OUT>; # $N = <OUT>; chomp $N; # for ($i=1; $i<=$N; $i++) { # ($natom,$x[$i],$y[$i],$z[$i]) = split ' ', <OUT>; # $atom[$i] = $ATOM[$natom]; # } # } # } # close OUT; # unlink 'symm.TINP'; # # @symm_atoms = symm_atoms($generator); # pp @symm_atoms; # SYMM_ATOMS: # foreach (@symm_atoms) { # #print "@{$_}\n"; # my $sum = 0; # foreach (@{$_}) { # next SYMM_ATOMS if $ppm[$_] eq ''; # $sum += $ppm[$_]; # } # my $aver = $sum/@{$_}; # $new_atom[$_->[0]] = join ',', map {$new_atom[$_]} @$_; # $ppm[shift @{$_}] = $aver; # foreach (@{$_}) { # $ppm[$_] = ''; # $new_atom[$_] = ''; # } # # } # } #&write_molden; } #sub symm_atoms { # my $generator = shift; # my @symm_atoms; # unless ($generator) { # @symm_atoms = map {[$_]} 1..$N; # return @symm_atoms; # } # my %condition = ( # 'z' => '$x[$i]==$x[$j] && $y[$i]==$y[$j] && $z[$i]==-$z[$j]', # 'x,y' => 'abs($x[$i])==abs($x[$j]) && # abs($y[$i])==abs($y[$j]) && # $z[$i]==$z[$j]', # 'xy' => '$x[$i]==-$x[$j] && $y[$i]==-$y[$j] && $z[$i]==$z[$j]', # 'x,y,z' => 'abs($x[$i])==abs($x[$j]) && # abs($y[$i])==abs($y[$j]) && # abs($z[$i])==abs($z[$j])', # 'xyz' => '$x[$i]==-$x[$j] && $y[$i]==-$y[$j] && $z[$i]==-$z[$j]', # 'xy,z' => '($x[$i]==-$x[$j] && $y[$i]==-$y[$j] && $z[$i]==$z[$j]) || # ($x[$i]==$x[$j] && $y[$i]==$y[$j] && $z[$i]==-$z[$j])', # 'xy,xz,yz' => '($x[$i]==-$x[$j] && $y[$i]==-$y[$j] && $z[$i]==$z[$j]) || # ($x[$i]==-$x[$j] && $y[$i]==$y[$j] && $z[$i]==-$z[$j]) || # ($x[$i]==$x[$j] && $y[$i]==-$y[$j] && $z[$i]==-$z[$j])' # ); # my @at = @atom; # for ($i=1; $i<=$N; $i++) { # next unless $at[$i]; # push @symm_atoms, [($i)]; # #$at[$i] = undef; # for ($j=$i+1; $j<=$N; $j++) { # next unless $at[$j]; # if (exists $condition{$generator} && eval "$condition{$generator}") { # $at[$j] = undef; # push @{$symm_atoms[-1]}, $j; # } # } # } # return @symm_atoms; #} sub read_molden { return undef if eof(); # 1-st string if (<>=~/^\s*(\d+)\s*$/) { ($N,$energy,@atom,@x,@y,@z,@ppm) = 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],undef) = split ' ', <>; return undef unless $atom[$i]=~/^\w\w?/ && $x[$i]=~/^-?\d+\.\d+/ && $y[$i]=~/^-?\d+\.\d+/ && $z[$i]=~/^-?\d+\.\d+/; } return $N; } sub write_molden { if (@atom && @x && @y && @z) { print " $#atom\n"; print " Energy $energy" if $energy; print "\n"; for ($i=1; $i<=$#atom; $i++) { if ($atom[$i] && $x[$i] && $y[$i] && $z[$i]) { printf " %-2s %12.8f %12.8f %12.8f",$atom[$i],$x[$i],$y[$i],$z[$i]; printf $atom[$i] eq 'H' ? " %10.3f" : " %9.2f" , $ppm[$i] if $ppm[$i]; print "\n"; } } } } |
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||