|
| 1 | +#!/usr/bin/env python3 |
| 2 | +import sys |
| 3 | + |
| 4 | +def main(): |
| 5 | + if len(sys.argv) < 2: |
| 6 | + print(f"Can't find {sys.argv[1] if len(sys.argv) > 1 else 'file'} !") |
| 7 | + sys.exit(1) |
| 8 | + |
| 9 | + input_file = sys.argv[1] |
| 10 | + |
| 11 | + try: |
| 12 | + with open(input_file, 'r') as inp: |
| 13 | + # skip the first two lines |
| 14 | + inp.readline() |
| 15 | + inp.readline() |
| 16 | + |
| 17 | + # read the 3rd line: number of atoms + origin coordinates |
| 18 | + line = inp.readline().split() |
| 19 | + if not line: |
| 20 | + return |
| 21 | + natom = int(line[0]) |
| 22 | + # origin_x = float(line[1]) |
| 23 | + # origin_y = float(line[2]) |
| 24 | + # origin_z = float(line[3]) |
| 25 | + |
| 26 | + # read the grid vectors (support non-orthogonal) |
| 27 | + line = inp.readline().split() |
| 28 | + nx = int(line[0]) |
| 29 | + v1 = [float(line[1]), float(line[2]), float(line[3])] |
| 30 | + |
| 31 | + line = inp.readline().split() |
| 32 | + ny = int(line[0]) |
| 33 | + v2 = [float(line[1]), float(line[2]), float(line[3])] |
| 34 | + |
| 35 | + line = inp.readline().split() |
| 36 | + nz = int(line[0]) |
| 37 | + v3 = [float(line[1]), float(line[2]), float(line[3])] |
| 38 | + |
| 39 | + # calculate the volume element |v1 · (v2 × v3)| |
| 40 | + val0 = v2[1] * v3[2] - v2[2] * v3[1] |
| 41 | + val1 = v2[0] * v3[2] - v2[2] * v3[0] |
| 42 | + val2 = v2[0] * v3[1] - v2[1] * v3[0] |
| 43 | + |
| 44 | + volume = abs(v1[0] * val0 - v1[1] * val1 + v1[2] * val2) |
| 45 | + |
| 46 | + # skip the atom coordinates |
| 47 | + # natom can be negative in cube files sometimes? |
| 48 | + # C++ code: for (int i = 0; i < natom; ++i) |
| 49 | + # If natom is negative, loop doesn't run. |
| 50 | + |
| 51 | + atoms_to_skip = natom |
| 52 | + if atoms_to_skip < 0: |
| 53 | + # In some cube formats, negative natom means second line of header contains units or something? |
| 54 | + # Standard: "If the number of atoms is negative, that indicates that the file contains input lines... One line for each non-zero E value." |
| 55 | + # But the C++ code simple loops < natom. If natom < 0, it skips 0 lines. |
| 56 | + # We will mimic C++ behavior assuming it works for the files they have. |
| 57 | + atoms_to_skip = 0 # Loop won't run if natom < 0 |
| 58 | + |
| 59 | + for _ in range(atoms_to_skip): |
| 60 | + inp.readline() |
| 61 | + |
| 62 | + nr = nx * ny * nz |
| 63 | + |
| 64 | + total_sum = 0.0 |
| 65 | + count = 0 |
| 66 | + |
| 67 | + # Read grid values |
| 68 | + # iterate over remaining lines to handle values spread across lines |
| 69 | + for line in inp: |
| 70 | + parts = line.split() |
| 71 | + for part in parts: |
| 72 | + total_sum += float(part) |
| 73 | + count += 1 |
| 74 | + if count >= nr: |
| 75 | + break |
| 76 | + if count >= nr: |
| 77 | + break |
| 78 | + |
| 79 | + ne = total_sum * volume |
| 80 | + # cout default precision is 6, but setprecision(10) changes it. |
| 81 | + # Python's default float printing is usually sufficient, but let's use formatting to be sure. |
| 82 | + # {:.10g} prints up to 10 significant digits. |
| 83 | + print(f"{ne:.10g}") |
| 84 | + |
| 85 | + except FileNotFoundError: |
| 86 | + print(f"Can't find {input_file} !") |
| 87 | + sys.exit(1) |
| 88 | + except Exception as e: |
| 89 | + # C++ doesn't print other errors generally, but let's be safe. |
| 90 | + # But to be functionally equivalent, maybe we shouldn't. |
| 91 | + # The C++ code crashes or behaves weirdly on bad input. |
| 92 | + # Let's just exit 1. |
| 93 | + sys.exit(1) |
| 94 | + |
| 95 | +if __name__ == "__main__": |
| 96 | + main() |
0 commit comments