qsym2/symmetry/symmetry_core/symmetry_core_asymmetric.rs
1//! Molecular symmetry element detection for asymmetric tops.
2
3use anyhow::{self, ensure, format_err};
4use itertools::Itertools;
5use log;
6use nalgebra::Vector3;
7
8use crate::rotsym::RotationalSymmetry;
9use crate::symmetry::symmetry_core::_search_proper_rotations;
10use crate::symmetry::symmetry_element::{SymmetryElement, ROT, SIG, TRROT, TRSIG};
11use crate::symmetry::symmetry_element_order::{ORDER_1, ORDER_2};
12
13use super::{PreSymmetry, Symmetry};
14
15impl Symmetry {
16 /// Performs point-group detection analysis for an asymmetric-top molecule.
17 ///
18 /// The possible asymmetric top point groups are:
19 ///
20 /// * $`\mathcal{C}_{1}`$ and $`\mathcal{C}_{2}`$,
21 /// * $`\mathcal{C}_{2v}`$,
22 /// * $`\mathcal{C}_{2h}`$,
23 /// * $`\mathcal{C}_{s}`$,
24 /// * $`\mathcal{D}_{2}`$,
25 /// * $`\mathcal{D}_{2h}`$, and
26 /// * $`\mathcal{C}_{i}`$.
27 ///
28 /// These are all Abelian groups.
29 ///
30 /// # Arguments
31 ///
32 /// * `presym` - A pre-symmetry-analysis structure containing information about the molecular
33 /// system.
34 /// * `tr` - A flag indicating if time reversal should also be considered. A time-reversed
35 /// symmetry element will only be considered if its non-time-reversed version turns out to be
36 /// not a symmetry element.
37 #[allow(clippy::too_many_lines)]
38 pub(super) fn analyse_asymmetric(
39 &mut self,
40 presym: &PreSymmetry,
41 tr: bool,
42 ) -> Result<(), anyhow::Error> {
43 let (_mois, principal_axes) = presym.recentred_molecule.calc_moi();
44
45 ensure!(
46 matches!(
47 presym.rotational_symmetry,
48 RotationalSymmetry::AsymmetricPlanar | RotationalSymmetry::AsymmetricNonPlanar
49 ),
50 "Unexpected rotational symmetry -- expected: {} or {}, actual: {}",
51 RotationalSymmetry::AsymmetricPlanar,
52 RotationalSymmetry::AsymmetricNonPlanar,
53 presym.rotational_symmetry
54 );
55
56 _search_proper_rotations(presym, self, true, tr)?;
57 log::debug!("Proper elements found: {:?}", self.get_elements(&ROT));
58 log::debug!(
59 "Time-reversed proper elements found: {:?}",
60 self.get_elements(&TRROT)
61 );
62
63 // Classify into point groups
64 let count_c2 = self
65 .get_proper(&ORDER_2)
66 .map_or(0, |proper_elements| proper_elements.len());
67 ensure!(
68 count_c2 == 0 || count_c2 == 1 || count_c2 == 3,
69 "Unexpected number of C2 axes: {count_c2}."
70 );
71
72 let max_ord = self.get_max_proper_order();
73
74 if count_c2 == 3 {
75 // Dihedral, either D2h or D2.
76 log::debug!("Dihedral family (asymmetric top).");
77 ensure!(
78 max_ord == ORDER_2,
79 "Unexpected principal order -- expected: 2, actual: {max_ord}."
80 );
81
82 // Principal axis, which is C2, is also a generator.
83 // If the group is a black-white magnetic group, then one C2 axis is non-time-reversed,
84 // while the other two are. Hence, one C2 generator is non-time-reversed, and the other
85 // must be. We take care of this by sorting `c2s` to put any non-time-reversed elements
86 // first.
87 let mut c2s = self
88 .get_proper(&ORDER_2)
89 .ok_or_else(|| format_err!("No C2 elements found."))?
90 .into_iter()
91 .cloned()
92 .collect_vec();
93 c2s.sort_by_key(SymmetryElement::contains_time_reversal);
94 let mut c2s = c2s.into_iter();
95 let c2 = c2s
96 .next()
97 .ok_or_else(|| format_err!("No C2 elements found."))?;
98 self.add_proper(
99 max_ord,
100 c2.raw_axis(),
101 true,
102 presym.dist_threshold,
103 c2.contains_time_reversal(),
104 );
105
106 // One other C2 axis is also a generator.
107 let another_c2 = c2s
108 .next()
109 .ok_or_else(|| format_err!("No more C2s found."))?;
110 self.add_proper(
111 max_ord,
112 another_c2.raw_axis(),
113 true,
114 presym.dist_threshold,
115 another_c2.contains_time_reversal(),
116 );
117
118 let z_vec = Vector3::new(0.0, 0.0, 1.0);
119 if let Some(improper_kind) = presym.check_improper(&ORDER_2, &z_vec, &SIG, tr) {
120 // Inversion centre, D2h
121 log::debug!("Located an inversion centre.");
122 self.set_group_name("D2h".to_owned());
123 self.add_improper(
124 ORDER_2,
125 &z_vec,
126 false,
127 SIG,
128 None,
129 presym.dist_threshold,
130 improper_kind.contains_time_reversal(),
131 );
132
133 // Add remaining mirror planes, each of which is
134 // perpendicular to a C2 axis.
135 let c2s = self
136 .get_proper(&ORDER_2)
137 .ok_or_else(|| format_err!("No C2 elements found."))?
138 .into_iter()
139 .cloned()
140 .collect_vec();
141 for c2 in &c2s {
142 let improper_check = presym.check_improper(&ORDER_1, c2.raw_axis(), &SIG, tr);
143 ensure!(
144 improper_check.is_some(),
145 "Expected improper element not found."
146 );
147 self.add_improper(
148 ORDER_1,
149 c2.raw_axis(),
150 false,
151 SIG,
152 None,
153 presym.dist_threshold,
154 improper_check
155 .ok_or_else(|| {
156 format_err!(
157 "Expected mirror plane perpendicular to `{}` not found.",
158 c2.raw_axis()
159 )
160 })?
161 .contains_time_reversal(),
162 );
163 }
164 let principal_element_axis = *self.get_proper_principal_element().raw_axis();
165 let improper_check =
166 presym.check_improper(&ORDER_1, &principal_element_axis, &SIG, tr);
167 ensure!(
168 improper_check.is_some(),
169 "Expected improper element not found."
170 );
171 self.add_improper(
172 ORDER_1,
173 &principal_element_axis,
174 true,
175 SIG,
176 None,
177 presym.dist_threshold,
178 improper_check
179 .ok_or_else(||
180 format_err!("Expected mirror plane perpendicular to the principal axis not found.")
181 )?
182 .contains_time_reversal(),
183 );
184 } else {
185 // Chiral, D2
186 self.set_group_name("D2".to_owned());
187 }
188 } else if count_c2 == 1 {
189 // Non-dihedral, either C2, C2v, or C2h
190 log::debug!("Non-dihedral family (asymmetric top).");
191 ensure!(
192 max_ord == ORDER_2,
193 "Unexpected principal order -- expected: 2, actual: {max_ord}."
194 );
195
196 // Principal axis, which is C2, is also a generator.
197 let c2s = self
198 .get_proper(&ORDER_2)
199 .ok_or_else(|| format_err!("No C2 elements found."))?;
200 let c2 = (*c2s
201 .iter()
202 .next()
203 .ok_or_else(|| format_err!("No C2 elements found."))?)
204 .clone();
205 self.add_proper(
206 max_ord,
207 c2.raw_axis(),
208 true,
209 presym.dist_threshold,
210 c2.contains_time_reversal(),
211 );
212
213 let z_vec = Vector3::new(0.0, 0.0, 1.0);
214 if let Some(improper_kind) = presym.check_improper(&ORDER_2, &z_vec, &SIG, tr) {
215 // Inversion centre, C2h
216 log::debug!("Located an inversion centre.");
217 self.add_improper(
218 ORDER_2,
219 &z_vec,
220 false,
221 SIG,
222 None,
223 presym.dist_threshold,
224 improper_kind.contains_time_reversal(),
225 );
226 self.set_group_name("C2h".to_owned());
227
228 // There is one σh.
229 let c2 = (*self
230 .get_proper(&ORDER_2)
231 .ok_or_else(|| format_err!("No C2 elements found."))?
232 .iter()
233 .next()
234 .ok_or_else(|| format_err!("No C2 elements found."))?)
235 .clone();
236
237 let improper_check = presym.check_improper(&ORDER_1, c2.raw_axis(), &SIG, tr);
238 ensure!(
239 improper_check.is_some(),
240 "Expected improper element not found."
241 );
242 self.add_improper(
243 ORDER_1,
244 c2.raw_axis(),
245 false,
246 SIG,
247 Some("h".to_owned()),
248 presym.dist_threshold,
249 improper_check
250 .as_ref()
251 .ok_or_else(|| {
252 format_err!(
253 "Expected mirror plane perpendicular to {} not found.",
254 c2.raw_axis()
255 )
256 })?
257 .contains_time_reversal(),
258 );
259 self.add_improper(
260 ORDER_1,
261 c2.raw_axis(),
262 true,
263 SIG,
264 Some("h".to_owned()),
265 presym.dist_threshold,
266 improper_check
267 .ok_or_else(|| {
268 format_err!(
269 "Expected mirror plane perpendicular to {} not found.",
270 c2.raw_axis()
271 )
272 })?
273 .contains_time_reversal(),
274 );
275 } else {
276 // No inversion centres.
277 // Locate σ planes
278 let mut count_sigma = 0;
279 if matches!(
280 presym.rotational_symmetry,
281 RotationalSymmetry::AsymmetricPlanar
282 ) {
283 // Planar system. The plane of the system (perpendicular to the highest-MoI
284 // principal axis) might be a symmetry element: time-reversed in the presence of
285 // a magnetic field (which must also lie in this plane), or both in the absence
286 // of a magnetic field.
287 if let Some(improper_kind) =
288 presym.check_improper(&ORDER_1, &principal_axes[2], &SIG, tr)
289 {
290 if presym.recentred_molecule.magnetic_atoms.is_some() {
291 ensure!(
292 improper_kind.contains_time_reversal(),
293 "Expected time-reversed improper element not found."
294 );
295 }
296 count_sigma += u32::from(self.add_improper(
297 ORDER_1,
298 &principal_axes[2],
299 false,
300 SIG,
301 Some("v".to_owned()),
302 presym.dist_threshold,
303 improper_kind.contains_time_reversal(),
304 ));
305 }
306 }
307
308 let sea_groups = &presym.sea_groups;
309 for sea_group in sea_groups.iter() {
310 if count_sigma == 2 {
311 break;
312 }
313 if sea_group.len() < 2 {
314 continue;
315 }
316 for atom2s in sea_group.iter().combinations(2) {
317 if count_sigma == 2 {
318 break;
319 }
320 let normal = (atom2s[0].coordinates.coords - atom2s[1].coordinates.coords)
321 .normalize();
322 if let Some(improper_kind) =
323 presym.check_improper(&ORDER_1, &normal, &SIG, tr)
324 {
325 if c2.contains_time_reversal()
326 && !improper_kind.contains_time_reversal()
327 {
328 log::debug!("The C2 axis is actually θ·C2. The non-time-reversed σv will be assigned as σh.");
329 count_sigma += u32::from(self.add_improper(
330 ORDER_1,
331 &normal,
332 false,
333 SIG,
334 Some("h".to_owned()),
335 presym.dist_threshold,
336 improper_kind.contains_time_reversal(),
337 ));
338 } else {
339 count_sigma += u32::from(self.add_improper(
340 ORDER_1,
341 &normal,
342 false,
343 SIG,
344 Some("v".to_owned()),
345 presym.dist_threshold,
346 improper_kind.contains_time_reversal(),
347 ));
348 }
349 }
350 }
351 }
352
353 log::debug!(
354 "Located {} σ ({} σv and {} σh).",
355 count_sigma,
356 self.get_sigma_elements("v")
357 .map_or(0, |sigmavs| sigmavs.len()),
358 self.get_sigma_elements("h")
359 .map_or(0, |sigmavs| sigmavs.len()),
360 );
361 if count_sigma == 2 {
362 self.set_group_name("C2v".to_owned());
363
364 // In C2v, one of the σ's is also a generator. We prioritise the
365 // non-time-reversed one as the generator.
366 let mut sigmas = self
367 .get_sigma_elements("v")
368 .or_else(|| {
369 log::debug!("No σv found. Searching for σh instead.");
370 self.get_sigma_elements("h")
371 })
372 .ok_or_else(|| format_err!("No σv nor σh found."))?
373 .into_iter()
374 .chain(self.get_sigma_elements("h").unwrap_or_default().into_iter())
375 .cloned()
376 .collect_vec();
377 sigmas.sort_by_key(SymmetryElement::contains_time_reversal);
378 let sigma = sigmas
379 .first()
380 .ok_or_else(|| format_err!("No σv or σh found."))?;
381 self.add_improper(
382 ORDER_1,
383 sigma.raw_axis(),
384 true,
385 SIG,
386 Some(sigma.additional_subscript.clone()),
387 presym.dist_threshold,
388 sigma.contains_time_reversal(),
389 );
390 } else {
391 ensure!(
392 count_sigma == 0,
393 "Unexpected number of σ mirror planes: {count_sigma}."
394 );
395 self.set_group_name("C2".to_owned());
396 }
397 }
398 } else {
399 // No C2 axes, so either C1, Ci, or Cs
400 log::debug!("No C2 axes found.");
401 let z_vec = Vector3::new(0.0, 0.0, 1.0);
402 if let Some(improper_kind) = presym.check_improper(&ORDER_2, &z_vec, &SIG, tr) {
403 log::debug!("Located an inversion centre.");
404 self.set_group_name("Ci".to_owned());
405 self.add_improper(
406 ORDER_2,
407 &z_vec,
408 false,
409 SIG,
410 None,
411 presym.dist_threshold,
412 improper_kind.contains_time_reversal(),
413 );
414 self.add_improper(
415 ORDER_2,
416 &z_vec,
417 true,
418 SIG,
419 None,
420 presym.dist_threshold,
421 improper_kind.contains_time_reversal(),
422 );
423 } else {
424 log::debug!("No inversion centres found.");
425 // Locate mirror planes
426 let sea_groups = &presym.sea_groups;
427 let mut count_sigma = 0;
428 for sea_group in sea_groups.iter() {
429 if count_sigma > 0 {
430 break;
431 }
432 if sea_group.len() < 2 {
433 continue;
434 }
435 for atom2s in sea_group.iter().combinations(2) {
436 let normal = (atom2s[0].coordinates.coords - atom2s[1].coordinates.coords)
437 .normalize();
438 if let Some(improper_kind) =
439 presym.check_improper(&ORDER_1, &normal, &SIG, tr)
440 {
441 count_sigma += u32::from(self.add_improper(
442 ORDER_1,
443 &normal,
444 false,
445 SIG,
446 None,
447 presym.dist_threshold,
448 improper_kind.contains_time_reversal(),
449 ));
450 }
451 }
452 }
453
454 if count_sigma == 0
455 && matches!(
456 presym.rotational_symmetry,
457 RotationalSymmetry::AsymmetricPlanar
458 )
459 {
460 log::debug!("Planar molecule based on MoIs but no σ found from SEA groups.");
461 log::debug!("Locating the planar mirror plane based on MoIs...");
462 let sigma_check = presym.check_improper(&ORDER_1, &principal_axes[2], &SIG, tr);
463 if sigma_check.is_some() {
464 ensure!(
465 self.add_improper(
466 ORDER_1,
467 &principal_axes[2],
468 false,
469 SIG,
470 None,
471 presym.dist_threshold,
472 sigma_check
473 .ok_or_else(|| format_err!(
474 "Expected {}mirror plane perpendicular to the highest-MoI principal axis not found.",
475 if tr { "time-reversed " } else { "" }
476 ))?
477 .contains_time_reversal(),
478 ),
479 "Failed to add {}mirror plane perpendicular to the highest-MoI principal axis.",
480 if tr { "time-reversed " } else { "" }
481 );
482 log::debug!(
483 "Located one planar {}mirror plane based on MoIs.",
484 if tr { "time-reversed " } else { "" }
485 );
486 count_sigma += 1;
487 } else {
488 assert!(!tr, "The only way for a planar molecule to not have a planar mirror plane is when a magnetic field is present but time reversal is not considered.");
489 log::debug!("No additional planar mirror planes found.");
490 }
491
492 // Old algorithm
493 // for atom3s in presym.recentred_molecule.atoms.iter().combinations(3) {
494 // let normal = (atom3s[1].coordinates.coords - atom3s[0].coordinates.coords)
495 // .cross(&(atom3s[2].coordinates.coords - atom3s[0].coordinates.coords));
496 // if normal.norm() < presym.dist_threshold {
497 // if let Some(e_atoms) = &presym.recentred_molecule.electric_atoms {
498 // let normal = (atom3s[1].coordinates.coords
499 // - atom3s[0].coordinates.coords)
500 // .cross(
501 // &(e_atoms[1].coordinates.coords
502 // - e_atoms[0].coordinates.coords),
503 // );
504 // } else {
505 // continue;
506 // }
507 // }
508 // let normal = normal.normalize();
509 // if presym.check_improper(&ORDER_1, &normal, &SIG) {
510 // count_sigma += self.add_improper(
511 // ORDER_1,
512 // normal,
513 // false,
514 // SIG.clone(),
515 // None,
516 // presym.dist_threshold,
517 // ) as u32;
518 // break;
519 // }
520 // }
521 // assert_eq!(count_sigma, 1);
522 // log::debug!("Located one planar mirror plane based on MoIs.");
523 }
524
525 log::debug!("Located {} σ.", count_sigma);
526 if count_sigma > 0 {
527 ensure!(
528 count_sigma == 1,
529 "Unexpected number of σ mirror planes: {count_sigma}."
530 );
531 let old_sigmas = self
532 .get_elements_mut(&SIG)
533 .and_then(|sigmas| sigmas.remove(&ORDER_1))
534 .or_else(|| {
535 self.get_elements_mut(&TRSIG)
536 .and_then(|tr_sigmas| tr_sigmas.remove(&ORDER_1))
537 })
538 .ok_or_else(|| {
539 format_err!("No normal or time-reversed mirror planes found.")
540 })?;
541 ensure!(
542 old_sigmas.len() == 1,
543 "Unexpected number of old σ mirror planes: {}.",
544 old_sigmas.len()
545 );
546 let old_sigma = old_sigmas
547 .into_iter()
548 .next()
549 .ok_or_else(|| format_err!("No σ found."))?;
550 self.add_improper(
551 ORDER_1,
552 old_sigma.raw_axis(),
553 false,
554 SIG,
555 Some("h".to_owned()),
556 presym.dist_threshold,
557 old_sigma.contains_time_reversal(),
558 );
559 self.add_improper(
560 ORDER_1,
561 old_sigma.raw_axis(),
562 true,
563 SIG,
564 Some("h".to_owned()),
565 presym.dist_threshold,
566 old_sigma.contains_time_reversal(),
567 );
568
569 self.set_group_name("Cs".to_owned());
570 } else {
571 let identity = (*self
572 .get_proper(&ORDER_1)
573 .ok_or_else(|| format_err!("No identity found."))?
574 .iter()
575 .next()
576 .ok_or_else(|| format_err!("No identity found."))?)
577 .clone();
578
579 self.add_proper(
580 ORDER_1,
581 identity.raw_axis(),
582 true,
583 presym.dist_threshold,
584 false,
585 );
586 self.set_group_name("C1".to_owned());
587 }
588 }
589 }
590
591 Ok(())
592 }
593}