spdkit_surface/
probe.rs

1// [[file:../spdkit-surface.note::db9e4081][db9e4081]]
2use crate::common::*;
3
4use gchemol::neighbors::{Neighbor, Neighborhood};
5use gchemol::Lattice;
6use vecfx::*;
7// db9e4081 ends here
8
9// [[file:../spdkit-surface.note::1fc3542c][1fc3542c]]
10fn guess_appropriate_zstep(lat: &Lattice, radius: f64) -> f64 {
11    assert!(radius.is_sign_positive(), "invalid radius: {}", radius);
12    let [_, _, zwidth] = lat.widths();
13    radius / zwidth
14}
15
16fn probe_surface_neighgbors_along_z(
17    probe: &Neighborhood,
18    lat: &gchemol::Lattice,
19    x: f64,
20    y: f64,
21    zhigh: f64,
22    zstep: f64,
23    r_cutoff: f64,
24) -> Result<Vec<Neighbor>> {
25    assert!(zhigh >= 0.0, "invalid zhigh: {}", zhigh);
26
27    // test if current probe is above or below the surface
28    let p = lat.to_cart([x, y, zhigh]);
29    let n = probe.search(p.into(), r_cutoff).collect_vec();
30    let mut i = 0;
31    if n.is_empty() {
32        // above the surface: move down step by step until touch the surface atoms
33        let mut z = zhigh;
34        loop {
35            i += 1;
36            z -= zstep;
37            let p = lat.to_cart([x, y, z]);
38            let n = probe.search(p.into(), r_cutoff).collect_vec();
39            if n.len() > 0 || z <= 0.0 {
40                trace!("probe surface atoms by moving down done in {} iterations", i);
41                break Ok(n);
42            }
43        }
44    } else {
45        // below the surface: move up step by step until move away the surface atoms
46        let mut z = zhigh;
47        let mut n_pre = n;
48        loop {
49            i += 1;
50            z += zstep;
51            let p = lat.to_cart([x, y, z]);
52            let n = probe.search(p.into(), r_cutoff).collect_vec();
53            if n.len() == 0 || z >= 1.0 {
54                trace!("probe surface atoms by moving up done in {} iterations", i);
55                break Ok(n_pre);
56            }
57            n_pre = n;
58        }
59    }
60}
61// 1fc3542c ends here
62
63// [[file:../spdkit-surface.note::c4191fe3][c4191fe3]]
64/// Probe surface atoms approaching slab surface from top in z-axis
65/// gradually
66///
67/// # Parameters
68/// * mol: the molecule to be probed (slab model)
69/// * r_cutoff: the cutoff radius for probing surface atoms
70/// * n_probes: the number of random probe atoms
71pub fn probe_surface_atoms(mol: &Molecule, r_cutoff: f64, n_probes: usize) -> Result<Vec<usize>> {
72    let n_probes = 800;
73    ensure!(mol.lattice.is_some(), "only work for periodic system");
74
75    let zhigh = mol.get_scaled_positions().unwrap().map(|[_, _, fz]| fz).float_max();
76    let probe = mol.create_neighbor_probe();
77    let lat = mol.get_lattice().ok_or(anyhow!("no lattice"))?;
78
79    let zstep = guess_appropriate_zstep(mol.lattice.as_ref().unwrap(), 0.5);
80    debug!("probe step size along z = {}", zstep);
81    let mut surface_nodes = std::collections::HashSet::new();
82    for _ in 0..n_probes {
83        let [x, y] = crate::sample::random_frac_xy();
84        let nn = probe_surface_neighgbors_along_z(&probe, lat, x, y, zhigh, zstep, r_cutoff)?;
85        for n in nn {
86            surface_nodes.insert(n.node);
87        }
88    }
89    let nodes = surface_nodes.into_iter().collect();
90    Ok(nodes)
91}
92// c4191fe3 ends here
93
94// [[file:../spdkit-surface.note::c797229e][c797229e]]
95#[test]
96fn test_probe_surface_atoms() -> Result<()> {
97    gut::cli::setup_logger_for_test();
98
99    let file = "./tests/files/Ni211.cif";
100    let mol = Molecule::from_file(file)?;
101    let surface_atoms = probe_surface_atoms(&mol, 1.8, 500)?;
102    let s = gut::utils::abbreviate_numbers_human_readable(&surface_atoms)?;
103    assert_eq!(
104        "2,4,8,10,12,16,18,20,24,26,28,32,34,36,40,42,44,48,50,52,56,58,60,64,66,68,72,74,76,80",
105        s
106    );
107
108    Ok(())
109}
110// c797229e ends here
spdkit_surface/probe.rs

spdkit_surface/
probe.rs
