#!/usr/bin/perl -ws

our ($h, $help, $NACP, $all_BCP);

if ($h or $help) {
(my $program = $0) =~ s/^.*[\/\\]//;
print "Usage: $program file.sum
Зависимости: perl

Окучивает для молдена результаты AIMAll (http://aim.tkgristmill.com).
Понимает sum- и mgp-файлы (в последних нет зарядов на атомах).
Создаются xyz-файлы с координатами атомов и критических точек. 
На экран печатается номер крит. точки в xyz-файле и информация о ней.
Обозначения критических точек можно задать свои, см. в начале скрипта, 
например, NACP и NNACP - XX; BCP, RCP и CCP - He. Сейчас все XX.
По умолчанию критические точки на всех ядрах (NACP) и между связями (BCP), 
которые нарисует молден, не включаются. С опциями -NACP и -all_BCP включаются.
Пятой колонкой в xyz-файле идут для атомов заряды q(A) и для критических точек 
элекронная плотность Rho (их можно посмотреть в молдене Label - atom+charge). 
";

exit;
}

# Обозначения критических точек (нужно фикт. атом XX или элемент)
# Элемент лучше брать какой-нибудь неходовой, например He, 
# и, возможно, уменьшить радиус этого элементы в файле .moldenrc
my %CP_name = (
  NACP  => 'XX',   # Nuclear Attractor Critical Point                                      
  NNACP => 'XX',   # Non-Nuclear Attractor Critical Point                                 
  BCP   => 'XX',   # Bond Critical Point                                                    
  RCP   => 'XX',   # Ring Critical Point                                                    
  CCP   => 'XX',   # Cage Critical Point                                                    
);

my $bohr = 0.529177249; # bohr => angstrom
my $num = qr/-?\d+\.\d+E[+-]\d+/; # Number

# Ковалентные радиусы из файла .moldenrc
my %radius = qw(
 H 0.43  Al 1.55  Mn 1.55  Rb 1.67  In 1.83  Pm 2.00  Ta 1.63  At 0.20  Bk 0.20
He 0.741 Si 1.40  Fe 1.54  Sr 1.32  Sn 1.66  Sm 2.00   W 1.57  Rn 0.20  Cf 1.73
Li 0.88   P 1.25  Co 1.53   Y 1.98  Sb 1.66  Eu 2.19  Re 1.55  Fr 0.20
Be 0.55   S 1.22  Ni 1.70  Zr 1.76  Te 1.67  Gd 1.99  Os 1.57  Ra 2.10
 B 1.03  Cl 1.19  Cu 1.72  Nb 1.68   I 1.60  Tb 1.96  Ir 1.52  Ac 2.08
 C 0.90  Ar 0.995 Zn 1.65  Mo 1.67  Xe 1.75  Dy 1.95  Pt 1.70  Th 1.99
 N 0.88   K 1.53  Ga 1.42  Tc 1.55  Cs 1.87  Ho 1.94  Au 1.70  Pa 1.81
 O 0.88  Ca 1.19  Ge 1.37  Ru 1.60  Ba 1.54  Er 1.93  Hg 1.90   U 1.78
 F 0.84  Sc 1.64  As 1.41  Rh 1.65  La 2.07  Tm 1.92  Tl 1.75  Np 1.75
Ne 0.815 Ti 1.67  Se 1.42  Pd 1.70  Ce 2.03  Yb 2.14  Pb 1.74  Pu 0.20
Na 1.17   V 1.53  Br 1.41  Ag 1.79  Pr 2.02  Lu 1.92  Bi 1.74  Am 1.71
Mg 1.30  Cr 1.55  Kr 1.069 Cd 1.89  Nd 2.01  Hf 1.77  Po 1.88  Cm 0.20
);



foreach my $file (@ARGV) {
  (my $file_name = $file) =~ s/\.[^.]*$//;
  open MGP, '<', $file or die "Can't open $file: $!\n";
  my $mol = [{}];
  my $N;
  print "\n$file\n";
  
  # Go
  while (<MGP>) {
    
    # Cartesian Coordinates
    if (m/^Nuclear Charges and Cartesian Coordinates:/) {
      <MGP>; <MGP>; <MGP>;
      while (<MGP>) {
        last if m/^\s*$/; # Termination by empty string
        my (undef,$Charge,$X,$Y,$Z) = split;
        push @$mol, [element(int $Charge),$X*$bohr,$Y*$bohr,$Z*$bohr];
      }
      $N = $#{$mol}; 
    }
    
    # Atomic Properties (in sum file only)
    if (m/^Some Atomic Properties:/) {
      while (<MGP>) {
        # Net Charge of Atoms
        if (m/^Atom A          q\(A\)/) {
          <MGP>;
          for (my $i=1; $i<=$N; $i++) {
            my $q_A = (split ' ', <MGP>)[1];
            $mol->[$i][4] = $q_A if $q_A =~ /^$num$/;
          }
        }
        # Delete bad charges
        if (m/^Warning!  Significant integration error for the following atom:\s+[A-Za-z]+(\d+)/) {
          undef $mol->[$1][4];
        }
        # Energy
        if (m/^Molecular energy E\(Mol\) from the wfn file:\s+($num)/) {
          $mol->[0]{Energy} = $1;
        }
        last if m/Electron Density Critical Point Analysis of Molecular Structure/; # Termination
      }
    }
    
    # Critical Points
    if (m/^CP#\s+\d+\s+Coords =\s*/) {
      my ($X,$Y,$Z) = split ' ', $';
      my ($Type,$Rho,$Type_str);
      while (<MGP>) {
        last if m/^\s*$/;
        if (m/^\s+Type = \(\S+\)\s+([A-Z]+)/ ) {
          $Type = $1;
          $Type_str = $_;
        }
        if (m/^\s+Rho =\s*(\S+)/) {
          $Rho = $1;
        }
      }
      
      my $need_CP = 1;
      $need_CP = 0 if 
        ( $Type eq 'NACP' && 
          ! $NACP
        ) ||
        ( $Type eq 'BCP' && 
          ! $all_BCP &&
          $Type_str =~ /[A-Za-z]+(\d+)\s+[A-Za-z]+(\d+)\s*$/ &&
          distance($mol,$1,$2) < $radius{$mol->[$1][0]}+$radius{$mol->[$2][0]}
        );
      
      if ($need_CP) {
        push @$mol, [$CP_name{$Type},$X*$bohr,$Y*$bohr,$Z*$bohr,$Rho];
        print $#{$mol}, $Type_str;
      }
    }
  }
  close MGP;
  
  write_molden($mol, "$file_name.xyz");
}

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 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";
    printf F " Energy %.8f", $mol->[0]{Energy} if $mol->[0]{Energy};
    print F "\n";
    for (my $i=1; $i<=$N; $i++) {
      printf F " %-2s %12.8f %12.8f %12.8f", @{$mol->[$i]}[0..3];
      printf F " %8.4f", $mol->[$i][4] if defined $mol->[$i][4];
      print F "\n";
    }
  }
  close F;
}

sub distance {
  my ($mol,$n1,$n2) = @_;
  my $dist = ($mol->[$n1][1] - $mol->[$n2][1])**2 +
             ($mol->[$n1][2] - $mol->[$n2][2])**2 +
             ($mol->[$n1][3] - $mol->[$n2][3])**2;
  return sqrt($dist);
}