gchemol_readwrite/formats/
gaussian_input.rs

1// [[file:../../gchemol-readwrite.note::*imports][imports:1]]
2use super::*;
3use super::parser::*;
4// imports:1 ends here
5
6// [[file:../../gchemol-readwrite.note::*header][header:1]]
7// Gaussian input file
8//
9// Reference
10// ---------
11// http://gaussian.com/input/?tabid=0
12// header:1 ends here
13
14// [[file:../../gchemol-readwrite.note::d068aeaf][d068aeaf]]
15// for reading/setting properties for `Atom` and `Molecule`
16const ATOM_KEY: &'static str = "gaussian-atom-private";
17const MOLE_KEY: &'static str = "gaussian-mole-private";
18
19/// Extra Atom properties pertinent to Gaussian input file
20#[derive(Debug, Clone, Serialize, Deserialize)]
21pub struct GaussianAtomInfo {
22    element_label: String,
23    mm_type: Option<String>,
24    mm_charge: Option<f64>,
25    frozen_code: Option<isize>,
26    position: [f64; 3],
27    properties: HashMap<String, String>,
28    oniom_info: OniomInfo,
29}
30
31impl Default for GaussianAtomInfo {
32    fn default() -> Self {
33        GaussianAtomInfo {
34            element_label: "C".to_owned(),
35            mm_type: None,
36            mm_charge: None,
37            frozen_code: None,
38            properties: HashMap::new(),
39            position: [0.0; 3],
40            oniom_info: OniomInfo::default(),
41        }
42    }
43}
44
45impl GaussianAtomInfo {
46    /// Set ONIOM layer for atom
47    pub fn set_oniom_layer(&mut self, layer: &str) {
48        self.oniom_info.layer = layer.to_owned().into();
49    }
50
51    /// Get ONIOM layer from atom
52    pub fn get_oniom_layer(&self) -> Option<&str> {
53        self.oniom_info.layer.as_ref().map(|x| x.as_str())
54    }
55
56    /// Set ONIOM link host atom
57    pub fn set_oniom_link_host(&mut self, link_host: usize) {
58        self.oniom_info.link_host = link_host.into();
59    }
60
61    /// Get ONIOM link host atom
62    pub fn get_oniom_link_host(&mut self) -> Option<usize> {
63        self.oniom_info.link_host
64    }
65
66    /// Set ONIOM link atom type
67    pub fn set_oniom_link_atom(&mut self, link_atom: &str) {
68        self.oniom_info.link_atom = link_atom.to_owned().into();
69    }
70
71    /// Get ONIOM link atom type
72    pub fn get_oniom_link_atom(&mut self) -> Option<&str> {
73        self.oniom_info.link_atom.as_ref().map(|x| x.as_str())
74    }
75
76    /// Attach/store extra properties to an `atom`.
77    pub fn attach(&self, atom: &mut Atom) {
78        let _ = atom.properties.store(ATOM_KEY, &self);
79    }
80
81    /// Extract extra properties from `atom`.
82    pub fn extract(atom: &Atom) -> Result<Self> {
83        let x = atom.properties.load(ATOM_KEY)?;
84        Ok(x)
85    }
86}
87// d068aeaf ends here
88
89// [[file:../../gchemol-readwrite.note::2d4f6dc2][2d4f6dc2]]
90/// Extra Molecule properties pertinent to Gaussian input file
91#[derive(Debug, Clone, Serialize, Deserialize, Default)]
92pub struct GaussianMoleculeInfo {
93    route: Option<String>,
94    link0: Option<String>,
95    title: Option<String>,
96    charge_and_multiplicity: Option<String>,
97}
98
99impl std::fmt::Display for GaussianMoleculeInfo {
100    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
101        let route = self.route.as_deref().unwrap_or("#p sp scf=tight HF/3-21G* geom=connect test").trim();
102        let link0 = self.link0.as_deref().unwrap_or("%nproc=1\n%mem=20MW").trim();
103        let title = self.title.as_deref().unwrap_or("Title Card Required").trim();
104        let csline = self.charge_and_multiplicity.as_deref().unwrap_or("0 1").trim();
105        write!(f, "{link0}\n{route}\n\n{title}\n\n{csline}")
106    }
107}
108
109impl GaussianMoleculeInfo {
110    /// Attach/store extra properties to `Molecule` `mol`.
111    pub fn attach(&self, mol: &mut Molecule) {
112        let _ = mol.properties.store(MOLE_KEY, &self);
113    }
114
115    /// Set route card for Gaussian input. NOTE: the "#" prefix must
116    /// be present in `route`.
117    pub fn set_route(&mut self, route: &str) {
118        self.route = route.to_owned().into();
119    }
120
121    /// Return the route in Gaussian input.
122    pub fn get_route(&self) -> Option<&str> {
123        self.route.as_deref()
124    }
125
126    /// Return the link0 commands in Gaussian input.
127    pub fn get_link0(&self) -> Option<&str> {
128        self.link0.as_deref()
129    }
130
131    /// Set route card for Gaussian input. NOTE: the "%" prefix must
132    /// be present in each link0 command as in `link0`.
133    pub fn set_link0(&mut self, link0: &str) {
134        self.route = link0.to_owned().into();
135    }
136
137    /// Extract extra properties from `Molecule`.
138    pub fn extract(mol: &Molecule) -> Result<Self> {
139        let x = mol.properties.load(MOLE_KEY)?;
140        Ok(x)
141    }
142}
143// 2d4f6dc2 ends here
144
145// [[file:../../gchemol-readwrite.note::576b287c][576b287c]]
146fn link0_cmd(s: &str) -> IResult<&str, &str> {
147    let prefix = tag("%");
148    do_parse!(s, prefix >> cmd: read_until_eol >> (cmd))
149}
150
151#[test]
152fn test_link0_cmd() {
153    let (_, cmd) = link0_cmd("%Mem=64MB\n").unwrap();
154    assert_eq!("Mem=64MB", cmd);
155
156    let (_, cmd) = link0_cmd("%save\n").unwrap();
157    assert_eq!("save", cmd);
158}
159
160fn link0_section(s: &str) -> IResult<&str, Vec<&str>> {
161    many0(link0_cmd)(s)
162}
163
164#[test]
165fn test_link0_section() {
166    let lines = "%chk=C5H12.chk
167%nproc=8
168%mem=5GB
169#p opt freq=noraman nosymm B3LYP/6-31+G** test geom=connect
170";
171
172    let (_, link0_cmds) = link0_section(lines).expect("gjf link0 section");
173    assert_eq!(3, link0_cmds.len());
174}
175// 576b287c ends here
176
177// [[file:../../gchemol-readwrite.note::*route section][route section:1]]
178fn blank_line(s: &str) -> IResult<&str, ()> {
179    do_parse!(s, space0 >> eol >> (()))
180}
181
182fn route_section(s: &str) -> IResult<&str, String> {
183    let mut pound = tag("#");
184    let mut print_level = opt(alt((tag_no_case("N"), tag_no_case("P"), tag_no_case("T"))));
185    let mut keywords = many_till(read_until_eol, blank_line);
186    do_parse!(
187        s,
188        pound >> print_level >>      // #P, #N, #T, #
189        space0 >> lines: keywords >> // opt B3LYP/6-31G* test
190        ({
191            lines.0.join(" ")
192        })
193    )
194}
195
196#[test]
197fn test_route_section() {
198    let lines = "#opt freq=noraman nosymm B3LYP/6-31+G** test geom=connect
199
200";
201    let x = route_section(lines).expect("gjf route section");
202
203    let lines = "#p opt freq=noraman nosymm
204B3LYP/6-31+G** test geom=connect
205
206";
207    let x = route_section(lines).expect("gjf route section multi-lines");
208}
209// route section:1 ends here
210
211// [[file:../../gchemol-readwrite.note::*title section][title section:1]]
212fn title_section(s: &str) -> IResult<&str, String> {
213    let mut title_lines = many_till(read_until_eol, blank_line);
214    do_parse!(
215        s,
216        lines: title_lines >> // xx
217        ({
218            lines.0.join(" ")
219        })
220    )
221}
222// title section:1 ends here
223
224// [[file:../../gchemol-readwrite.note::766dccd9][766dccd9]]
225// Specifies the net electric charge (a signed integer) and the spin
226// multiplicity (usually a positive integer)
227fn read_charge_and_spin_list(s: &str) -> IResult<&str, Vec<isize>> {
228    let mut values = separated_list1(space1, signed_digit);
229    do_parse!(s, space0 >> v: values >> eol >> (v))
230}
231
232#[test]
233fn test_charge_and_spin() {
234    let (_, x) = read_charge_and_spin_list("0	1\n").expect("gjf charge & spin");
235    assert_eq!(2, x.len());
236
237    let line = " 0 1 \t 0  1 -3 2 \n";
238    let (r, x) = read_charge_and_spin_list(line).expect("gjf charge & spin");
239    assert_eq!(6, x.len());
240}
241// 766dccd9 ends here
242
243// [[file:../../gchemol-readwrite.note::*element info][element info:1]]
244// single atom parameter entry
245// Fragment=1, or Iso=13, or Spin=3
246fn atom_param(s: &str) -> IResult<&str, (&str, &str)> {
247    separated_pair(alphanumeric1, tag("="), alphanumeric1)(s)
248}
249
250#[test]
251fn test_gjf_atom_param() {
252    let (_, (param, value)) = atom_param("fragment=1 ").unwrap();
253    assert_eq!("fragment", param);
254    assert_eq!("1", value);
255}
256
257// multiple property entries
258// (fragment=1,iso=13,spin=3)
259fn atom_params(s: &str) -> IResult<&str, Vec<(&str, &str)>> {
260    let params = separated_list1(tag(","), atom_param);
261    delimited(tag("("), params, tag(")"))(s)
262}
263
264#[test]
265fn test_gjf_atom_params() {
266    let (_, d) = atom_params("(fragment=1,iso=13,spin=3) ").expect("gjf atom properties");
267    assert_eq!(3, d.len())
268}
269
270// MM parameters, such as atom type and partial charge
271// -CA--0.25
272fn atom_mm_info(s: &str) -> IResult<&str, (&str, Option<f64>)> {
273    let mm_type = preceded(tag("-"), is_not("- \t"));
274    let mm_charge = preceded(tag("-"), double);
275    pair(mm_type, opt(mm_charge))(s)
276}
277
278#[test]
279fn test_gjf_atom_mm_info() {
280    let (_, (mm_type, mm_charge)) = atom_mm_info("-CA--0.25").unwrap();
281    assert_eq!("CA", mm_type);
282    assert_eq!(Some(-0.25), mm_charge);
283
284    let (_, (mm_type, mm_charge)) = atom_mm_info("-CA").unwrap();
285    assert_eq!("CA", mm_type);
286    assert_eq!(None, mm_charge);
287}
288// element info:1 ends here
289
290// [[file:../../gchemol-readwrite.note::0d3a42da][0d3a42da]]
291#[derive(Debug, Clone, Serialize, Deserialize, Default)]
292struct OniomInfo {
293    layer: Option<String>,
294    link_atom: Option<String>,
295    link_host: Option<usize>,
296    // TODO: scale factor
297}
298
299impl std::fmt::Display for OniomInfo {
300    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
301        if let Some(layer) = &self.layer {
302            write!(f, "{layer}")?;
303            if let Some(link_atom) = &self.link_atom {
304                write!(f, " {link_atom}")?;
305                if let Some(link_host) = &self.link_host {
306                    write!(f, " {link_host}")?;
307                }
308            }
309        }
310        Ok(())
311    }
312}
313
314impl OniomInfo {
315    fn from_str(s: &str) -> Self {
316        let mut oniom = Self::default();
317        let parts = s.split_whitespace().collect_vec();
318        if parts.is_empty() {
319            return oniom;
320        }
321        oniom.layer = parts[0].to_owned().into();
322        if parts.len() >= 2 {
323            oniom.link_atom = parts[1].to_owned().into();
324            if parts.len() >= 3 {
325                oniom.link_host = parts[2].parse().ok();
326            }
327        }
328        return oniom;
329    }
330}
331
332// C-CA--0.25   0   -4.703834   -1.841116   -0.779093 L
333// C-CA--0.25   0   -3.331033   -1.841116   -0.779093 L H-HA-0.1  3
334fn atom_oniom_info(s: &str) -> IResult<&str, OniomInfo> {
335    do_parse!(
336        s,
337        s: read_line >> // ONIOM layer: H, M, L
338        (OniomInfo::from_str(s))
339    )
340}
341
342#[test]
343fn test_gjf_atom_oniom_params() {
344    let (_, oniom) = atom_oniom_info(" L H-Ha-0.1 3\n").unwrap();
345    assert_eq!("L", oniom.layer.unwrap());
346    assert_eq!("H-Ha-0.1", oniom.link_atom.unwrap());
347    assert_eq!(Some(3), oniom.link_host);
348
349    let (_, oniom) = atom_oniom_info(" L\n").unwrap();
350    assert_eq!("L", oniom.layer.unwrap());
351}
352// 0d3a42da ends here
353
354// [[file:../../gchemol-readwrite.note::dfcbceb2][dfcbceb2]]
355use std::collections::HashMap;
356use std::iter::FromIterator;
357
358fn frozen_code(s: &str) -> IResult<&str, isize> {
359    terminated(signed_digit, space1)(s)
360}
361
362// How about this: C-CA--0.25(fragment=1,iso=13,spin=3) 0 0.0 1.2 3.4 H H-H_
363fn atom_line(s: &str) -> IResult<&str, GaussianAtomInfo> {
364    let mut mm_info = opt(atom_mm_info);
365    let mut params = opt(atom_params);
366    let mut frozen = opt(frozen_code);
367    do_parse!(
368        s,
369        space0 >> e: alphanumeric1      >> // element label
370        m: mm_info                      >> // MM type and MM partial charge
371        p: params                       >> // Optional atom info, such as fragment, nuclei props
372        space1                          >> // xx
373        f: frozen                       >> // frozen code: 0, -1
374        c: xyz_array                    >> // x, y, z coords
375        o: atom_oniom_info              >> // H H-H_
376        ({
377            GaussianAtomInfo {
378                element_label: e.to_owned(),
379                mm_type: m.and_then(|x| Some(x.0.to_owned())),
380                mm_charge: m.and_then(|x| x.1),
381                frozen_code: f,
382                position: c,
383                properties: if p.is_some() {p.unwrap().into_iter().map(|(k, v)| (k.to_owned(), v.to_owned())).collect()} else {HashMap::new()},
384                oniom_info: o,
385                ..Default::default()
386            }
387        })
388    )
389}
390// dfcbceb2 ends here
391
392// [[file:../../gchemol-readwrite.note::58144795][58144795]]
393#[test]
394fn test_gjf_atom_line() {
395    let line = "C-CA--0.25(fragment=1,iso=13,spin=3)  -1 0.00   0.00   0.00 L H-HA-0.1  3\n";
396    let (_, ga) = atom_line(line).expect("full oniom atom line");
397    assert_eq!(ga.oniom_info.layer.unwrap(), "L");
398
399    let line = "C-CA--0.25(fragment=1,iso=13,spin=3) -1 0.00   0.00   0.00 L \n";
400    let (_, ga) = atom_line(line).expect("oniom without link atom");
401    assert_eq!(ga.oniom_info.layer.unwrap(), "L");
402    assert_eq!(ga.frozen_code, Some(-1));
403
404    let line = "C-CA--0.25(fragment=1,iso=13,spin=3) -1 0.00   0.00   0.00\n";
405    let (_, ga) = atom_line(line).expect("frozen + xyz");
406    assert_eq!(ga.frozen_code, Some(-1));
407
408    let line = "C-CA--0.25(fragment=1,iso=13,spin=3) 0.00   0.00   0.00\n";
409    let (_, ga) = atom_line(line).expect("missing frozen code");
410    assert_eq!(ga.mm_charge, Some(-0.25));
411    assert_eq!(ga.mm_type.unwrap(), "CA");
412    assert_eq!(ga.properties["fragment"], "1");
413
414    let (_, ga) = atom_line("C-CA--0.25  0.00   0.00   0.00\n").expect("no key-value params");
415    assert_eq!(ga.mm_charge, Some(-0.25));
416    assert_eq!(ga.mm_type.unwrap(), "CA");
417
418    let (_, ga) = atom_line("C-C_3 0.00 0.00 0.00\n").expect("no mm charge");
419    assert_eq!(ga.mm_type.unwrap(), "C_3");
420    assert_eq!(ga.mm_charge, None);
421
422    let line = " C12(fragment=1)  0.00   0.00   0.00\n";
423    let (_, ga) = atom_line(line).expect("key-value params only");
424    assert_eq!(ga.properties["fragment"], "1");
425    assert_eq!(ga.position[0], 0.0);
426
427    let line = " C12  0.00   0.00   0.00\n";
428    let (_, ga) = atom_line(line).expect("simple");
429    assert_eq!(ga.position[0], 0.0);
430}
431// 58144795 ends here
432
433// [[file:../../gchemol-readwrite.note::*connectivity specs][connectivity specs:1]]
434// Connectivity specs example:
435//
436// 1 2 1.0 3 1.0 4 1.0 5 1.0
437// 2
438// 3
439fn bond_item(s: &str) -> IResult<&str, (usize, f64)> {
440    do_parse!(
441        s,
442        space1 >> n: unsigned_digit >> // bond atom number: 2
443        space1 >> o: double         >> // bond order 1.0
444        (n, o)
445    )
446}
447
448fn read_connect_line(s: &str) -> IResult<&str, Vec<(usize, usize, Bond)>> {
449    let mut read_bond_items = many0(bond_item);
450    do_parse!(
451        s,
452        space0 >> n: unsigned_digit    >> // host atom number
453        others: read_bond_items >> eol >> // bond atoms and bond orders
454        ({
455            let mut bonds = vec![];
456            for (index2, order) in others {
457                bonds.push((n, index2, Bond::single()));
458            }
459
460            bonds
461        })
462    )
463}
464
465#[test]
466fn test_gjf_connectivity() {
467    let (_, x) = read_connect_line(" 1 2 1.0 3 1.0 4 1.0 5 1.0\n").expect("gjf connectivity");
468    assert_eq!(4, x.len());
469}
470// connectivity specs:1 ends here
471
472// [[file:../../gchemol-readwrite.note::fc56d173][fc56d173]]
473fn parse_molecule_extra(s: &str) -> IResult<&str, GaussianMoleculeInfo> {
474    let mut link0 = opt(link0_section);
475    do_parse!(
476        s,
477        l: link0         >> // Link0 commands, which is optional section
478        r: route_section >> // route section contains job keywords
479        t: title_section >> // molecular title
480        ({
481            let mut extra  = GaussianMoleculeInfo::default();
482            // add prefix "%" for each link0 cmd
483            extra.link0 = l.map(|x| x.iter().map(|s| format!("%{s}")).join("\n"));
484            extra.route = format!("#P {r}").to_owned().into();
485            extra.title = t.to_owned().into();
486            extra
487        })
488    )
489}
490
491fn parse_molecule_specs<'a>(s: &'a str, extra: &mut GaussianMoleculeInfo) -> IResult<&'a str, Molecule> {
492    let mut read_atoms = many1(atom_line);
493    let mut read_bonds = opt(many1(read_connect_line));
494    do_parse!(
495        s,
496        cslist: read_charge_and_spin_list >> // charge and spin multipy
497        atoms: read_atoms                 >> // atom symbol, coordinates, atom type, ...
498        blank_line                        >>
499        bonds: read_bonds                 >> // connectivity section
500        ({
501            // set charge and multiplicity
502            extra.charge_and_multiplicity = cslist.iter().map(|x| x.to_string()).join(" ").into();
503            let mut mol = Molecule::default();
504            let mut lat_vectors = vec![];
505            let atoms = atoms.into_iter().filter_map(|x| {
506                // Handle dummy TV atoms (transitional vector)
507                if &x.element_label.to_uppercase() == "TV" {
508                    debug!("found TV dummy atom");
509                    lat_vectors.push(x.position.clone());
510                    None
511                } else {
512                    // set atom freezing flag
513                    let mut a = Atom::new(x.element_label.clone(), x.position.clone());
514                    if let Some(f) = x.frozen_code {
515                        if f == -1 {
516                            a.set_freezing([true; 3]);
517                        }
518                    }
519                    // store extra atom properties such as ONIOM layer or fragment index.
520                    x.attach(&mut a);
521                    Some(a)
522                }
523            });
524            mol.add_atoms_from((1..).zip(atoms));
525            if lat_vectors.len() == 3 {
526                let lat = Lattice::new([lat_vectors[0], lat_vectors[1], lat_vectors[2]]);
527                mol.set_lattice(lat);
528            } else if !lat_vectors.is_empty() {
529                error!("Expect 3, but found {} TV atoms.", lat_vectors.len());
530            }
531
532            // handle bonds
533            if let Some(bonds) = bonds {
534                for (u, v, b) in bonds.into_iter().flatten() {
535                    mol.add_bond(u, v, b);
536                }
537            }
538            mol
539        })
540    )
541}
542
543// FIXME: how about gaussian extra input
544pub fn parse_molecule(s: &str) -> Result<Molecule> {
545    let (r, mut mol_extra) = parse_molecule_extra(s).map_err(|e| format_err!("{}", e))?;
546    // We replace comma with space in molecular specification part for easy
547    // parsing.
548    let r = r.replace(",", " ");
549    let (r, mut mol) = parse_molecule_specs(&r, &mut mol_extra)
550        .map_err(|e| format_err!("parsing gaussian input molecule specification failure: {:}.\n Input stream: {}", e, r))?;
551    // NOTE: title card could be verbose
552    // if let Some(title) = &mol_extra.title {
553    //     mol.set_title(&title);
554    // }
555    mol_extra.attach(&mut mol);
556
557    Ok(mol)
558}
559
560#[test]
561fn test_parse_gaussian_molecule() {
562    let txt = "%chk=C5H12.chk
563%nproc=8
564%mem=5GB
565#p opt freq=noraman nosymm B3LYP/6-31+G** test geom=connect
566
567Title Card
568Required
569
5700 1
571 C(Fragment=1)   -1         -1.29639700         -0.54790000         -0.04565800 L
572 H(Fragment=1)    0         -0.94903500         -1.58509500         -0.09306500 L
573 H(Fragment=1)    0         -0.93491200         -0.03582400         -0.94211800 L
574 H(Fragment=1)    0         -2.39017900         -0.56423400         -0.09090100 L
575 C(Fragment=2)    0         -0.80594100          0.14211700          1.23074100 L
576 H(Fragment=2)    0         -1.13863100          1.18750400          1.23095600 L
577 H(Fragment=2)    0          0.29065200          0.17299300          1.23044100 L
578 C(Fragment=3)    0         -1.29047900         -0.54051400          2.51480900 L
579 H(Fragment=3)    0         -0.95681000         -1.58684300          2.51564600 L
580 H(Fragment=3)    0         -2.38822700         -0.57344700          2.51397600 L
581 C(Fragment=4)    0         -0.80793200          0.14352700          3.79887800 L
582 H(Fragment=4)    0         -1.14052400          1.18894500          3.79694800 L
583 H(Fragment=4)    0          0.28866200          0.17430200          3.80089400 L
584 C(Fragment=5)    0         -1.30048800         -0.54500600          5.07527000 L
585 H(Fragment=5)    0         -0.95334000         -1.58219500          5.12438500 L
586 H(Fragment=5)    0         -0.94034400         -0.03198400          5.97172800 L
587 H(Fragment=5)    0         -2.39434000         -0.56114800          5.11880700 L
588
589 1 2 1.0 3 1.0 4 1.0 5 1.0
590 2
591 3
592 4
593 5 6 1.0 7 1.0 8 1.0
594 6
595 7
596 8 9 1.0 10 1.0 11 1.0
597 9
598 10
599 11 12 1.0 13 1.0 14 1.0
600 12
601 13
602 14 15 1.0 16 1.0 17 1.0
603 15
604 16
605 17
606";
607
608    let mol = parse_molecule(txt).expect("gjf molecule");
609    assert_eq!(17, mol.natoms());
610    assert_eq!(16, mol.nbonds());
611}
612// fc56d173 ends here
613
614// [[file:../../gchemol-readwrite.note::5605d45c][5605d45c]]
615// TODO: atom properties
616fn format_atom(a: &Atom) -> String {
617    let [x, y, z] = a.position();
618    let symbol = a.symbol();
619    let fcode = if a.freezing() == [true; 3] { -1 } else { 0 };
620    let part = format!(" {symbol:15} {fcode:2} {x:14.8} {y:14.8} {z:14.8}");
621
622    // format ONIOM layer, link atom, link host
623    if let Ok(extra) = GaussianAtomInfo::extract(&a) {
624        let oniom = extra.oniom_info;
625        format!("{part} {oniom}\n")
626    } else {
627        format!("{part}\n")
628    }
629}
630
631// string representation in gaussian input file format
632fn format_molecule(mol: &Molecule) -> String {
633    let mut lines = String::new();
634
635    let mut extra = GaussianMoleculeInfo::extract(mol).unwrap_or_default();
636    if extra.title.is_none() {
637        extra.title = mol.title().into();
638    }
639
640    lines.push_str(&format!("{extra}\n"));
641
642    for (_, a) in mol.atoms() {
643        let line = format_atom(&a);
644        lines.push_str(&line);
645    }
646
647    // crystal vectors
648    if let Some(lattice) = mol.lattice {
649        // let va = lattice.vector_a();
650        // let vb = lattice.vector_b();
651        // let vc = lattice.vector_c();
652        for l in lattice.vectors().iter() {
653            lines.push_str(&format!(" TV              {:14.8}{:14.8}{:14.8}\n", l.x, l.y, l.z));
654        }
655    }
656
657    // connectivity
658    lines.push_str("\n");
659    let mut map = HashMap::new();
660    for (i, j, b) in mol.bonds() {
661        let mut neighbors = map.entry(i).or_insert(vec![]);
662        neighbors.push((j, b.order()));
663    }
664    for (i, a) in mol.atoms() {
665        let mut line = format!("{:<5}", i);
666        if let Some(neighbors) = map.get(&i) {
667            for (j, o) in neighbors {
668                line.push_str(&format!(" {:<} {:<.1}", j, o));
669            }
670        }
671        lines.push_str(&format!("{}\n", line));
672    }
673
674    lines.push_str("\n");
675    lines
676}
677// 5605d45c ends here
678
679// [[file:../../gchemol-readwrite.note::d13a6041][d13a6041]]
680#[derive(Clone, Copy, Debug)]
681/// Gaussian input file
682pub struct GaussianInputFile();
683
684impl GaussianInputFile {
685    /// Extract extra atom information from `atom`, which can be
686    /// attached to `Atom` later.
687    pub fn extra_atom_info(atom: &Atom) -> GaussianAtomInfo {
688        GaussianAtomInfo::extract(atom).ok().unwrap_or_default()
689    }
690
691    /// Extract extra molecule information from `mol`, which can be
692    /// attached to `Molecule` later.
693    pub fn extra_molecule_info(mol: &Molecule) -> GaussianMoleculeInfo {
694        GaussianMoleculeInfo::extract(mol).ok().unwrap_or_default()
695    }
696}
697
698/// References
699/// http://gaussian.com/input/?tabid=0
700impl ChemicalFile for GaussianInputFile {
701    fn ftype(&self) -> &str {
702        "gaussian/input"
703    }
704
705    fn possible_extensions(&self) -> Vec<&str> {
706        vec![".gjf", ".com", ".gau"]
707    }
708
709    fn format_molecule(&self, mol: &Molecule) -> Result<String> {
710        Ok(format_molecule(mol))
711    }
712}
713
714impl ParseMolecule for GaussianInputFile {
715    fn parse_molecule(&self, input: &str) -> Result<Molecule> {
716        parse_molecule(input)
717    }
718}
719// d13a6041 ends here
720
721// [[file:../../gchemol-readwrite.note::2530a669][2530a669]]
722impl ReadPart for GaussianInputFile {
723    fn read_next(&self, context: ReadContext) -> ReadAction {
724        Terminated(|line: &str| {
725            let link1 = "--link1--\n";
726            line.to_lowercase() == link1
727        })
728        .read_next(context)
729    }
730}
731
732impl GaussianInputFile {
733    pub fn partitions<R: BufRead + Seek>(&self, mut r: TextReader<R>) -> Result<impl Iterator<Item = String>> {
734        Ok(r.partitions(*self))
735    }
736}
737// 2530a669 ends here
738
739// [[file:../../gchemol-readwrite.note::457b015d][457b015d]]
740#[test]
741fn test_gaussian_input_file() -> Result<()> {
742    let s = gut::fs::read_file("tests/files/gaussian/test0769.com")?;
743    let mut mol = GaussianInputFile().parse_molecule(&s)?;
744    let mut atoms_high_layer = vec![];
745    for (i, a) in mol.atoms() {
746        let extra = GaussianInputFile::extra_atom_info(&a);
747        match extra.get_oniom_layer() {
748            Some("H") => atoms_high_layer.push(i),
749            _ => {}
750        }
751    }
752    assert_eq!(atoms_high_layer.len(), 5);
753
754    // set molecule level attributes
755    let mut extra = GaussianInputFile::extra_molecule_info(&mol);
756    extra.set_route("# opt B3LYP/6-31G**");
757    extra.attach(&mut mol);
758    let fs = GaussianInputFile().format_molecule(&mol)?;
759    let mol_ = GaussianInputFile().parse_molecule(&fs)?;
760    assert_eq!(mol.natoms(), mol_.natoms());
761    let extra_ = GaussianInputFile::extra_molecule_info(&mol);
762    assert_eq!(extra.get_route(), extra_.get_route());
763
764    Ok(())
765}
766// 457b015d ends here
gchemol_readwrite/formats/gaussian_input.rs

gchemol_readwrite/formats/
gaussian_input.rs
