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

addXX

#!/usr/bin/perl -ws

use vars qw($h $Z);
if (!@ARGV or $h) {
  (my $program = $0) =~ s/.*[\/\\]//;
  print "
Usage: $program [-Z] 'i,j,k... l,m... ...' file.xyz ...

Зависимости: perl

Добавляет фиктивные атомы (XX) в геометрический центр каждого из наборов
атомов, перечисленных в первом параметре. i,j,k,l,m - номера атомов.
Если после набора атомов в круглых скобках даны расстояния (все без пробелов),
то XX добавляются на этих расстояниях над центром кольца.
Если дана опция -Z и хотя бы тройка атомов, выполняется переориентация 
молекулы: центр координат помещается в центр перечисленных атомов, 
ось Z направлена перпендикулярно плоскости, задаваемой этими атомами,
a ось X - в сторону первого из атомов тройки.
Нужно все это, в частности, для определения точек, в которых считать NICS.
file.xyz - файл в молденовском xyz-формате. 
Примеры:
$program 1-6 file.xyz - добавить XX в центр атомов с 1 по 6-й 
  (например, в центр бензольного кольца)
cat file.xyz | $program '1,2 3-5' - добавить два XX, между 1-м и 2-м и
  между 3-5 атомами
$program '1-6(0,-1,1)' file.xyz - добавить XX в центр кольца, заданного атомами
  1-6, а также на расстоянии 1 ангстрем под и над кольцом.
$program -Z '1-6()' file.xyz - просто переориентировать молекулу, не добавляя XX
  \n";
  exit;
}


my ($N,$energy,@atom,@x,@y,@z,@charge);

@at_num = split ' ', shift @ARGV;

while (&read_molden()) {
  $OK = 1;
  $first_at_num = 1;
  AT_NUM:
  foreach (@at_num) {
    if ( m/([^(]+)\(([^)]*)\)/ ) {
      @nums = get_atom_numbers($1);
      if (@nums<3) {
        warn "Not enough numbers: @nums\n";
        undef $OK;
        next AT_NUM;
      }
      
      @distances = split /,/, $2;
      foreach (@distances) {
        warn "Invalid: @distances\n" and next AT_NUM unless m/^-?\d+(\.\d+)?$/;
      }
      #warn "@nums (@distances)\n";
      
      if ($Z && $first_at_num) {
        XXorientation(@nums);
        $first_at_num = undef;
      }
      
      ($x,$y,$z) = centre_coordinates(@nums);
      undef $OK and next AT_NUM unless defined $x;
      
      ($A,$B,$C) = plane_normal(@nums);
      undef $OK and next AT_NUM unless defined $A;
      
      foreach (@distances) {
        add_XX($x + $A*$_, $y + $B*$_, $z + $C*$_);
      }
      
    } else {
      @nums = get_atom_numbers($_);
      next AT_NUM unless @nums;
      (warn("Only one number: @nums\n"), undef $OK, next AT_NUM) if @nums==1;

      if ($Z && $first_at_num) {
        XXorientation(@nums);
        $first_at_num = undef;
      }
      
      ($x,$y,$z) = centre_coordinates(@nums);
      undef $OK and next AT_NUM unless defined $x;

      add_XX($x, $y, $z);
    }
  }
  &write_molden if $OK;
}

