mopac_parser.py

import re
import os

class MOPACParser:
    """
    A parser for MOPAC output files (.out) to extract key data such as energy, charges, dipole moment, etc.
    """

    def __init__(self, file_path):
        """
        Initialize the parser with the path to the MOPAC output file.

        :param file_path: Path to the .out file
        """
        self.file_path = file_path
        self.content = None
        self.data = {}

    def read_file(self):
        """
        Read the content of the file.
        """
        if not os.path.exists(self.file_path):
            raise FileNotFoundError(f"File {self.file_path} not found.")
        with open(self.file_path, 'r') as f:
            self.content = f.read()

    def extract_energy(self):
        """
        Extract the final heat of formation or total energy.
        """
        # Look for "FINAL HEAT OF FORMATION =" or "TOTAL ENERGY ="
        match = re.search(r'FINAL HEAT OF FORMATION\s*=\s*([-\d.]+)', self.content, re.IGNORECASE)
        if not match:
            match = re.search(r'TOTAL ENERGY\s*=\s*([-\d.]+)', self.content, re.IGNORECASE)
        if match:
            self.data['energy'] = float(match.group(1))
        else:
            self.data['energy'] = None

    def extract_charges(self):
        """
        Extract atomic charges from the NET ATOMIC CHARGES section.
        """
        charges = []
        # Find the start of the charges section
        start_match = re.search(r'NET ATOMIC CHARGES', self.content, re.IGNORECASE)
        if start_match:
            # Extract lines after the header
            lines = self.content[start_match.end():].split('\n')
            for line in lines:
                line = line.strip()
                if not line:
                    continue
                parts = line.split()
                if len(parts) == 6 and parts[0].isdigit():
                    try:
                        charge = float(parts[2])
                        charges.append({'atom': int(parts[0]), 'symbol': parts[1], 'charge': charge})
                    except ValueError:
                        continue
                elif len(parts) > 0 and not parts[0].isdigit():
                    continue  # skip headers
                else:
                    break
        self.data['charges'] = charges

    def extract_dipole_moment(self):
        """
        Extract the dipole moment components and total.
        """
        # Look for "DIPOLE" or "ELECTRIC DIPOLE MOMENT"
        match = re.search(r'DIPOLE\s+X\s+([-\d.]+)\s+Y\s+([-\d.]+)\s+Z\s+([-\d.]+)', self.content, re.IGNORECASE)
        if match:
            dx = float(match.group(1))
            dy = float(match.group(2))
            dz = float(match.group(3))
            total = (dx**2 + dy**2 + dz**2)**0.5
            self.data['dipole'] = {'x': dx, 'y': dy, 'z': dz, 'total': total}
        else:
            self.data['dipole'] = None

    def extract_geometry(self):
        """
        Extract the last geometry if present.
        """
        geometry = []
        # Try multiple patterns for geometry sections
        patterns = [
            r'CARTESIAN COORDINATES',
            r'FINAL GEOMETRY',
            r'OPTIMIZED GEOMETRY',
            r'GEOMETRY AT END OF JOB'
        ]
        start_pos = -1
        for pattern in patterns:
            matches = list(re.finditer(pattern, self.content, re.IGNORECASE))
            if matches:
                start_pos = max(start_pos, matches[-1].end())  # Take the last occurrence

        if start_pos != -1:
            lines = self.content[start_pos:].split('\n')
            for line in lines:
                line = line.strip()
                if not line:
                    continue
                parts = line.split()
                if len(parts) == 5 and parts[0].isdigit():
                    try:
                        atom_num = int(parts[0])
                        symbol = parts[1]
                        x = float(parts[2])
                        y = float(parts[3])
                        z = float(parts[4]) if len(parts) > 4 else 0.0
                        geometry.append({'atom': atom_num, 'symbol': symbol, 'x': x, 'y': y, 'z': z})
                    except ValueError:
                        continue
                elif len(parts) > 0 and not parts[0].isdigit():
                    continue  # skip headers like "NO. ATOM X Y Z"
                else:
                    break
        self.data['geometry'] = geometry

    def extract_molecular_formula(self):
        """
        Extract the molecular formula from the Empirical Formula line or geometry.
        """
        match = re.search(r'Empirical Formula:\s*(.*)', self.content, re.IGNORECASE)
        if match:
            self.data['formula'] = match.group(1).strip()
        else:
            # Fallback to deriving from geometry
            if not self.data.get('geometry'):
                self.data['formula'] = None
                return
            from collections import Counter
            symbols = [atom['symbol'] for atom in self.data['geometry']]
            counts = Counter(symbols)
            # Sort by atomic number or alphabetically
            formula = ''.join(f"{sym}{count if count > 1 else ''}" for sym, count in sorted(counts.items()))
            self.data['formula'] = formula

    def parse(self):
        """
        Parse the file and extract all data.
        """
        self.read_file()
        self.extract_energy()
        self.extract_charges()
        self.extract_dipole_moment()
        self.extract_geometry()
        self.extract_molecular_formula()
        return self.data

import json

# Main script to process all .out files in the current directory
if __name__ == "__main__":
    current_dir = '.'
    out_files = [f for f in os.listdir(current_dir) if f.endswith('.out')]

    if not out_files:
        print("No .out files found in the current directory.")
    else:
        results = {}
        for file in out_files:
            parser = MOPACParser(file)
            try:
                data = parser.parse()
                results[file] = data
                print(f"Parsed {file}: Formula = {data.get('formula')}, Energy = {data.get('energy')}, Dipole total = {data.get('dipole', {}).get('total') if data.get('dipole') else 'N/A'}")

                # Save geometry to GJF file if available
                if data.get('geometry'):
                    gjf_filename = file.replace('.out', '.gjf')
                    with open(gjf_filename, 'w') as f:
                        f.write(f"%chk={file.replace('.out', '.chk')}\n")
                        f.write("# pm7\n")
                        f.write("\n")
                        f.write(f"Geometry from {file}\n")
                        f.write("\n")
                        f.write("0 1\n")
                        for atom in data['geometry']:
                            f.write(f"{atom['symbol']} {atom['x']:.6f} {atom['y']:.6f} {atom['z']:.6f}\n")
                        f.write("\n")
                    print(f"Geometry saved to {gjf_filename}")
                    # Remove geometry from JSON data
                    del data['geometry']

            except Exception as e:
                print(f"Error parsing {file}: {e}")

        # Save results to JSON file
        with open('parsed_results.json', 'w') as f:
            json.dump(results, f, indent=4)
        print("All results saved to parsed_results.json")