Новосибирский институт органической химии им. Н.Н. Ворожцова СО РАН

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

tenzor


#!/usr/bin/perl -ws

use List::Util qw(sum);
#use Data::Dump qw(pp);

use vars '$h';
if ($h) {
  (my $program = $0) =~ s/.*[\/\\]//;
  print <<HELP;

   Usage: $program filenames|stdin
   
   Зависимости: Perl

   Переориентация молекул в соответствии с осями эллипсоида инерции.
Осью x становится наибольшая ось эллипсоида, осью z - наименьшая.
   Если различия между соседними молекулами не очень большие, то
каждая последующая молекула как бы подстраивается под предыдущую.
Это бывает нужно для получения гладкого движения по координате реакции.
С этой задачей вообще то лучше справляется программа smooth, но для нее
(в отличие от $program) необходима согласованная нумерация атомов в молекулах,
что не всегда бывает, особенно при стыковке двух ветвей координаты реакции. 
Согласованную нумерацию можно сделать вручную, с помощью программы renum, 
а можно автоматически, с помощью $program -sort.
   Warning! На высокосимметричных молекулах $program глючит.

Опции:
-sort   перенумерация 2-й и последующих молекул с подстройкой под предыдущую
-sort=2 1-я молекула перенумеровывается по Хиллу, затем как при -sort
HELP
  exit;
}

## Определение наличия модулей.
#BEGIN { 
#  my($found, @DBs, $mod);
#  $found = 0; 
#  @DBs = qw(Math::Cephes::Matrix Math::MatrixReal); 
#  for $mod (@DBs) { 
#    if (eval "require $mod") { 
#      $mod->import();
#      $module = $mod; # $module содержит имя первого найденного модуля.
#      $found = 1; 
#      last;
#    }
#  }
#  # Если нет модулей, то печатаем файл/поток и выход
#  unless ($found) {
#    print while (<>);
#    #die "None of @DBs loaded";
#    exit;
#  }
#}

# Массы атомов
%massa = qw(
H  1.007947    He 4.0026022   Li 6.9412      Be 9.0121823   B  10.8117
C  12.01078    N  14.00672    O  15.99943    F  18.99840325 Ne 20.17976
Na 22.9897702  Mg 24.30506    Al 26.9815382  Si 28.08553    P  30.9737612
S  32.0655     Cl 35.4532     Ar 39.9481     K  39.09831    Ca 40.0784
Sc 44.9559108  Ti 47.8671     V  50.94151    Cr 51.99616    Mn 54.9380499
Fe 55.8452     Co 58.9332009  Ni 58.69342    Cu 63.5463     Zn 65.4094
Ga 69.7231     Ge 72.641      As 74.921602   Se 78.963      Br 79.9041
Kr 83.7982     Rb 85.46783    Sr 87.621      Y  88.905852   Zr 91.2242
Nb 92.906382   Mo 95.942      Ru 101.072     Rh 102.905502  Pd 106.421
Ag 107.86822   Cd 112.4118    In 114.8183    Sn 118.7107    Sb 121.7601
Te 127.603     I  126.904473  Xe 131.2936    Cs 132.905452  Ba 137.3277
La 138.90552   Ce 140.1161    Pr 140.907652  Nd 144.243     Sm 150.363
Eu 151.9641    Gd 157.253     Tb 158.925342  Dy 162.5001    Ho 164.930322
Er 167.2593    Tm 168.934212  Yb 173.043     Lu 174.9671    Hf 178.492
Ta 180.94791   W  183.841     Re 186.2071    Os 190.233     Ir 192.2173
Pt 195.0782    Au 196.966552  Hg 200.592     Tl 204.38332   Pb 207.21
Bi 208.980382  Th 232.03811   Pa 231.035882  U  238.028913
);

$count =0; # Счетчик геометрий

