Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/tools/msi2lmp/test/data-compare.pl
d1bbe9546c02lammm-devel

Last Change
Raw File

data-compare.pl
View Options

#!/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;
c4science · Help