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

mskxyz

#!/usr/bin/perl -ws

use vars '$h';
if ($h) {
  (my $program = $0) =~ s/.*[\/\\]//;
  print "
Usage: $program file.xyzppm | msketch -

Зависимости: Perl, msketch из пакета Marvin Beans 

Показывает структуру(ы) с хим.сдвигами в MarvinSketch
Строго говаря, msketch для работы скрипта не нужен, 
но скрипт без него не имеет смысла.
  \n";
  exit;
}

use List::Util qw(sum max min);

#$molconvert = 'molconvert';

$N0 = 0;
$endX = 0;

while (&read_molden) {
  push @elementType, @atom[1..$N];
  
  foreach (1..$N) {
    push(@mrvExtraLabel, '0'), next if $ppm[$_] eq '';
    if ($atom[$_] eq 'H') {
      push @mrvExtraLabel, sprintf("%.2f", $ppm[$_]);
    } else {
      push @mrvExtraLabel, sprintf("%.1f", $ppm[$_]);
    }
  }
  
  # Перемещаем геометрию
  $minX = min(@x[1..$N]);
  $minY = min(@y[1..$N]);
  $centrZ = sum(@z[1..$N]) / $N;
  $_ -= $minX for @x[1..$N];
  $_ -= $minY for @y[1..$N];
  $_ -= $centrZ for @z[1..$N];
  
  # Смещаем по x
  push @x3, map {$_+$endX} @x[1..$N];
  push @y3, @y[1..$N];
  push @z3, @z[1..$N];
  #print "@x3\n\n";
  $endX += max(@x[1..$N])+5;
  
  for ($i=1; $i<=$N; $i++) {
    for ($j=$i+1; $j<=$N; $j++) {
      if (($x[$i]-$x[$j])**2+($y[$i]-$y[$j])**2+($z[$i]-$z[$j])**2 < 2.9) {
        push @bondArray, 'a'.($i+$N0).' a'.($j+$N0);
      }
    }
  }

#  open XYZ, "| $molconvert mrv - > TMP.mrv" or die "Can't run $molconvert: $!\n";
#  write_molden('XYZ');
#  close XYZ;
#  
#  open MRV, 'TMP.mrv' or die "Can't open temp mrv: $!\n";
#  while (<MRV>) {
#    while (m/<bond atomRefs2="a(\d+)\s+a(\d+)/g) {
#      push @{$bondArray{$1+$N0}}, $2+$N0;
#    }
#  }
#  close MRV;
  
  
  $N0 += $N;
}
@atomID = map {"a$_"} 1..$N0;

print qq(<?xml version="1.0" ?>
<MDocument>
  <MChemicalStruct>
    <molecule molID="m1">
      <atomArray
          atomID="@atomID"
          elementType="@elementType"
          mrvExtraLabel="@mrvExtraLabel"
          x3="@x3"
          y3="@y3"
          z3="@z3"
          />
      <bondArray>\n);
foreach (@bondArray) {
  print qq(        <bond atomRefs2="$_" order="1" />\n);
}
print qq(      </bondArray>
    </molecule>   
  </MChemicalStruct>
</MDocument>
);

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 {
  my $fh = shift || 'STDOUT';
  if (@atom && @x && @y && @z) {
    print $fh " $#atom\n";
    print $fh " Energy $energy" if $energy;
    print $fh "\n";
    for ($i=1; $i<=$#atom; $i++) {
      if ($atom[$i] && $x[$i] && $y[$i] && $z[$i]) {
        printf $fh " %-2s %12.8f %12.8f %12.8f",$atom[$i],$x[$i],$y[$i],$z[$i];
        #printf $fh $atom[$i] eq 'H' ? " %10.3f" : " %9.2f"  , $ppm[$i] if $ppm[$i];
        print $fh "\n";
      }
    }
  }
}