#!/usr/bin/perl -ws

# Полное имя программы symmetry. Симметрия молекулы будет автоматически 
# определена и правильно задана в группе $DATA.
my $symmetry_x = `which symmetry`;

# Полное имя программы babel. Для задания молекулы z-матрицей
my $babel_x = `which babel`;

# Полное имя файла Природовских базисов basis.in
# Для -basis=priSET
my $basis_in = '/usr/local/lib/priroda/basis/basis.in';

#use Data::Dump 'pp';

# Переменные опций
our ($h,$help,$run,$dft,$mp,$basis,$charge,$mult,$dlc,$zmt,$nosymm);

if ($h || $help) {
  (my $program = $0) =~ s/^.*[\/\\]//;
  print <<HELP;
Prepare input files for GAMESS from xyz geometry. 

Usage: $program [-option[=string]] file.xyz [file1.xyz ...]
Зависимости: perl, [symmetry, babel].


-run=str  Default optimize. For -run=SADPOINT sets hess=calc in \$statpt.
          For run=IRC puts \$irc group with evaluated IRC step.
-dft=str  Functional for DFT. -dft is -dft=B3LYP
-mp=N     MPLEVL. -mp is -mp=2
-basis='str'  Default 6-31G*. Pople's N21/N31/N311 is parsed.
              F.e. -basis='6-311++G**'
              If -basis=priSET, then it is Priroda basis SET from 
              basis.in (after each atom in\$DATA group).
              f.e. -basis=priL1 
-mult=N    Multiplicity. scftyp=UHF for -mult=2 
-charge=N  Charge. If dosn't specified, then calculated (0 or 1).
-zmt       Transform xyz to z-matrice. Babel needs. Force -nodlc. 
-nodlc     Force no DLC. Default for -zmt.
-nosymm    Force C1 symmetry.

HELP
  exit;
}

$run ||= 'optimize';
$dft = 'B3LYP' if defined($dft) && $dft=~/^1$/;
$mp ||= 0;
$mp = 2 if $mp == 1;
$basis ||= '6-31G*';
$mult ||= 1;

chomp $babel_x;
if (-x $babel_x) {
  $nodlc = 1 if $zmt;
} else {
  warn "No babel for xyz to zmat transformation. Force coord=unique.\n" if $zmt;
  $zmt = 0;
}

chomp $symmetry_x;
unless (-x $symmetry_x) {
  warn "No symmetry program. Force -nosymm=1.\n" unless $nosymm;
  $nosymm = 1;
}

while (@ARGV) {
  my $file = shift;
  my $mol = (read_molden($file))[-1];
  next if @$mol < 2;
  
  xyz2zmt($mol) if $zmt;
  
  symmetr($mol) unless $nosymm;
  my $symm = $mol->[0]{Symmetry} || 'C1';
  
  #pp $mol; exit;
  (my $name = $file) =~ s/\.xyz(?:ppm)?$//;
  my $charge = calc_charge($mol, $mult) unless $charge;
  my $coord = $zmt ? 'zmt' : 'unique';
  #pp $mol;
  my $out_file = write_gamess_input(file     => $name,
                                    mol      => $mol,
                                    run      => $run, 
                                    mult     => $mult,
                                    symmerty => $symm,   
                                    charge   => $charge,
                                    basis    => $basis,
                                    dft      => $dft,
                                    mp       => $mp,
                                    dlc      => $nodlc ? 0 : 1,
                                    coord    => $coord,
                                    );
  
  if ($basis && $basis =~ /^pri(.+)/) {
    my $set = $1;
    my %pri_bas = pribas($mol,$basis_in,$set);
    #pp %pri_bas;
    $out_file = include_pribas($out_file, $set, %pri_bas);
  }
  print "$out_file\n";
}

