chematic_mol/
gaussian.rs

1//! Gaussian 16 / Gaussian 09 file format support.
2//!
3//! ## Supported formats
4//!
5//! | Format | Read | Write | Description |
6//! |--------|------|-------|-------------|
7//! | `.gjf` / `.com` | ✅ | ✅ | Gaussian input file |
8//! | `.log` / `.out` | ✅ | — | Gaussian output file (geometry + energy) |
9//!
10//! ## GJF structure
11//!
12//! ```text
13//! %chk=ethanol.chk
14//! # B3LYP/6-31G* opt
15//!
16//! Ethanol geometry optimization
17//!
18//! 0 1
19//! C   0.000000  0.000000  0.000000
20//! C   1.531000  0.000000  0.000000
21//! O   2.058000  1.198000  0.000000
22//! H   ...
23//! ```
24//!
25//! Coordinates are returned as `Vec<(f64, f64, f64)>` (Ångströms, same
26//! convention as `chematic-mol`'s other parsers).
27
28use chematic_core::{Atom, AtomIdx, Element, Molecule, MoleculeBuilder};
29
30/// Return type shared by the GJF parser and its coordinate helper:
31/// `(Molecule, coords, charge, multiplicity)`.
32type GjfResult = (Molecule, Vec<(f64, f64, f64)>, i32, u32);
33
34// ---------------------------------------------------------------------------
35// Error type
36// ---------------------------------------------------------------------------
37
38/// Errors that can occur when parsing Gaussian files.
39#[derive(Debug, Clone, PartialEq)]
40pub enum GaussianError {
41    /// The charge/multiplicity line was missing or malformed.
42    MissingChargeMultiplicity,
43    /// No atom coordinates were found in the file.
44    NoAtoms,
45    /// An element symbol or atomic number was not recognised.
46    UnknownElement(String),
47    /// A coordinate value could not be parsed as a float.
48    InvalidCoordinate(String),
49    /// The log file contained no `Standard orientation:` block.
50    NoStandardOrientation,
51}
52
53impl core::fmt::Display for GaussianError {
54    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
55        match self {
56            Self::MissingChargeMultiplicity => {
57                write!(f, "charge/multiplicity line not found in GJF")
58            }
59            Self::NoAtoms => write!(f, "no atom coordinates found in file"),
60            Self::UnknownElement(s) => write!(f, "unknown element '{s}'"),
61            Self::InvalidCoordinate(s) => write!(f, "invalid coordinate value '{s}'"),
62            Self::NoStandardOrientation => {
63                write!(f, "no 'Standard orientation' block found in Gaussian log")
64            }
65        }
66    }
67}
68
69impl std::error::Error for GaussianError {}
70
71// ---------------------------------------------------------------------------
72// Result type for LOG files
73// ---------------------------------------------------------------------------
74
75/// Data extracted from a Gaussian output (.log / .out) file.
76pub struct GaussianLogResult {
77    /// Molecular topology (atoms only; no bonds).
78    pub mol: Molecule,
79    /// Atomic coordinates in Ångströms `(x, y, z)`.
80    pub coords: Vec<(f64, f64, f64)>,
81    /// Final SCF energy in hartree, if present.
82    pub scf_energy: Option<f64>,
83}
84
85// ---------------------------------------------------------------------------
86// GJF parser
87// ---------------------------------------------------------------------------
88
89/// Parse a Gaussian input file (`.gjf` / `.com`).
90///
91/// Returns `(Molecule, coords, charge, multiplicity)`.
92/// `coords` is a `Vec<(f64, f64, f64)>` in Ångströms.
93/// The `Molecule` has no bonds; use `determine_bonds` to infer connectivity.
94pub fn parse_gjf(input: &str) -> Result<GjfResult, GaussianError> {
95    // Split into sections at blank lines.
96    let sections: Vec<Vec<&str>> = {
97        let mut secs: Vec<Vec<&str>> = Vec::new();
98        let mut current: Vec<&str> = Vec::new();
99        for line in input.lines() {
100            let trimmed = line.trim();
101            if trimmed.is_empty() {
102                if !current.is_empty() {
103                    secs.push(current.clone());
104                    current.clear();
105                }
106            } else {
107                current.push(trimmed);
108            }
109        }
110        if !current.is_empty() {
111            secs.push(current);
112        }
113        secs
114    };
115
116    // GJF structure: [%directives +] #route [blank] title [blank] charge/mult coords...
117    // Find the route section (a section containing a line that starts with '#'),
118    // then the charge/multiplicity section is two sections later.
119    // This avoids falsely matching a title line that happens to be "0 1".
120    let route_idx = sections
121        .iter()
122        .position(|sec| sec.iter().any(|l| l.starts_with('#')))
123        .ok_or(GaussianError::MissingChargeMultiplicity)?;
124    let cm_idx = route_idx + 2;
125    if cm_idx >= sections.len() {
126        return Err(GaussianError::MissingChargeMultiplicity);
127    }
128    // Verify the section looks like a charge/multiplicity line.
129    {
130        let parts: Vec<&str> = sections[cm_idx][0].split_whitespace().collect();
131        if parts.len() < 2 || parts[0].parse::<i32>().is_err() || parts[1].parse::<u32>().is_err() {
132            return Err(GaussianError::MissingChargeMultiplicity);
133        }
134    }
135
136    let parts: Vec<&str> = sections[cm_idx][0].split_whitespace().collect();
137    let charge: i32 = parts[0].parse().unwrap();
138    let multiplicity: u32 = parts[1].parse().unwrap();
139
140    // Coordinates follow the charge line in the same or next section.
141    let coord_lines: &[&str] = if sections[cm_idx].len() > 1 {
142        &sections[cm_idx][1..]
143    } else if let Some(next) = sections.get(cm_idx + 1) {
144        next.as_slice()
145    } else {
146        return Err(GaussianError::NoAtoms);
147    };
148
149    parse_atom_coords(coord_lines, charge, multiplicity)
150}
151
152fn parse_atom_coords(
153    lines: &[&str],
154    charge: i32,
155    multiplicity: u32,
156) -> Result<GjfResult, GaussianError> {
157    let mut builder = MoleculeBuilder::new();
158    let mut coords: Vec<(f64, f64, f64)> = Vec::new();
159
160    for line in lines {
161        let line = line.trim();
162        if line.is_empty() || line.starts_with('!') {
163            continue;
164        }
165        if line.starts_with("Variables:") || line.starts_with("Constants:") {
166            break;
167        }
168        let parts: Vec<&str> = line.split_whitespace().collect();
169        if parts.len() < 4 {
170            continue;
171        }
172        let raw_sym = parts[0].trim_end_matches(|c: char| c.is_ascii_digit());
173        let elem = if raw_sym.is_empty() {
174            // Bare atomic number (e.g. "6" for carbon — valid Gaussian input).
175            let atomic_num: u8 = parts[0]
176                .parse()
177                .map_err(|_| GaussianError::UnknownElement(parts[0].to_string()))?;
178            Element::from_atomic_number(atomic_num)
179                .ok_or_else(|| GaussianError::UnknownElement(parts[0].to_string()))?
180        } else {
181            Element::from_symbol(raw_sym)
182                .ok_or_else(|| GaussianError::UnknownElement(raw_sym.to_string()))?
183        };
184
185        let x: f64 = parts[1]
186            .parse()
187            .map_err(|_| GaussianError::InvalidCoordinate(parts[1].to_string()))?;
188        let y: f64 = parts[2]
189            .parse()
190            .map_err(|_| GaussianError::InvalidCoordinate(parts[2].to_string()))?;
191        let z: f64 = parts[3]
192            .parse()
193            .map_err(|_| GaussianError::InvalidCoordinate(parts[3].to_string()))?;
194
195        builder.add_atom(Atom::new(elem));
196        coords.push((x, y, z));
197    }
198
199    if coords.is_empty() {
200        return Err(GaussianError::NoAtoms);
201    }
202
203    Ok((builder.build(), coords, charge, multiplicity))
204}
205
206// ---------------------------------------------------------------------------
207// GJF writer
208// ---------------------------------------------------------------------------
209
210/// Write a Gaussian input file (`.gjf`) string.
211///
212/// - `method`: route section keywords, e.g. `"B3LYP/6-31G* opt"`.
213/// - `title`: job title comment (uses `"chematic"` if empty).
214pub fn write_gjf(
215    mol: &Molecule,
216    coords: &[(f64, f64, f64)],
217    charge: i32,
218    multiplicity: u32,
219    method: &str,
220    title: &str,
221) -> String {
222    let method = if method.is_empty() {
223        "B3LYP/6-31G* opt"
224    } else {
225        method
226    };
227    let title = if title.is_empty() { "chematic" } else { title };
228
229    let mut out = format!("# {method}\n\n{title}\n\n{charge} {multiplicity}\n");
230
231    for i in 0..mol.atom_count() {
232        let idx = AtomIdx(i as u32);
233        let sym = mol.atom(idx).element.symbol();
234        let (x, y, z) = coords.get(i).copied().unwrap_or((0.0, 0.0, 0.0));
235        out.push_str(&format!("{sym:<3}  {:12.6}  {:12.6}  {:12.6}\n", x, y, z));
236    }
237    out.push('\n');
238    out
239}
240
241// ---------------------------------------------------------------------------
242// Gaussian LOG parser
243// ---------------------------------------------------------------------------
244
245/// Parse a Gaussian output file (`.log` / `.out`).
246///
247/// Extracts:
248/// - The **last** `Standard orientation:` geometry block.
249/// - The **last** `SCF Done:` energy in hartree (if present).
250///
251/// Returns [`GaussianLogResult`] with no bond information.
252pub fn parse_gaussian_log(input: &str) -> Result<GaussianLogResult, GaussianError> {
253    let lines: Vec<&str> = input.lines().collect();
254
255    // Locate last Standard orientation: block.
256    let last_orient_start = lines
257        .iter()
258        .rposition(|l| l.contains("Standard orientation:"))
259        .ok_or(GaussianError::NoStandardOrientation)?;
260
261    let mut builder = MoleculeBuilder::new();
262    let mut coords: Vec<(f64, f64, f64)> = Vec::new();
263    let mut dashes_seen = 0usize;
264    let mut in_table = false;
265
266    for line in &lines[last_orient_start + 1..] {
267        let trimmed = line.trim();
268        if trimmed.starts_with("---") {
269            dashes_seen += 1;
270            if dashes_seen == 2 {
271                in_table = true; // after 2nd dashes line
272            } else if in_table {
273                break; // closing dashes = end of table
274            }
275            continue;
276        }
277        if !in_table {
278            continue;
279        }
280        // Gaussian 09/16: Center_Num Atomic_Num Atomic_Type X Y Z  (6 cols)
281        // Gaussian 03:    Center_Num Atomic_Num X Y Z             (5 cols)
282        let parts: Vec<&str> = trimmed.split_whitespace().collect();
283        let (an_col, x_col) = match parts.len() {
284            n if n >= 6 => (1, 3), // standard 6-column format
285            5 => (1, 2),           // Gaussian 03 5-column format (no Atomic_Type)
286            _ => continue,
287        };
288        let atomic_num: u8 = parts[an_col]
289            .parse()
290            .map_err(|_| GaussianError::UnknownElement(parts[an_col].to_string()))?;
291        let elem = Element::from_atomic_number(atomic_num)
292            .ok_or_else(|| GaussianError::UnknownElement(parts[an_col].to_string()))?;
293        let x: f64 = parts[x_col]
294            .parse()
295            .map_err(|_| GaussianError::InvalidCoordinate(parts[x_col].to_string()))?;
296        let y: f64 = parts[x_col + 1]
297            .parse()
298            .map_err(|_| GaussianError::InvalidCoordinate(parts[x_col + 1].to_string()))?;
299        let z: f64 = parts[x_col + 2]
300            .parse()
301            .map_err(|_| GaussianError::InvalidCoordinate(parts[x_col + 2].to_string()))?;
302
303        builder.add_atom(Atom::new(elem));
304        coords.push((x, y, z));
305    }
306
307    if coords.is_empty() {
308        return Err(GaussianError::NoAtoms);
309    }
310
311    // Extract last SCF Done energy.
312    let scf_energy = lines
313        .iter()
314        .rfind(|l| l.contains("SCF Done:"))
315        .and_then(|l| {
316            let after = l[l.find('=')? + 1..].trim();
317            after.split_whitespace().next()?.parse::<f64>().ok()
318        });
319
320    Ok(GaussianLogResult {
321        mol: builder.build(),
322        coords,
323        scf_energy,
324    })
325}
326
327// ---------------------------------------------------------------------------
328// Tests
329// ---------------------------------------------------------------------------
330
331#[cfg(test)]
332mod tests {
333    use super::*;
334
335    const ETHANOL_GJF: &str = r#"# B3LYP/6-31G* opt
336
337Ethanol
338
3390 1
340C   0.000000   0.000000   0.000000
341C   1.531000   0.000000   0.000000
342O   2.058000   1.198000   0.000000
343H  -0.390000   1.020000   0.000000
344H  -0.390000  -0.510000   0.884000
345H  -0.390000  -0.510000  -0.884000
346H   1.921000  -1.020000   0.000000
347H   1.921000   0.510000  -0.884000
348H   2.028000   1.604000   0.886000
349
350"#;
351
352    const LOG_SNIPPET: &str = r#"
353 Standard orientation:
354 ---------------------------------------------------------------------
355 Center     Atomic      Atomic             Coordinates (Angstroms)
356 Number     Number       Type             X           Y           Z
357 ---------------------------------------------------------------------
358      1          6           0        0.000000    0.000000    0.000000
359      2          6           0        1.531000    0.000000    0.000000
360      3          8           0        2.058000    1.198000    0.000000
361 ---------------------------------------------------------------------
362 SCF Done:  E(RB3LYP) =  -154.0987765432  A.U. after   12 cycles
363"#;
364
365    #[test]
366    fn parse_gjf_atom_count() {
367        let (mol, coords, charge, mult) = parse_gjf(ETHANOL_GJF).unwrap();
368        assert_eq!(mol.atom_count(), 9);
369        assert_eq!(coords.len(), 9);
370        assert_eq!(charge, 0);
371        assert_eq!(mult, 1);
372    }
373
374    #[test]
375    fn parse_gjf_first_coord() {
376        let (_, coords, _, _) = parse_gjf(ETHANOL_GJF).unwrap();
377        let (x, y, z) = coords[0];
378        assert!(x.abs() < 1e-6 && y.abs() < 1e-6 && z.abs() < 1e-6);
379    }
380
381    #[test]
382    fn write_gjf_roundtrip() {
383        let (mol, coords, charge, mult) = parse_gjf(ETHANOL_GJF).unwrap();
384        let out = write_gjf(&mol, &coords, charge, mult, "B3LYP/6-31G* opt", "Ethanol");
385        assert!(out.contains("# B3LYP/6-31G* opt"));
386        assert!(out.contains("0 1"));
387        let (mol2, coords2, _, _) = parse_gjf(&out).unwrap();
388        assert_eq!(mol2.atom_count(), 9);
389        assert_eq!(coords2.len(), 9);
390    }
391
392    #[test]
393    fn parse_gjf_charged_molecule() {
394        let gjf =
395            "# HF/STO-3G\n\nwater cation\n\n1 2\nO 0.0 0.0 0.0\nH 0.96 0.0 0.0\nH -0.24 0.93 0.0\n";
396        let (mol, _, charge, mult) = parse_gjf(gjf).unwrap();
397        assert_eq!(mol.atom_count(), 3);
398        assert_eq!(charge, 1);
399        assert_eq!(mult, 2);
400    }
401
402    #[test]
403    fn parse_log_geometry() {
404        let r = parse_gaussian_log(LOG_SNIPPET).unwrap();
405        assert_eq!(r.mol.atom_count(), 3);
406        assert_eq!(r.coords.len(), 3);
407        let (x, y, z) = r.coords[0];
408        assert!(x.abs() < 1e-6 && y.abs() < 1e-6 && z.abs() < 1e-6);
409    }
410
411    #[test]
412    fn parse_log_energy() {
413        let r = parse_gaussian_log(LOG_SNIPPET).unwrap();
414        let e = r.scf_energy.expect("energy should be present");
415        assert!((e - (-154.0987765432)).abs() < 1e-8);
416    }
417
418    #[test]
419    fn parse_log_last_orientation() {
420        // When multiple Standard orientation blocks exist, last one wins.
421        let input = format!(
422            "{}\n Standard orientation:\n ---\n ---\n      1          6           0        9.0    9.0    9.0\n ---\n{}",
423            LOG_SNIPPET, ""
424        );
425        let r = parse_gaussian_log(&input).unwrap();
426        // Last block has atom at (9, 9, 9)
427        assert!((r.coords[0].0 - 9.0).abs() < 1e-6);
428    }
429}
chematic_mol/gaussian.rs

chematic_mol/
gaussian.rs