# Для каждой геометрии
while (&read_molden) {
  # Создаем массив атомных масс
  $m[$_] = $massa{$atom[$_]} for 1..$N;

  # Вычисляем молекулярную массу
  $M = sum(@m[1..$N]);
  #print "$M\n";

  # Перемещаем геометрию в центр масс
  foreach my $X (qw(x y z)) {
    my $sum = sum(map { ${$X}[$_]*$m[$_] } 1..$N);
    $sum /= $M;                                        
    $_ -= $sum for @{$X}[1..$N];
  }
  #&write_molden;
  
  #  Суть: молекулу представляем в виде эллипсоида инерции.
  #  Длина главных осей этого эллипсоида определяются собственными значениями 
  #  тензора инерции, а их направление - соответствующими собственными векторами. 
  #  Г. Голдстейн. Классическая механика. Глава 5.
  
  # Вычисляем тензор инерции
  # | Ixx Ixy Ixz |
  # | Ixy Iyy Iyz |
  # | Ixz Iyz Izz |
  $r2[$_] = $x[$_]**2+$y[$_]**2+$z[$_]**2 for 1..$N;
  $Ixx = sum(map {$m[$_]*($r2[$_]-$x[$_]**2)} 1..$N);
  $Iyy = sum(map {$m[$_]*($r2[$_]-$y[$_]**2)} 1..$N);
  $Izz = sum(map {$m[$_]*($r2[$_]-$z[$_]**2)} 1..$N);
  $Ixy = -sum(map {$m[$_]*$x[$_]*$y[$_]} 1..$N);
  $Ixz = -sum(map {$m[$_]*$x[$_]*$z[$_]} 1..$N);
  $Iyz = -sum(map {$m[$_]*$y[$_]*$z[$_]} 1..$N);

  #($Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz) = (6,-2,2,5,0,7);
  
  #print "Тензор инерции:\n";
  #printf "%12.8f %12.8f %12.8f\n%12.8f %12.8f %12.8f\n%12.8f %12.8f %12.8f\n",
  #         $Ixx,  $Ixy,  $Ixz,   $Ixy,  $Iyy,  $Iyz,   $Ixz,  $Iyz,  $Izz;
  
#  ############# Вычисления с использованием Math::MatrixReal
#  if ($module eq 'Math::MatrixReal') {
#    # Создаем матрицу-тензор инерции
#    my $matrix = Math::MatrixReal->new_from_rows([[$Ixx,$Ixy,$Ixz],[$Ixy,$Iyy,$Iyz],[$Ixz,$Iyz,$Izz]]);
#
#    # Получаем собственные значения и собственные вектора тензора инерции
#    # Собственные вектора (в столбцах $V) ортонормированы
#    my ($l, $V) = $matrix->sym_diagonalize();
#    #pp to_aref($l); pp to_aref($V);
#    #exit;
#
#    # Преобразуем собственные числа из обьекта Math::MatrixReal в обычный массив
#    @l = to_aref($l);
#    @l = map {$_->[0]} @l;
#
#    # Создаем матрицу-геометрию
#    my $geom = Math::MatrixReal->new_from_cols([[@x[1..$N]],[@y[1..$N]],[@z[1..$N]]]);
#    #print $geom;
#
#    # Преобразуем геометрию к новой системе координат (собственные вектора)
#    my $new_geom = ~($geom * $V);
#    @new_geom = to_aref($new_geom);
#    #pp @new_geom;
#  }
#
#  sub to_aref {
#    (my $Vstr = sprintf($_[0])) =~ s/\n/,/g;
#    $Vstr =~ s/(\d)\s+([\-\d])/$1, $2/g;
#    return eval $Vstr;
#  }
#  ############# Конец использования Math::MatrixReal
#  
#  ############# Вычисления с использованием Math::Cephes::Matrix
#  if ($module eq 'Math::Cephes::Matrix') {
#    # Создаем матрицу-тензор инерции
#    my $matrix = Math::Cephes::Matrix->new([[$Ixx,$Ixy,$Ixz],[$Ixy,$Iyy,$Iyz],[$Ixz,$Iyz,$Izz]]);
#
#    # Получаем собственные значения и собственные вектора тензора инерции
#    # Собственные вектора (в столбцах $V) ортонормированы
#    my ($l, $V1) = $matrix->eigens();
#    my $V = ($V1->transp)->coef;
#    #pp $l; pp $V;
#    # Массив собственных чисел
#    @l = @$l;
#
#    # Преобразуем геометрию к новой системе координат (собственные вектора)
#    for ($i=0; $i<$N; $i++) {
#      for (0..2) {
#        $new_geom[$_][$i] = $x[$i+1]*$V->[0][$_] + 
#                            $y[$i+1]*$V->[1][$_] + 
#                            $z[$i+1]*$V->[2][$_];
#      }
#    }
#    #write_molden();
#  }
#  ############# Конец использования Math::Cephes::Matrix
  
  my @eigen = jacobi($Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz) or die "Eigen error on geom $count\n";
  #pp @eigen;
  # Преобразуем геометрию к новой системе координат (собственные вектора)
  for ($i=0; $i<$N; $i++) {
    for (0..2) {
      $new_geom[$_][$i] = $x[$i+1]*$eigen[$_][1][0] + 
                          $y[$i+1]*$eigen[$_][1][1] + 
                          $z[$i+1]*$eigen[$_][1][2];
    }
  }
  #write_molden();
  
#  # Сортируем индексы массива собственных значений 
#  # в порядке убывания значений массива
#  my @a;
#  $a[$_] = [$l[$_],$_] for 0..$#l;
#  @a = sort {$a->[0] <=> $b->[0]} @a;
#  @ind = map {$_->[1]} @a;
  
  # Переставляем оси эллипсоида, чтобы хорошо смотрелось в молдене
  @x[1..$N] = @{$new_geom[1]};
  @y[1..$N] = @{$new_geom[2]};
  @z[1..$N] = @{$new_geom[0]};
  
  # Дальше - самое сложное. Нужно согласовать направления осей текущей
  # геометрии с направлениями предыдущей. Эти направления после уже 
  # проделанных преобразования могут отличаться только знаками. Сложность
  # же заключается в том, что порядок нумерации атомов в двух геометриях
  # может быть различен. Но предполагаем, что две соседние геометрии 
  # не сильно отличаются друг от друга
  
  # Создаем структуру данных "молекула", 
  # объединяющую символы атомов, их порядок в геометрии и координаты.
  # Это хэш (атомных символов) хэшей (их номера от 1 до N) массивов (координат)
  my %mol;
  foreach (1..$N) {
    $mol{$atom[$_]}->{$_} = [$x[$_],$y[$_],$z[$_],$ppm[$_]];
    #$mol{$atom[$_]}->{$_} = [sprintf('%.12f',$x[$_]),sprintf('%.12f',$y[$_]),sprintf('%.12f',$z[$_]),$ppm[$_]];
  }
  #warn "___________________________\nИсходная геометрия\n";
  #pp $mol{C};
  
  # Если опция сортировки = 2
  # Сортировка порядка атомов 1-й молекулы по системе Хилла
  if ($sort && $sort == 2 && $count == 0) {
    # Создаем вспомогательный массив массивов [атом,индекс]
    my @at;
    $at[$_-1] = [$atom[$_],$_] for 1..$#atom;
    #pp @at;

    # Сортируем его по атомам с приоритетом C и H
    @at = sort {  
      return -1 if uc($a->[0]) eq 'C' and uc($b->[0]) ne 'C';
      return  1 if uc($a->[0]) ne 'C' and uc($b->[0]) eq 'C';
      return -1 if uc($a->[0]) eq 'H' and uc($b->[0]) ne 'H';
      return  1 if uc($a->[0]) ne 'H' and uc($b->[0]) eq 'H';
      uc($a->[0]) cmp uc($b->[0]);
    } @at;
    # (["N", 1], ["H", 2], ["F", 3], ["Cl", 4]) @at до сортировки
    # (["H", 2], ["Cl", 4], ["F", 3], ["N", 1]) @at после сортировки
    # (2, 4, 3, 1)                              map {$_->[1]} @at
    #pp @at;
    
    # Переносим полученные координаты из "молекулы" %mol в массивы @x,@y,@z,@ppm
    my $i = 1;
    foreach $j (map {$_->[1]} @at) {
      foreach my $atom (keys %mol) {
        if (exists $mol{$atom}{$j}) {
          ($x[$i],$y[$i],$z[$i],$ppm[$i]) = @{$mol{$atom}{$j}};
          $atom[$i++] = $atom;
        }
      }
    }
    # Пересоздаем структуру данных "молекула"
    undef %mol;
    foreach (1..$N) {
      $mol{$atom[$_]}->{$_} = [$x[$_],$y[$_],$z[$_],$ppm[$_]];
    }
    
    $sort = 1; # Остальные молекулы будут сортироваться по 1-й
  }

  # Если не первая геометрия
  if ($count != 0) {
    # Вычисляем "кинетич. энергию" между текущей и предыдущей молекулами
    # В возвращаемом массиве 0-й элемент - это "кинетич. энергия", остальные - 
    # пермутационный вектор, определяющий порядок нумерации текущей молекулы,
    # при котором она лучше всего накладывается на предыдущую.
    my @jmin = kinen(\%mol,\%mol0);
    
    # Далее будем последовательно менять знаки координат текущей молекулы
    # (отражение) и смотреть, при каком отражении будет мин. кин. энергия
    # Последовательность отражений:
    #  0 +x +y +z    
    #  1 +x -y +z    
    #  2 -x -y +z    
    #  3 -x +y +z    
    #  4 -x +y -z    
    #  5 +x +y -z    
    #  6 +x -y -z    
    #  7 -x -y -z    
    
    my $reflect = 0;    # Счетчик отражений
    my $reflectmin = 0; # Номер отражения с минимальн. кинетич. энергией
    
    # Для последовательности отражений, подобранной так, чтобы меняя каждый раз 
    # знак только одной координаты (0 - x, 1 - y, 2 - z), охватить все варианты
    foreach (1,0,1,2,0,1,0) {
      # Меняем знак соответствующей координаты
      znak(\%mol, $_);
      # Вычисляем кин. энергию
      @j = kinen(\%mol,\%mol0);
      #print "@j\n";
      $reflect++;  # Счетчик отражений
      # Смотрим, не лучше ли очередная кин. энергия
      if ($j[0] < $jmin[0]) {
        @jmin = @j;
        $reflectmin = $reflect;
      }
    }
    #warn "$reflectmin @jmin\n";

    # Вспомогательный массив, определяющий, как нужно поменять знаки координат,
    # чтобы из молекулы, полученной после последнего отражения (-x -y -z), 
    # восстановить молекулу после i-го отражения
    my @ar = ([0,1,2],[0,2],[2],[1,2],[1],[0,1],[0],[]);
    
    # Меняем знаки, чтобы получить наименьшую кин. энергию 
    znak(\%mol, @{$ar[$reflectmin]});
    
    # Переносим полученные координаты из "молекулы" %mol в массивы @x,@y,@z
    # в зависимости от значения опции -sort
    if ($sort) {       # Сортировка по предыдущей геометрии
      my $i = 1;
      #print "@jmin\n";
      foreach $j (@jmin[1..$#jmin]) {
        foreach my $atom (keys %mol) {
          if (exists $mol{$atom}{$j}) {
            ($x[$i],$y[$i],$z[$i],$ppm[$i]) = @{$mol{$atom}{$j}};
            $atom[$i++] = $atom;
          }
        }
      }
    }
    else {                           # Без сортировки
      foreach my $atom (keys %mol) {
        foreach (keys %{$mol{$atom}}) {
          ($x[$_],$y[$_],$z[$_]) = @{$mol{$atom}{$_}};
        }
      }
    }
  }
  
  # Делаем текущую молекулу предыдущей
  foreach (1..$N) {
    $mol0{$atom[$_]}->{$_} = [$x[$_],$y[$_],$z[$_],$ppm[$_]];
    #$mol0{$atom[$_]}->{$_} = [sprintf('%.12f',$x[$_]),sprintf('%.12f',$y[$_]),sprintf('%.12f',$z[$_]),$ppm[$_]];
  }
  #warn "Преобразованная геометрия\n";
  #pp $mol{C};
  
  $count++; # Счетчик геометрий
  
  # Печатаем
  &write_molden;
}

# Менять знаки координат в структуре $h{$atom}{$index}[$xyz]
# Usage: znak(\%h, 0|1|2); x - 0, y - 1, z - 2
sub znak {
  my %h = %{$_[0]} and shift;
  foreach my $xyz (@_) {
    foreach my $atom (keys %h) {
      for (keys %{$h{$atom}}) {
        $h{$atom}{$_}[$xyz] *= -1;
      }
    }
  }
}

sub copymol {
  my $molref = $_[0];
  my $newmolref;
  foreach my $atom (keys %{$molref}) {
    while (my ($ind,$aref) = each %{$molref->{$atom}}) {
      $newmolref->{$atom}->{$ind} = [@{$aref}];
    }
  }
  return $newmolref;
}


sub kinen {
  # Делаем независимую копию 1-й (текущей) молекулы
  my $molref  = copymol($_[0]);
  my %mol  = %{$molref};
  # 2-я (предыдущая) молекула
  my %mol0 = %{$_[1]};
  
  my $E = 0;
  
  # Для углеродов, водородов и т.д.
  foreach my $atom (sort keys %mol0) {
    
    # Если в молекулах разное кол-во, например, углеродов, то выход
    if (scalar(keys %{$mol{$atom}}) != scalar(keys %{$mol0{$atom}})) {
      warn "Not isomers ($atom)\n";
      return;
    }
    
    # Для каждого из атомов данного сорта 2-й молекулы
    foreach my $i (sort {$a<=>$b} keys %{$mol0{$atom}}) {
      my $rmin = 1e10;
      # перебираем атомы того же сорта 1-й молекулы
      foreach my $j (keys %{$mol{$atom}}) {
        # вычисляем расстояние между атомами двух молекул
        $r = sum(map {($mol{$atom}{$j}[$_] - $mol0{$atom}{$i}[$_])**2} 0..2);
        # запоминая ближайший атом 1-й молекулы
        if ($r<$rmin) {
          $rmin = $r;
          $jmin = $j;
        }
      }
      # Сумимируем кинетическую энергию
      $E += $rmin*$massa{$atom};
      # i-му атому 2-й молекулы ставим в соответствие ближайший атом из 1-й
      $j[$i] = $jmin;
      # и удаляем его из 1-й молекулы
      delete $mol{$atom}{$jmin};
    }
    # Записываем накопленную кин. энергию в начало массива
    $j[0] = $E;
  }
  # Возвращаем массив
  return @j;
}

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";
      }
    }
  }
}