sub write_gamess_input {
  my %par = (run      => 'optimize',
             charge   => 0,
             mult     => 1,
             symmerty => 'C1',
             basis    => '6-31G*',
             dft      => '',
             mp       => 0,
             dlc      => 1,
             file     => '',
             coord    => 'unique',
             @_);
  
#  print "$par{dlc}\n";
  $par{run} = uc $par{run};

  my $symm = $par{symmerty};
  if ($symm ne 'C1') {
    $symm =~ s/^([CD])(\d)([vhd]?)$/${1}n$3 $2/;
    $symm =~ s/^S(\d)$/'S2n '.int($1\/2)/e;
    $symm .= "\n";
  }
  
  my $file = $par{file};
  $file .= '.'.uc($par{dft}) if $par{dft};
  $file .= ".MP$par{mp}" if $par{mp};
  
  my $basis = $par{basis};
  #print "$basis\n";
  if ($basis =~ /^(\d)G?-?(\d+)G?(\+*)G?(\**)$/) {
    $basis = "gbasis=n$2 ngauss=$1 ";
    $basis .= 'ndfunc=1 '   if length($4) >= 1;
    $basis .= 'npfunc=1 '   if length($4) >= 2;
    $basis .= 'diffsp=.t. ' if length($3) >= 1;
    $basis .= 'diffs=.t. '  if length($3) >= 2;
    $basis =~ s/\s+$//;
    
    my $to_name = $par{basis};
    $to_name =~ s/[G\-]//g;
    $to_name =~ s/\+/p/g;
    $to_name =~ s/\*/x/g;
    $to_name =~ s/\W/_/g;
    $file .= ".$to_name";
  }
  if ($basis =~ /^pri(.+)/) {
    $file .= ".$1";
  }
  
  my $scftyp = $par{mult}%2 ? 'RHF' : 'UHF';
  my $hess = $par{run} eq 'SADPOINT' ? 'calc' : 'guess';
  
  my $mol = $par{mol};
  my $nzvar = ($par{dlc} || $zmt) ? 3*@$mol-9 : 0;
  
  $file .= '.inp';
  my $oldfh;
  if ($par{file}) {
    #warn "$file\n";
    open F, '>', $file or do {warn "Can't write $file: $!\n"; return 0};
    $oldfh = select F;
  }
  
  my $irc_step = IRC_step($mol) if $par{run} eq 'IRC';
  
  print " \$contrl coord=$par{coord} icharg=$par{charge} mult=$par{mult} scftyp=$scftyp maxit=100\n";
  print "  nzvar=$nzvar ispher=+1 mplevl=$par{mp} runtyp=$run exetyp=run \$end\n";  
  print " \$system timlim=3000 mwords=32 memddi=8 \$end\n";
  print " \$statpt hess=$hess nstep=300 \$end\n";
  print " \$basis  $basis \$end\n";
  print " \$dft dfttyp=$par{dft} \$end\n" if $par{dft};
  print " \$zmat dlc=.t. auto=.t. \$end\n" if $par{dlc};
  print " \$irc pace=gs2 stride=$irc_step npoint=100 saddle=.t. forwrd=.f. nprt=-2 npun=-2 \$end\n" if $par{run} eq 'IRC';
  print " \$scf dirscf=.t. \$end\n";
  print " \$guess guess=huckel \$end\n";
  print " \$data\n";
  print " \n";
  print " $symm\n";
  
  if (exists $mol->[0]{Zmt}) {
    print @{$mol->[0]{Zmt}};
  }
  else {
    if (exists($mol->[0]{Unique}) && ($mol->[0]{Symmetry} ne 'C1') ) {
      foreach (@{$mol->[0]{Unique}}) {
        my ($n,$x,$y,$z) = split;
        printf " %-2s%6.1f%8.3f%8.3f%8.3f\n", 
                element($n),$n,$x,$y,$z; 
      }
    }
    else {
    #pp $mol;
    #print @{$mol->[0]{Zmt}};
      for (my $i=1; $i<@$mol; $i++) {
        printf " %-2s%6.1f%13.8f%13.8f%13.8f\n", 
         $mol->[$i][0], element($mol->[$i][0]), 
         $mol->[$i][1], $mol->[$i][2], $mol->[$i][3]; 
      }
    }
  }
  
  print " \$end\n";
  
  do {select $oldfh; close F} if $par{file};
  return $file;
}


sub calc_charge {
 ############################################################################
 ## Принимает молекулу и мультиплетность
 ## Возвращает заряд (0 или 1)
 ############################################################################
  my ($mol, $mult) = @_;
  my $sum = $mult-1;
  for (my $i=1; $i<@$mol; $i++) {
    $sum += element($mol->[$i][0]);
  }
  return $sum%2;
}


