forked from Angel-Jia/VASP-script
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathchgdiff.pl
88 lines (69 loc) · 1.76 KB
/
chgdiff.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#!/usr/bin/env perl
#;-*- Perl -*-
@args = @ARGV;
@args == 2 || die "usage: chgdiff.pl <reference CHGCAR> <CHGCAR2>\n";
open (IN1,$args[0]) || die ("Can't open file $!");
open (IN2,$args[1]) || die ("Can't open file $!");
open (OUT,">CHGCAR_diff");
for ($i=0; $i<5; $i++) {
$line1 = <IN1>;
$line2 = <IN2>;
$header1 .= $line1;
}
# check whether it is a vasp5 format
$line1 = <IN1>;
$header1 .= $line1;
$line1 =~ s/^\s+//;
@line1 = split(/\s+/,$line1);
if ($line1[0] =~ /^\d+$/) {
@atoms1 = @line1;
} else {
$atoms1 = <IN1>;
$header1 .= $atoms1;
@atoms1 = split(/\s+/,$atoms1);
}
$line2 = <IN2>;
$line2 =~ s/^\s+//;
@line2 = split(/\s+/,$line2);
if ($line2[0] =~ /^\d+$/) {
@atoms2 = @line2;
} else {
$atoms2 = <IN2>;
@atoms2 = split(/\s+/,$atoms2);
}
$sum1 += $_ for @atoms1;
$sum2 += $_ for @atoms2;
print "Atoms in file1: ".$sum1.", Atoms in file2: ".$sum2."\n";
for ($i=0; $i<$sum1+2; $i++) {
$header1 .= <IN1>;
}
for ($i=0; $i<$sum2+2; $i++) {
$dummy = <IN2>;
}
$points1 = <IN1>;
$header1 .= $points1;
$points2 = <IN2>;
@points1 = split(/\s+/,$points1);
@points2 = split(/\s+/,$points2);
$psum1 = 1;
$psum2 = 1;
for ($i=1; $i<@points1; $i++) {
$psum1 *= $points1[$i];
$psum2 *= $points2[$i];
}
print "Points in file1: ".$psum1.", Points in file2: ".$psum2."\n";
if ($psum1 != $psum2) {die ("Number of points not same in two files!");}
print OUT $header1;
for ($i=0; $i<$psum1/5; $i++) {
$line1 = <IN1>;
$line2 = <IN2>;
@line1 = split(/\s+/,$line1);
@line2 = split(/\s+/,$line2);
for ($j=1; $j<@line1; $j++) {
$line1[$j] = $line2[$j]-$line1[$j];
}
printf OUT "%18.11E %18.11E %18.11E %18.11E %18.11E\n",$line1[1],$line1[2],$line1[3],$line1[4],$line1[5];
}
close(OUT);
close(IN2);
close(IN1);