Parsing Output Files¶

from pathlib import Path
from calcflow import parse_qchem_multi_job_output

results = parse_qchem_multi_job_output(Path("mom_calc.out").read_text())
gs_result  = results[0]  # ground state job
es_result  = results[1]  # excited state job

if result.termination_status != "NORMAL":
    raise RuntimeError(f"Calculation failed: {result.termination_status}")

result.final_energy              # float | None — total energy in Hartree
result.nuclear_repulsion_energy  # float | None — Hartree
result.dispersion                # DispersionCorrection | None
result.smd                       # SmdResults | None

if result.dispersion:
    print(result.dispersion.energy_hartree)   # DFT-D correction in Hartree
    print(result.dispersion.method)           # e.g. "D3BJ"

if result.smd:
    print(result.smd.total_energy_kcal_mol)
    print(result.smd.electrostatic_kcal_mol)
    print(result.smd.non_electrostatic_kcal_mol)

scf = result.scf
if scf:
    print(scf.converged)       # bool
    print(scf.energy)          # float — final SCF energy in Hartree
    print(scf.n_iterations)    # int

    for i, it in enumerate(scf.iterations):
        print(f"  iter {i+1}: {it.energy:.10f} Hartree  delta={it.delta_energy:.2e}")

if scf.components:
    c = scf.components
    print(c.nuclear_repulsion)    # Hartree
    print(c.one_electron)         # Hartree
    print(c.coulomb)              # Hartree
    print(c.exchange_correlation) # Hartree

orbs = result.orbitals
if orbs:
    for orb in orbs.alpha_orbitals:
        occ = "occ" if orb.occupied else "virt"
        print(f"  MO {orb.index}: {orb.energy:.4f} Hartree  [{occ}]")

    homo = next(o for o in reversed(orbs.alpha_orbitals) if o.occupied)
    lumo = next(o for o in orbs.alpha_orbitals if not o.occupied)
    gap  = lumo.energy - homo.energy
    print(f"HOMO-LUMO gap: {gap * 27.211:.2f} eV")

    # Unrestricted: beta orbitals
    if orbs.beta_orbitals:
        for orb in orbs.beta_orbitals:
            print(f"  beta MO {orb.index}: {orb.energy:.4f} Hartree")

mulliken  = result.get_charges("Mulliken")
hirshfeld = result.get_charges("Hirshfeld")
cm5       = result.get_charges("CM5")
loewdin   = result.get_charges("Loewdin")   # ORCA

if mulliken:
    # charges is a dict: atom_index (0-based) -> charge value
    for idx, charge in mulliken.charges.items():
        print(f"  atom {idx}: {charge:+.4f}")

    # spin densities (unrestricted calculations)
    if mulliken.spins:
        for idx, spin in mulliken.spins.items():
            print(f"  atom {idx} spin: {spin:+.4f}")

mp = result.multipole
if mp:
    dx, dy, dz = mp.dipole.x, mp.dipole.y, mp.dipole.z
    total = (dx**2 + dy**2 + dz**2) ** 0.5
    print(f"Dipole moment: {total:.4f} Debye")

tddft = result.tddft
if tddft:
    # TDA states (Tamm-Dancoff approximation)
    for state in tddft.tda_states:
        print(f"TDA S{state.state_number}: {state.excitation_energy_ev:.3f} eV"
              f"  f={state.oscillator_strength:.4f}")

    # Full TDDFT states
    for state in tddft.tddft_states:
        ev = state.excitation_energy_ev
        nm = 1239.8 / ev
        print(f"S{state.state_number}: {ev:.3f} eV ({nm:.0f} nm)"
              f"  f={state.oscillator_strength:.4f}")

    # Orbital transitions for a specific state
    state = tddft.tddft_states[0]
    for trans in state.transitions:
        print(f"  {trans.from_orbital} -> {trans.to_orbital}  coeff={trans.coefficient:.4f}")

if tddft.nto_analyses:
    for nto in tddft.nto_analyses:
        print(f"State {nto.state_number}: NTO analysis available")

if tddft.transition_density_matrices:
    for tdm in tddft.transition_density_matrices:
        print(f"State {tdm.state_number}:")
        print(f"  hole-electron separation: {tdm.exciton.rhe:.3f} Å")
        print(f"  CT number: {tdm.exciton.ct_number:.4f}")

if tddft.unrelaxed_density_matrices:
    for udm in tddft.unrelaxed_density_matrices:
        print(f"State {udm.state_number}:")
        if udm.mulliken:
            print(f"  hole/electron Mulliken charges available")

adc = result.adc
if adc:
    print(f"Method: {adc.method}")  # e.g. "ADC(2)"

    gs = adc.ground_state
    print(f"HF energy:  {gs.hf_energy:.6f} Hartree")
    print(f"MP2 energy: {gs.mp2_energy:.6f} Hartree")

    for state in adc.excited_states:
        print(f"State {state.state_number} ({state.multiplicity}):")
        print(f"  {state.excitation_energy_ev:.3f} eV"
              f"  OPA f={state.opa_oscillator_strength:.4f}")
        if state.tpa_cross_section is not None:
            print(f"  TPA cross section: {state.tpa_cross_section:.2f} GM")

if result.input_geometry:
    print(result.input_geometry.to_xyz_str())

if result.final_geometry:
    result.final_geometry.to_xyz_file("optimized.xyz")
    print(f"Optimized energy: {result.final_geometry.energy}")

if result.timing:
    for module, duration in result.timing.modules.items():
        print(f"{module}: {duration:.1f} s")

# Serialize — raw_output excluded, schema_version included
json_str = result.to_json()

# Reload — works across calcflow versions (with migrations)
from calcflow.common.results import CalculationResult
result2 = CalculationResult.from_json(json_str)

# Dict form
data = result.to_dict()
result3 = CalculationResult.from_dict(data)

from pathlib import Path
from calcflow import parse_qchem_output, spectrum_from_excited_states, make_energy_grid
import numpy as np

result = parse_qchem_output(Path("tddft.out").read_text())
states = result.tddft.tddft_states

grid     = make_energy_grid(
    energies=[s.excitation_energy_ev for s in states],
    padding=1.0, n_points=2000
)
spectrum = spectrum_from_excited_states(states, grid, fwhm=0.3)  # eV FWHM

np.savetxt("spectrum.dat", np.column_stack([grid, spectrum]))