#!/usr/bin/env python3 """Calculate paramagnetic NMR shifts from complete ORCA output files. by Dr. Robert Gericke, Adrian Näder and Dr. Michael Patzschke This implements the susceptibility formalism in Pell, Pintacuda and Grey, Prog. Nucl. Magn. Reson. Spectrosc. 111 (2019) 1-271, https://doi.org/10.1016/j.pnmrs.2018.05.001. Inputs are (1) a DFT output containing HFC A matrices in MHz (and optionally a g-matrix), and (2) a CASSCF/NEVPT2 output containing temperature-dependent chi*T in cm^3 K mol^-1 and effective-Hamiltonian total g-matrices. Choose the susceptibility level, g source, temperature and nuclei at the prompts or with CLI flags. ORCA atom indices are retained verbatim. Examples: python ORCA_2FCSPCS.py dft.out casscf_nevpt2.out python ORCA_2FCSPCS.py dft.out casscf_nevpt2.out \ --level NEVPT2 --g-source DFT --temperature 298 \ --nuclei 1H 13C 19F 31P 195Pt --output shifts.txt --csv shifts.csv The program requires NumPy. When ORCA prints only a total A matrix, the total shift remains available and unprinted decomposition columns read nan. The two calculations must have the same atomic coordinates and Cartesian axes. For a given matrix A, the total paramagnetic shift (TPMS) in ppm is evaluated via: delta_para = 2*pi*1e-18 / (3*mu0*muB*gamma_I) * Tr[chi * inv(g.T) * A] A can be decomposed into isotropic and anisotropic parts, corresponding to contact shift (CS) and pseudo-contact shift (PCS), respectively: A_iso = I * Tr(A)/3 A_aniso = A - A_iso CS = delta_para(A_iso) PCS = delta_para(A_aniso) TPMS = CS + PCS = delta_para(A) According to the contributions outputted by Orca (printlevel 3 in eprnmr module) A_iso and A_aniso can be further split into: A_iso = A_FC + Tr(A_SD)/3 + Tr(A_SOC)/3 + Tr(A_GR)/3 A_aniso = (A_SD - I*Tr(A_SD)/3) + (A_SOC - I*Tr(A_SOC)/3) + (A_GR - I*Tr(A_GR)/3) The resulting output table contains the TPMS and the following contributions (for each requested nucleus) defined by their contribution to the HFC A-tensor: abbreviation A-contribution definition ------------------------------------------------------------------------------ TPMS full A-tensor total paramagnetic shift CS A_iso contact shift PCS A_aniso pseudo-contact shift FC A_iso(FC) Fermi-contact shift (due to spin-density at nucleus) PSO_iso A_iso(SOC) isotropic paramagnetic spin-orbit part SD A_aniso(SD) spin-dipolar shift (SD_iso) A_iso(SD) spin-dipolar contrib. to CS --> should be negligible PSO_aniso A_aniso(SOC) anisotropic paramagnetic spin-orbit part GR A(GR) Gauge correction (anisotropic part can become siginifcant) The text table contains: Index Nucleus TPMS CS PCS FC PSO_iso SD PSO_aniso PSO_total GR PSO_total and GR are aggregates already included in TPMS, not extra terms to add. The CSV uses Total_Paramagnetic_Shift (TPMS) as its descriptive TPMS header and includes the rank-1 correction and symmetric-only PCS for comparison. """ from __future__ import annotations import argparse import bisect import csv import io import math import re import sys from dataclasses import dataclass from pathlib import Path from typing import Iterable, Optional import numpy as np # constants MU0 = 4.0 * math.pi * 1e-7 # T^2 m^3 J^-1 MU_B = 9.2740100783e-24 # J T^-1 N_A = 6.02214129e23 # mol^-1 # Gyromagnetic ratios in rad s^-1 T^-1 # workflow and the ORCA P-factor ratios in the supplied HFC output. GAMMA_RAD_S_T: dict[tuple[str, int], float] = { ("H", 1): 26.7522e7, #CODATA 2018 ("C", 13): 6.7282e7, #NMR Enc. 1996 ("F", 19): 25.1662e7, #NMR Enc. 1996 ("P", 31): 10.8291e7, #NMR Enc. 1996 ("Pt", 195): 5.8385e7, #NMR Enc. 1996 ("N", 15): -2.7116e7, #NMR Enc. 1996 ("Si", 29): -5.3190e7, #NMR Enc. 1996 ("O", 17): -3.6264e7, #NMR Enc. 1996 ("Li", 7): 10.3962e7, #NMR Enc. 1996 ("B", 11): 8.5847e7, #NMR Enc. 1996 ("Na", 23): 7.0761e7, #NMR Enc. 1996 ("Al", 27): 6.9763e7, #NMR Enc. 1996 } # ORCA's printed proton P factor. The P factors carry the sign of gamma and # permit other NMR-active isotopes without an ever-growing hard-coded table. PROTON_P_MHZ_BOHR3 = 533.5514 FLOAT = r"[-+]?\d+(?:\.\d*)?(?:[EeDd][-+]?\d+)?" @dataclass(frozen=True) class HFCData: """One ORCA HFC block; ``atom_index`` is exactly ORCA's index.""" atom_index: int element: str isotope: int p_mhz_bohr3: Optional[float] a_total_mhz: np.ndarray a_sd_mhz: Optional[np.ndarray] a_soc_mhz: Optional[np.ndarray] a_gauge_mhz: Optional[np.ndarray] a_fc_printed_mhz: Optional[float] @property def atom_label(self) -> str: return f"{self.atom_index}{self.element}" @property def isotope_label(self) -> str: return f"{self.isotope}{self.element}" @property def a_fc_residual_mhz(self) -> np.ndarray: if any(x is None for x in (self.a_sd_mhz, self.a_soc_mhz, self.a_gauge_mhz)): raise ValueError(f"Missing HFC decomposition for {self.atom_label}") return self.a_total_mhz - self.a_sd_mhz - self.a_soc_mhz - self.a_gauge_mhz @dataclass(frozen=True) class MagneticData: temperature_k: float chi_t_cgs_cm3_k_mol: np.ndarray g_matrix: np.ndarray @dataclass class MagneticSeries: """One level's zero-field chi*T grid and effective-Hamiltonian total g.""" temperatures: list[float] chi_tensors: list[np.ndarray] g_matrix: Optional[np.ndarray] = None @dataclass(frozen=True) class ShiftResult: index: int nucleus: str atom_label: str tpms: float cs: float pcs: float fc: float pso_iso: float sd: float sd_iso: float pso_aniso: float pso_total: float gr: float gr_iso: float gr_aniso: float fc_aniso: float rank1_in_pcs: float pcs_symmetric_only: float component_closure: float fc_printed_difference_mhz: Optional[float] def _float(value: str) -> float: return float(value.replace("D", "E").replace("d", "e")) def read_text(path: Path) -> str: return path.read_text(errors="replace") def _matrix_rows(lines: Iterable[str], marker: str) -> np.ndarray: """Parse exactly three ORCA matrix rows, with an optional 0/1/2 header.""" rows: list[list[float]] = [] for count, line in enumerate(lines): if count >= 10: break values = line.split() if values == ["0", "1", "2"]: continue if len(values) == 4 and values[0] == str(len(rows)): values = values[1:] if len(values) == 3 and all(re.fullmatch(FLOAT, x) for x in values): rows.append([_float(x) for x in values]) if len(rows) == 3: result = np.asarray(rows, dtype=float) if not np.all(np.isfinite(result)): break return result elif rows: break raise ValueError(f"Incomplete or invalid 3x3 matrix after {marker!r}") def parse_3x3_after(text: str, marker: str, start: int = 0) -> np.ndarray: """Read a plain or row-indexed ORCA 3x3 matrix after ``marker``.""" position = text.find(marker, start) if position < 0: raise ValueError(f"Could not find matrix block {marker!r}") return _matrix_rows( iter(text[position + len(marker):].splitlines()), marker ) def parse_magnetic_series(lines: Iterable[str]) -> dict[str, MagneticSeries]: """Scan a full ORCA file once, associating chi grids and g with QDPT levels. Only the first total g-matrix following an EFFECTIVE HAMILTONIAN heading is retained; the later Kramers-pair, S and L contributions are not substitutes. A small extracted block without QDPT headings is marked UNLABELLED. """ series: dict[str, MagneticSeries] = {} level = "UNLABELLED" in_chi_table = False new_chi_table = False pending_temperature: Optional[float] = None awaiting_g = False iterator = iter(lines) for line in iterator: level_match = re.search( r"QDPT WITH\s+(CASSCF|NEVPT2)\s+DIAGONAL ENERGIES", line, re.I ) if level_match: level = level_match.group(1).upper() in_chi_table = False awaiting_g = False if "TEMPERATURE DEPENDENT MOLAR MAGNETIC SUSCEPTIBILITY TENSOR" in line: in_chi_table = True new_chi_table = True pending_temperature = None elif "TEMPERATURE DEPENDENT MAGNETIC SUSCEPTIBILITY" in line: in_chi_table = False if in_chi_table: if "Magnetic field (Gauss):" in line: field_values = re.findall(FLOAT, line.split(":", 1)[1]) if len(field_values) >= 3 and any(abs(_float(x)) > 1e-8 for x in field_values[:3]): in_chi_table = False # Kurland–McGarvey input: zero field continue temp_match = re.search(rf"\bTEMPERATURE/K:\s*({FLOAT})", line) if temp_match: pending_temperature = _float(temp_match.group(1)) elif "Tensor in molecular frame (cm3*K/mol):" in line: if pending_temperature is None or pending_temperature <= 0: raise ValueError("Susceptibility tensor lacks a positive temperature") chi_t = _matrix_rows(iterator, "susceptibility tensor") if not np.allclose(chi_t, chi_t.T, rtol=0, atol=2e-6): raise ValueError(f"Non-symmetric chi*T at {pending_temperature:g} K") if new_chi_table: # A later zero-field table at this level supersedes the old. series[level] = MagneticSeries([], []) new_chi_table = False section = series[level] section.temperatures.append(pending_temperature) section.chi_tensors.append(chi_t) pending_temperature = None if "ELECTRONIC G-MATRIX FROM EFFECTIVE HAMILTONIAN" in line: in_chi_table = False awaiting_g = True elif awaiting_g and re.match(r"^\s*total g-matrix:\s*$", line, re.I): section = series.setdefault(level, MagneticSeries([], [])) if section.g_matrix is None: section.g_matrix = _matrix_rows(iterator, "effective-Hamiltonian total g") awaiting_g = False if not series: raise ValueError("No zero-field susceptibility or effective-Hamiltonian g found") for name, section in series.items(): if section.temperatures and len(set(section.temperatures)) != len(section.temperatures): raise ValueError(f"Repeated printed temperatures in {name} susceptibility grid") return series def read_magnetic_file(path: Path) -> dict[str, MagneticSeries]: with path.open(errors="replace") as handle: return parse_magnetic_series(handle) def susceptibility_at_temperature( section: MagneticSeries, temperature: float ) -> np.ndarray: """Use the printed tensor, or interpolate chi*T between two printed rows.""" if not math.isfinite(temperature) or temperature <= 0: raise ValueError("Temperature must be a finite, positive number in K") if not section.temperatures: raise ValueError("The selected level contains no susceptibility tensors") order = sorted(range(len(section.temperatures)), key=section.temperatures.__getitem__) temperatures = [section.temperatures[i] for i in order] position = bisect.bisect_left(temperatures, temperature) if position < len(temperatures) and math.isclose(temperatures[position], temperature, abs_tol=1e-8): return section.chi_tensors[order[position]] if position == 0 or position == len(temperatures): raise ValueError( f"{temperature:g} K is outside the printed chi*T range " f"{temperatures[0]:g}–{temperatures[-1]:g} K; no extrapolation is performed" ) lower, upper = temperatures[position - 1:position + 1] weight = (temperature - lower) / (upper - lower) return ((1.0 - weight) * section.chi_tensors[order[position - 1]] + weight * section.chi_tensors[order[position]]) def parse_magnetic_data(text: str) -> MagneticData: """Compatibility helper for a single extracted chi/g block.""" sections = parse_magnetic_series(io.StringIO(text)) section = next((s for s in sections.values() if s.temperatures and s.g_matrix is not None), None) if section is None: raise ValueError("No matching susceptibility tensor and g-matrix found") return MagneticData(section.temperatures[0], section.chi_tensors[0], section.g_matrix) def read_geometry(path: Path) -> Optional[list[tuple[str, np.ndarray]]]: """Read the first molecular-frame Cartesian geometry for an orientation check.""" with path.open(errors="replace") as handle: for line in handle: if "CARTESIAN COORDINATES (ANGSTROEM)" not in line: continue geometry: list[tuple[str, np.ndarray]] = [] for row in handle: fields = row.split() if len(fields) == 4 and re.fullmatch(r"[A-Za-z]+", fields[0]) and all( re.fullmatch(FLOAT, field) for field in fields[1:] ): geometry.append((fields[0].capitalize(), np.array([_float(x) for x in fields[1:]]))) elif geometry: return geometry return geometry or None return None def validate_same_geometry(dft_path: Path, mr_path: Path) -> None: """Tensor contractions require matching molecular axes and atom order.""" left, right = read_geometry(dft_path), read_geometry(mr_path) if left is None or right is None: # historical truncated extracts return if len(left) != len(right) or any(a[0] != b[0] for a, b in zip(left, right)): raise ValueError("DFT and multireference atom lists differ; tensor frames cannot be combined") difference = max(float(np.linalg.norm(a[1] - b[1])) for a, b in zip(left, right)) if difference > 1e-3: raise ValueError( f"DFT and multireference coordinates differ by up to {difference:.4g} Å; " "use outputs in the same molecular Cartesian frame" ) def _optional_matrix(block: str, markers: Iterable[str]) -> Optional[np.ndarray]: for marker in markers: if marker in block: return parse_3x3_after(block, marker) return None HFC_HEADER = re.compile( r"Nucleus\s+(\d+)\s*([A-Za-z]+)\s*:\s*A\s*:\s*" r"(?:Isotope|ISTP)\s*=\s*(\d+)\b[^\n]*", re.I ) def _parse_hfc_block(block: str) -> HFCData: match = HFC_HEADER.search(block) if match is None: raise ValueError("Invalid ORCA HFC nucleus heading") atom_index = int(match.group(1)) raw_element = match.group(2) element = raw_element[0].upper() + raw_element[1:].lower() isotope = int(match.group(3)) p_match = re.search(rf"\bP\s*=\s*({FLOAT})\b", match.group(0), re.I) p_value = _float(p_match.group(1)) if p_match else None total = _optional_matrix(block, ( "Total HFC matrix (all values in MHz):", "Raw HFC matrix (all values in MHz):", )) if total is None: raise ValueError(f"No HFC A matrix in ORCA block for {atom_index}{element}") fc_match = re.search(rf"^\s*A\(FC\)\s+({FLOAT})", block, re.M) return HFCData( atom_index=atom_index, element=element, isotope=isotope, p_mhz_bohr3=p_value, a_total_mhz=total, a_sd_mhz=_optional_matrix(block, ("SD contribution to HFC matrix:",)), a_soc_mhz=_optional_matrix(block, ( "NOC/SOC contribution to HFC matrix:", "Orbital contribution to HFC matrix:", )), a_gauge_mhz=_optional_matrix(block, ( "Gauge correction (diamagnetic) contribution to HFC matrix:", "Gauge correction contribution to HFC matrix:", )), a_fc_printed_mhz=_float(fc_match.group(1)) if fc_match else None, ) def parse_all_hfc(text: str) -> list[HFCData]: """Compatibility helper for extracted HFC text.""" matches = list(HFC_HEADER.finditer(text)) if not matches: raise ValueError("No ORCA HFC nucleus blocks were found") return [_parse_hfc_block(text[m.start():matches[i + 1].start() if i + 1 < len(matches) else len(text)]) for i, m in enumerate(matches)] def read_dft_file(path: Path) -> tuple[list[HFCData], Optional[np.ndarray]]: """Stream a full ORCA DFT output, retaining only small HFC blocks and g.""" nuclei: list[HFCData] = [] g_matrix: Optional[np.ndarray] = None block: list[str] = [] with path.open(errors="replace") as handle: for line in handle: if HFC_HEADER.search(line): if block: nuclei.append(_parse_hfc_block("".join(block))) block = [line] elif block: block.append(line) # The A-tensor note closes the block; in old outputs the next # Nucleus heading or EOF also closes it. if "The A matrix conforms" in line: nuclei.append(_parse_hfc_block("".join(block))) block = [] elif len(block) > 400: raise ValueError("HFC nucleus block exceeds 400 lines; check ORCA format") if re.match(r"^\s*The g-matrix:\s*$", line, re.I): g_matrix = _matrix_rows(handle, "DFT g-matrix") if block: nuclei.append(_parse_hfc_block("".join(block))) if not nuclei: raise ValueError(f"No HFC A matrices in {path}; a DFT A calculation is required") if len({(x.atom_index, x.element, x.isotope) for x in nuclei}) != len(nuclei): raise ValueError("Duplicate HFC nucleus blocks; supply a single ORCA job") return nuclei, g_matrix def susceptibility_angstrom3(data: MagneticData) -> np.ndarray: """Convert ORCA chi*T (cgs cm^3 K mol^-1) to molecular chi in A^3.""" return 4.0 * math.pi * 1e24 * data.chi_t_cgs_cm3_k_mol / (N_A * data.temperature_k) def rank0(matrix: np.ndarray) -> np.ndarray: """Spherical rank-0 part""" return np.eye(3) * (np.trace(matrix) / 3.0) def rank1(matrix: np.ndarray) -> np.ndarray: """Antisymmetric rank-1 part""" return 0.5 * (matrix - matrix.T) def rank2(matrix: np.ndarray) -> np.ndarray: """Symmetric traceless rank-2 part""" symmetric = 0.5 * (matrix + matrix.T) return symmetric - np.eye(3) * (np.trace(symmetric) / 3.0) def traceless(matrix: np.ndarray) -> np.ndarray: """Full traceless part of a matrix.""" return matrix - rank0(matrix) def shift_ppm( matrix_mhz: np.ndarray, chi_angstrom3: np.ndarray, g_matrix: np.ndarray, gamma_rad_s_t: float, ) -> float: """Directly expressed as chemical shift.""" prefactor = 2.0 * math.pi * 1e-18 / (3.0 * MU0 * MU_B * gamma_rad_s_t) contraction = np.trace(chi_angstrom3 @ np.linalg.inv(g_matrix.T) @ matrix_mhz) return float(prefactor * contraction) def nuclear_gamma(hfc: HFCData) -> float: """Use legacy calibrated values or ORCA's isotope-specific signed P factor.""" key = (hfc.element, hfc.isotope) if key in GAMMA_RAD_S_T: return GAMMA_RAD_S_T[key] if hfc.p_mhz_bohr3 is not None and hfc.p_mhz_bohr3 != 0: return GAMMA_RAD_S_T[("H", 1)] * hfc.p_mhz_bohr3 / PROTON_P_MHZ_BOHR3 raise ValueError( f"No nonzero gyromagnetic ratio for {hfc.isotope_label}: " "the ORCA nucleus heading must contain a nonzero P factor" ) def calculate_shift(hfc: HFCData, magnetic: MagneticData) -> ShiftResult: gamma = nuclear_gamma(hfc) chi = susceptibility_angstrom3(magnetic) contract = lambda matrix: shift_ppm(matrix, chi, magnetic.g_matrix, gamma) total = hfc.a_total_mhz # trace without symmetrizing A. cs = contract(rank0(total)) pcs = contract(traceless(total)) tpms = cs + pcs rank1_in_pcs = contract(rank1(total)) pcs_symmetric_only = contract(rank2(total)) # ORCA can print only the total A tensor. CS/PCS/TPMS remain valid, but # unpublished components must not be misrepresented as zero. if any(x is None for x in (hfc.a_sd_mhz, hfc.a_soc_mhz, hfc.a_gauge_mhz)): missing = math.nan return ShiftResult( index=hfc.atom_index, nucleus=hfc.isotope_label, atom_label=hfc.atom_label, tpms=tpms, cs=cs, pcs=pcs, fc=missing, pso_iso=missing, sd=missing, sd_iso=missing, pso_aniso=missing, pso_total=missing, gr=missing, gr_iso=missing, gr_aniso=missing, fc_aniso=missing, rank1_in_pcs=rank1_in_pcs, pcs_symmetric_only=pcs_symmetric_only, component_closure=missing, fc_printed_difference_mhz=None, ) fc = hfc.a_fc_residual_mhz sd = hfc.a_sd_mhz soc = hfc.a_soc_mhz gauge_matrix = hfc.a_gauge_mhz fc_shift = contract(rank0(fc)) fc_aniso = contract(traceless(fc)) sd_iso = contract(rank0(sd)) sd_shift = contract(traceless(sd)) pso_iso = contract(rank0(soc)) pso_aniso = contract(traceless(soc)) gr_iso = contract(rank0(gauge_matrix)) gr_aniso = contract(traceless(gauge_matrix)) pso_total = pso_iso + pso_aniso gr = gr_iso + gr_aniso # Every named contribution is included exactly once in this check. component_sum = (fc_shift + fc_aniso + sd_iso + sd_shift + pso_iso + pso_aniso + gr_iso + gr_aniso) fc_printed_difference = None if hfc.a_fc_printed_mhz is not None: fc_printed_difference = np.trace(fc) / 3.0 - hfc.a_fc_printed_mhz return ShiftResult( index=hfc.atom_index, nucleus=hfc.isotope_label, atom_label=hfc.atom_label, tpms=tpms, cs=cs, pcs=pcs, fc=fc_shift, pso_iso=pso_iso, sd=sd_shift, sd_iso=sd_iso, pso_aniso=pso_aniso, pso_total=pso_total, gr=gr, gr_iso=gr_iso, gr_aniso=gr_aniso, fc_aniso=fc_aniso, rank1_in_pcs=rank1_in_pcs, pcs_symmetric_only=pcs_symmetric_only, component_closure=tpms - component_sum, fc_printed_difference_mhz=fc_printed_difference, ) TABLE_COLUMNS = ( "Index", "Nucleus", "TPMS", "CS", "PCS", "FC", "PSO_iso", "SD", "PSO_aniso", "PSO_total", "GR", ) CSV_COLUMNS = (*TABLE_COLUMNS[:2], "Total_Paramagnetic_Shift", *TABLE_COLUMNS[3:]) def format_table(results: list[ShiftResult]) -> str: header = (f"{'Index':<8s} {'Nucleus':<10s} " + " ".join(f"{name:>12s}" for name in TABLE_COLUMNS[2:])) lines = [header] for item in results: values = (item.tpms, item.cs, item.pcs, item.fc, item.pso_iso, item.sd, item.pso_aniso, item.pso_total, item.gr) lines.append( f"{item.index:<8d} {item.nucleus:<10s} " + " ".join(f"{value:12.2f}" for value in values) ) return "\n".join(lines) + "\n" def write_csv(results: list[ShiftResult], path: Path) -> None: fields = [ *CSV_COLUMNS, "SD_iso", "GR_iso", "GR_aniso", "FC_aniso", "Rank1_in_PCS", "PCS_symmetric_only", "Component_closure", "Atom_label", "FC_printed_difference_MHz", ] with path.open("w", newline="") as handle: writer = csv.DictWriter(handle, fieldnames=fields) writer.writeheader() for result in results: writer.writerow( { "Index": result.index, "Nucleus": result.nucleus, "Total_Paramagnetic_Shift": f"{result.tpms:.10f}", "CS": f"{result.cs:.10f}", "PCS": f"{result.pcs:.10f}", "FC": f"{result.fc:.10f}", "PSO_iso": f"{result.pso_iso:.10f}", "SD": f"{result.sd:.10f}", "SD_iso": f"{result.sd_iso:.10f}", "PSO_aniso": f"{result.pso_aniso:.10f}", "PSO_total": f"{result.pso_total:.10f}", "GR": f"{result.gr:.10f}", "GR_iso": f"{result.gr_iso:.10f}", "GR_aniso": f"{result.gr_aniso:.10f}", "FC_aniso": f"{result.fc_aniso:.10f}", "Rank1_in_PCS": f"{result.rank1_in_pcs:.10f}", "PCS_symmetric_only": f"{result.pcs_symmetric_only:.10f}", "Component_closure": f"{result.component_closure:.3e}", "Atom_label": result.atom_label, "FC_printed_difference_MHz": ( "" if result.fc_printed_difference_mhz is None else f"{result.fc_printed_difference_mhz:.10f}" ), } ) def _canonical_atom_label(value: str) -> str: match = re.fullmatch(r"(\d+)([A-Za-z]+)", value) if match is None: raise ValueError(f"Invalid atom label {value!r}; use e.g. 1Pt or 12H") element = match.group(2)[0].upper() + match.group(2)[1:].lower() return f"{int(match.group(1))}{element}" def _canonical_isotope_label(value: str) -> str: match = re.fullmatch(r"(\d+)([A-Za-z]+)", value.strip()) if match is None: raise ValueError(f"Invalid isotope {value!r}; use e.g. 31P or 195Pt") element = match.group(2).capitalize() return f"{int(match.group(1))}{element}" def select_nuclei( all_hfc: list[HFCData], isotope_labels: list[str], atom_labels: Optional[list[str]] ) -> list[HFCData]: available = {item.isotope_label for item in all_hfc} requested = {_canonical_isotope_label(label) for label in isotope_labels} missing_isotopes = requested - available if missing_isotopes: raise ValueError( "No HFC A tensor for " + ", ".join(sorted(missing_isotopes)) + "; available: " + ", ".join(sorted(available)) ) selected = [item for item in all_hfc if item.isotope_label in requested] if atom_labels: requested = {_canonical_atom_label(label) for label in atom_labels} selected = [item for item in selected if item.atom_label in requested] missing = requested - {item.atom_label for item in selected} if missing: raise ValueError("Requested atoms not found: " + ", ".join(sorted(missing))) if not selected: raise ValueError("No requested HFC nuclei were found in the DFT output") return selected def print_tensors(hfc_items: list[HFCData], magnetic: MagneticData) -> None: print(f"Temperature: {magnetic.temperature_k:.2f} K") print("chi*T [cm^3 K mol^-1]:") print(np.array2string(magnetic.chi_t_cgs_cm3_k_mol, precision=9)) print("chi [Angstrom^3 per molecule]:") print(np.array2string(susceptibility_angstrom3(magnetic), precision=9)) print("g matrix (used without symmetrization):") print(np.array2string(magnetic.g_matrix, precision=9)) for hfc in hfc_items: print(f"\nORCA atom {hfc.atom_index}, {hfc.isotope_label}") for name, matrix in ( ("A_total / MHz", hfc.a_total_mhz), ("A_FC residual / MHz", hfc.a_fc_residual_mhz if all( x is not None for x in (hfc.a_sd_mhz, hfc.a_soc_mhz, hfc.a_gauge_mhz) ) else None), ("A_SD / MHz", hfc.a_sd_mhz), ("A_NOC_SOC / MHz", hfc.a_soc_mhz), ("A_gauge / MHz", hfc.a_gauge_mhz), ): print(name + ":") print(np.array2string(matrix, precision=9) if matrix is not None else "not printed by ORCA") def _ask(prompt: str, default: Optional[str] = None) -> str: suffix = f" [{default}]" if default is not None else "" try: value = input(f"{prompt}{suffix}: ").strip() except EOFError as exc: raise ValueError( f"Missing {prompt.lower()}; supply the corresponding command-line option" ) from exc return value or (default if default is not None else "") def _choice(prompt: str, allowed: list[str], provided: Optional[str], default: str) -> str: if not allowed: raise ValueError(f"No available choices for {prompt.lower()}") if provided is not None: candidate = provided.upper() if candidate == "NEVPT": candidate = "NEVPT2" if candidate not in allowed: raise ValueError(f"{prompt}: {candidate} unavailable; choose {', '.join(allowed)}") return candidate while True: candidate = _ask(f"{prompt} ({', '.join(allowed)})", default).upper() if candidate == "NEVPT": candidate = "NEVPT2" if candidate in allowed: return candidate print("Please choose one of: " + ", ".join(allowed)) def _selected_series(sections: dict[str, MagneticSeries], level: str) -> MagneticSeries: if level in sections: return sections[level] if set(sections) == {"UNLABELLED"}: print(f"Note: the extracted magnetic block has no QDPT level label; assuming {level}.", file=sys.stderr) return sections["UNLABELLED"] raise ValueError(f"No {level} section in multireference output; present: {', '.join(sections)}") def build_parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser( description=( "Predict CS, PCS, SD and PSO contributions from ORCA " "DFT HFC and CASSCF/NEVPT2 susceptibility/g outputs." ) ) parser.add_argument("dft_out", type=Path, nargs="?", help="complete ORCA DFT HFC output") parser.add_argument( "nevpt2_out", type=Path, nargs="?", help="complete ORCA CASSCF/NEVPT2 chi*T and g output" ) parser.add_argument( "--nuclei", nargs="+", default=None, help="isotopes with A tensors, e.g. --nuclei 31P 1H 13C 19F 195Pt; or 'all'", ) parser.add_argument("--g-source", choices=("DFT", "CASSCF", "NEVPT2"), type=lambda s: "NEVPT2" if s.upper() == "NEVPT" else s.upper(), help="g-matrix source (prompted if omitted)") parser.add_argument("--level", choices=("CASSCF", "NEVPT2"), type=lambda s: "NEVPT2" if s.upper() == "NEVPT" else s.upper(), help="susceptibility level (prompted if omitted)") parser.add_argument( "--atoms", nargs="+", default=None, help="optional ORCA atom labels, e.g. --atoms 1Pt 12H", ) parser.add_argument( "--temperature", type=float, default=None, help="temperature in K within the printed chi*T grid (interpolated if needed)", ) parser.add_argument( "--output", type=Path, default=Path("paramagnetic_shifts.txt"), help="text table (default: paramagnetic_shifts.txt)", ) parser.add_argument( "--csv", type=Path, default=None, help="optional full-precision CSV with diagnostic columns", ) parser.add_argument( "--show-tensors", action="store_true", help="print parsed input tensors" ) return parser def main() -> None: parser = build_parser() args = parser.parse_args() try: dft_path = args.dft_out or Path(_ask("DFT ORCA output path")) mr_path = args.nevpt2_out or Path(_ask("CASSCF/NEVPT2 ORCA output path")) validate_same_geometry(dft_path, mr_path) hfc_items, dft_g = read_dft_file(dft_path) sections = read_magnetic_file(mr_path) levels = [name for name in ("CASSCF", "NEVPT2") if name in sections and sections[name].temperatures] if not levels and "UNLABELLED" in sections and sections["UNLABELLED"].temperatures: levels = ["CASSCF", "NEVPT2"] # explicitly assumed for an excerpt if not levels: raise ValueError("No zero-field chi*T tensor at CASSCF or NEVPT2 level") level = _choice("Susceptibility level", levels, args.level, "NEVPT2" if "NEVPT2" in levels else levels[0]) section = _selected_series(sections, level) if not section.temperatures: raise ValueError(f"No chi*T tensors at {level}") range_text = f"{min(section.temperatures):g}–{max(section.temperatures):g} K" default_temperature = 298.0 if min(section.temperatures) <= 298 <= max(section.temperatures) else section.temperatures[0] temp_text = (str(args.temperature) if args.temperature is not None else _ask(f"Temperature in K ({range_text})", f"{default_temperature:g}")) temperature = float(temp_text) chi_t = susceptibility_at_temperature(section, temperature) g_options = (["DFT"] if dft_g is not None else []) + [ name for name in ("CASSCF", "NEVPT2") if name in sections and sections[name].g_matrix is not None ] if set(sections) == {"UNLABELLED"} and sections["UNLABELLED"].g_matrix is not None: g_options += ["CASSCF", "NEVPT2"] g_source = _choice("g-matrix source", g_options, args.g_source, level if level in g_options else g_options[0]) g_matrix = dft_g if g_source == "DFT" else _selected_series(sections, g_source).g_matrix if g_matrix is None or np.linalg.cond(g_matrix) > 1e12: raise ValueError(f"The {g_source} g-matrix is missing or singular") magnetic = MagneticData(temperature, chi_t, g_matrix) available = sorted({item.isotope_label for item in hfc_items}, key=lambda s: (int(re.match(r"\d+", s).group()), s)) nuclei_text = (args.nuclei if args.nuclei is not None else _ask("NMR nuclei (comma-separated, or all; A available for " + ", ".join(available) + ")", "all").replace(",", " ").split()) labels = available if [x.lower() for x in nuclei_text] == ["all"] else [ part for value in nuclei_text for part in value.split(",") if part ] selected = select_nuclei(hfc_items, labels, args.atoms) results = [calculate_shift(item, magnetic) for item in selected] inputs = {dft_path.resolve(), mr_path.resolve()} if args.output.resolve() in inputs or (args.csv is not None and args.csv.resolve() in inputs): raise ValueError("Output path matches an ORCA input file") if args.csv is not None and args.csv.resolve() == args.output.resolve(): raise ValueError("--csv and --output must name different files") except (ValueError, OSError) as exc: parser.error(str(exc)) if args.show_tensors: print_tensors(selected, magnetic) print() table = format_table(results) print(f"chi*T level: {level}; temperature: {temperature:g} K; g source: {g_source}") if temperature not in section.temperatures: print("chi*T was linearly interpolated between printed temperatures.") if any(math.isnan(item.pso_total) for item in results): print("Some A decompositions were not printed; their component columns are nan.") print(table, end="") args.output.parent.mkdir(parents=True, exist_ok=True) args.output.write_text(table) if args.csv is not None: args.csv.parent.mkdir(parents=True, exist_ok=True) write_csv(results, args.csv) print(f"\nWrote {args.output}") if args.csv is not None: print(f"Wrote {args.csv}") print("All shifts are in ppm; ORCA atom indices are unchanged.") if __name__ == "__main__": main()