#!/usr/bin/perl -w
# Tool to validate and compare two LAMMPS data files
# with "inexact" floating point comparisons
# July 2013 by Axel Kohlmeyer <akohlmey@gmail.com>

use strict;
use warnings;

my $version = 'v0.3';

# delta for floating point comparisons.
my $small = 1.0e-4;
# two hashes for storing system information
my %data1;
my %data2;

# simple checks after reading a section header
sub section_check {
  my ($fh,$data,$section) = @_;
  $_ = <$fh>;
  # skip empty and whitespace-only lines
  die "Line following '$section' is not empty" unless (/^\s*$/);
  die "Incomplete or incorrect header" unless ($data->{natoms} > 0);
  die "Incomplete or incorrect header" unless ($data->{natomtypes} > 0);
}

sub get_next {
  my ($fh) = @_;

  while (<$fh>) {
    chomp;
    # trim off comments
    $_ =~ s/#.*$//;
    # skip empty and whitespace-only lines
    next if (/^\s*$/);
    last;
  }
  return $_;
}

# fill hash with default data.
sub data_defaults {
  my ($data) = @_;

  $data->{natoms} = 0;
  $data->{nbonds} = 0;
  $data->{nangles} = 0;
  $data->{ndihedrals} = 0;
  $data->{nimpropers} = 0;

  $data->{natomtypes} = 0;
  $data->{nbondtypes} = 0;
  $data->{nangletypes} = 0;
  $data->{ndihedraltypes} = 0;
  $data->{nimpropertypes} = 0;

  $data->{xlo} =  0.5;
  $data->{xhi} = -0.5;
  $data->{ylo} =  0.5;
  $data->{yhi} = -0.5;
  $data->{zlo} =  0.5;
  $data->{zhi} = -0.5;
  $data->{xy} =   0.0;
  $data->{xz} =   0.0;
  $data->{yz} =   0.0;
  $data->{triclinic} = 0;
}