sub XXorientation {
  my @nums = @_;
  warn "Not enough numbers (@nums) for reorientation\n" and return undef if @nums<3;
  
  my ($x,$y,$z) = centre_coordinates(@nums);
  warn "Can't centre_coordinates(@nums)\n" and return undef unless defined $x;

  my ($A,$B,$C) = plane_normal(@nums);
  warn "Can't plane_normal(@nums)\n" and return undef unless defined $A;
  
  add_XX($x, $y, $z);
  add_XX($x+$A, $y+$B, $z+$C);
  
  re_orientation($#atom-1,$#atom,$nums[0]);
  
  $#atom -= 2; $#x -= 2; $#y -= 2; $#z -= 2;
  my @ttt = @x; @x = @y; @y = @z; @z = @ttt;
}

sub centre_coordinates {
  my $x = 0; my $y = 0; my $z = 0;
  my @nums = @_;
  return undef unless @nums; 
  foreach (@nums) {
    warn "Invalid atom number: $_" and return undef if $_>$#atom;
    $x += $x[$_]; $y += $y[$_]; $z += $z[$_]; 
  }
  $x /= $#nums+1; $y /= $#nums+1; $z /= $#nums+1;
  return ($x,$y,$z);
}

sub plane_normal {
  my @nums = @_;
  my ($A,$B,$C);
  return undef if @nums<3;
  
  my ($x,$y,$z) = centre_coordinates(@nums);
  my (@x1,@y1,@z1);
  foreach (@nums) {
    $x1[$_] = $x[$_]-$x; $y1[$_] = $y[$_]-$y; $z1[$_] = $z[$_]-$z;
  }
  my ($XX,$XY,$XZ,$YY,$YZ,$ZZ);
  foreach (@nums) {
    #$X  += $x[$_];        $Y  += $y[$_];        $Z  += $z[$_];
    $XX += $x1[$_]*$x1[$_]; $YY += $y1[$_]*$y1[$_]; $ZZ += $z1[$_]*$z1[$_];
    $XY += $x1[$_]*$y1[$_]; $YZ += $y1[$_]*$z1[$_]; $XZ += $z1[$_]*$x1[$_];
  }
  #print STDERR "$XX  $XY  $XZ\n$YY  $YZ\n$ZZ\n";
  my @eigen = jacobi($XX,$XY,$XZ,$YY,$YZ,$ZZ) or die "Eigen error\n";
  ($A,$B,$C) = @{$eigen[2][1]};
#  my $d = 1e-3;
#  if (@nums == 3) {
#    my ($i,$j,$k) = @nums;
#    $A = ($y[$j]-$y[$i])*($z[$k]-$z[$i])-($z[$j]-$z[$i])*($y[$k]-$y[$i]);
#    $B = ($z[$j]-$z[$i])*($x[$k]-$x[$i])-($x[$j]-$x[$i])*($z[$k]-$z[$i]);
#    $C = ($x[$j]-$x[$i])*($y[$k]-$y[$i])-($y[$j]-$y[$i])*($x[$k]-$x[$i]);
#    #print "@nums  $A $B $C\n";
#  } else {
#    my ($X,$Y,$Z,$XX,$YY,$ZZ,$XY,$YZ,$XZ);
#    foreach (@nums) {
#      $X  += $x[$_];        $Y  += $y[$_];        $Z  += $z[$_];
#      $XX += $x[$_]*$x[$_]; $YY += $y[$_]*$y[$_]; $ZZ += $z[$_]*$z[$_];
#      $XY += $x[$_]*$y[$_]; $YZ += $y[$_]*$z[$_]; $XZ += $z[$_]*$x[$_];
#    }
#    print STDERR "$X,$Y,$Z,$XX,$YY,$ZZ,$XY,$YZ,$XZ\n";
#    
#    if (abs($XX)<$d && abs($YY)<$d or 
#        abs($XX)<$d && abs($ZZ)<$d or 
#        abs($YY)<$d && abs($ZZ)<$d) {
#      return undef;
#    }
#    
#    if    (abs($ZZ)<$d) {
#      $A = 0;
#      $B = 0;
#      $C = 1;
#    }
#    elsif (abs($YY)<$d) {
#      $A = 0;
#      $B = 1;
#      $C = 0;
#    }
#    elsif (abs($XX)<$d) {
#      $A = 1;
#      $B = 0;
#      $C = 0;
#    }
#    else {
#      # A*XX + B*XY + C*XZ = -X   XX XY XZ | -X    -X XY XZ   XX -X XZ   XX XY -X
#      # A*XY + B*YY + C*YZ = -Y   XY YY YZ | -Y    -Y YY YZ   XY -Y YZ   XY YY -Y
#      # A*XZ + B*YZ + C*ZZ = -Z   XZ YZ ZZ | -Z    -Z YZ ZZ   XZ -Z ZZ   XZ YZ -Z
#      $A = ($X*($YY*$ZZ-$YZ**2)+$XY*($YZ*$Z-$Y*$ZZ)+$XZ*($Y*$YZ-$YY*$Z));
#      $B = -($XX*($YZ*$Z-$Y*$ZZ)+$X*$XY*$ZZ+$XZ*(-$XY*$Z-$X*$YZ)+$XZ**2*$Y);
#      $C = ($XX*($YY*$Z-$Y*$YZ)-$XY**2*$Z+$X*$XY*$YZ+$XZ*($XY*$Y-$X*$YY));
##      $A = $X*($YZ*$YZ-$YY*$ZZ) + $Y*($XY*$ZZ-$YZ*$XZ) + $Z*($YY*$XZ-$XY*$YZ);
##      $B = $X*($XY*$ZZ-$XZ*$YZ) + $Y*($XZ*$XZ-$XX*$ZZ) + $Z*($XX*$YZ-$XY*$XZ);
##      $C = $X*($YY*$XZ-$YZ*$XY) + $Y*($YZ*$XX-$XY*$XZ) + $Z*($XY*$XY-$YY*$XX);
#    }
#  }
  my $normal = sqrt($A**2 + $B**2 +$C**2);
  #print STDERR "@nums  $A $B $C\n";
  $A /= $normal; $B /= $normal; $C /= $normal;
  return ($A,$B,$C);
}

#sub plane_norm {
#  my @nums = @_;
#  return undef if @nums<3;
#  my $A = 0; my $B = 0; my $C = 0;
#  for ($num=0; $num<=$#nums-2; $num++) {
#    $i = $nums[$num]; $j = $nums[$num+1]; $k = $nums[$num+2]; 
#    $A += ($y[$j]-$y[$i])*($z[$k]-$z[$i])-($z[$j]-$z[$i])*($y[$k]-$y[$i]);
#    $B += ($z[$j]-$z[$i])*($x[$k]-$x[$i])-($x[$j]-$x[$i])*($z[$k]-$z[$i]);
#    $C += ($x[$j]-$x[$i])*($y[$k]-$y[$i])-($y[$j]-$y[$i])*($x[$k]-$x[$i]);
#  }
#  $A /= $#nums-1; $B /= $#nums-1; $C /= $#nums-1;
#  my $normal = sqrt($A**2 + $B**2 +$C**2);
#  $A /= $normal; $B /= $normal; $C /= $normal;
#  #print "@nums  $A $B $C\n";
#  return ($A,$B,$C);
#}

sub add_XX {
  push @x, $_[0];
  push @y, $_[1];
  push @z, $_[2];
  push @atom, 'XX';
}

sub centre_to_atom {
  my $i = shift;
  return undef if $i > @atom;
  my ($xi,$yi,$zi) = ($x[$i],$y[$i],$z[$i]);
  for ($n=1; $n<=$#atom; $n++) {
    $x[$n] -= $xi; $y[$n] -= $yi; $z[$n] -= $zi; 
  }
}

sub re_orientation {
  my ($i,$j,$k) = @_;
  centre_to_atom($i);
  #&write_molden;
  my $xj = sqrt($x[$j]**2 + $y[$j]**2 + $z[$j]**2);
  my $xk = ($x[$k]*$x[$j] + $y[$k]*$y[$j] + $z[$k]*$z[$j]) / $xj;
  my $yk = sqrt($x[$k]**2 + $y[$k]**2 + $z[$k]**2 - $xk**2);
  my ($A1,$B1,$C1) = plane_normal($i,$j,$k);
  
  for ($n=1; $n<@atom; $n++) {
    if ($n != $i && $n != $j && $n != $k) {
      my $xn = ($x[$n]*$x[$j] + $y[$n]*$y[$j] + $z[$n]*$z[$j])/$xj;
      my $yn = ($x[$n]*$x[$k] + $y[$n]*$y[$k] + $z[$n]*$z[$k] - $xn*$xk) / $yk;
      my $zn = sqrt($x[$n]**2 + $y[$n]**2 + $z[$n]**2 - $xn**2 - $yn**2);
      my ($A2,$B2,$C2) = plane_normal($j,$k,$n);
      my $Vec = ($B1*$C2-$C1*$B2)*($x[$k]-$x[$j])+
                ($C1*$A2-$A1*$C2)*($y[$k]-$y[$j])+
                ($A1*$B2-$B1*$A2)*($z[$k]-$z[$j]);
      $zn = -$zn if $Vec < 0;
      #print "$Vec\n";
      ($x[$n],$y[$n],$z[$n]) =  ($xn,$yn,$zn);
    }
  }
  ($x[$j],$y[$j],$z[$j],$x[$k],$y[$k],$z[$k]) =  ($xj,0,0,$xk,$yk,0);
  #&write_molden;
}

sub get_atom_numbers {
  my ($vector,@result);
  $vector = shift;
  $vector =~ s/^\s+//;   # Убираем пробелы в начале
  $vector =~ s/\s+$//;   # Убираем пробелы в конце
  $vector =~ s/\s+/,/g;  # Заменяем пробелы на запятые

  unless ($vector =~ /^(\d+(-\d+)?,)*\d+(-\d+)?$/g) {
    warn "Not valid atom numbers\n";
    return ();
  }

  foreach (split /,/, $vector) {
    push @result, /(\d+)-(\d+)/ ? ($1<=$2 ? $1..$2 : reverse($2..$1)) : $_;
  }
  return @result;
}

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 (defined($atom[$i]) && defined($x[$i]) && defined($y[$i]) && defined($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";
      }
    }
  }
}

