utils.py

import re
import math

# ===========================================================================
# NBO E(2) Parser
# ===========================================================================

_DONOR_RE = (
    r"\s*(\d+)\.\s+"
    r"(BD\*?|LP\*?|CR|RY\*?)\s*"
    r"\(\s*(\d+)\)\s+"
    r"([A-Za-z]+)\s+(\d+)"
    r"(?:\s*-\s*([A-Za-z]+)\s+(\d+))?"
)
_ACCEPTOR_RE = (
    r"\s*/\s*(\d+)\.\s+"
    r"(BD\*?|LP\*?|CR|RY\*?)\s*"
    r"\(\s*(\d+)\)\s+"
    r"([A-Za-z]+)\s+(\d+)"
    r"(?:\s*-\s*([A-Za-z]+)\s+(\d+))?"
)
_VALUES_RE = r"\s+([\d.]+)\s+([\d.]+)\s+([\d.]+)\s*$"
_E2_LINE_RE = re.compile(_DONOR_RE + _ACCEPTOR_RE + _VALUES_RE)
_NBO_SUMMARY_RE = re.compile(
    r"\s*(\d+)\.\s+"
    r"(BD\*?|LP\*?|CR|RY\*?)\s*"
    r"\(\s*(\d+)\)\s+"
    r"([A-Za-z]+)\s+(\d+)"
    r"(?:\s*-\s*([A-Za-z]+)\s+(\d+))?"
    r"\s+([\d.]+)"
    r"\s+(-?[\d.]+)"
)

def parse_e2_section(lines):
    records = []
    in_section = False
    section_label = ""
    for line in lines:
        if "Second Order Perturbation Theory Analysis" in line:
            in_section = True; continue
        if not in_section: continue
        if "Natural Bond Orbitals (Summary)" in line: break
        stripped = line.strip()
        if stripped.startswith("within unit") or stripped.startswith("from unit"):
            section_label = stripped; continue
        m = _E2_LINE_RE.search(line)
        if m is None: continue
        g = m.groups()
        d_atoms = f"{g[3]} {g[4]}" + (f" - {g[5]} {g[6]}" if g[5] else "")
        a_atoms = f"{g[10]} {g[11]}" + (f" - {g[12]} {g[13]}" if g[12] else "")
        records.append({
            "section": section_label,
            "donor_idx": int(g[0]), "donor_type": g[1], "donor_type_idx": int(g[2]),
            "donor_atoms": d_atoms,
            "acceptor_idx": int(g[7]), "acceptor_type": g[8], "acceptor_type_idx": int(g[9]),
            "acceptor_atoms": a_atoms,
            "E2_kcal_mol": float(g[14]), "Ej_Ei_au": float(g[15]), "Fij_au": float(g[16]),
        })
    return records

def parse_nbo_summary(lines):
    nbo_map = {}
    in_section = False
    for line in lines:
        if "Natural Bond Orbitals (Summary)" in line: in_section = True; continue
        if not in_section: continue
        s = line.strip()
        if s.startswith("NHO Directionality") or s.startswith("NBO state"): break
        if s.startswith("$END") or s.startswith("Effective longest"): break
        m = _NBO_SUMMARY_RE.search(line)
        if m is None: continue
        idx = int(m.group(1))
        nbo_map[idx] = {"occupancy": float(m.group(8)), "energy_au": float(m.group(9))}
    return nbo_map

def merge_nbo_data(e2_records, nbo_map):
    rows = []
    for rec in e2_records:
        d = nbo_map.get(rec["donor_idx"], {})
        a = nbo_map.get(rec["acceptor_idx"], {})
        rows.append({**rec,
            "donor_occupancy": d.get("occupancy"),
            "donor_orbital_energy_au": d.get("energy_au"),
            "acceptor_occupancy": a.get("occupancy"),
            "acceptor_orbital_energy_au": a.get("energy_au"),
        })
    return rows

def _find_halogen_atoms(rows, halogen):
    atoms = set()
    pat = re.compile(rf"\b{re.escape(halogen)}\s+\d+")
    for r in rows:
        for col in ("donor_atoms", "acceptor_atoms"):
            for m in pat.finditer(r[col]):
                atoms.add(m.group())
    return atoms

def _auto_detect_halogen(rows):
    for sym in ("I", "Br", "Cl", "At"):
        if _find_halogen_atoms(rows, sym): return sym
    return None