# read/parse lammps data file
sub read_data {
  my ($fh,$data) = @_;
  my $section;

  # read header. first line is already chopped off
  while (get_next($fh)) {

    if (/^\s*([0-9]+)\s+atoms\s*$/)     { $data->{natoms} = $1;     next; }
    if (/^\s*([0-9]+)\s+bonds\s*$/)     { $data->{nbonds} = $1;     next; }
    if (/^\s*([0-9]+)\s+angles\s*$/)    { $data->{nangles} = $1;    next; }
    if (/^\s*([0-9]+)\s+dihedrals\s*$/) { $data->{ndihedrals} = $1; next; }
    if (/^\s*([0-9]+)\s+impropers\s*$/) { $data->{nimpropers} = $1; next; }

    if (/^\s*([0-9]+)\s+atom types\s*$/)   { $data->{natomtypes} = $1;  next; }
    if (/^\s*([0-9]+)\s+bond types\s*$/)   { $data->{nbondtypes} = $1;  next; }
    if (/^\s*([0-9]+)\s+angle types\s*$/)  { $data->{nangletypes} = $1; next; }
    if (/^\s*([0-9]+)\s+dihedral types\s*$/)
      { $data->{ndihedraltypes} = $1; next; }
    if (/^\s*([0-9]+)\s+improper types\s*$/)
      { $data->{nimpropertypes} =$1 ; next; }

    if (/^\s*([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+xlo xhi\s*$/) {
      $data->{xlo}=$1;
      $data->{xhi}=$2;
      next;
    }
    if (/^\s*([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+ylo yhi\s*$/) {
      $data->{ylo}=$1;
      $data->{yhi}=$2;
      next;
    }
    if (/^\s*([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+zlo zhi\s*$/) {
      $data->{zlo}=$1;
      $data->{zhi}=$2;
      next;
    }
    if (/^\s*([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+xy xz yz\s*$/) {
      $data->{xy}=$1;
      $data->{xz}=$2;
      $data->{yz}=$3;
      $data->{triclinic} = 1;
      next;
    }

    # if we reach this point, the header has ended;
    last;
  }

  $a = $data->{natoms};
  $b = $data->{natomtypes};
  die "Invalid number of atoms: $a" unless ($a > 0);
  die "Invalid number of atom types: $b" unless ($b > 0);

  my ($i,$j,$k);
  while (1) {
    if (/^\s*(\S+|\S+ Coeffs)\s*$/) {
      if ($1 eq "Masses") {
        $data->{mass} = [];
        $i = 0;
        section_check($fh,$data,"Masses");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9]+)\s*$/) {
            $j = $1 - 1;
            die "Atom type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{natomtypes}));
            ++$i;
            $data->{mass}[$j] = $2;
            next;
          }

          die "Too many entries in Masses section"
            if ($i > $data->{natomtypes});
          die "Too few entries in Masses section"
            if ($i < $data->{natomtypes});
          die "Multiple mass assignments to the same atom type"
            if (scalar @{$data->{mass}} != $data->{natomtypes});

          last;
        }
      } elsif ($1 eq "Pair Coeffs") {
        $data->{paircoeff} = [];
        $i = 0;
        section_check($fh,$data,"Pair Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "Atom type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{natomtypes}));
            ++$i;
            $data->{paircoeff}[$j] = $2;
            next;
          }

          die "Too many entries in Pair Coeffs section"
            if ($i > $data->{natomtypes});
          die "Too few entries in Pair Coeffs section"
            if ($i < $data->{natomtypes});
          die "Multiple pair coefficient assignments to the same atom type"
            if (scalar @{$data->{paircoeff}} != $data->{natomtypes});

          last;
        }
      } elsif ($1 eq "Bond Coeffs") {
        $data->{bondcoeff} = [];
        $i = 0;
        section_check($fh,$data,"Bond Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "Bond type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{nbondtypes}));
            ++$i;
            $data->{bondcoeff}[$j] = $2;
            next;
          }

          die "Too many entries in Bond Coeffs section"
            if ($i > $data->{nbondtypes});
          die "Too few entries in Bond Coeffs section"
            if ($i < $data->{nbondtypes});
          die "Multiple bond coefficient assignments to the same bond type"
            if (scalar @{$data->{bondcoeff}} != $data->{nbondtypes});

          last;
        }
      } elsif ($1 eq "Angle Coeffs") {
        $data->{anglecoeff} = [];
        $i = 0;
        section_check($fh,$data,"Angle Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "Angle type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{nangletypes}));
            ++$i;
            $data->{anglecoeff}[$j] = $2;
            next;
          }

          die "Too many entries in Angle Coeffs section"
            if ($i > $data->{nangletypes});
          die "Too few entries in Angle Coeffs section"
            if ($i < $data->{nangletypes});
          die "Multiple angle coefficient assignments to the same angle type"
            if (scalar @{$data->{anglecoeff}} != $data->{nangletypes});

          last;
        }
      } elsif ($1 eq "BondBond Coeffs") {
        $data->{bondbondcoeff} = [];
        $i = 0;
        section_check($fh,$data,"BondBond Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "Angle type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{nangletypes}));
            ++$i;
            $data->{bondbondcoeff}[$j] = $2;
            next;
          }

          die "Too many entries in BondBond Coeffs section"
            if ($i > $data->{nangletypes});
          die "Too few entries in BondBond Coeffs section"
            if ($i < $data->{nangletypes});
          die "Multiple angle coefficient assignments to the same angle type"
            if (scalar @{$data->{bondbondcoeff}} != $data->{nangletypes});

          last;
        }
      } elsif ($1 eq "BondAngle Coeffs") {
        $data->{bondanglecoeff} = [];
        $i = 0;
        section_check($fh,$data,"BondAngle Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "Angle type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{nangletypes}));
            ++$i;
            $data->{bondanglecoeff}[$j] = $2;
            next;
          }

          die "Too many entries in BondAngle Coeffs section"
            if ($i > $data->{nangletypes});
          die "Too few entries in BondAngle Coeffs section"
            if ($i < $data->{nangletypes});
          die "Multiple bondangle coefficient assignments to the same angle type"
            if (scalar @{$data->{bondanglecoeff}} != $data->{nangletypes});

          last;
        }
      } elsif ($1 eq "Dihedral Coeffs") {
        $data->{dihedralcoeff} = [];
        $i = 0;
        section_check($fh,$data,"Dihedral Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "Dihedral type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
            ++$i;
            $data->{dihedralcoeff}[$j] = $2;
            next;
          }

          die "Too many entries in Dihedral Coeffs section"
            if ($i > $data->{ndihedraltypes});
          die "Too few entries in Dihedral Coeffs section"
            if ($i < $data->{ndihedraltypes});
          die "Multiple dihedral coefficient assignments to the same dihedral type"
            if (scalar @{$data->{dihedralcoeff}} != $data->{ndihedraltypes});

          last;
        }
      } elsif ($1 eq "AngleAngleTorsion Coeffs") {
        $data->{angleangletorsioncoeff} = [];
        $i = 0;
        section_check($fh,$data,"AngleAngleTorsion Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "AngleAngleTorsion type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
            ++$i;
            $data->{angleangletorsioncoeff}[$j] = $2;
            next;
          }

          die "Too many entries in AngleAngleTorsion Coeffs section"
            if ($i > $data->{ndihedraltypes});
          die "Too few entries in AngleAngleTorsion Coeffs section"
            if ($i < $data->{ndihedraltypes});
          die "Multiple dihedral coefficient assignments to the same dihedral type"
            if (scalar @{$data->{angleangletorsioncoeff}} != $data->{ndihedraltypes});

          last;
        }
      } elsif ($1 eq "EndBondTorsion Coeffs") {
        $data->{endbondtorsioncoeff} = [];
        $i = 0;
        section_check($fh,$data,"EndBondTorsion Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "EndBondTorsion type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
            ++$i;
            $data->{endbondtorsioncoeff}[$j] = $2;
            next;
          }

          die "Too many entries in EndBondTorsion Coeffs section"
            if ($i > $data->{ndihedraltypes});
          die "Too few entries in EndBondTorsion Coeffs section"
            if ($i < $data->{ndihedraltypes});
          die "Multiple dihedral coefficient assignments to the same dihedral type"
            if (scalar @{$data->{endbondtorsioncoeff}} != $data->{ndihedraltypes});

          last;
        }
      } elsif ($1 eq "MiddleBondTorsion Coeffs") {
        $data->{middlebondtorsioncoeff} = [];
        $i = 0;
        section_check($fh,$data,"MiddleBondTorsion Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "MiddleBondTorsion type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
            ++$i;
            $data->{middlebondtorsioncoeff}[$j] = $2;
            next;
          }

          die "Too many entries in MiddleBondTorsion Coeffs section"
            if ($i > $data->{ndihedraltypes});
          die "Too few entries in MiddleBondTorsion Coeffs section"
            if ($i < $data->{ndihedraltypes});
          die "Multiple dihedral coefficient assignments to the same dihedral type"
            if (scalar @{$data->{middlebondtorsioncoeff}} != $data->{ndihedraltypes});

          last;
        }
      } elsif ($1 eq "BondBond13 Coeffs") {
        $data->{bondbond13coeff} = [];
        $i = 0;
        section_check($fh,$data,"BondBond13 Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "BondBond13 type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
            ++$i;
            $data->{bondbond13coeff}[$j] = $2;
            next;
          }

          die "Too many entries in BondBond13 Coeffs section"
            if ($i > $data->{ndihedraltypes});
          die "Too few entries in BondBond13 Coeffs section"
            if ($i < $data->{ndihedraltypes});
          die "Multiple dihedral coefficient assignments to the same dihedral type"
            if (scalar @{$data->{bondbond13coeff}} != $data->{ndihedraltypes});

          last;
        }
      } elsif ($1 eq "AngleTorsion Coeffs") {
        $data->{angletorsioncoeff} = [];
        $i = 0;
        section_check($fh,$data,"AngleTorsion Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "AngleTorsion type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
            ++$i;
            $data->{angletorsioncoeff}[$j] = $2;
            next;
          }

          die "Too many entries in AngleTorsion Coeffs section"
            if ($i > $data->{ndihedraltypes});
          die "Too few entries in AngleTorsion Coeffs section"
            if ($i < $data->{ndihedraltypes});
          die "Multiple dihedral coefficient assignments to the same dihedral type"
            if (scalar @{$data->{angletorsioncoeff}} != $data->{ndihedraltypes});

          last;
        }
      } elsif ($1 eq "Improper Coeffs") {
        $data->{impropercoeff} = [];
        $i = 0;
        section_check($fh,$data,"Improper Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "Improper type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{nimpropertypes}));
            ++$i;
            $data->{impropercoeff}[$j] = $2;
            next;
          }

          die "Too many entries in Improper Coeffs section"
            if ($i > $data->{nimpropertypes});
          die "Too few entries in Improper Coeffs section"
            if ($i < $data->{nimpropertypes});
          die "Multiple improper coefficient assignments to the same improper type"
            if (scalar @{$data->{impropercoeff}} != $data->{nimpropertypes});

          last;
        }
      } elsif ($1 eq "AngleAngle Coeffs") {
        $data->{angleanglecoeff} = [];
        $i = 0;
        section_check($fh,$data,"AngleAngle Coeffs");

        while (get_next($fh)) {

          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
            $j = $1 - 1;
            die "AngleAngle type $1 is out of range" 
              if (($1 < 1) || ($1 > $data->{nimpropertypes}));
            ++$i;
            $data->{angleanglecoeff}[$j] = $2;
            next;
          }

          die "Too many entries in AngleAngle Coeffs section"
            if ($i > $data->{nimpropertypes});
          die "Too few entries in AngleAngle Coeffs section"
            if ($i < $data->{nimpropertypes});
          die "Multiple angleangle coefficient assignments to the same angle type"
            if (scalar @{$data->{angleanglecoeff}} != $data->{nimpropertypes});

          last;
        }
      } elsif ($1 eq "Atoms") {
        $data->{tag} = [];
        $data->{type} = [];
        $data->{molid} = [];
        $data->{charge} = [];
        $data->{posx} = [];
        $data->{posy} = [];
        $data->{posz} = [];
        $data->{imgx} = [];
        $data->{imgy} = [];
        $data->{imgz} = [];
        $i = 0;
        section_check($fh,$data,"Atoms");

        while (get_next($fh)) {
          if (/^\s*([0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+([-+.eE0-9]+)(|\s+(-?[0-9]+)\s+(-?[0-9]+)\s+(-?[0-9]+))\s*$/) {

            $k = $1 - 1;
            die "Atom id $1 is out of range"
              if (($1 < 1) || ($1 > $data->{natoms}));

            $j = $3 - 1;
            die "Atom type $2 is out of range"
              if (($3 < 1) || ($3 > $data->{natomtypes}));

            ++$i;
            $data->{tag}[$k] = $1;
            $data->{molid}[$k] = $2;
            $data->{type}[$k] = $3;
            $data->{charge}[$k] = $4;
            $data->{posx}[$k] = $5;
            $data->{posy}[$k] = $6;
            $data->{posz}[$k] = $7;
            $data->{imgx}[$k] = 0;
            $data->{imgy}[$k] = 0;
            $data->{imgz}[$k] = 0;
            if (! $8 eq "") {
              $data->{imgx}[$k] = $9;
              $data->{imgy}[$k] = $10;
              $data->{imgz}[$k] = $11;
            }
            next;
#          } else {
#            print "Atoms: $_\n";
          }

          die "Too many entries in Atoms section: $i vs. $data->{natoms}"
            if ($i > $data->{natoms});
          die "Too few entries in Atoms section: $i vs. $data->{natoms}"
            if ($i < $data->{natoms});
          die "Multiple atoms assigned to the same atom ID"
            if (scalar @{$data->{tag}} != $data->{natoms});

          last;
        }
      } elsif ($1 eq "Velocities") {
        $data->{velx} = [];
        $data->{vely} = [];
        $data->{velz} = [];
        $i = 0;
        section_check($fh,$data,"Velocities");

        while (get_next($fh)) {
          if (/^\s*([0-9]+)\s+([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+([-+.eE0-9]+)\s*$/) {

            $k = $1 - 1;
            die "Atom id $1 is out of range"
              if (($1 < 1) || ($1 > $data->{natoms}));

            ++$i;
            $data->{velx}[$k] = $2;
            $data->{vely}[$k] = $3;
            $data->{velz}[$k] = $4;
            next;
          }

          die "Too many entries in Velocities section"
            if ($i > $data->{natoms});
          die "Too few entries in Velocities section"
            if ($i < $data->{natoms});
          die "Multiple velocities assigned to the same atom ID"
            if (scalar @{$data->{velx}} != $data->{natoms});

          last;
        }
      } elsif ($1 eq "Bonds") {
        $data->{bondt} = [];
        $data->{bond1} = [];
        $data->{bond2} = [];
        $i = 0;
        section_check($fh,$data,"Bonds");

        while (get_next($fh)) {
          if (/^\s*([0-9]+)\s+(-?[0-9]+)\s+([0-9]+)\s+([0-9]+)\s*$/) {

            $k = $1 - 1;
            die "Bond id $1 is out of range"
              if (($1 < 1) || ($1 > $data->{nbonds}));

            die "Bond type $2 is out of range"
              if (($2 == 0) || ($2 > $data->{nbondtypes}));

            die "Bond atom 1 ID $3 is out of range"
              if (($3 < 1) || ($3 > $data->{natoms}));
            die "Bond atom 2 ID $4 is out of range"
              if (($4 < 1) || ($4 > $data->{natoms}));

            ++$i;
            $data->{bondt}[$k] = $2;
            $data->{bond1}[$k] = $3;
            $data->{bond2}[$k] = $4;
            next;
          }

          die "Too many entries in Bonds section"
            if ($i > $data->{nbonds});
          die "Too few entries in Bonds section"
            if ($i < $data->{nbonds});
          die "Multiple bonds assigned to the same bond ID"
            if (scalar @{$data->{bondt}} != $data->{nbonds});

          last;
        }
      } elsif ($1 eq "Angles") {
        $data->{anglet} = [];
        $data->{angle1} = [];
        $data->{angle2} = [];
        $data->{angle3} = [];
        $i = 0;
        section_check($fh,$data,"Angles");

        while (get_next($fh)) {
          if (/^\s*([0-9]+)\s+(-?[0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s*$/) {

            $k = $1 - 1;
            die "Angle id $1 is out of range"
              if (($1 < 1) || ($1 > $data->{nangles}));

            die "Angle type $2 is out of range"
              if (($2 == 0) || ($2 > $data->{nangletypes}));

            die "Angle atom 1 ID $3 is out of range"
              if (($3 < 1) || ($3 > $data->{natoms}));
            die "Angle atom 2 ID $4 is out of range"
              if (($4 < 1) || ($4 > $data->{natoms}));
            die "Angle atom 3 ID $5 is out of range"
              if (($5 < 1) || ($5 > $data->{natoms}));

            ++$i;
            $data->{anglet}[$k] = $2;
            $data->{angle1}[$k] = $3;
            $data->{angle2}[$k] = $4;
            $data->{angle3}[$k] = $5;
            next;
          }

          die "Too many entries in Angles section"
            if ($i > $data->{nangles});
          die "Too few entries in Angles section"
            if ($i < $data->{nangles});
          die "Multiple angles assigned to the same angle ID"
            if (scalar @{$data->{anglet}} != $data->{nangles});

          last;
        }
      } elsif ($1 eq "Dihedrals") {
        $data->{dihedralt} = [];
        $data->{dihedral1} = [];
        $data->{dihedral2} = [];
        $data->{dihedral3} = [];
        $data->{dihedral4} = [];
        $i = 0;
        section_check($fh,$data,"Dihedrals");

        while (get_next($fh)) {
          if (/^\s*([0-9]+)\s+(-?[0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s*$/) {

            $k = $1 - 1;
            die "Dihedral id $1 is out of range"
              if (($1 < 1) || ($1 > $data->{ndihedrals}));

            die "Dihedral type $2 is out of range"
              if (($2 == 0) || ($2 > $data->{ndihedraltypes}));

            die "Dihedral atom 1 ID $3 is out of range"
              if (($3 < 1) || ($3 > $data->{natoms}));
            die "Dihedral atom 2 ID $4 is out of range"
              if (($4 < 1) || ($4 > $data->{natoms}));
            die "Dihedral atom 3 ID $5 is out of range"
              if (($5 < 1) || ($5 > $data->{natoms}));
            die "Dihedral atom 4 ID $6 is out of range"
              if (($6 < 1) || ($6 > $data->{natoms}));

            ++$i;
            $data->{dihedralt}[$k] = $2;
            $data->{dihedral1}[$k] = $3;
            $data->{dihedral2}[$k] = $4;
            $data->{dihedral3}[$k] = $5;
            $data->{dihedral4}[$k] = $6;
            next;
          }

          die "Too many entries in Dihedrals section"
            if ($i > $data->{ndihedrals});
          die "Too few entries in Dihedrals section"
            if ($i < $data->{ndihedrals});
          die "Multiple dihedrals assigned to the same dihedral ID"
            if (scalar @{$data->{dihedralt}} != $data->{ndihedrals});

          last;
        }
      } elsif ($1 eq "Impropers") {
        $data->{impropert} = [];
        $data->{improper1} = [];
        $data->{improper2} = [];
        $data->{improper3} = [];
        $data->{improper4} = [];
        $i = 0;
        section_check($fh,$data,"Impropers");

        while (get_next($fh)) {
          if (/^\s*([0-9]+)\s+(-?[0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s*$/) {

            $k = $1 - 1;
            die "Improper id $1 is out of range"
              if (($1 < 1) || ($1 > $data->{nimpropers}));

            die "Improper type $2 is out of range"
              if (($2 == 0) || ($2 > $data->{nimpropertypes}));

            die "Improper atom 1 ID $3 is out of range"
              if (($3 < 1) || ($3 > $data->{natoms}));
            die "Improper atom 2 ID $4 is out of range"
              if (($4 < 1) || ($4 > $data->{natoms}));
            die "Improper atom 3 ID $5 is out of range"
              if (($5 < 1) || ($5 > $data->{natoms}));
            die "Improper atom 4 ID $6 is out of range"
              if (($6 < 1) || ($6 > $data->{natoms}));

            ++$i;
            $data->{impropert}[$k] = $2;
            $data->{improper1}[$k] = $3;
            $data->{improper2}[$k] = $4;
            $data->{improper3}[$k] = $5;
            $data->{improper4}[$k] = $6;
            next;
          }

          die "Too many entries in Impropers section"
            if ($i > $data->{nimpropers});
          die "Too few entries in Impropers section"
            if ($i < $data->{nimpropers});
          die "Multiple impropers assigned to the same improper ID"
            if (scalar @{$data->{impropert}} != $data->{nimpropers});

          last;
        }
      } else {
        die "Bad data: $_";
      }

      last unless ($_);
    } else {
      die "Bad data: $_";
    }
  }
}

sub floatdiff {
  my ($n1,$n2,$rel) = @_;

  my $diff = abs($n1-$n2);
  my $avg = (abs($n1)+abs($n2))*0.5;
  return 0 if ($avg == 0.0);
  if ($rel) {
#    print "relative difference: ",$diff/$avg," vs. $small\n";
    return 0 if ($diff/$avg < $small);
  } else {
#    print "absolute difference: ",$diff," vs. $small\n";
    return 0 if ($diff < $small);
  }
  return 1;
}

sub coeffcompare {
  my ($d1,$d2,$coeff,$type) = @_;
  my (@c1,@c2,$a,$b);
  my ($field,$count,$i,$j,$t) = ($coeff . 'coeff', 'n' . $type . 'types', 0,0,0);

  if (exists $d1->{$field} && exists $d2->{$field}) {
    for ($i=0; $i < $d1->{$count}; ++$i) {
      $t = $i+1;
      @c1 = split /\s+/, ${$$d1{$field}}[$i];
      @c2 = split /\s+/, ${$$d2{$field}}[$i];
      die "Inconsistent number of $coeff coefficients for $type type $t: $#c1 vs $#c2\n"
        if ($#c1 != $#c2);

      for ($j = 0; $j <= $#c1; ++$j) {
        $a = $c1[$j]; $b = $c2[$j];
        die "Inconsistent $coeff coefficient ", $j+1,
          " for $type type $t: $a vs. $b" if (floatdiff($a,$b,1));
      }
    }
  } else {
    die "Field $field only exists in data file 1" if (exists $d1->{$field});
    die "Field $field only exists in data file 2" if (exists $d2->{$field});
  }
}

sub topocompare {
  my ($d1,$d2,$type,$count) = @_;
  my ($num,$a,$b,$field,$i,$j,$t);

  $field = 'n' . $type . 's';
  $num = $d1->{$field};

  for ($i=0; $i < $num; ++$i) {
    $t = $i+1;
    $field = $type . 't';
    $a = $d1->{$field}[$i]; $b = $d2->{$field}[$i];
    die "Inconsistent $type types for $type $t: $a vs. $b" if ($a != $b);
    for ($j=1; $j <= $count; ++$j) {
      $field = $type . $j;
      $a = $d1->{$field}[$i]; $b = $d2->{$field}[$i];
      die "Inconsistent $type atom $j for $type $t: $a vs. $b" if ($a != $b);
    }
  }
}

sub syscompare {
  my ($d1,$d2) = @_;
  my ($i,$j,$t,$a,$b,@l);

  # check atoms.
  die "Number of atoms does not match"
    if ($d1->{natoms} != $d2->{natoms});
  die "Number of atom types does not match"
    if ($d1->{natomtypes} != $d2->{natomtypes});

  # check bonded interactions
  @l = ('bond','angle','dihedral','improper');
  foreach $i (@l) {
    $t = 'n' . $i . 's';
    $a = $d1->{$t};
    $b = $d2->{$t};
    die "Number of ",$i,"s does not match: $a vs $b" unless ($a == $b);

    $t = 'n' . $i . 'types';
    $a = $d1->{$t};
    $b = $d2->{$t};
    die "Number of ",$i," types does not match: $a vs $b" unless ($a == $b);
  }

  topocompare($d1,$d2,'bond',2);
  topocompare($d1,$d2,'angle',3);
  topocompare($d1,$d2,'dihedral',4);
  topocompare($d1,$d2,'improper',4);

  coeffcompare($d1,$d2,'pair','atom');
  coeffcompare($d1,$d2,'bond','bond');
  coeffcompare($d1,$d2,'angle','angle');
  coeffcompare($d1,$d2,'bondbond','angle');
  coeffcompare($d1,$d2,'bondangle','angle');
  coeffcompare($d1,$d2,'dihedral','dihedral');
  coeffcompare($d1,$d2,'angleangletorsion','dihedral');
  coeffcompare($d1,$d2,'bondbond13','dihedral');
  coeffcompare($d1,$d2,'endbondtorsion','dihedral');
  coeffcompare($d1,$d2,'middlebondtorsion','dihedral');
  coeffcompare($d1,$d2,'improper','improper');
  coeffcompare($d1,$d2,'angleangle','improper');

  for ($i=0; $i < $d1->{natomtypes}; ++$i) {
    $j = $i+1;
    if (exists $d1->{mass}[$i]) {
      $a = $d1->{mass}[$i];
    } else {
      die "No mass for atom type $j in data file 1";
    }
    if (exists $d2->{mass}[$i]) {
      $a = $d2->{mass}[$i];
    } else {
      die "No mass for atom type $j in data file 2";
    }
  }

  # check box information
  die "Inconsistent box shape" if ($d1->{triclinic} != $d2->{triclinic});

  @l = ('xlo','xhi','ylo','yhi','zlo','zhi');
  if ($d1->{triclinic}) { push @l, ('xy','xz','yz'); }
  for $i (@l) {
    $a = $d1->{$i};
    $b = $d2->{$i};
    die "Box data for $i does not match: $a $b" if (floatdiff($a,$b,0));
  }

  for ($i=0; $i < $d1->{natoms}; ++$i) {
    $j = $i+1;
    for $t ('tag','molid','type','imgx','imgy','imgz') {
      if (exists $d1->{$t}[$i]) {
        $a = $d1->{$t}[$i];
      } else {
        $a = 0;
      }
      if (exists $d2->{$t}[$i]) {
        $b = $d2->{$t}[$i];
      } else {
        $b = 0;
      }
      die "Inconsistent data for $t, atom $j: $a vs. $b" if ($a != $b);
    }

    for $t ('charge','posx','posy','posz') {
      if (exists $d1->{$t}[$i]) {
        $a = $d1->{$t}[$i];
      } else {
        $a = 0;
      }
      if (exists $d2->{$t}[$i]) {
        $b = $d2->{$t}[$i];
      } else {
        $b = 0;
      }
      die "Inconsistent data for $t, atom $j: $a vs. $b" if (floatdiff($a,$b,0));
    }
  }

}

########################################################################
# main program

my $fp;

if ($#ARGV < 1) {
  die "usage $0 <file 1> <file 2>";
}

print "\nLAMMPS data file validation tool. $version\n\n";

data_defaults(\%data1);
data_defaults(\%data2);

# read in first data file
open($fp, '<', $ARGV[0]) or die $!;
print "opened data file 1: $ARGV[0]\n";
$_=<$fp>;
print;
read_data($fp,\%data1);
print "done reading data file 1\n\n";
close $fp;

# read in second data file
open($fp, '<', $ARGV[1]) or die $!;
print "opened data file 2: $ARGV[1]\n";
$_=<$fp>;
print;
read_data($fp,\%data2);
print "done reading data file 2\n\n";
close $fp;

# compare data sets
syscompare(\%data1,\%data2);

print "File $ARGV[0] and $ARGV[1] match\n\n";

exit 0;