#!/usr/bin/perl -ws

our($asymm,$debug,$all);

$debug ||= 0;
$all ||= 0;
#use Data::Dump 'pp';

# Если нет параметров, то печатаем справку
if (! @ARGV or $ARGV[0] eq '-h') {
  (my $program = $0) =~ s/.*[\/\\]//;
  print "\n Usage: $program file [files]\n
Зависимости: perl\n
Этот перл-скрипт читает выдачи 
квантово-химической программы ПРИРОДА (автор Д.Н. Лайков) и делает файлы, 
понятные molden'у и заготовки для следующего расчета.
Проверенные версии ПРИРОДЫ: p210, p407, p6, p11 (linux), p202 (Windows).\n
Создаются следующие файлы:
file.xyz    декартовы координаты с зарядами для всех шагов оптимизации
            или для каждой точки irc;
file.mos    собственные вектора (орбитали) в формате molden'а
            (почему-то обрабатываются molden'ом очень медленно);
file.freq   колебания в формате molden'а
            (при task=hessian);
file.MOL    оптимизированная геометрия и гессиан, 
            которые можно вставить в след. расчет (не для molden'а);
file.ppm    хим. сдвиги (sigma и, если посчитан стандарт, то и delta)
            (не для molden'а);
file.xyzppm геометрия для molden'а, с хим. сдвигами (delta) вместо зарядов
            (molden желательно немного переделать, чтобы он понимал 
            трехзначные \"заряды\" (знаю, как)).\n\n";
  exit;
}

# Вгоняем таблицу Менделеева
@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
           );

%HoF = (
PBE => {L1=>{C=>-38.0843685585,
             H=>-0.5872542349,
             N=>-54.7271355293,
             O=>-75.1115133857},
       },
);
my $num = qr/-?\d+\.\d+/;


NEW_FILE:
foreach $file (@ARGV) {
  my (@xXX,@yXX,@zXX);
  my %formula;
  my (@thermo,@intensity);
  my @crit_points;
  undef $hessian;
  # Открываем файлы
  $ppm    = "$file.ppm";
  $mos    = "$file.mos";
  $freq   = "$file.freq";
  $mol    = "$file.MOL";
  $xyz    = "$file.xyz";
  $xyzppm = "$file.xyzppm";
  open (INP,         $file) || die "Cannot open input file $file";
  open (MOS,       ">$mos") || die "Cannot open input file $mos";
  open (FREQ,     ">$freq") || die "Cannot open input file $freq";
  open (MOL,       ">$mol") || die "Cannot open input file $mol";
  open (XYZ,       ">$xyz") || die "Cannot open input file $xyz";
  open (PPM,       ">$ppm") || die "Cannot open input file $ppm";
  open (XYZPPM, ">$xyzppm") || die "Cannot open input file $xyzppm";
  
  # по дефолту вид расчета и
  # ограничитель, после нахождения которого нужно печатать геометрию
  # (для IRC и SCAN печатаем только оптимизированные на каждом шаге геометрии) 
  $task = 'optimize';
  $delimiter = 'eng>\$end';

  # scan input file
  my $opt_conv = 0;
  while (<INP>) {
    
    # program  Priroda  version 4.07 (14.09.2004)
    # program  Priroda  version 2.10 (14.01.2002)
    # Priroda 6 (2006.08.20)
    if (m/\bPriroda\b\D*(\d+)/) {
      $pri_ver = $1;
    }
    
    if (/THERMOSTAT/) {
      close MOS; close FREQ; close MOL; close PPM;  close XYZPPM;
      unlink($mos,$freq,$mol,$ppm,$xyzppm);
      seek (INP, 0, 0);
      thermostat(); 
      close INP; close XYZ;
      next NEW_FILE;
    }
    if (/Intrinsic Reaction Coordinate/) {$task = 'irc'};
    if (/SCAN\>/ | /task=scan/)          {$task = 'scan'};
    if (/OPTIMIZATION CONVERGED/)        {$opt_conv++};
    if (/NMR Shielding Tensors/)         {$task = 'nmr'};
    if (/COMPUTATION OF SECOND DERIVATIVES/)         {$task = 'hessian'};
    #print "$task\n";
    #if (/Calculation of NMR Chemical Shifts/)         {$task = 'nmr'};

    # Считываем название основного базиса
    $basis  = $1 if m/^ Basis set input: '(\S+)'/i;
    $basis =~ s/.*\/// if $basis;
    # Дополнительный базис
    $basis2 = $1 if m/^ Basis set input: '(\S+)'/i;
    $basis2 =~ s/.*\/// if $basis;

    # Считываем число атомов
    $natoms = $1 if m/^\s*(\d+)\s+atoms,\s*\d+\s+electrons/;
    $natoms = $1 if m/Number of atoms\s+(\d+)/;
    
    # Заряд молекулы и ее мультиплетность
    $Charge = $1 if /Charge of molecule\s+(-?\d+)/;
    $Mult = $1 if /Spin multiplicity\s+(\d+)/;

#    if (/molecule input:/) { 
#      $_=<INP>; 
#      ($natoms, undef)=split; 
#    }
    # Функционал
    $functional = $1 if m/Approximation to E\(xc\)\s+(\S+)/;
    
    if (m/formula: (\w+)/) {
      my $brutto = $1;
      $brutto =~ s/([A-Z][a-z]?)(?!\d)/${1}1/g;
      %formula = split /(\d+)/, $brutto;
    }
    
    $ZPVE = $1 if /ZPVE =\s*(\S+)/;
    if (m/^\s*T =\s*(\S+)\s+K/) {
      my $T = $1;
      $T =~ s/0+$//; $T =~ s/\.$//;
      push @thermo, {T => $T};
      <INP>;<INP>;<INP>;<INP>;<INP>;
      $_ = <INP>;
      my @a = split;
      $_ = sprintf("%.2f",$_) for @a[3..5];
      $thermo[-1]{lnQ} = $a[1];
      $thermo[-1]{S} = $a[2];
      $thermo[-1]{U} = $a[3];
      $thermo[-1]{H} = $a[4];
      $thermo[-1]{G} = $a[5];
    }
    
    if (m/qm file:/) {
      $qm = 1;
      $delimiter = 'mol>\$end';
    }
  }
  
  if ($task ne 'irc' && $opt_conv > 1) {
    $task = 'scan'
  }
  
  if ($task eq 'irc' || $task eq 'scan') {
    $delimiter = 'OPTIMIZATION CONVERGED';
  }
  if ($task eq 'nmr') {
    $delimiter = 'Atomic Coordinates:';
  }
  
  if ($all and $task eq 'irc' || $task eq 'scan') {
    $delimiter = 'eng>\$end';
    $delimiter = 'mol>\$end' if $qm;
  }

  print "\$natoms = $natoms\n" if ($debug >=1);
  print "\$task = $task\n" if ($debug >=1);
  print "\$theory = QM\n" if ($debug >=1 and $qm);
  print "\$delimiter = $delimiter\n" if ($debug >=1);
  
  # Переходим в начало файла
  seek (INP, 0, 0);
    
  # (флаг, найдены ли заряды и спины)
  $mulOK=0;

  # Вновь сканируем файл  
  while (<INP>) {
    # Читаем координаты
    if (/Atomic Coordinates:/) { 
      for ($i=1; $i<=$natoms; $i++) {
        $_=<INP>;
        ($charge[$i], $x[$i], $y[$i], $z[$i]) = split; 
      }
    }
    #print "@charge\n" if ($debug >=1); 

    # Считываем энергию
    $energy=$1 if $qm && /^ Energy:\s*(\S+)/;
    $energy=0+$1 if m/^  E  =\s*(\S+)/;
    $energy=0+$1 if m/^eng> E=(\S+)/;
    $dipole = $1 if m/ dipole = .+? (\S+) D/;
    if ($task eq 'irc') {
      if (m/^IRC> point.+? g = (\d+\.\d+)/) {
        $gradient = $1;
      }
    }

    # Читаем заряды и спины
    if (m/atom    charge    spin\s*/) { 
      $mulOK=1;
      for ($i=1; $i<=$natoms; $i++) {
        $_=<INP>;
        (undef, undef, $mulcharge[$i], $spin[$i]) = split; 
      }
    }
    
    if (m/Critical Points of the Density/) {
      my @tmp;
      while (<INP>) {
        last if m/Atomic Coordinates:/;
        if (m/r=\s*($num),\s*($num),\s*($num)\s*f=($num)/) {
          push @tmp, ['XX', $1,$2,$3,$4];
        }
      }
      push @crit_points, \@tmp;
    }
    
    if (m/Geometry Optimization Step (\d+)/) {
      if (defined($step) && $step > $1) {
        #print "$step > $1\n";
        $new_optimization = 1;
      }
      $step = $1;
    }
    
    $geo_opt = 1 if m/^ OPTIMIZATION CONVERGED/;
    
    # Печатаем число атомов, энергию и координаты 
    # (а также зарядамы и спины, если они найдены)
#    if (/$delimiter/ or ($new_optimization and !$geo_opt)) {
    if (/$delimiter/) {
      #print "($new_optimization and !$geo_opt)\n";
      $new_optimization = 0;
      $geo_opt = 0;
      printf XYZ " %d\n", @crit_points ? $natoms+@{$crit_points[-1]} : $natoms;
      print XYZ " Energy $energy";
      print XYZ "  Charge $Charge" if $Charge;
      print XYZ "  Mult $Mult" if $Mult && $Mult != 1;
      print XYZ "  Dipole $dipole" if $dipole; undef $dipole;
      if ($task eq 'irc' and ! $all and $gradient) {
        print XYZ " Gradient $gradient";
      }
      if ($ZPVE) {
        print XYZ "  ZPE $ZPVE";
      }
      if (@thermo) {
        foreach (@thermo) {
          print XYZ "  G($_->{T}) $_->{G}";
        }
      }
      if ($energy && $functional) {
        if (! %formula) {
          foreach my $at (@charge) {
            next unless $at;
            #print "$at\n";
            $formula{$ATOM[$at]}++;
          }
        }
        my $HoF = $energy;
        foreach my $element (keys %formula) {
          if (! exists $HoF{$functional}{$basis}{$element}) {
            undef $HoF;
            last;
          }
          $HoF -= $HoF{$functional}{$basis}{$element}*$formula{$element};
        }
        if (defined $HoF) {
          printf XYZ "  HoF %.2f kcal", $HoF*627.5;
        }
      }
      print XYZ "\n";
      #print XYZ ($mulOK ? "                         Charge    Spin\n" : "\n");
      if ($mulOK) {
        for ($i=1; $i<=$natoms; $i++)  { 
          printf XYZ "%2.4s %12.8f %12.8f %12.8f %10.4f %8.4f\n", $ATOM[$charge[$i]], $x[$i], $y[$i], $z[$i], $mulcharge[$i], $spin[$i]; 
        } 
      } 
      else {  
        for ($i=1; $i<=$natoms; $i++)  { 
          printf XYZ "%2.4s %12.8f %12.8f %12.8f\n", $ATOM[$charge[$i]], $x[$i], $y[$i], $z[$i]; 
        } 
      }
      if (@crit_points) {
        my $aref = pop @crit_points;
        #pp $aref;
        for ($i=0; $i<@$aref; $i++)  { 
          printf XYZ "%2.4s %12.8f %12.8f %12.8f %10.4f\n", @{$aref->[$i]}; 
        } 
      }
    }
    
    # При простой оптимизации геометрии (в т.ч. saddle)
    if ($task eq 'optimize') {
      print XYZ if (/OPTIMIZATION CONVERGED/);
    }
    
    # Создаем файлы орбиталей и частот колебаний и геометрия+гессиан
    print MOS  if s/mos>(.*$)/$1/;
    print MOL  if s/MOL>(.*$)/$1/;
    $hessian=1 if (/SECOND DERIVATIVES/);
    if ($hessian) {
      print MOL  if s/eng>(.*$)/$1/;
    }
    $vibration=1 if (/^ VIBRATIONAL ANALYSIS and THERMOCHEMISTRY/);
    if ($vibration && $pri_ver) {
      if ($pri_ver == 2) {
        push @intensity, split if s/^ ir i\.//;
      }
#      if ($pri_ver == 4 or $pri_ver == 6) {
      if ($pri_ver >= 4) {
        if (/^ \| Mode \| Freq\.   \| Mass\. \| IR Int\./) {
          <INP>;
          while (($_=<INP>) =~ /^\s*\|\s*\d+\s*\|/) {
            push @intensity, (split /\|/,$_)[4];
          }
        }
      }
      print FREQ if s/freq>//;
    }

    # Разбираемся с хим. сдвигами
    if (/Calculation of NMR Chemical Shifts/) {
      <INP>;
      if (<INP> =~ /additional points/) {
        while (my($x,$y,$z) = split ' ', <INP>) {
          push @xXX, $x;
          push @yXX, $y;
          push @zXX, $z;
        }
      }
    }
    if (/NMR Shielding Tensors/) {
      
      %TMS = (
        'PBE' => {
          'iglo.bas'     =>{'C'=>181.7934,'H'=>31.567, 'Si'=>332.2421},
          'II-ccpv3z.bas'=>{'C'=>181.7934,'H'=>31.567, 'Si'=>332.2421},
          'iglo3z.bas'   =>{'C'=>182.088, 'H'=>31.599, 'Si'=>334.1568},
          'II-3z.bas'    =>{'C'=>182.088, 'H'=>31.599, 'Si'=>334.1568},
          'L2.bas'       =>{'C'=>182.94-7.16,  'H'=>31.179-0.031, 'Si'=>352.9499},
          'L1.bas'       =>{'C'=>182.94-7.16,  'H'=>31.179-0.031, 'Si'=>352.9499},
          'L22.bas'      =>{'C'=>182.94-10,  'H'=>31.179-0.031,  # TMS-empir.
                            'Si'=>352.9499, # TMS
#                            'N'=>263.2636,  # NH3 L1 geom.
                            'N'=>228.74,    # THEOCHEM
                            'O'=>186.57+77.3, # PhOH [OMR 1993, 31, 489-494]
                            },
          '3z.bas'       =>{'C'=>182.000, 'H'=>31.392, 'Si'=>328.8347, 
                            'F'=>318.91,  'N'=>-134.594,'O'=>327.9538,
                            'Al'=>555.98, 'Ti'=>-1085, 'Cl'=>-20.2,
                            #'B'=> 94.05,       # sigma(BF3-Et2O)
                            'B'=>[-0.935433,92.888],
                            'P'=>277.08,       # sigma(H3PO4)
                            'P'=>353.96+(-62), # sigma(PMe3)+delta(PMe3)
                             },
          '3z_sbk.bas'   =>{'C'=>182.000, 'H'=>31.392, 'Si'=>328.8347, 
                            'F'=>318.91,  'N'=>-134.594,'O'=>327.9538,
                            'Al'=>555.98, 'Ti'=>-1085, 'Cl'=>-20.2,
                            'B'=> 93.90,       # sigma(BF3-Et2O)
                            'P'=>277.08,       # sigma(H3PO4)
                            'P'=>353.96+(-62), # sigma(PMe3)+delta(PMe3)
                             },
          'ccpv3z.bas'   =>{'C'=>183.089, 'H'=>31.199, 'Si'=>335.1516},
          'ccpv4z.bas'   =>{'C'=>180.54,  'H'=>31.116, 'Si'=>344.3229},
          'III-3z.bas'   =>{'C'=>188.0585+2.1, 'H'=>31.288+0.233},
          'IV-3z.bas'    =>{'C'=>188.1792+2.1, 'H'=>31.1894+0.233}
        },
        'B3LYP' => {
          '3z.bas' => {'C'=>166.214, 'H'=>29.602, 'Si'=>289.219}
        }
      );
      
      foreach my $f (keys %TMS) {
        foreach my $b (keys %{$TMS{$f}}) {
          if ($b =~ /(.+)\.bas$/) {
            $TMS{$f}{$1} = $TMS{$f}{$b};
          }
        }
      }
      
      # print header
      $basis  =~ s/.*[\/\\]//;
      $basis2 =~ s/.*[\/\\]//;
      
      print PPM " basis set '$basis/$basis2' \n";
      print PPM "  Atom      SIGMA  DELTA\n";
      print XYZPPM '  ',$natoms+@xXX,"\n";
      print XYZPPM " Energy $energy\n" if $energy;

      # Печатаем хим. сдвиги
      for ($i=1; $i<=$natoms; $i++) {
        $_=<INP>; # читаем строку
        #(undef, $number, $atom, undef, undef, $sigma) = split;
        ($number,$atom,$sigma,$principal) = m/\s*Atom\s*(\d+)\s+(\w+).+?(-?\d+\.\d+)(.+)/;
        if ($asymm) {
          # Journal of Molecular Structure: THEOCHEM 953 (2010) 83-90
          my ($s11,$s22,$s33) = (split ' ', $principal)[1..3];
          # Int. J. Mol. Sci. 2002, 3, ~890
          ($s33,$s11,$s22) = sort {abs($b-$sigma)<=>abs($a-$sigma)} ($s11,$s22,$s33);
          # после сортировки нмчего не изменилось
          #print "$s33,$s11,$s22 - $sigma\n";
          $sigma += ($s22-$s11)/($s33-$sigma);
        }
        $atom =~ s/(.+):/$1/; # удаляем двоеточие в названии атома
        printf PPM "%4d %-4.4s %8.4f", $number, $atom, $sigma;
        printf XYZPPM "%2.4s %12.8f %12.8f %12.8f", $atom, $x[$i], $y[$i], $z[$i]; 
        if (exists $TMS{$functional}{$basis}{$atom}) {
          my $tms = $TMS{$functional}{$basis}{$atom};
          my $ppm;
          if (ref($tms) eq 'ARRAY') {
            $ppm = $sigma*$tms->[0]+$tms->[1];
          }
          else {
            $ppm = $tms - $sigma;
          }
          printf PPM " %10.4f", $ppm;
          printf XYZPPM "%10.4f", $ppm; 
        }
        printf PPM "\n";
        printf XYZPPM "\n";
        <INP>; <INP>; <INP>; #Пропускаем три строки
      }
      # Печатаем NICS
      for ($i=0; $i<@xXX; $i++) {
        $_=<INP>; # читаем строку
        ($sigma) = /isotropic=\s+(\S+)/;
        printf PPM "%4d %-4.4s %8.4f %10.4f\n", $natoms+$i+1, 'XX', $sigma, -$sigma;
        printf XYZPPM "%2.4s %12.8f %12.8f %12.8f %10.4f\n", 
                      'XX', $xXX[$i], $yXX[$i], $zXX[$i], -$sigma; 
        <INP>; <INP>; <INP>; #Пропускаем три строки
      }
    }
  }
    
  # Если не закрыть файлы, то они могут не (сразу) записаться на диск
  close XYZ; close MOS; close FREQ; close MOL; close PPM;  close XYZPPM;

  # Пустые файлы удаляем
  if (-z $xyz)    {unlink $xyz}
  if (-z $ppm)    {unlink $ppm}
  if (-z $xyzppm) {unlink $xyzppm}
  if (-z $freq)   {unlink $freq}
  if (-z $mos)    {unlink $mos}
  if (-z $mol)    {unlink $mol}
  
    # Дописываем интенсивности ИК 
  if (@intensity) {
    if (-f $freq) {
      open (FREQ,">>$freq") || die "Cannot open $freq to add";
      print FREQ "[INT]\n";
      foreach (@intensity) {
        print FREQ " $_\n";
      }
      close FREQ;
    }
  }
#  if (@crit_points) {
#    pp @crit_points;
#  }
}
 sub thermostat {
   my (@mols,@trj);
   while (<INP>) {
     if (/^ point /) {
       my @t = /(-?\d+(?:\.\d+)?)/g;
       push @trj, \@t;
     }
     if (/Atomic Coordinates:/) {
       my @mol;
       $mol[0] = $trj[-1] if @trj;
       while (<INP>) {
         last if /#/;
         s/\s*(\d+)\s/ $ATOM[$1] /;
         push @mol, $_;
       }
       push @mols, \@mol;
     }
   }
   unshift @{$mols[0]}, $trj[0];
   
   my $N = $#{$mols[0]};
   foreach my $mol (@mols) {
     print XYZ "$N\n";
     print XYZ "  Energy $mol->[0][2]  point $mol->[0][0]  s $mol->[0][1]  g $mol->[0][3]\n";
     shift @$mol;
     print XYZ @$mol;
   }
   
   open TRJ, '>', "$file.trj" or die "Can't open $file.trj: $!\n";
   foreach (@trj) {
     printf TRJ "%4d %10f %14.8f %12.8f\n", $_->[0], $_->[1], $_->[2], $_->[3];
   }
   close TRJ;
 }