def halogen_bond_analysis(rows, halogen=None):
    if halogen is None:
        halogen = _auto_detect_halogen(rows)
        if halogen is None:
            return {"error": "ハロゲン原子 (I/Br/Cl) が見つかりません。"}
    X = halogen
    x_atoms = _find_halogen_atoms(rows, X)
    if not x_atoms:
        return {"error": f"{X} 原子が見つかりません。"}

    inter = [r for r in rows if not r["section"].startswith("within unit")]
    bdstar = set()
    for r in inter:
        if r["acceptor_type"] == "BD*" and any(xa in r["acceptor_atoms"] for xa in x_atoms):
            bdstar.add(r["acceptor_idx"])
        if r["donor_type"] == "BD*" and any(xa in r["donor_atoms"] for xa in x_atoms):
            bdstar.add(r["donor_idx"])
    ry_lp = set()
    for r in inter:
        if r["acceptor_type"] in ("RY*", "LP*") and r["acceptor_atoms"] in x_atoms:
            ry_lp.add(r["acceptor_idx"])

    cats = []
    global_element_stats = {}
    for label, title, filt in [
        (f"1_All→BD*(C-{X})", f"ハロゲン結合: All → BD*(C-{X})",
         lambda r: r["acceptor_idx"] in bdstar),
        (f"2_All→RY*/LP*({X})", f"d軌道関与: All → RY*/LP*({X})",
         lambda r: r["acceptor_idx"] in ry_lp),
        (f"3_LP({X})→Lewis_base", f"バックドネーション: LP({X}) → Lewis base",
         lambda r: r["donor_type"] == "LP"
                   and any(xa in r["donor_atoms"] for xa in x_atoms)
                   and r["acceptor_idx"] not in bdstar
                   and r["acceptor_idx"] not in ry_lp),
    ]:
        cat_rows = sorted([{**r, "category": label} for r in inter if filt(r)],
                          key=lambda r: r["E2_kcal_mol"], reverse=True)
        
        cat_stats = {}
        for r in cat_rows:
            target = r["acceptor_atoms"] if label.startswith("3_") else r["donor_atoms"]
            parts = target.split()
            if len(parts) >= 3 and parts[1] == "-":
                lb_elem = f"{parts[0]}-{parts[2]}"
            elif len(parts) >= 1:
                lb_elem = parts[0]
            else:
                lb_elem = "Unknown"
            
            r["lewis_base_element"] = lb_elem
            
            for s_dict in (cat_stats, global_element_stats):
                if lb_elem not in s_dict:
                    s_dict[lb_elem] = {"count": 0, "sum_e2": 0.0, "max_e2": 0.0}
                s_dict[lb_elem]["count"] += 1
                s_dict[lb_elem]["sum_e2"] += r["E2_kcal_mol"]
                s_dict[lb_elem]["max_e2"] = max(s_dict[lb_elem]["max_e2"], r["E2_kcal_mol"])
        
        stat_list = []
        for elem, s in sorted(cat_stats.items(), key=lambda x: -x[1]["count"]):
            stat_list.append({
                "element": elem,
                "count": s["count"],
                "avg_e2": round(s["sum_e2"] / s["count"], 2),
                "max_e2": round(s["max_e2"], 2)
            })

        secs = list(set(r["section"] for r in cat_rows))
        cats.append({"label": label, "title": title, "count": len(cat_rows),
                     "sections": " / ".join(sorted(secs)), "rows": cat_rows,
                     "stats": stat_list})

    global_stat_list = []
    for elem, s in sorted(global_element_stats.items(), key=lambda x: -x[1]["count"]):
        global_stat_list.append({
            "element": elem,
            "count": s["count"],
            "avg_e2": round(s["sum_e2"] / s["count"], 2),
            "max_e2": round(s["max_e2"], 2)
        })

    return {"halogen": X, "x_atoms": sorted(x_atoms),
            "total": sum(c["count"] for c in cats), "categories": cats,
            "global_stats": global_stat_list}


# ===========================================================================
# Gaussian Complete Extractor
# ===========================================================================

def _safe_floats(text):
    return [float(x) for x in re.findall(r"[-+]?\d*\.\d+", text)]

def _dist(a1, a2):
    return math.sqrt((a1['x']-a2['x'])**2 + (a1['y']-a2['y'])**2 + (a1['z']-a2['z'])**2)

