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