1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

//! Molecular symmetry element detection for linear systems.

use anyhow::{self, ensure, format_err};

use super::{PreSymmetry, Symmetry};
use crate::rotsym::RotationalSymmetry;
use crate::symmetry::symmetry_element::SIG;
use crate::symmetry::symmetry_element_order::{ElementOrder, ORDER_1, ORDER_2, ORDER_I};
use approx;
use log;
use nalgebra::Vector3;

impl Symmetry {
    /// Performs point-group detection analysis for a linear molecule.
    ///
    /// # Arguments
    ///
    /// * `presym` - A pre-symmetry-analysis structure containing information about the molecular
    /// system.
    /// * `tr` - A flag indicating if time reversal should also be considered. A time-reversed
    /// symmetry element will only be considered if its non-time-reversed version turns out to be
    /// not a symmetry element.
    pub(super) fn analyse_linear(
        &mut self,
        presym: &PreSymmetry,
        tr: bool,
    ) -> Result<(), anyhow::Error> {
        let (mois, principal_axes) = presym.recentred_molecule.calc_moi();

        ensure!(
            matches!(
                presym.rotational_symmetry,
                RotationalSymmetry::ProlateLinear
            ),
            "Unexpected rotational symmetry -- expected: {}, actual: {}",
            RotationalSymmetry::ProlateLinear,
            presym.rotational_symmetry
        );
        ensure!(
            presym
                .sea_groups
                .iter()
                .all(|sea_group| sea_group.len() <= 2),
            "Unexpected SEA groups of more than two atoms found for a linear molecule."
        );

        // C∞
        ensure!(
            approx::relative_eq!(
                mois[0],
                0.0,
                epsilon = presym.dist_threshold,
                max_relative = presym.dist_threshold
            ),
            "Unexpected non-zero smallest principal moment of inertia."
        );
        ensure!(
            self.add_proper(
                ORDER_I,
                &principal_axes[0],
                true,
                presym.dist_threshold,
                false
            ),
            "Expected C∞ axis not added."
        );
        if let Some(improper_kind) = presym.check_improper(
            &ElementOrder::Int(2),
            &Vector3::new(0.0, 0.0, 1.0),
            &SIG,
            tr,
        ) {
            // i
            log::debug!("Located an inversion centre.");
            ensure!(
                self.add_improper(
                    ORDER_2,
                    &Vector3::new(0.0, 0.0, 1.0),
                    false,
                    SIG,
                    None,
                    presym.dist_threshold,
                    improper_kind.contains_time_reversal()
                ),
                "Expected inversion centre not added."
            );

            // σh must exist if C∞ and i both exist.
            log::debug!("σh implied from C∞ and i.");
            let sigma_check = presym.check_improper(&ORDER_1, &principal_axes[0], &SIG, tr);
            ensure!(
                sigma_check.is_some(),
                "Expected σh implied by C∞ and i not found."
            );
            ensure!(
                self.add_improper(
                    ORDER_1,
                    &principal_axes[0],
                    true,
                    SIG,
                    Some("h".to_owned()),
                    presym.dist_threshold,
                    sigma_check
                        .ok_or_else(|| format_err!(
                            "Expected mirror plane implied by C∞ and i not found."
                        ))?
                        .contains_time_reversal(),
                ),
                "Expected σh implied by C∞ and i not added."
            );

            if let Some(proper_kind) = presym.check_proper(&ORDER_2, &principal_axes[1], tr) {
                // C2
                log::debug!("Located a C2 axis perpendicular to C∞.");
                self.add_proper(
                    ORDER_2,
                    &principal_axes[1],
                    true,
                    presym.dist_threshold,
                    proper_kind.contains_time_reversal(),
                );
                self.set_group_name("D∞h".to_owned());
            } else {
                // No C2
                log::debug!("No C2 axes perpendicular to C∞ found.");
                self.set_group_name("C∞h".to_owned());
            }
        } else {
            // No i
            log::debug!("No inversion centres found.");
            if let Some(improper_kind) =
                presym.check_improper(&ORDER_1, &principal_axes[1], &SIG, tr)
            {
                // σv
                log::debug!("Located a σv plane.");
                self.add_improper(
                    ORDER_1,
                    &principal_axes[1],
                    true,
                    SIG,
                    Some("v".to_owned()),
                    presym.dist_threshold,
                    improper_kind.contains_time_reversal(),
                );
                self.set_group_name("C∞v".to_owned());
            } else {
                // No σv
                log::debug!("No σv planes found.");
                self.set_group_name("C∞".to_owned());
            }
        }

        Ok(())
    }
}