def parse_gaussian_complete(text, file_name="unknown"):
    lines = text.splitlines()
    results = []
    state = {
        "step": 0, "method": "Unknown", "basis": "Unknown", "solvent": "None",
        "scf_energy": None, "homo": None, "lumo": None, "dipole": None,
        "s2_gs": None, "geometry": [], "nbo_charges": [], "mulliken_charges": [],
        "imaginary_freqs": 0, "zpe": None, "therm_g": None, "sum_g": None,
        "td_data": {}, "is_td": False,
    }

    reading_geo = reading_mulliken = reading_nbo = next_line_dipole = False
    reading_route = False
    route_buffer = []
    temp_occ, temp_virt = [], []
    current_td_state = None

    # Patterns
    route_start = re.compile(r"^\s*#\s?p?")
    dash_line = re.compile(r"^\s*-{5,}")
    solvent_body = re.compile(r"Solvent\s*:\s+([A-Za-z0-9\-\_]+)", re.I)
    scf_pat = re.compile(r"SCF Done:\s+E\(\w+\)\s+=\s+([-+]?\d*\.\d+)")
    freq_pat = re.compile(r"Frequencies\s+--\s+(.+)")
    zp_pat = re.compile(r"Zero-point correction=\s+([-.\d]+)")
    therm_g_pat = re.compile(r"Thermal correction to Gibbs Free Energy=\s+([-.\d]+)")
    sum_g_pat = re.compile(r"Sum of electronic and thermal Free Energies=\s+([-.\d]+)")
    dipole_head = re.compile(r"Dipole moment")
    dipole_val = re.compile(r"Tot=\s+([\d\.]+)")
    mul_start = re.compile(r"Mulliken charges:")
    mul_row = re.compile(r"^\s+(\d+)\s+([A-Z][a-z]?)\s+([-.\d]+)")
    nbo_start = re.compile(r"Summary of Natural Population Analysis")
    nbo_row = re.compile(r"^\s+([A-Z][a-z]?)\s+(\d+)\s+([-.\d]+)")
    occ_pat = re.compile(r"Alpha\s+occ\.\s+eigenvalues\s+--\s+(.+)")
    virt_pat = re.compile(r"Alpha\s+virt\.\s+eigenvalues\s+--\s+(.+)")
    geo_start = re.compile(r"(Standard|Input)\s+orientation:")
    geo_row = re.compile(r"^\s+(\d+)\s+(\d+)\s+\d+\s+([-.\d]+)\s+([-.\d]+)\s+([-.\d]+)")
    geo_end = re.compile(r"^\s*-{5,}")
    exc_pat = re.compile(r"Excited State\s+(\d+):.*?\s+([\d\.]+)\s+eV\s+([\d\.]+)\s+nm\s+f=([\d\.]+)")
    exc_s2 = re.compile(r"<S\*\*2>=([-.\d]+)")
    trans_pat = re.compile(r"^\s+(\d+)\s+->\s+(\d+)\s+([-.\d]+)")
    td_total = re.compile(r"Total Energy, E\(TD.*?\)\s*=\s*([-+]?\d+\.\d+)")
    term_pat = re.compile(r"Normal termination of Gaussian")

    for line in lines:
        if route_start.search(line):
            reading_route = True; route_buffer = []
            state["is_td"] = False; state["td_data"] = {}
            temp_occ, temp_virt = [], []
            state["imaginary_freqs"] = 0
            state["method"] = "Unknown"; state["basis"] = "Unknown"; state["solvent"] = "None"
        if reading_route:
            if dash_line.search(line):
                reading_route = False
                full = " ".join(route_buffer).lower()
                if 'td' in full: state["is_td"] = True
                for t in full.split():
                    if '/' in t:
                        try:
                            m, b = t.split('/')[:2]
                            state["method"] = m.strip(); state["basis"] = b.strip()
                        except: pass
                if 'solvent' in full:
                    ms = re.search(r"solvent\s*=\s*([a-z0-9\-_]+)", full)
                    state["solvent"] = ms.group(1) if ms else "Yes(Check)"
            else:
                c = line.strip()
                if c: route_buffer.append(c)
        if state["solvent"] in ("None", "Unknown", "Yes(Check)"):
            ms = solvent_body.search(line)
            if ms: state["solvent"] = ms.group(1)
        if geo_start.search(line):
            reading_geo = True; state["geometry"] = []; continue
        if reading_geo:
            if geo_end.search(line) and state["geometry"]:
                reading_geo = False
            else:
                m = geo_row.search(line)
                if m:
                    state["geometry"].append({
                        'id': int(m.group(1)), 'atom': int(m.group(2)),
                        'x': float(m.group(3)), 'y': float(m.group(4)), 'z': float(m.group(5))
                    })
        ms = scf_pat.search(line)
        if ms: state["scf_energy"] = float(ms.group(1))
        ms = freq_pat.search(line)
        if ms:
            for f in _safe_floats(ms.group(1)):
                if f < 0: state["imaginary_freqs"] += 1
        ms = zp_pat.search(line)
        if ms: state["zpe"] = float(ms.group(1))
        ms = therm_g_pat.search(line)
        if ms: state["therm_g"] = float(ms.group(1))
        ms = sum_g_pat.search(line)
        if ms: state["sum_g"] = float(ms.group(1))
        ms = occ_pat.search(line)
        if ms: temp_occ.extend(_safe_floats(ms.group(1)))
        ms = virt_pat.search(line)
        if ms: temp_virt.extend(_safe_floats(ms.group(1)))
        if dipole_head.search(line):
            ms = dipole_val.search(line)
            if ms: state["dipole"] = float(ms.group(1))
            else: next_line_dipole = True
        elif next_line_dipole:
            ms = dipole_val.search(line)
            if ms: state["dipole"] = float(ms.group(1)); next_line_dipole = False
            elif "Quadrupole" in line: next_line_dipole = False
        if "<S**2>" in line and "Excited State" not in line:
            try: state["s2_gs"] = float(line.split("=")[1].split()[0])
            except: pass
        if mul_start.search(line):
            reading_mulliken = True; state["mulliken_charges"] = []; continue
        if reading_mulliken:
            if "Sum of Mulliken" in line: reading_mulliken = False
            else:
                m = mul_row.search(line)
                if m and m.group(2) in ("I", "Br", "Cl"):
                    state["mulliken_charges"].append(f"{m.group(2)}{m.group(1)}:{float(m.group(3)):.3f}")
        if nbo_start.search(line):
            reading_nbo = True; state["nbo_charges"] = []; continue
        if reading_nbo:
            if "=======" in line: reading_nbo = False
            else:
                m = nbo_row.search(line)
                if m and m.group(1) in ("I", "Br", "Cl"):
                    state["nbo_charges"].append(f"{m.group(1)}{m.group(2)}:{float(m.group(3)):.3f}")
        if state["is_td"]:
            ms = exc_pat.search(line)
            if ms:
                sid = int(ms.group(1)); current_td_state = sid
                if sid not in state["td_data"]: state["td_data"][sid] = {"trans": []}
                es2 = None
                ms2 = exc_s2.search(line)
                if ms2: es2 = float(ms2.group(1))
                state["td_data"][sid].update({"ev": float(ms.group(2)), "nm": float(ms.group(3)),
                                              "f": float(ms.group(4)), "s2": es2})
            ms = trans_pat.search(line)
            if ms and current_td_state:
                state["td_data"][current_td_state]["trans"].append(
                    f"{ms.group(1)}->{ms.group(2)}({ms.group(3)})")
            ms = td_total.search(line)
            if ms and current_td_state:
                state["td_data"][current_td_state]["total_e"] = float(ms.group(1))
        
        if term_pat.search(line):
            state["step"] += 1
            if temp_occ: state["homo"] = temp_occ[-1]
            if temp_virt: state["lumo"] = temp_virt[0]
            homo_ev = state["homo"] * 27.2114 if state["homo"] else None
            lumo_ev = state["lumo"] * 27.2114 if state["lumo"] else None
            gap_ev = (lumo_ev - homo_ev) if (homo_ev and lumo_ev) else None

            yx_d, xc_d, cy_d = [], [], []
            if state["geometry"]:
                Xs = [a for a in state["geometry"] if a['atom'] in (53, 35, 17)]
                Ys = [a for a in state["geometry"] if a['atom'] in (8, 16, 34)]
                Cs = [a for a in state["geometry"] if a['atom'] == 6]
                
                atom_sym = {6: 'C', 8: 'O', 16: 'S', 17: 'Cl', 34: 'Se', 35: 'Br', 53: 'I'}
                
                for x in Xs:
                    xsym = atom_sym.get(x['atom'], 'X')
                    ds = sorted([(round(_dist(x, y), 3), y['id'], atom_sym.get(y['atom'], 'Y')) for y in Ys if _dist(x, y) < 3.5])
                    if ds: yx_d.append("; ".join([f"{d[2]}{d[1]}-{xsym}{x['id']}({d[0]})" for d in ds]))
                    dsc = sorted([(round(_dist(x, c), 3), c['id']) for c in Cs if _dist(x, c) < 2.5])
                    if dsc: xc_d.append("; ".join([f"C{d[1]}-{xsym}{x['id']}({d[0]})" for d in dsc]))
                
                for c in Cs:
                    for y in Ys:
                        d = _dist(c, y)
                        ysym = atom_sym.get(y['atom'], 'Y')
                        if d < 2.5: cy_d.append(f"C{c['id']}-{ysym}{y['id']}({d:.3f})")

            base = {
                "Filename": file_name, "Step": state["step"],
                "Method": state["method"], "Basis": state["basis"], "CPCM": state["solvent"],
                "SCF Energy (Ha)": state["scf_energy"],
                "ZPE Corr (Ha)": state["zpe"], "Therm Corr to G (Ha)": state["therm_g"],
                "Sum Elec + G (Ha)": state["sum_g"],
                "Dipole (Debye)": state["dipole"], "GS <S**2>": state["s2_gs"],
                "Imaginary Freqs": state["imaginary_freqs"],
                "HOMO (eV)": round(homo_ev, 4) if homo_ev else None,
                "LUMO (eV)": round(lumo_ev, 4) if lumo_ev else None,
                "Gap (eV)": round(gap_ev, 4) if gap_ev else None,
                "NBO Charge": " / ".join(state["nbo_charges"]),
                "Mulliken Charge": " / ".join(state["mulliken_charges"]),
                "C-X Dist (<2.5A)": " / ".join(xc_d),
                "Y-X Dist (<3.5A)": " / ".join(yx_d),
                "C-Y Dist (<2.5A)": " / ".join(cy_d) if cy_d else "",
            }
            if state["is_td"] and state["td_data"]:
                for sid, data in state["td_data"].items():
                    row = {**base, "State": sid,
                           "Exc. Energy (eV)": data.get("ev"),
                           "Wavelength (nm)": data.get("nm"),
                           "f": data.get("f"), "ES <S**2>": data.get("s2"),
                           "S1 Total Energy (Ha)": data.get("total_e"),
                           "Transitions": "; ".join(data.get("trans", []))}
                    results.append(row)
            else:
                results.append(base)
    return results