sub jacobi {
 ###########################################################################
 ## Собственные числа и векторы симметричной матрицы методом Якоби.
 ## Код из Жориной программы dte (C).
 ## Принимает массив ($Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz).
 ## Возвращает отсортированный (по убыванию собств. чисел) массив двухэлементных
 ## массивов. 1-й элемент - собств. число, 2-й - соответствующий ему вектор.
 ###########################################################################
  my @a = ([ $_[0],$_[1],$_[2] ],
           [ $_[1],$_[3],$_[4] ],
           [ $_[2],$_[4],$_[5] ]);
  my ($j, $iq, $ip, $i, $n);
  my ($tresh, $theta, $tau, $t, $sm, $s, $h, $g, $c);
  my $maxit = 50;

  $n = 3;

  my @v = ([1,0,0], [0,1,0], [0,0,1]);
  my @b = my @d = ($a[0][0],$a[1][1],$a[2][2]);
  my @z = (0,0,0);
  
  my $nrot = 0;
  for ($i=0; $i<$maxit; $i++) {
    
    $sm = abs($a[0][1]) + abs($a[0][2]) + abs($a[1][2]);
    
    if ($sm == 0) {
      return sort {$b->[0] <=> $a->[0]}
              ([$d[0], [$v[0][0],$v[1][0],$v[2][0]]],
               [$d[1], [$v[0][1],$v[1][1],$v[2][1]]],
               [$d[2], [$v[0][2],$v[1][2],$v[2][2]]]);
    }
    # The normal return, which relies on 
    # quadratic convergence to machine underflow. 
    
    $tresh = ($i < 4) ? 0.2*$sm/($n*$n) : 0;
    # ...on the first three sweeps. 
    # ...thereafter. 
    
    for ($ip=0; $ip<$n-1; $ip++) {
      for ($iq=$ip+1; $iq<$n; $iq++) {
        $g = 100*abs($a[$ip][$iq]);
        # After four sweeps, skip the rotation if the off-diagonal element is small. 
        if ($i>4 && abs($d[$ip])+$g == abs($d[$ip]) && abs($d[$iq])+$g == abs($d[$iq])) {
          $a[$ip][$iq] = 0;
        }
        elsif (abs ($a[$ip][$iq]) > $tresh) {
          $h = $d[$iq]-$d[$ip];
          if (abs($h)+$g == abs($h)) {
            $t = $a[$ip][$iq]/$h;
          }  # $t = 1/(2theta) 
          else
          {
            $theta = 0.5*$h/$a[$ip][$iq];
            $t = 1/(abs($theta)+sqrt(1+$theta*$theta));
            if ($theta < 0) {$t = -$t;}
          }
          $c = 1/sqrt(1+$t*$t);
          $s = $t*$c;
          $tau = $s/(1+$c);
          $h = $t*$a[$ip][$iq];
          $z[$ip] -= $h;
          $z[$iq] += $h;
          $d[$ip] -= $h;
          $d[$iq] += $h;
          $a[$ip][$iq] = 0;
          for ($j=0; $j<$ip; $j++) {  # Case of rotations 0 <= $j < p. 
            $g = $a[$j][$ip]; 
            $h = $a[$j][$iq];
            $a[$j][$ip] = $g-$s*($h+$g*$tau);
            $a[$j][$iq] = $h+$s*($g-$h*$tau);
          }
          for ($j=$ip+1; $j<$iq; $j++) {  # Case of rotations p < $j < q. 
            $g = $a[$ip][$j]; 
            $h = $a[$j][$iq];
            $a[$ip][$j] = $g-$s*($h+$g*$tau);
            $a[$j][$iq] = $h+$s*($g-$h*$tau);
          }
          for ($j=$iq+1; $j<$n; $j++) {  # Case of rotations q < $j < $n. 
            $g = $a[$ip][$j]; 
            $h = $a[$iq][$j];
            $a[$ip][$j] = $g-$s*($h+$g*$tau);
            $a[$iq][$j] = $h+$s*($g-$h*$tau);
          }
          for ($j=0; $j<$n; $j++) {
            $g = $v[$j][$ip]; 
            $h = $v[$j][$iq];
            $v[$j][$ip] = $g-$s*($h+$g*$tau);
            $v[$j][$iq] = $h+$s*($g-$h*$tau);
          }
          ++$nrot;
        }
      }
    }
    for ($ip=0; $ip<$n; $ip++) {
      $b[$ip] += $z[$ip];
      $d[$ip] = $b[$ip];      # Update $d with the sum of tapq, 
      $z[$ip] = 0;        # and reinitialize $z. 
    }
  }
  return undef;
}