#!/usr/bin/perl -w

if (! @ARGV or $ARGV[0]=~/^--?h(elp)?/) {
  print "\n Usage: symm file.xyz
Зависимости: perl, symmetry\n
Проводит симметризацию геометрии с помощью программы symmetry, 
постепенно увеличивая параметр -final последней, вплоть до 0.2
(и подстраивая параметр -primary).

Читает последнюю геометрию из file.xyz, если удалось обнаружить
симметрию, то результат в формате xyz печатается на stdout.
\n";
  exit;
}

$symmetry = '/usr/local/bin/symmetry';
die "Program symmetry not found\n" unless (-x $symmetry);

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

foreach (@ARGV) {
  open INP, '<', $_ || warn "Can't open $_: $!\n" , next;
  while (read_molden(\*INP)) {};
  close INP;
  open TINP, ">", "symm.TINP" || warn "Can't write temporary file: $!\n" , next;
  print TINP " $N\n";
  for ($i=1; $i<=$N; $i++) {
    print TINP " $ATOM{$atom[$i]}   $x[$i]  $y[$i]  $z[$i]\n";
  }
  close TINP;

  $group = 'C1';
  $prim = 0.05;
  foreach my $ext (reverse (0.7, 1..7)) {
    $final = sprintf '%.8f', 10**(-$ext) - 1e-8;
    $prim = $final if $prim < $final;

#    $final = 10**(-$_);
    open OUT, '-|', "$symmetry -primary $prim -final $final symm.TINP";
    while (<OUT>) {
      if (m/^It seems to be the (\w*) point group$/) {
        if ($1 ne $group) {
          $group = $1;
          while (<OUT>) {
            if (m/^Symmetrized coordinates of all atoms/) {
              <OUT>; <OUT>;
              for ($i=1; $i<=$N; $i++) {
                (undef,$x[$i],$y[$i],$z[$i]) = split ' ', <OUT>;
              }
              #$energy .= "   Group $group";
              &write_molden(-without_charges);
              last;
            }
          }
        }
        last;
      }
    }
    close OUT;
  }

  unlink "symm.TINP";
#  print @output;
}

#symmetry [option value ...] [filename]
#Valid options are:
#  -verbose      (  0) Determines verbosity level
#                      All values above 0 are intended for debugging purposes
#  -maxaxisorder ( 20) Maximum order of rotation axis to look for
#  -maxoptcycles (200) Maximum allowed number of cycles in symmetry element optimization
#  --                  Terminates option processing
#Defaults should be Ok for these:
#  -same         (   0.001) Atoms are colliding if distance falls below this value
#  -primary      (    0.05) Initial loose criterion for atom equivalence
#  -final        (  0.0001) Final criterion for atom equivalence
#  -maxoptstep   (     0.5) Largest step allowed in symmetry element optimization
#  -minoptstep   (   1e-07) Termination criterion in symmetry element optimization
#  -gradstep     (   1e-07) Finite step used in numeric gradient evaluation
#  -minchange    (   1e-10) Minimum allowed change in target function
#  -minchgcycles (       5) Number of minchange cycles before optimization stops
#  -maxuniq                 Reorient molecule preferring more unique atoms
#  -minuniq      (default)  Reorient molecule preferring less unique atoms
#  -symmplane      (all)    Symmetrize atoms according to this plane number only
#  -symmaxis       (all)    Symmetrize atoms according to this axis number only
#  -symmimproper   (all)    Symmetrize atoms according to this improper axis number only
#  -symmcenter     (all)    Symmetrize atoms according to this inversion center number only
#  -na                      Prefer non-abelian symmetry elements for reorientation
#  -minaxis                 Prefer lower order axes for reorientation
#                           (use together with -na for tetrahedral groups)
#
#Input is expected in the following format:
#number_of_atoms
#AtomicNumber X Y Z
#...
#
#Note that only primitive rotations will be reported
#This is version $Revision: 1.15 $ ($Date: 2000/01/25 16:47:17 $)

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
  ($energy) = <$inp>=~/(-?\d+\.\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;
}

sub write_molden {
# &write_molden(\*FILEHANDLER [, -without_charges]);
# &write_molden([-without_charges]);   # To STDOUT
  if (@atom && @x && @y && @z) {
    my ($out, $i, $opts);
    if (@_) {
      $out = shift;
      if ($out=~/^-/) {
        $opts = $out;
        $out = 'STDOUT';
      } else {
        $opts = shift;
      }
    } else {
      $out = 'STDOUT';
    }
    my %opts = &get_options($opts);
    print $out " $#atom\n";
    print $out " Energy $energy" if $energy;
    print $out "   Symmetry $group" if $group;
    print $out "\n";
    for ($i=1; $i<=$#atom; $i++) {
      if ($atom[$i] && $x[$i] && $y[$i] && $z[$i]) {
        printf $out " %-2s %12.8f %12.8f %12.8f",$atom[$i],$x[$i],$y[$i],$z[$i];
        print $out "   $charge[$i]" if $charge[$i] && ! exists $opts{'without_charges'};
        print $out "\n";
      }
    }
  }
}

sub get_options {
#%h = &get_options('a b -c d -e \*A -f e');
#e => undef
#c => d
#f => e
  return unless $_[0];
   return $_[0]=~/-(\w+)(?:\s+)?(\w+)?/g;
}