#!/usr/bin/perl -ws

#use Data::Dump 'pp';
our ($atom,$no_del,$gen,$file,$h,$help);
$| = 1;
(my $script = $0) =~ s/.*\///;

if ($h || $help) {
  print <<HELP;

Генерирование базисов (включая auxiliary) для Природы с помощью orca или g09

Usage: $script [-atom=1,4..9] basis_name
Зависимости: perl, orca или g09

Базисы серии L достаточны для решения многих задач, но они нигде, кроме как в Природе
не используются. Между тем, изредка возникает потребность в других, более стандартных
базисах. Например, чтобы воспроизвести расчет, сделанный в другой программе.

Этот скрипт генерирует для Природы базисы, доступные как стандартные для программ 
orca или g09.  Эти программы запускается для каждого из элементов с ключевыми словами 
ORCA -     '! PBE basis_name PrintBasis NoIter'
гауссиан - '# test BLYP/basis_name DenFit GFInput Guess=Only Pop=None'

Опции:
-gen=orca|g09  Программа-генератор базисов. Default -gen=orca 
               (генерируемые ею дополнительные базисы оказались лучше для Природы).
-atom=1..104   Номера элементов, для которых нужны базисы. Default -atom='1,4..9'
-no_del        Не удалять инпуты/аутпуты orca и g09 (для тестирования)
-file=file_name  Названия базисов берутся из файла (бьется по пробелам/строкам).
HELP
  exit;
}

$atom ||= '1,4..9';
$gen ||= 'orca';
if ($gen !~ /^(g09|orca)$/) {
  die "-gen=g09|orca\n";
}

if ($file && -f $file) {
  open F, '<', $file or die "Can't open $file: $!\n";
  my @l = <F>;
  close F;
  foreach (@l) {
    push @ARGV, split;
  }
}

die "usage $script [-atom=1,4..9] basis_name\n" unless @ARGV;

my $g09 = 'rung09';
die "No $g09\n" if $gen eq 'go9' && ! -x $g09;
my $orca = '/usr/local/bin/orca';
die "No $orca\n" if $gen eq 'orca' && ! -x $orca;

# Вгоняем таблицу Менделеева
my @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 Rf Db Sg Bh Hs Mt Ds Rg
           );

my %SPD = ('S'=>0,'P'=>1,'D'=>2,'F'=>3,'G'=>4,'H'=>5,
           'I'=>6,'J'=>7,'K'=>8,'L'=>8,'M'=>10,'N'=>11);

