source: util/src/BOSSMatcher.pl.in@ 0ab37b

#!@PERL@ -w
#
# In the course of numerous fragment calculations with breaks and interrupts in between, there might
# arise some hickhack and mess with the sequence of the calculations. BOSSMatcher compares the nuclear
# coordinates given in the config files of each fragment with those written out in the forces file of
# ESPACK and prints which fragments do not match with their config file coordinates.

use strict;

use Getopt::Long;
use File::Copy;

my ($Help, $Prefix, $Inputdir);

GetOptions("help"    => \$Help,
        "prefix=s",             =>      \$Prefix, 
        "inputdir=s"    =>      \$Inputdir);

if ($Help) { usage(); }

sub usage {
    print STDERR <<"EOF";
Usage: $0 [OPTIONS]
 compares the nuclear
 coordinates given in the config files of each fragment with those written out in the forces file of
 ESPACK and prints which fragments do not match with their config file coordinates.
  --prefix <prefix> prefix in front of .forces.all and .energy.all
  --inputdir <inputdir> read data from this directory
  --help                        this help
EOF
    exit 1;
}

if ((!defined $Inputdir) || (!defined $Prefix)) {
        print "Inputdir and/or Prefix not specified! See --help.\n";
        exit 1;
}

# 1. find the digits of the fragment number
print "Finding the correct number of digits.\n";
my $digits="";
my $digit=0;
my $file;
do {
        $digits=$digits."0";    # add a zero
        $digit++;
        $file=$Inputdir."/BondFragment".$digits.".conf";
        print "Looking for $file ... \n";
} while ((! -e $file) && ($digit < 10));

if (!($digit < 10)) {
        print "In the directory specified no initial fragment config file was found, typo?\n";
        exit 1;
}

# 2. loop over all fragments and parse config and forces file, comparing the coordinates

my $number=0;
my $failures=0;
my $failflag=0;
my @lines;
my $i;
my @configcoordinates;
my @forcecoordinates;
my @line;
my $IsAngstroem=0;
do {
        @configcoordinates = ();
        @forcecoordinates = ();
        $failflag=0;
        $file = sprintf("%s%0${digit}d%s", $Inputdir."/BondFragment",$number,".conf");
        if (-e $file) {
                print "Scanning $file ... ";
                # 2a. parse the config file
                open(CONFIG, $file);
                @lines = <CONFIG>;
                for ($i=0; $i < scalar(@lines); $i++)   {
                        if (${lines[$i]}  =~ /Ion_Type._/) {
                                @line = split(/[ \t]+/,${lines[$i]});
                                my @coord = (${line[1]}, ${line[2]}, ${line[3]});
                                push(@configcoordinates, \@coord);
                        }
                        if (${lines[$i]}  =~ /IsAngstroem/) {
                                @line = split(/[ \t]+/,${lines[$i]});
                                if (${line[1]} == 0) {
                                        $IsAngstroem = 1;
                                } else {
                                        $IsAngstroem = 1.8897261;
                                }
                                #print "Found IsAngstroem ".${line[1]}.", setting to $IsAngstroem.\n";
                        }
                }
                close(CONFIG);
                # 2b. parse the pcp.forces.file
                $file = sprintf("%s%0${digit}d%s", $Inputdir."/BondFragment",$number,"/${Prefix}.forces.all");
                if (-e $file) {
                        open(FORCES, $file);
                        @lines = <FORCES>;              # read file into buffer
                        for ($i=1; $i < (scalar(@lines)-1); $i++)       {               # skip first and last line
                                @line = split(/[ \t]+/,${lines[$i]});
                                my @coord = (${line[2]}, ${line[3]}, ${line[4]});
                                push(@forcecoordinates, \@coord);
                        }
                        close(FORCES);
                } else {
                        $failflag = 1;
                }
                # 2c. compare both coordinates and sum up the deviation
                if (scalar(@configcoordinates) != scalar(@forcecoordinates)     ) {
                        $failflag = 1;
                        print "mismatch in number of ions ".scalar(@configcoordinates)." vs. ".scalar(@forcecoordinates)." ... ";
                } else {
                        #print "match in number of ions ".scalar(@configcoordinates)." vs. ".scalar(@forcecoordinates)." ... ";
                }
                my $forces;
                my $configs;
                my $error = 0;
                do {
                        $forces = pop(@forcecoordinates);
                        $configs = pop(@configcoordinates);
                        my $currenterror = 0;
                        for ($i=0;$i<3;$i++) {
                                @$configs[$i] *= $IsAngstroem;
                                $currenterror += (@$configs[$i]-@$forces[$i])*(@$configs[$i]-@$forces[$i]);
                        }
                        $error += sqrt($currenterror);
                        #print "$currenterror ";
                } while (scalar(@configcoordinates) > 0);
                print "deviation: ";
                if ($error > 1e-4) {
                        $failflag=1;
                        print "$error ";
                } else {
                        printf "<1e-4 ";
                }

                # 2d. print final status
                if ($failflag) {
                        $failures++;
                        print "failed.\n";      
                } else {
                        print "done.\n";        
                }
                $number++;
        } else {
                print "Last config file nr. $number reached, done.\n";
        }
} while (-e $file);

# 3. finish
print "Done checking, I noticed $failures failure(s).\n";
exit 0