1#!/usr/bin/perl -w
2# Tool to validate and compare two LAMMPS data files
3# with "inexact" floating point comparisons
4# July 2013 by Axel Kohlmeyer <akohlmey@gmail.com>
5
6use strict;
7use warnings;
8
9my $version = 'v0.3';
10
11# delta for floating point comparisons.
12my $small = 1.0e-4;
13# two hashes for storing system information
14my %data1;
15my %data2;
16
17# simple checks after reading a section header
18sub section_check {
19  my ($fh,$data,$section) = @_;
20  $_ = <$fh>;
21  # skip empty and whitespace-only lines
22  die "Line following '$section' is not empty" unless (/^\s*$/);
23  die "Incomplete or incorrect header" unless ($data->{natoms} > 0);
24  die "Incomplete or incorrect header" unless ($data->{natomtypes} > 0);
25}
26
27sub get_next {
28  my ($fh) = @_;
29
30  while (<$fh>) {
31    chomp;
32    # trim off comments
33    $_ =~ s/#.*$//;
34    # skip empty and whitespace-only lines
35    next if (/^\s*$/);
36    last;
37  }
38  return $_;
39}
40
41# fill hash with default data.
42sub data_defaults {
43  my ($data) = @_;
44
45  $data->{natoms} = 0;
46  $data->{nbonds} = 0;
47  $data->{nangles} = 0;
48  $data->{ndihedrals} = 0;
49  $data->{nimpropers} = 0;
50
51  $data->{natomtypes} = 0;
52  $data->{nbondtypes} = 0;
53  $data->{nangletypes} = 0;
54  $data->{ndihedraltypes} = 0;
55  $data->{nimpropertypes} = 0;
56
57  $data->{xlo} =  0.5;
58  $data->{xhi} = -0.5;
59  $data->{ylo} =  0.5;
60  $data->{yhi} = -0.5;
61  $data->{zlo} =  0.5;
62  $data->{zhi} = -0.5;
63  $data->{xy} =   0.0;
64  $data->{xz} =   0.0;
65  $data->{yz} =   0.0;
66  $data->{triclinic} = 0;
67}
68
69# read/parse lammps data file
70sub read_data {
71  my ($fh,$data) = @_;
72  my $section;
73
74  # read header. first line is already chopped off
75  while (get_next($fh)) {
76
77    if (/^\s*([0-9]+)\s+atoms\s*$/)     { $data->{natoms} = $1;     next; }
78    if (/^\s*([0-9]+)\s+bonds\s*$/)     { $data->{nbonds} = $1;     next; }
79    if (/^\s*([0-9]+)\s+angles\s*$/)    { $data->{nangles} = $1;    next; }
80    if (/^\s*([0-9]+)\s+dihedrals\s*$/) { $data->{ndihedrals} = $1; next; }
81    if (/^\s*([0-9]+)\s+impropers\s*$/) { $data->{nimpropers} = $1; next; }
82
83    if (/^\s*([0-9]+)\s+atom types\s*$/)   { $data->{natomtypes} = $1;  next; }
84    if (/^\s*([0-9]+)\s+bond types\s*$/)   { $data->{nbondtypes} = $1;  next; }
85    if (/^\s*([0-9]+)\s+angle types\s*$/)  { $data->{nangletypes} = $1; next; }
86    if (/^\s*([0-9]+)\s+dihedral types\s*$/)
87      { $data->{ndihedraltypes} = $1; next; }
88    if (/^\s*([0-9]+)\s+improper types\s*$/)
89      { $data->{nimpropertypes} =$1 ; next; }
90
91    if (/^\s*([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+xlo xhi\s*$/) {
92      $data->{xlo}=$1;
93      $data->{xhi}=$2;
94      next;
95    }
96    if (/^\s*([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+ylo yhi\s*$/) {
97      $data->{ylo}=$1;
98      $data->{yhi}=$2;
99      next;
100    }
101    if (/^\s*([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+zlo zhi\s*$/) {
102      $data->{zlo}=$1;
103      $data->{zhi}=$2;
104      next;
105    }
106    if (/^\s*([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+xy xz yz\s*$/) {
107      $data->{xy}=$1;
108      $data->{xz}=$2;
109      $data->{yz}=$3;
110      $data->{triclinic} = 1;
111      next;
112    }
113
114    # if we reach this point, the header has ended;
115    last;
116  }
117
118  $a = $data->{natoms};
119  $b = $data->{natomtypes};
120  die "Invalid number of atoms: $a" unless ($a > 0);
121  die "Invalid number of atom types: $b" unless ($b > 0);
122
123  my ($i,$j,$k);
124  while (1) {
125    if (/^\s*(\S+|\S+ Coeffs)\s*$/) {
126      if ($1 eq "Masses") {
127        $data->{mass} = [];
128        $i = 0;
129        section_check($fh,$data,"Masses");
130
131        while (get_next($fh)) {
132
133          if (/^\s*([0-9]+)\s+([-+.eE0-9]+)\s*$/) {
134            $j = $1 - 1;
135            die "Atom type $1 is out of range"
136              if (($1 < 1) || ($1 > $data->{natomtypes}));
137            ++$i;
138            $data->{mass}[$j] = $2;
139            next;
140          }
141
142          die "Too many entries in Masses section"
143            if ($i > $data->{natomtypes});
144          die "Too few entries in Masses section"
145            if ($i < $data->{natomtypes});
146          die "Multiple mass assignments to the same atom type"
147            if (scalar @{$data->{mass}} != $data->{natomtypes});
148
149          last;
150        }
151      } elsif ($1 eq "Pair Coeffs") {
152        $data->{paircoeff} = [];
153        $i = 0;
154        section_check($fh,$data,"Pair Coeffs");
155
156        while (get_next($fh)) {
157
158          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
159            $j = $1 - 1;
160            die "Atom type $1 is out of range"
161              if (($1 < 1) || ($1 > $data->{natomtypes}));
162            ++$i;
163            $data->{paircoeff}[$j] = $2;
164            next;
165          }
166
167          die "Too many entries in Pair Coeffs section"
168            if ($i > $data->{natomtypes});
169          die "Too few entries in Pair Coeffs section"
170            if ($i < $data->{natomtypes});
171          die "Multiple pair coefficient assignments to the same atom type"
172            if (scalar @{$data->{paircoeff}} != $data->{natomtypes});
173
174          last;
175        }
176      } elsif ($1 eq "Bond Coeffs") {
177        $data->{bondcoeff} = [];
178        $i = 0;
179        section_check($fh,$data,"Bond Coeffs");
180
181        while (get_next($fh)) {
182
183          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
184            $j = $1 - 1;
185            die "Bond type $1 is out of range"
186              if (($1 < 1) || ($1 > $data->{nbondtypes}));
187            ++$i;
188            $data->{bondcoeff}[$j] = $2;
189            next;
190          }
191
192          die "Too many entries in Bond Coeffs section"
193            if ($i > $data->{nbondtypes});
194          die "Too few entries in Bond Coeffs section"
195            if ($i < $data->{nbondtypes});
196          die "Multiple bond coefficient assignments to the same bond type"
197            if (scalar @{$data->{bondcoeff}} != $data->{nbondtypes});
198
199          last;
200        }
201      } elsif ($1 eq "Angle Coeffs") {
202        $data->{anglecoeff} = [];
203        $i = 0;
204        section_check($fh,$data,"Angle Coeffs");
205
206        while (get_next($fh)) {
207
208          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
209            $j = $1 - 1;
210            die "Angle type $1 is out of range"
211              if (($1 < 1) || ($1 > $data->{nangletypes}));
212            ++$i;
213            $data->{anglecoeff}[$j] = $2;
214            next;
215          }
216
217          die "Too many entries in Angle Coeffs section"
218            if ($i > $data->{nangletypes});
219          die "Too few entries in Angle Coeffs section"
220            if ($i < $data->{nangletypes});
221          die "Multiple angle coefficient assignments to the same angle type"
222            if (scalar @{$data->{anglecoeff}} != $data->{nangletypes});
223
224          last;
225        }
226      } elsif ($1 eq "BondBond Coeffs") {
227        $data->{bondbondcoeff} = [];
228        $i = 0;
229        section_check($fh,$data,"BondBond Coeffs");
230
231        while (get_next($fh)) {
232
233          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
234            $j = $1 - 1;
235            die "Angle type $1 is out of range"
236              if (($1 < 1) || ($1 > $data->{nangletypes}));
237            ++$i;
238            $data->{bondbondcoeff}[$j] = $2;
239            next;
240          }
241
242          die "Too many entries in BondBond Coeffs section"
243            if ($i > $data->{nangletypes});
244          die "Too few entries in BondBond Coeffs section"
245            if ($i < $data->{nangletypes});
246          die "Multiple angle coefficient assignments to the same angle type"
247            if (scalar @{$data->{bondbondcoeff}} != $data->{nangletypes});
248
249          last;
250        }
251      } elsif ($1 eq "BondAngle Coeffs") {
252        $data->{bondanglecoeff} = [];
253        $i = 0;
254        section_check($fh,$data,"BondAngle Coeffs");
255
256        while (get_next($fh)) {
257
258          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
259            $j = $1 - 1;
260            die "Angle type $1 is out of range"
261              if (($1 < 1) || ($1 > $data->{nangletypes}));
262            ++$i;
263            $data->{bondanglecoeff}[$j] = $2;
264            next;
265          }
266
267          die "Too many entries in BondAngle Coeffs section"
268            if ($i > $data->{nangletypes});
269          die "Too few entries in BondAngle Coeffs section"
270            if ($i < $data->{nangletypes});
271          die "Multiple bondangle coefficient assignments to the same angle type"
272            if (scalar @{$data->{bondanglecoeff}} != $data->{nangletypes});
273
274          last;
275        }
276      } elsif ($1 eq "Dihedral Coeffs") {
277        $data->{dihedralcoeff} = [];
278        $i = 0;
279        section_check($fh,$data,"Dihedral Coeffs");
280
281        while (get_next($fh)) {
282
283          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
284            $j = $1 - 1;
285            die "Dihedral type $1 is out of range"
286              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
287            ++$i;
288            $data->{dihedralcoeff}[$j] = $2;
289            next;
290          }
291
292          die "Too many entries in Dihedral Coeffs section"
293            if ($i > $data->{ndihedraltypes});
294          die "Too few entries in Dihedral Coeffs section"
295            if ($i < $data->{ndihedraltypes});
296          die "Multiple dihedral coefficient assignments to the same dihedral type"
297            if (scalar @{$data->{dihedralcoeff}} != $data->{ndihedraltypes});
298
299          last;
300        }
301      } elsif ($1 eq "AngleAngleTorsion Coeffs") {
302        $data->{angleangletorsioncoeff} = [];
303        $i = 0;
304        section_check($fh,$data,"AngleAngleTorsion Coeffs");
305
306        while (get_next($fh)) {
307
308          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
309            $j = $1 - 1;
310            die "AngleAngleTorsion type $1 is out of range"
311              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
312            ++$i;
313            $data->{angleangletorsioncoeff}[$j] = $2;
314            next;
315          }
316
317          die "Too many entries in AngleAngleTorsion Coeffs section"
318            if ($i > $data->{ndihedraltypes});
319          die "Too few entries in AngleAngleTorsion Coeffs section"
320            if ($i < $data->{ndihedraltypes});
321          die "Multiple dihedral coefficient assignments to the same dihedral type"
322            if (scalar @{$data->{angleangletorsioncoeff}} != $data->{ndihedraltypes});
323
324          last;
325        }
326      } elsif ($1 eq "EndBondTorsion Coeffs") {
327        $data->{endbondtorsioncoeff} = [];
328        $i = 0;
329        section_check($fh,$data,"EndBondTorsion Coeffs");
330
331        while (get_next($fh)) {
332
333          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
334            $j = $1 - 1;
335            die "EndBondTorsion type $1 is out of range"
336              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
337            ++$i;
338            $data->{endbondtorsioncoeff}[$j] = $2;
339            next;
340          }
341
342          die "Too many entries in EndBondTorsion Coeffs section"
343            if ($i > $data->{ndihedraltypes});
344          die "Too few entries in EndBondTorsion Coeffs section"
345            if ($i < $data->{ndihedraltypes});
346          die "Multiple dihedral coefficient assignments to the same dihedral type"
347            if (scalar @{$data->{endbondtorsioncoeff}} != $data->{ndihedraltypes});
348
349          last;
350        }
351      } elsif ($1 eq "MiddleBondTorsion Coeffs") {
352        $data->{middlebondtorsioncoeff} = [];
353        $i = 0;
354        section_check($fh,$data,"MiddleBondTorsion Coeffs");
355
356        while (get_next($fh)) {
357
358          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
359            $j = $1 - 1;
360            die "MiddleBondTorsion type $1 is out of range"
361              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
362            ++$i;
363            $data->{middlebondtorsioncoeff}[$j] = $2;
364            next;
365          }
366
367          die "Too many entries in MiddleBondTorsion Coeffs section"
368            if ($i > $data->{ndihedraltypes});
369          die "Too few entries in MiddleBondTorsion Coeffs section"
370            if ($i < $data->{ndihedraltypes});
371          die "Multiple dihedral coefficient assignments to the same dihedral type"
372            if (scalar @{$data->{middlebondtorsioncoeff}} != $data->{ndihedraltypes});
373
374          last;
375        }
376      } elsif ($1 eq "BondBond13 Coeffs") {
377        $data->{bondbond13coeff} = [];
378        $i = 0;
379        section_check($fh,$data,"BondBond13 Coeffs");
380
381        while (get_next($fh)) {
382
383          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
384            $j = $1 - 1;
385            die "BondBond13 type $1 is out of range"
386              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
387            ++$i;
388            $data->{bondbond13coeff}[$j] = $2;
389            next;
390          }
391
392          die "Too many entries in BondBond13 Coeffs section"
393            if ($i > $data->{ndihedraltypes});
394          die "Too few entries in BondBond13 Coeffs section"
395            if ($i < $data->{ndihedraltypes});
396          die "Multiple dihedral coefficient assignments to the same dihedral type"
397            if (scalar @{$data->{bondbond13coeff}} != $data->{ndihedraltypes});
398
399          last;
400        }
401      } elsif ($1 eq "AngleTorsion Coeffs") {
402        $data->{angletorsioncoeff} = [];
403        $i = 0;
404        section_check($fh,$data,"AngleTorsion Coeffs");
405
406        while (get_next($fh)) {
407
408          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
409            $j = $1 - 1;
410            die "AngleTorsion type $1 is out of range"
411              if (($1 < 1) || ($1 > $data->{ndihedraltypes}));
412            ++$i;
413            $data->{angletorsioncoeff}[$j] = $2;
414            next;
415          }
416
417          die "Too many entries in AngleTorsion Coeffs section"
418            if ($i > $data->{ndihedraltypes});
419          die "Too few entries in AngleTorsion Coeffs section"
420            if ($i < $data->{ndihedraltypes});
421          die "Multiple dihedral coefficient assignments to the same dihedral type"
422            if (scalar @{$data->{angletorsioncoeff}} != $data->{ndihedraltypes});
423
424          last;
425        }
426      } elsif ($1 eq "Improper Coeffs") {
427        $data->{impropercoeff} = [];
428        $i = 0;
429        section_check($fh,$data,"Improper Coeffs");
430
431        while (get_next($fh)) {
432
433          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
434            $j = $1 - 1;
435            die "Improper type $1 is out of range"
436              if (($1 < 1) || ($1 > $data->{nimpropertypes}));
437            ++$i;
438            $data->{impropercoeff}[$j] = $2;
439            next;
440          }
441
442          die "Too many entries in Improper Coeffs section"
443            if ($i > $data->{nimpropertypes});
444          die "Too few entries in Improper Coeffs section"
445            if ($i < $data->{nimpropertypes});
446          die "Multiple improper coefficient assignments to the same improper type"
447            if (scalar @{$data->{impropercoeff}} != $data->{nimpropertypes});
448
449          last;
450        }
451      } elsif ($1 eq "AngleAngle Coeffs") {
452        $data->{angleanglecoeff} = [];
453        $i = 0;
454        section_check($fh,$data,"AngleAngle Coeffs");
455
456        while (get_next($fh)) {
457
458          if (/^\s*([0-9]+)\s+([-+.eE0-9 ]+)\s*$/) {
459            $j = $1 - 1;
460            die "AngleAngle type $1 is out of range"
461              if (($1 < 1) || ($1 > $data->{nimpropertypes}));
462            ++$i;
463            $data->{angleanglecoeff}[$j] = $2;
464            next;
465          }
466
467          die "Too many entries in AngleAngle Coeffs section"
468            if ($i > $data->{nimpropertypes});
469          die "Too few entries in AngleAngle Coeffs section"
470            if ($i < $data->{nimpropertypes});
471          die "Multiple angleangle coefficient assignments to the same angle type"
472            if (scalar @{$data->{angleanglecoeff}} != $data->{nimpropertypes});
473
474          last;
475        }
476      } elsif ($1 eq "Atoms") {
477        $data->{tag} = [];
478        $data->{type} = [];
479        $data->{molid} = [];
480        $data->{charge} = [];
481        $data->{posx} = [];
482        $data->{posy} = [];
483        $data->{posz} = [];
484        $data->{imgx} = [];
485        $data->{imgy} = [];
486        $data->{imgz} = [];
487        $i = 0;
488        section_check($fh,$data,"Atoms");
489
490        while (get_next($fh)) {
491          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*$/) {
492
493            $k = $1 - 1;
494            die "Atom id $1 is out of range"
495              if (($1 < 1) || ($1 > $data->{natoms}));
496
497            $j = $3 - 1;
498            die "Atom type $2 is out of range"
499              if (($3 < 1) || ($3 > $data->{natomtypes}));
500
501            ++$i;
502            $data->{tag}[$k] = $1;
503            $data->{molid}[$k] = $2;
504            $data->{type}[$k] = $3;
505            $data->{charge}[$k] = $4;
506            $data->{posx}[$k] = $5;
507            $data->{posy}[$k] = $6;
508            $data->{posz}[$k] = $7;
509            $data->{imgx}[$k] = 0;
510            $data->{imgy}[$k] = 0;
511            $data->{imgz}[$k] = 0;
512            if (! $8 eq "") {
513              $data->{imgx}[$k] = $9;
514              $data->{imgy}[$k] = $10;
515              $data->{imgz}[$k] = $11;
516            }
517            next;
518#          } else {
519#            print "Atoms: $_\n";
520          }
521
522          die "Too many entries in Atoms section: $i vs. $data->{natoms}"
523            if ($i > $data->{natoms});
524          die "Too few entries in Atoms section: $i vs. $data->{natoms}"
525            if ($i < $data->{natoms});
526          die "Multiple atoms assigned to the same atom ID"
527            if (scalar @{$data->{tag}} != $data->{natoms});
528
529          last;
530        }
531      } elsif ($1 eq "Velocities") {
532        $data->{velx} = [];
533        $data->{vely} = [];
534        $data->{velz} = [];
535        $i = 0;
536        section_check($fh,$data,"Velocities");
537
538        while (get_next($fh)) {
539          if (/^\s*([0-9]+)\s+([-+.eE0-9]+)\s+([-+.eE0-9]+)\s+([-+.eE0-9]+)\s*$/) {
540
541            $k = $1 - 1;
542            die "Atom id $1 is out of range"
543              if (($1 < 1) || ($1 > $data->{natoms}));
544
545            ++$i;
546            $data->{velx}[$k] = $2;
547            $data->{vely}[$k] = $3;
548            $data->{velz}[$k] = $4;
549            next;
550          }
551
552          die "Too many entries in Velocities section"
553            if ($i > $data->{natoms});
554          die "Too few entries in Velocities section"
555            if ($i < $data->{natoms});
556          die "Multiple velocities assigned to the same atom ID"
557            if (scalar @{$data->{velx}} != $data->{natoms});
558
559          last;
560        }
561      } elsif ($1 eq "Bonds") {
562        $data->{bondt} = [];
563        $data->{bond1} = [];
564        $data->{bond2} = [];
565        $i = 0;
566        section_check($fh,$data,"Bonds");
567
568        while (get_next($fh)) {
569          if (/^\s*([0-9]+)\s+(-?[0-9]+)\s+([0-9]+)\s+([0-9]+)\s*$/) {
570
571            $k = $1 - 1;
572            die "Bond id $1 is out of range"
573              if (($1 < 1) || ($1 > $data->{nbonds}));
574
575            die "Bond type $2 is out of range"
576              if (($2 == 0) || ($2 > $data->{nbondtypes}));
577
578            die "Bond atom 1 ID $3 is out of range"
579              if (($3 < 1) || ($3 > $data->{natoms}));
580            die "Bond atom 2 ID $4 is out of range"
581              if (($4 < 1) || ($4 > $data->{natoms}));
582
583            ++$i;
584            $data->{bondt}[$k] = $2;
585            $data->{bond1}[$k] = $3;
586            $data->{bond2}[$k] = $4;
587            next;
588          }
589
590          die "Too many entries in Bonds section"
591            if ($i > $data->{nbonds});
592          die "Too few entries in Bonds section"
593            if ($i < $data->{nbonds});
594          die "Multiple bonds assigned to the same bond ID"
595            if (scalar @{$data->{bondt}} != $data->{nbonds});
596
597          last;
598        }
599      } elsif ($1 eq "Angles") {
600        $data->{anglet} = [];
601        $data->{angle1} = [];
602        $data->{angle2} = [];
603        $data->{angle3} = [];
604        $i = 0;
605        section_check($fh,$data,"Angles");
606
607        while (get_next($fh)) {
608          if (/^\s*([0-9]+)\s+(-?[0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s*$/) {
609
610            $k = $1 - 1;
611            die "Angle id $1 is out of range"
612              if (($1 < 1) || ($1 > $data->{nangles}));
613
614            die "Angle type $2 is out of range"
615              if (($2 == 0) || ($2 > $data->{nangletypes}));
616
617            die "Angle atom 1 ID $3 is out of range"
618              if (($3 < 1) || ($3 > $data->{natoms}));
619            die "Angle atom 2 ID $4 is out of range"
620              if (($4 < 1) || ($4 > $data->{natoms}));
621            die "Angle atom 3 ID $5 is out of range"
622              if (($5 < 1) || ($5 > $data->{natoms}));
623
624            ++$i;
625            $data->{anglet}[$k] = $2;
626            $data->{angle1}[$k] = $3;
627            $data->{angle2}[$k] = $4;
628            $data->{angle3}[$k] = $5;
629            next;
630          }
631
632          die "Too many entries in Angles section"
633            if ($i > $data->{nangles});
634          die "Too few entries in Angles section"
635            if ($i < $data->{nangles});
636          die "Multiple angles assigned to the same angle ID"
637            if (scalar @{$data->{anglet}} != $data->{nangles});
638
639          last;
640        }
641      } elsif ($1 eq "Dihedrals") {
642        $data->{dihedralt} = [];
643        $data->{dihedral1} = [];
644        $data->{dihedral2} = [];
645        $data->{dihedral3} = [];
646        $data->{dihedral4} = [];
647        $i = 0;
648        section_check($fh,$data,"Dihedrals");
649
650        while (get_next($fh)) {
651          if (/^\s*([0-9]+)\s+(-?[0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s*$/) {
652
653            $k = $1 - 1;
654            die "Dihedral id $1 is out of range"
655              if (($1 < 1) || ($1 > $data->{ndihedrals}));
656
657            die "Dihedral type $2 is out of range"
658              if (($2 == 0) || ($2 > $data->{ndihedraltypes}));
659
660            die "Dihedral atom 1 ID $3 is out of range"
661              if (($3 < 1) || ($3 > $data->{natoms}));
662            die "Dihedral atom 2 ID $4 is out of range"
663              if (($4 < 1) || ($4 > $data->{natoms}));
664            die "Dihedral atom 3 ID $5 is out of range"
665              if (($5 < 1) || ($5 > $data->{natoms}));
666            die "Dihedral atom 4 ID $6 is out of range"
667              if (($6 < 1) || ($6 > $data->{natoms}));
668
669            ++$i;
670            $data->{dihedralt}[$k] = $2;
671            $data->{dihedral1}[$k] = $3;
672            $data->{dihedral2}[$k] = $4;
673            $data->{dihedral3}[$k] = $5;
674            $data->{dihedral4}[$k] = $6;
675            next;
676          }
677
678          die "Too many entries in Dihedrals section"
679            if ($i > $data->{ndihedrals});
680          die "Too few entries in Dihedrals section"
681            if ($i < $data->{ndihedrals});
682          die "Multiple dihedrals assigned to the same dihedral ID"
683            if (scalar @{$data->{dihedralt}} != $data->{ndihedrals});
684
685          last;
686        }
687      } elsif ($1 eq "Impropers") {
688        $data->{impropert} = [];
689        $data->{improper1} = [];
690        $data->{improper2} = [];
691        $data->{improper3} = [];
692        $data->{improper4} = [];
693        $i = 0;
694        section_check($fh,$data,"Impropers");
695
696        while (get_next($fh)) {
697          if (/^\s*([0-9]+)\s+(-?[0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s*$/) {
698
699            $k = $1 - 1;
700            die "Improper id $1 is out of range"
701              if (($1 < 1) || ($1 > $data->{nimpropers}));
702
703            die "Improper type $2 is out of range"
704              if (($2 == 0) || ($2 > $data->{nimpropertypes}));
705
706            die "Improper atom 1 ID $3 is out of range"
707              if (($3 < 1) || ($3 > $data->{natoms}));
708            die "Improper atom 2 ID $4 is out of range"
709              if (($4 < 1) || ($4 > $data->{natoms}));
710            die "Improper atom 3 ID $5 is out of range"
711              if (($5 < 1) || ($5 > $data->{natoms}));
712            die "Improper atom 4 ID $6 is out of range"
713              if (($6 < 1) || ($6 > $data->{natoms}));
714
715            ++$i;
716            $data->{impropert}[$k] = $2;
717            $data->{improper1}[$k] = $3;
718            $data->{improper2}[$k] = $4;
719            $data->{improper3}[$k] = $5;
720            $data->{improper4}[$k] = $6;
721            next;
722          }
723
724          die "Too many entries in Impropers section"
725            if ($i > $data->{nimpropers});
726          die "Too few entries in Impropers section"
727            if ($i < $data->{nimpropers});
728          die "Multiple impropers assigned to the same improper ID"
729            if (scalar @{$data->{impropert}} != $data->{nimpropers});
730
731          last;
732        }
733      } else {
734        die "Bad data: $_";
735      }
736
737      last unless ($_);
738    } else {
739      die "Bad data: $_";
740    }
741  }
742}
743
744sub floatdiff {
745  my ($n1,$n2,$rel) = @_;
746
747  my $diff = abs($n1-$n2);
748  my $avg = (abs($n1)+abs($n2))*0.5;
749  return 0 if ($avg == 0.0);
750  if ($rel) {
751#    print "relative difference: ",$diff/$avg," vs. $small\n";
752    return 0 if ($diff/$avg < $small);
753  } else {
754#    print "absolute difference: ",$diff," vs. $small\n";
755    return 0 if ($diff < $small);
756  }
757  return 1;
758}
759
760sub coeffcompare {
761  my ($d1,$d2,$coeff,$type) = @_;
762  my (@c1,@c2,$a,$b);
763  my ($field,$count,$i,$j,$t) = ($coeff . 'coeff', 'n' . $type . 'types', 0,0,0);
764
765  if (exists $d1->{$field} && exists $d2->{$field}) {
766    for ($i=0; $i < $d1->{$count}; ++$i) {
767      $t = $i+1;
768      @c1 = split /\s+/, ${$$d1{$field}}[$i];
769      @c2 = split /\s+/, ${$$d2{$field}}[$i];
770      die "Inconsistent number of $coeff coefficients for $type type $t: $#c1 vs $#c2\n"
771        if ($#c1 != $#c2);
772
773      for ($j = 0; $j <= $#c1; ++$j) {
774        $a = $c1[$j]; $b = $c2[$j];
775        die "Inconsistent $coeff coefficient ", $j+1,
776          " for $type type $t: $a vs. $b" if (floatdiff($a,$b,1));
777      }
778    }
779  } else {
780    die "Field $field only exists in data file 1" if (exists $d1->{$field});
781    die "Field $field only exists in data file 2" if (exists $d2->{$field});
782  }
783}
784
785sub topocompare {
786  my ($d1,$d2,$type,$count) = @_;
787  my ($num,$a,$b,$field,$i,$j,$t);
788
789  $field = 'n' . $type . 's';
790  $num = $d1->{$field};
791
792  for ($i=0; $i < $num; ++$i) {
793    $t = $i+1;
794    $field = $type . 't';
795    $a = $d1->{$field}[$i]; $b = $d2->{$field}[$i];
796    die "Inconsistent $type types for $type $t: $a vs. $b" if ($a != $b);
797    for ($j=1; $j <= $count; ++$j) {
798      $field = $type . $j;
799      $a = $d1->{$field}[$i]; $b = $d2->{$field}[$i];
800      die "Inconsistent $type atom $j for $type $t: $a vs. $b" if ($a != $b);
801    }
802  }
803}
804
805sub syscompare {
806  my ($d1,$d2) = @_;
807  my ($i,$j,$t,$a,$b,@l);
808
809  # check atoms.
810  die "Number of atoms does not match"
811    if ($d1->{natoms} != $d2->{natoms});
812  die "Number of atom types does not match"
813    if ($d1->{natomtypes} != $d2->{natomtypes});
814
815  # check bonded interactions
816  @l = ('bond','angle','dihedral','improper');
817  foreach $i (@l) {
818    $t = 'n' . $i . 's';
819    $a = $d1->{$t};
820    $b = $d2->{$t};
821    die "Number of ",$i,"s does not match: $a vs $b" unless ($a == $b);
822
823    $t = 'n' . $i . 'types';
824    $a = $d1->{$t};
825    $b = $d2->{$t};
826    die "Number of ",$i," types does not match: $a vs $b" unless ($a == $b);
827  }
828
829  topocompare($d1,$d2,'bond',2);
830  topocompare($d1,$d2,'angle',3);
831  topocompare($d1,$d2,'dihedral',4);
832  topocompare($d1,$d2,'improper',4);
833
834  coeffcompare($d1,$d2,'pair','atom');
835  coeffcompare($d1,$d2,'bond','bond');
836  coeffcompare($d1,$d2,'angle','angle');
837  coeffcompare($d1,$d2,'bondbond','angle');
838  coeffcompare($d1,$d2,'bondangle','angle');
839  coeffcompare($d1,$d2,'dihedral','dihedral');
840  coeffcompare($d1,$d2,'angleangletorsion','dihedral');
841  coeffcompare($d1,$d2,'bondbond13','dihedral');
842  coeffcompare($d1,$d2,'endbondtorsion','dihedral');
843  coeffcompare($d1,$d2,'middlebondtorsion','dihedral');
844  coeffcompare($d1,$d2,'improper','improper');
845  coeffcompare($d1,$d2,'angleangle','improper');
846
847  for ($i=0; $i < $d1->{natomtypes}; ++$i) {
848    $j = $i+1;
849    if (exists $d1->{mass}[$i]) {
850      $a = $d1->{mass}[$i];
851    } else {
852      die "No mass for atom type $j in data file 1";
853    }
854    if (exists $d2->{mass}[$i]) {
855      $a = $d2->{mass}[$i];
856    } else {
857      die "No mass for atom type $j in data file 2";
858    }
859  }
860
861  # check box information
862  die "Inconsistent box shape" if ($d1->{triclinic} != $d2->{triclinic});
863
864  @l = ('xlo','xhi','ylo','yhi','zlo','zhi');
865  if ($d1->{triclinic}) { push @l, ('xy','xz','yz'); }
866  for $i (@l) {
867    $a = $d1->{$i};
868    $b = $d2->{$i};
869    die "Box data for $i does not match: $a $b" if (floatdiff($a,$b,0));
870  }
871
872  for ($i=0; $i < $d1->{natoms}; ++$i) {
873    $j = $i+1;
874    for $t ('tag','molid','type','imgx','imgy','imgz') {
875      if (exists $d1->{$t}[$i]) {
876        $a = $d1->{$t}[$i];
877      } else {
878        $a = 0;
879      }
880      if (exists $d2->{$t}[$i]) {
881        $b = $d2->{$t}[$i];
882      } else {
883        $b = 0;
884      }
885      die "Inconsistent data for $t, atom $j: $a vs. $b" if ($a != $b);
886    }
887
888    for $t ('charge','posx','posy','posz') {
889      if (exists $d1->{$t}[$i]) {
890        $a = $d1->{$t}[$i];
891      } else {
892        $a = 0;
893      }
894      if (exists $d2->{$t}[$i]) {
895        $b = $d2->{$t}[$i];
896      } else {
897        $b = 0;
898      }
899      die "Inconsistent data for $t, atom $j: $a vs. $b" if (floatdiff($a,$b,0));
900    }
901  }
902
903}
904
905########################################################################
906# main program
907
908my $fp;
909
910if ($#ARGV < 1) {
911  die "usage $0 <file 1> <file 2>";
912}
913
914print "\nLAMMPS data file validation tool. $version\n\n";
915
916data_defaults(\%data1);
917data_defaults(\%data2);
918
919# read in first data file
920open($fp, '<', $ARGV[0]) or die $!;
921print "opened data file 1: $ARGV[0]\n";
922$_=<$fp>;
923print;
924read_data($fp,\%data1);
925print "done reading data file 1\n\n";
926close $fp;
927
928# read in second data file
929open($fp, '<', $ARGV[1]) or die $!;
930print "opened data file 2: $ARGV[1]\n";
931$_=<$fp>;
932print;
933read_data($fp,\%data2);
934print "done reading data file 2\n\n";
935close $fp;
936
937# compare data sets
938syscompare(\%data1,\%data2);
939
940print "File $ARGV[0] and $ARGV[1] match\n\n";
941
942exit 0;
943