BASIS:
foreach my $basis_name (@ARGV) {
  print "$basis_name: ";
  my $info = "  Basis and auxiliary basis sets for PRIRODA 
  generated by $script script from $gen program outputs\n";
  $info .= sprintf(" %10s%-30s%s\n",'','basis','basis2') if $gen eq 'orca';
  
  my $basis = "\$basis\ntype=general\ndefault=$basis_name\n\nset=$basis_name\n";
  my $basis2 = "\$basis2\ntype=general\ndefault=$basis_name\n\nset=$basis_name\n";
  ELEMENT:
  foreach my $element_number (eval $atom) {
    my $charge = $element_number % 2;
    $charge = 1-$charge if $gen eq 'orca';
    my $element = $ATOM[$element_number];
    print "$element ";
    if ($gen eq 'g09') {
      open L, '>', "$element.com" or die "Can't wtite $element.com: $!\n";
      print L "# test BLYP/$basis_name DenFit GFInput Guess=Only Pop=None\n\n$element\n\n$charge 1\n$element\n\n\n";
      close L;

      system($g09, "$element.com");
  #    my @args = ($g09, "$element.com");
  #    system(@args) == 0 or die "\nSystem\n@args\nfailed: $?\n";

      open L, '<', "$element.log" or die "Can't open $element.log: $!\n";
      my (@g09_basis, @g09_fit_basis);
      while (<L>) {
        if (m/ Atomic number out of range/) {
          print "\b"x(length($element)+1);
          unlink "$element.com", "$element.log" unless $no_del;
          next ELEMENT;
        }
        if (m/ No stored fitting set matched to this AO basis/) {
          print "No stored fitting set matched to this AO basis\n";
          unlink "$element.com", "$element.log" unless $no_del;
          last ELEMENT;
        }
        if (m/ AO basis set in the form of general basis input/ .. m/ \*\*\*\*/) {
          push @g09_basis, $_;
        }
        if (m/ Standard density basis:\s+(\S+)/) {
          $info .= "\t\t\t\tbasis2: $1\n";
        }
        if (m/ Density basis set in the form of general basis input/ .. m/ \*\*\*\*/) {
          push @g09_fit_basis, $_;
        }
      }
      close L;

      if (@g09_basis > 2) {
        shift @g09_basis; pop @g09_basis;
      } else {
        print STDERR "AO basis set $basis_name for $element is wrong\n";
        next;
      }
      if (@g09_fit_basis > 2) {
        shift @g09_fit_basis; pop @g09_fit_basis;
      } else {
        print STDERR "Density basis $basis_name for $element set is wrong\n";
        next;
      }

      unlink "$element.com", "$element.log" unless $no_del;

      $info .= "$element\t ";

      my @b = g2pribas(@g09_basis);
      my $sum = 0;
      foreach (@b) {
        $sum += @$_;
      }
      $basis .= sprintf "%9d%4d\n", $element_number, $sum;
      for (my $i=0; $i<@b; $i++) {
        foreach (@{$b[$i]}) {
          $basis .= sprintf "%d%3d\n", $i, 0+@$_;
          $info .= 0+@$_;
          foreach (@$_) {
            $basis .= join(' ', map {Lprintf($_)} @$_) . "\n";
          }
        }
        $info .= '/';
      }

      $info =~ s/\/$/\n/;
      #pp @b; 
      my @b2 = g2pribas(@g09_fit_basis);
      my $sum2 = 0;
      foreach (@b2) {
        $sum2 += @$_;
      }
      $basis2 .= sprintf "%9d%4d\n", $element_number, $sum2;
      for (my $i=0; $i<@b2; $i++) {
        foreach (@{$b2[$i]}) {
          $basis2 .= sprintf "%d%3d\n", $i, 0+@$_;
          foreach (@$_) {
            $basis2 .= join(' ', map {Lprintf($_)} @$_) . "\n";
          }
        }
      }
  #    print "@g09_fit_basis\n";
      #pp @b2;
    }
    elsif ($gen eq 'orca') {
      open L, '>', "$element.inp" or die "Can't wtite $element.inp: $!\n";
      print L "! PBE $basis_name PrintBasis NoIter MiniPrint\n\n* xyz $charge 1\n $element 0 0 0\n*\n";
      close L;

      open F, '-|', "$orca $element.inp 2>/dev/null" or die "Can't run $orca: $!\n";
      my (@b,@b2,$inf,$inf2);
      my $count = 0;
      my $count2 = 0;
      while (<F>) {
        last if m/One Electron integrals/;
        if (m/There is no basis function on atom/) {
          print "\b"x(length($element)+1);
          unlink glob("$element.*") unless $no_del;
          close F;
          next ELEMENT;
        }
        
        if (m/(\S+) contracted to (\S+) pattern/) {
          if (! $inf) {
            $inf = "{$1}/[$2]";
          } else {
            $inf2 = "{$1}/[$2]";
          }
        }
        
        if (m/^ # Basis set for element/..m/end;/) {
          if (m/^ ([D-S]) (\d+)/) {
            push @b, sprintf "%d%3d\n", $SPD{$1}, $2;
            $count++;
            for (my $i=1; $i<=$2; $i++) {
              my ($alpha,$d) = map {Lprintf($_)} (split ' ', <F>)[1,2];
              push @b, "$alpha $d\n";
            }
          }
        }
        
        if (m/^ # Auxiliary basis set for element/..m/end;/) {
          if (m/^ ([D-S]) (\d+)/) {
            push @b2, sprintf "%d%3d\n", $SPD{$1}, $2;
            $count2++;
            for (my $i=1; $i<=$2; $i++) {
              my ($alpha,$d) = map {Lprintf($_)} (split ' ', <F>)[1,2];
              push @b2, "$alpha $d\n";
            }
          }
        }
      }
      close F;
      unless (@b) {
        print STDERR "No basis\n";
        unlink glob("$element.*") unless $no_del;
        next BASIS if $element eq 'H';
        next ELEMENT;
      }
      unless (@b2) {
        print STDERR "No auxiliary basis\n";
        next ELEMENT;
      }
      unshift @b, sprintf "%9d%4d\n", $element_number, $count;
      unshift @b2, sprintf "%9d%4d\n", $element_number, $count2;
      $basis .= join '', @b;
      $basis2 .= join '', @b2;
      $info .= sprintf(" %-5s%-30s%s\n",$element,$inf,$inf2) if $inf;
      unlink glob("$element.*") unless $no_del;
    }
  }

  my %h;
  $info = join "\n", grep {! $h{$_}++} split "\n", $info if $gen eq 'g09';
  open L, '>', $basis_name or die "Can't open $basis_name: $!\n";
  print L "$info\n\n";
  print L "$basis\n\$end\n";
  print L "$basis2\n\$end\n";
  close L;
  print "\n";
}

sub g2pribas {
  my @pribas;
  while (@_) {
    $_ = shift @_;
    if (m/^\s*[A-Z]/) {
      my ($shells, $contractations, $scale_factor, $unknown) = split;
      if ($scale_factor != 1 or $unknown != 0) {
        print;
        return undef;
      }
      my (@alpha,@d);
      for (my $i=0; $i<$contractations; $i++) {
        ($alpha[$i],@{$d[$i]}) = split ' ', shift @_;
        #print "$alpha[$i],@{$d[$i]}\n";
        if (length($shells) != @{$d[$i]}) {
          print STDERR "$shells but @{$d[$i]}\n";
          return undef;
        }
      }
      #pp @d;
      my $count=0;
      foreach (split '', $shells) {
        my @bas;
        for (my $i=0; $i<$contractations; $i++) {
          $bas[$i] = [$alpha[$i], $d[$i][$count]];
        }
        #pp @bas;
        push @{$pribas[$SPD{$_}]}, \@bas;
        $count++;
      }
    }
  }
  return @pribas;
}

sub Lprintf {
  my $num = shift;
  $num =~ s/d/e/i;
  my ($m,$e) = split 'e', sprintf("%+.12e", $num);
  $m =~ s/\.//;
  $m =~ s/^([+-])/$1./;
  $e = sprintf "%+03d", $e+1;
  return $m.'00000e'.$e;
}