sub read_molden {
 ############################################################################
 ## Молекула -- ссылка на массив, 0-й элемент -- свойства, следующие - атомы.
 ## Свойства -- ссылка на хэш с ключами Energy, Symmetry
 ## Атом -- ссылка на массив [atom, x, y, z, ppm] (ppm может не быть).
 ##
 ## Читает xyz. Параметры - имена xyz-файлов. Если параметров нет, то <>.
 ## Возвращает массив найденных молекул.
 ############################################################################
  local @ARGV = @_ ? @_ : @ARGV;
  my $num = qr/-?\d+(?:\.\d+)?/;
  my @mols;
  my $line;
  LOOP:
  while ($line || defined($line = <>)) {
    #print $line;
    if ($line =~/^\s*(\d+)\s*$/) {
      my @mol;
      my $N = $1;
      last LOOP if eof();
      next LOOP if eof(ARGV);
      $line = <>;
      ($mol[0]{Energy}) = $line =~ /($num)/o;
      ($mol[0]{Symmetry}) = $line =~ /symm\S*\s+(\S+)/i;
      for (my $i=1; $i<=$N; $i++) {
        last LOOP if eof();
        next LOOP if eof(ARGV);
        $line = <>;
        #print $line;
        if ($line =~ /^\s*([A-Z]{1,2})\s+($num)\s+($num)\s+($num)\s*(.*)/io) {
          $mol[$i] = [$1,$2,$3,$4,$5];
        } else {
          next LOOP;
        }
      }
      push @mols, \@mol;
      last LOOP if eof();
    } else {
      undef $line;
    }
  }
  return @mols;
}

sub element {
 ############################################################################
 ## Возвращает атомный номер по символу либо символ по номеру
 ## или undef при отсутствии элемента в таблице Менделеева
 ############################################################################
  # Таблица Менделеева
  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 $atom = shift;
  if ($atom=~/^\d+$/) {
    return undef if $atom > $#ATOM;
    return $ATOM[$atom];
  }
  $atom = ucfirst $atom;
  my %ATOM;
  @ATOM{@ATOM} = 0..$#ATOM;
  return $ATOM{$atom} if exists $ATOM{$atom};
  return undef;
}

sub sum {
  my $sum = 0;
  foreach (@_) {
    $sum += $_;
  }
  return $sum;
}

sub write_molden {
  my $file = pop @_;
  open F, '>', $file or do {warn "Can't write to $file: $!\n"; return};
  foreach my $mol (@_) {;
    my $N = $#{$mol};
    print F " $N\n";
    print F " Energy $mol->[0]{Energy}" if $mol->[0]{Energy};
    print F "  Point $mol->[0]{Point}" if $mol->[0]{Point};
    print F "  Symmetry $mol->[0]{Symmetry}" if $mol->[0]{Symmetry};
    print F "  Ellips $mol->[0]{Ellips}" if $mol->[0]{Ellips};
    print F "\n";
    for (my $i=1; $i<=$N; $i++) {
      printf F " %-2s %12.8f %12.8f %12.8f\n", @{$mol->[$i]};
    }
  }
  close F;
}

sub xyz2zmt {
 ############################################################################
 ## Принимает ссылку на молекулу.
 ## В свойства дописывает элемент хэша Zmt => 'строки z-матрицы'.
 ## Требуется babel.
 ############################################################################

  my $mol = shift;
  write_molden($mol, 'tmp_TMP_tmp.xyz');
  my @output = `$babel_x -ixyz tmp_TMP_tmp.xyz -ofh 2>/dev/null`;
  unlink 'tmp_TMP_tmp.xyz';
  shift @output;
  do {warn "Bad babel output\n"; return} if @output != @$mol;
  shift @output;
  $_ = " $_" foreach @output;
  $mol->[0]{Zmt} = \@output;
}

sub IRC_step {
 ############################################################################
 ## Принимает ссылку на молекулу. Возвращает шаг IRC в bohr.(amu)^0.5
 ############################################################################

  # Массы изотопов
  my %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
  );

  my @atoms = @{$_[0]};
  shift @atoms;
  my $N = @atoms;
  #print "$N atoms\n";

  my %formula;
  $formula{ucfirst lc $_->[0]}++ foreach @atoms;

  my $M;
  $M += $massa{$_}*$formula{$_} foreach keys %formula;
  #print "MW $M\n";
  
  # Примем за "нормальный" шаг 0.5 ангстр.(amu)^0.5 для N=40 M=240
  # Будем нормировать шаг на (N*M)^0.5
  my $norm_step = 0.5;
  my $norm_N = 40;
  my $norm_M = 240;
  
  my $step = $norm_step*sqrt($M*$N/$norm_N/$norm_M);
  $step /= 0.53;  # angstrem to bohr
  #print "$step\n";
  # Округляем до 1 значащей цифры
  $step = 0 + sprintf("%.0e", $step);
  return $step;
}