#Собственные числа и векторы симметричной матрицы методом Якоби.
#Код из Жориной программы dte (C).
#Принимает массив ($Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz).
#Возвращает отсортированный (по убыванию собств. чисел) массив двухэлементных
#массивов. 1-й элемент - собств. число, 2-й - соответствующий ему вектор.
sub jacobi 
{
  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]);
#  for ($ip=0; $ip<$n; $ip++)  # Initialize to the identity matrix.
#  {
#    for ($iq=0; $iq<$n; $iq++) {$v[$ip][$iq] = 0;}
#    $v[$ip][$ip] = 1;
#  }
  
  my @b = my @d = ($a[0][0],$a[1][1],$a[2][2]);
  my @z = (0,0,0);
#  for ($ip=0; $ip<$n; $ip++)
#  {
#    $b[$ip] = $d[$ip] = $a[$ip][$ip];  # Initialize $b and $d to the diagonal of $a.
#    $z[$ip] = 0;         # This vector will accumulate terms of the form tapq 
#  }
  
  my $nrot = 0;
  for ($i=0; $i<$maxit; $i++) {
    
    $sm = abs($a[0][1]) + abs($a[0][2]) + abs($a[1][2]);
#    $sm = 0;        # Sum off-diagonal elements. 
#    for ($ip=0; $ip<$n-1; $ip++) {
#      for ($iq=$ip+1; $iq<$n; $iq++) {
#        $sm += abs ($a[$ip][$iq]);
#      }
#    }
    
    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]]]);
    }
#    if ($sm == 0) {return 0;}  # The normal return, which relies on 
#      # quadratic convergence to machine underflow. 
    
    $tresh = ($i < 4) ? 0.2*$sm/($n*$n) : 0;
#    if ($i < 4) {$tresh = 0.2*$sm/($n*$n);}  # ...on the first three sweeps. 
#    else {$tresh = 0;}    # ...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;
}