#!/usr/bin/perl -ws

unless (@ARGV) {
  print "
 Usage: priin [-gjc|-mop|-xyz|-mrv] file [files]
 
Simple converter into Priroda input files (Perl)
From: Gaussian Input (*.GJ[CF]), created by Chem3D (-gjc option)
      Mopac *.MOP or *.dat or *.arc (-mop option)
      Molden *.xyz or *.xyzppm (-xyz option)
      MarvinSketch *.mrv (-mrv option)
To:   Input (*.in) for PRIRODA
    Если задана опция, то для всех файлов будет предполагаться
    заданный опцией формат. В противном случае формат будет 
    определяться  из расширения (case-insensitive).

-qm  theory=qm_n3  basis=qm.in

-md  task=thermostat + \$dynamic group
     Опции -md_e -md_d -md_v -md_h -md_n -md_m -md_i управляют соотв. параметрами
     Формат -md_i=1-9 делает несколько инпутов с разными i
-vol Использовать сферическую полость. 
     -md_v рассчитывается по плотности \"жидкости\" (типа -vol=0.89)

Если в мопаковских файлах некоторые геометрические параметры заморожены 
(0 а не 1 после значения параметра), то в in-файле будут соответствующие
fix=...

Если в исходном xyz-файле присутствуют фиктивные атомы XX (см. addXX),
то будет task=nmr и группа \$nmr для расчета экранирования в положениях
фиктивных атомов.\n";
  exit;
}

our($mem,$disk,$basis,$rel,$task,$theory,$H,$qm);

$mem ||= 2000;
$disk ||= -4000;
$basis ||= 'L1';
$rel ||= '';
$basis = 'qm.in' if $qm;
$basis = '3z' if $z;
$task ||= 'optimize';
$task = 'thermostat' if $md;
$task = 'hessian' if $H;
$theory ||= 'dft';
$theory = 'qm_n3' if $qm;

$md_e ||= 0.001;
$md_d ||= 0.03;
if ($vol || $md_v || $md_h) {
  $md_v ||= 500;
  $md_h ||= 0.1;
  $vol  ||= 1;
}
$md_n ||= '1000,100';
$md_m ||= 10;
$md_i ||= 1; # -md_i=1-10 will be 10 in-files
if ($md_T) {
  $md_e = sprintf "%.6f", $md_T/315775;
}

my $num = qr/-?\d+\.\d+/;

# Вгоняем таблицу Менделеева
@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
           );
for ($i=0; $i<@ATOM; $i++) { $ATOM{$ATOM[$i]} = $i };
#foreach $k (keys %ATOM) {print "$k\t$ATOM{$k}\n"};

our %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
);

#@MW = qw (0 1.00794 4.002602 6.941 9.012182 10.811 12.0107 14.0067 15.9994 
#18.9984032 20.1797 22.989770 24.3050 26.981538 28.0855 30.973761 32.065 
#35.453 39.948 39.0983 40.078 44.955910 47.867 50.9415 51.9961 54.938049 55.845 
#58.933200 58.6934 63.546 65.409 69.723 72.64 74.92160 78.96 79.904 83.798 
#85.4678 87.62 88.90585 91.224 92.90638 95.94 98 101.07 102.90550 106.42 
#107.8682 112.411 114.818 118.710 121.760 127.60 126.90447 131.293 132.90545 
#137.327 138.9055 140.116 140.90765 144.24 145 150.36 151.964 157.25 158.92534 
#162.500 164.93032 167.259 168.93421 173.04 174.967 178.49 180.9479 183.84 
#186.207 190.23 192.217 195.078 196.96655 200.59 204.3833 207.2 208.98038 209 
#210 222 223 226 227 232.0381 231.03588 238.02891 237 244 243 247 247 251 252 
#257 258 259 262 261 262 266 264 277 268);