sub symmetr {
 ############################################################################
 ## Принимает ссылку на молекулу.
 ## В свойства дописывает элементы хэша Symmetry => 'группа симметрии'
 ## и Unique => массив строк с уникальными координатами для GAMESS.
 ## Требуется symmetry.
 ############################################################################
  
   my $mol = shift;
  my $N = @$mol - 1;
  
  open TINP, ">", "symm.TINP" or warn "Can't write temporary file: $!\n",
                                 return undef;
  print TINP " $N\n";
  for (my $i=1; $i<=$N; $i++) {
    print TINP element($mol->[$i][0]), 
               "  $mol->[$i][1]  $mol->[$i][2]  $mol->[$i][3]\n";
  }
  close TINP;
  
  my $prim = 0.05;
  my $final = 0.2;
  open OUT, '-|', "$symmetry_x -primary $prim -final $final symm.TINP" or
                  warn "Can't run $symmetry_x: $!\n", return undef;
  while (<OUT>) {
    if (m/^It seems to be the (\w+) point group$/) {
      $mol->[0]{Symmetry} = $1;
    }
    if (m/^Coordinates of abelian symmetry unique atoms/) {
      <OUT>; <OUT>;
      while (<OUT>) {
        push @{$mol->[0]{Unique}}, $_;
      }
    }
  }
  close OUT;
  unlink 'symm.TINP';
}

sub pribas {
 ############################################################################
 ## Принимает ссылку на молекулу, файл базиса Природы, set.
 ## Для каждого атома пятым элементом дописывает базис в формате Гамесс
 ############################################################################
  
   my ($mol,$basfile,$set) = @_;
  
  # Hash of unique elements
  my %elements;
  for (my $i=1; $i<@$mol; $i++) {
    $elements{$mol->[$i][0]} = [];
  }
  
  open F, '<', $basfile or die "Can't open $basfile: $!\n";
  my $is_set = '';
  my $atom = '';
  while (<F>) {
    last if m/\Send/;
    if (m/set=(.*)/) {
      $is_set = grep {/^$set$/} split /\|/, $1;
    }
    if (m/\s+(\d+)\s+\d+/) {
      $atom = element($1);
    }
    if ($is_set && exists($elements{$atom}) && m/^[+\d]/) {
      push @{$elements{$atom}}, $_;
    }
  }
  close F;
  
  my @s = map {uc} qw(s p d f g h i j k l m n);
  foreach my $atom (keys %elements) {
    #print "$atom\n";
    my $str = '';
    my $contr;
    my $count;
    foreach (@{$elements{$atom}}) {
      if (m/^(\d+)\s+(\d+)/) {
        $str .= sprintf "%4s%8d\n", $s[$1],$2;
        $contr = $2;
        $count = 1;
        next;
      }
      #s/([+-])\./${1}0./g;
      $str .= sprintf "%6d%26.10f%14.10f\n", $count++, split;
    }
    
    if ($str) {
      $str .= "\n";
    } else {
      warn "No basis $set for $atom in $basfile\n";
    }
    $elements{$atom} = $str;
    #print $str;
  }
#  for (my $i=1; $i<@$mol; $i++) {
#    my $str = $elements{$mol->[$i][0]};
#    unless ($str) {
#      warn "No basis $set for $mol->[$i][0] in $basfile\n";
#      $mol->[$i][5] = '';
#      next;
#    }
#    $mol->[$i][5] = "$str\n" if $str;
#  }
  return %elements;
}

sub include_pribas {
  my $file = shift;
  my $set = shift;
  %pri_bas = @_;
#  (my $out = $file) =~ s/[^.]+\.inp/$set.inp/;
#  $out = "$file.pri" if $out eq $file;
  
  open O, '>', "$file.TMP" or die "Can't write to $file.TMP: $!\n";
  open F, '<', $file or die "Can't open $file: $!\n";
  while (<F>) {
    if (m/\$basis/i) {
      print O "  ! Priroda basis $set in \$data group\n";
      next;
    }
    print O;
    if (m/\$data/i) {
      $_ = <F>; print O;
      while (<F>) {
        print O;
        if (m/^\s*([A-Z][a-z]?)\s+/) {
          #print "$1\n";
          print O $pri_bas{$1} if $pri_bas{$1};
        }
        last if m/\$end/i;
      }
    }
  }
  close F;
  close O;
  rename "$file.TMP", $file;
  return $file;
}