# ===========================================================================
# XYZ Coordinates Extractor
# ===========================================================================

def parse_gaussian_xyz(text, file_name="unknown"):
    lines = text.splitlines()
    atom_sym = {1: 'H', 5: 'B', 6: 'C', 7: 'N', 8: 'O', 9: 'F', 14: 'Si', 15: 'P', 16: 'S', 17: 'Cl', 34: 'Se', 35: 'Br', 53: 'I'}
    
    geo_start = re.compile(r"(Standard|Input)\s+orientation:")
    geo_row = re.compile(r"^\s+(\d+)\s+(\d+)\s+\d+\s+([-.\d]+)\s+([-.\d]+)\s+([-.\d]+)")
    geo_end = re.compile(r"^\s*-{5,}")
    
    best_geo = []
    reading_geo = False
    temp_geo = []

    for line in lines:
        if geo_start.search(line):
            reading_geo = True
            temp_geo = []
            continue
        
        if reading_geo:
            if geo_end.search(line) and len(temp_geo) > 0:
                reading_geo = False
                best_geo = temp_geo
                continue
            
            m = geo_row.search(line)
            if m:
                atom_type = int(m.group(2))
                x, y, z = m.group(3), m.group(4), m.group(5)
                sym = atom_sym.get(atom_type, str(atom_type))
                temp_geo.append(f"{sym:<2} {x:>12} {y:>12} {z:>12}")

    if not best_geo:
        return None
        
    xyz_text = f"{len(best_geo)}\n{file_name}\n" + "\n".join(best_geo)
    return {"filename": file_name, "atoms": len(best_geo), "xyz_text": xyz_text}