if ($xyz || $mop || $gjc || $mrv) {
  $predefined_inp_format = 1;
  if    ($gjc) {$inp_format = 'gjc'}
  elsif ($mop) {$inp_format = 'mop'}
  elsif ($xyz) {$inp_format = 'xyz'}
  elsif ($mrv) {$inp_format = 'mrv'}
  else {print "Unknown option\n"; exit};
}

foreach ("/disk2/tmp/$ENV{USER}", '/disk2/tmp', '/tmp') {
  -w $_ or next;
  $path = $_;
  last;
}
$path ||= '.';

foreach (@ARGV) {
  $inp = $_;
  next if -d $inp;
  (my $name = $inp) =~ s/\.(\w+$)//i;
  unless ($predefined_inp_format) {
    if    ($1 =~ /gj[cf]/i) {$inp_format = 'gjc'}
    elsif ($1 =~ /mop|dat|arc/i) {$inp_format = 'mop'}
    elsif ($1 =~ /xyz|xyzppm/i) {$inp_format = 'xyz'}
    elsif ($1 =~ /mrv/i) {$inp_format = 'mrv'}
    else {print "Unknown extension .$1 ($name.$1)\n"; next};
  }

  open (INP, $inp ) || die "Cannot open input file $inp: $!\n";
  
  undef $charge; undef $mult; $sigma = 1;
  undef @atom; undef @x; undef @y; undef @z;
  undef $heat; undef $method;
  undef @f; undef $fix;
  my (@nXX,@points);
  
  
  if ($inp_format eq 'xyz') {
    1 while &read_molden(\*INP);
  }
  elsif ($inp_format eq 'gjc') {
    $n = 0;
    while (<INP>) {
      if (m/^\s*(\d+)\s+(\d+)\s*$/) {
        ($charge, $mult) = ($1,$2);
        last;
      }
    }
    while (<INP>) {
      if (m/^\s*(\w\w?)\s+\d\s+($num)\s+($num)\s+($num)\s*$/) {
        ($atom[$n], $x[$n], $y[$n], $z[$n]) = ($1,$2,$3,$4);
        $atom[$n] = ucfirst lc $atom[$n];
        if (! exists $ATOM{$atom[$n]}) {
          warn "Do not exists element $ATOM{$atom[$n]}\n";
          last;
        }
        $n++;
      }
      else {
        last;
      }
    }
  }
  elsif ($inp_format eq 'mop') {
    while (<INP>) {
      if (/HEAT +OF +FORMATION += +(-?\d+\.\d+)/i) {$heat = $1};
      if (/CHARGE=([+\-]?\d+)/i) {$charge = $1};
      if (/MULT=([+\-]?\d+)/i) {$mult = $1};
      if (/(PM3|AM1|MNDO|MINDO\S*)/i) {$method = $1};
      if (/([a-z]{1,2})( +0\.0+ +0){3}/i) {
        $f[0][0]=$f[0][1]=$f[0][2]=$f[1][1]=$f[1][2]=$f[2][2]=1;
        $atom[0] = $1;
        $_=<INP>; 
        ($atom[1], $x[1], $f[1][0], undef) = split;
        $k[1] = 1;
        $_=<INP>; 
        ($atom[2],$x[2],$f[2][0],$y[2],$f[2][1],undef,undef,$k[2],$l[2],undef) = split;
        $i = 2;
        while (<INP>) {
          last if (/^\s*$/);
          $i++;
          #    print "$_";
                ($atom[$i],$x[$i],$f[$i][0],$y[$i],$f[$i][1],$z[$i],$f[$i][2],$k[$i],$l[$i],$m[$i],undef) = split;
        }
        $n = $i;
        #  print "n=$n\n";
      }
    }
    foreach $i (0..$n) {
      foreach $j (0..2) {
        if ($f[$i][$j]==0) {
          $ii = $i+1;
          if    ($j==0) {$fix .= "  1,$ii,$k[$i],0,0"}
          elsif ($j==1) {$fix .= "  2,$ii,$k[$i],$l[$i],0"}
          elsif ($j==2) {$fix .= "  3,$ii,$k[$i],$l[$i],$m[$i]"}
        }
      }
    }
#  print "$heat\n$charge\n$mult\n$method\n";
  }
  elsif ($inp_format eq 'mrv') {
    my ($atom,$x,$y,$z);
    while (<INP>) {
      $atom = $1 if m/\selementType="(.+?)"/;  
      $x = $1    if m/\sx3="(.+?)"/;           
      $y = $1    if m/\sy3="(.+?)"/;           
      $z = $1    if m/\sz3="(.+?)"/;           
    }
    print "$atom\n";
    @atom = map { ucfirst } split ' ', $atom;
    @x = split ' ', $x;
    @y = split ' ', $y;
    @z = split ' ', $z;
    $N = @atom;
    if (!$N || @x!=$N || @y!=$N || @z!=$N) {
      warn "Error in $inp\n";
      next;
    }
  }
  close INP;
  
#  if ($heat) {
#    print IN " $method" if ($method);
#    print IN " HEAT OF FORMATION = $heat KCAL/MOLE\n";
#  }
  
  undef $SUM;
  foreach (@atom) {
    next unless defined $_;
    next if m/xx/i;
    $_ = ucfirst(lc($_));
    $SUM += $ATOM{$_};
    $rel = 'r' if ($basis=~/^L/ && $ATOM{$_}>36);
    
#    print "$_ $ATOM{$_}\n";
  }
  unless (defined $charge) {
    $charge = ($SUM%2) ? 1 : 0;
  }
#  print "$charge\n"; exit;

  if ($inp_format eq 'xyz') {
    #print "inp_format eq 'xyz'\n";
    for ($i=1; $i<=$N; $i++) {
      unshift @nXX, $i if $atom[$i] =~ /^xx?$/i;
    }
    #print "@nXX\n";
    foreach (@nXX) {
#      delete $atom[$_];
      splice @atom, $_, 1;
      unshift @points, sprintf("%14.8f%14.8f%14.8f\n",splice(@x,$_,1),splice(@y,$_,1),splice(@z,$_,1));
      $N--;
    }
  }
  $task = 'nmr' if @nXX;
  
  $in = "$name.$rel$basis.in";
  $in = "$name.qm.in" if $qm;
  print "Input file $inp  Format $inp_format  Output file $in\n";
  open (IN, ">:unix", $in) || die "Cannot write in output file $in: $!\n";

  # $system
  print IN "\$system mem=$mem disk=$disk";
  print IN " path=$path" if $disk > 0;
  print IN " \$end\n";
  
  # $control
  print IN "\$control
 task=$task
 theory=$theory
 basis=$rel$basis\n";
  #print IN " print=+charges\n" if $task eq 'optimize';
  print IN "\$end\n";
  
  # $grid
  print IN "\$grid accur=1e-8 \$end\n" if ! $md;
  
  # $optimize
  if (! $md && $task eq 'optimize') {
    print IN "\$optimize\n";
      if ($fix) {
        $fix =~ s/^  //;
        print IN " fix=$fix\n";
      }
      print IN " steps=300 tol=1e-5\n\$end\n";
  }
  
  # $dynamic
  if ($md) {
    print IN "\$dynamic\n e=$md_e d=$md_d\n";
    if ($vol) {
      my $mw = 0;
      for (my $i=1; $i<=$N; $i++) {
        $mw += $massa{$atom[$i]};
      }
      $md_v = sprintf "%d", ($mw/6.02e23)/$vol*1e24/0.52917721**3/$N;
      print IN " v=$md_v h=$md_h\n";
    }
    print IN " n=$md_n m=$md_m\n";
    $md_i =~ /^(\d+)/ && print IN " i=$1\n";
    print IN "\$end\n";
  }
  
  # hessian
  if ($task eq 'hessian') {
    print IN "\$thermo sigma=$sigma \$end\n" if $sigma;
  }
  # $nmr
  if (@nXX) {
    print IN "\$nmr\n points=\n";
    print IN @points;
    print IN "\$end\n";
  }
  
  # $molecule
  print IN "\$molecule\n";
  print IN " charge=$charge";
  (defined $mult) ? (print IN " mult=$mult\n") : (print IN " mult=1\n");
  if ($inp_format eq 'xyz') {
    print IN " cartesian\n";
    for ($i=1; $i<=$N; $i++) {
      printf IN "%4d%14.8f%14.8f%14.8f\n",$ATOM{$atom[$i]}, $x[$i], $y[$i], $z[$i] if $atom[$i] ne 'XX';
    }
  }
  if ($inp_format eq 'mrv') {
    print IN " cartesian\n";
    for ($i=0; $i<$N; $i++) {
      printf IN "%4d%14.8f%14.8f%14.8f\n",$ATOM{$atom[$i]}, $x[$i], $y[$i], $z[$i] if $atom[$i] ne 'XX';
    }
  }
  if ($inp_format eq 'gjc') {
    print IN " cartesian\n";
    for ($i=0; $i<$n; $i++) {
      printf IN "%4d%12.6f%12.6f%12.6f\n",$ATOM{$atom[$i]}, $x[$i], $y[$i], $z[$i];
    }
  }
  if ($inp_format eq 'mop') {
    print IN " z-matrix\n";
    printf IN "%4d\n",$ATOM{$atom[0]};
    printf IN "%4d%6d%10.6f\n",$ATOM{$atom[1]},$k[1],$x[1];
    printf IN "%4d%6d%10.6f%6d%10.4f\n",$ATOM{$atom[2]},$k[2],$x[2],$l[2],$y[2];
    for ($i=3; $i<=$n; $i++) {
      printf IN "%4d%6d%10.6f%6d%10.4f%6d%10.4f\n",$ATOM{$atom[$i]},$k[$i],$x[$i],$l[$i],$y[$i],$m[$i],$z[$i];
    }
  }
  print IN "\$end\n";  
  close IN;
  
  if ($md_i && $md_i=~/(\d+)-(\d+)/) {
    my ($i1,$i2) = ($1,$2);
    open L, '<', $in or die "Can't open $in: $!\n";
    my $in_content = join '', <L>;
    close L;
    foreach my $i ($i1..$i2) {
      my $format = '%0'.length($i2).'d';
      my $in_new;
      if ($in =~ /^([^.]+)\.(.+)$/) {
        $in_new = $1 . '-' . sprintf($format,$i) . '.' . $2;
      }
      else {
        $in_new = $in . '-' . sprintf($format,$i);
      } 
      $in_content =~ s/ i=\d+/ i=$i/;
      open L, '>', $in_new or die "Can't write to $in_new: $!\n";
      print L $in_content;
      close L;
    }
    unlink $in;
  }
}

sub read_molden {
  my $inp = shift;
  $inp = 'STDIN' unless $inp;
  return undef if eof($inp);
  # 1-st string
  if (<$inp>=~/^\s*(\d+)\s*$/) {
    ($N,$energy,@atom,@x,@y,@z,@charge) = undef;
    $N = $1;
  } else {
    return undef;
  }
  # 2-nd string
  my $line2 = <$inp>;
  ($energy) = $line2=~/(-?\d+\.\d+)/;
  $line2=~/Charge (-?\d+)/ and $charge = $1;
  $line2=~/Mult (\d+)/ and $mult = $1;
  $mult = 1 if (defined($charge) && ! $mult);
  $charge = 0 if (! defined($charge) && $mult);
  $sigma = $1 if $line2=~/sigma (\d+)/;
  # Geometry
  for ($i=1;$i<=$N;$i++) {
    ($atom[$i],$x[$i],$y[$i],$z[$i],$charge[$i],undef) = split ' ', <$inp>;
    return undef unless $atom[$i]=~/^\w\w?/ &&
                        $x[$i]=~/^-?\d+\.\d+/ &&
                        $y[$i]=~/^-?\d+\.\d+/ &&
                        $z[$i]=~/^-?\d+\.\d+/;
  }
  return $N;
}