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

// import

// [[file:~/Workspace/Programming/gchemol-rs/neighbors/neighbors.note::*import][import:1]]

// import:1 ends here

// mods

// [[file:~/Workspace/Programming/gchemol-rs/neighbors/neighbors.note::*mods][mods:1]]
mod aperiodic;
mod periodic;
// mods:1 ends here

// base

// [[file:~/Workspace/Programming/gchemol-rs/neighbors/neighbors.note::*base][base:1]]
mod base {
    use lattice::Lattice;
    use indexmap::IndexMap;
    use octree::Octree;
    use vecfx::Vector3f;

    pub type Point = [f64; 3];

    /// Helper struct for neighbors search result.
    #[derive(Debug, Clone, Copy)]
    pub struct Neighbor {
        /// The node connected to the host point.
        pub node: usize,
        /// The distance to the host point.
        pub distance: f64,
        /// Scaled displacment vector relative to origin cell if PBC enabled.
        pub image: Option<Vector3f>,
    }

    /// Neighborhood is a neighboring nodes detector, for given cutoff distance.
    #[derive(Debug, Clone, Default)]
    pub struct Neighborhood {
        /// particle coordinates
        pub(crate) points: IndexMap<usize, Point>,

        /// Octree object
        pub(crate) tree: Option<Octree>,

        /// Periodic lattice.
        pub(crate) lattice: Option<Lattice>,
    }
}
// base:1 ends here

// api

// [[file:~/Workspace/Programming/gchemol-rs/neighbors/neighbors.note::*api][api:1]]
mod api {
    use crate::base::*;
    use lattice::Lattice;
    use octree::Octree;

    impl Neighborhood {
        /// Constructs a neighborhood detector using the given `cutoff` distance.
        pub fn new() -> Self {
            Self { ..Default::default() }
        }

        /// Automatically build and update Neighborhood with contents of an
        /// iterator.
        ///
        /// The position of a point is associated with a permanent key in type
        /// of `usize`.
        pub fn update<I>(&mut self, iter: I)
        where
            I: IntoIterator<Item = (usize, Point)>,
        {
            // update data points
            for (k, v) in iter {
                self.points.insert(k, v);
            }
            // FIXME: how to reduce octree building
            let points: Vec<_> = self.points.values().copied().collect();
            let mut tree = Octree::new(points);
            let bucket_size = 100;
            tree.build(bucket_size);
            self.tree = Some(tree);
        }

        #[cfg(feature = "adhoc")]
        /// Reset internal data.
        pub fn clear(&mut self) {
            self.points.clear();
            self.lattice = None;
            self.tree = None;
        }

        /// Return a list of the nodes connected to the node `n`.
        ///
        /// Parameters
        /// ----------
        /// * n: the key of host node for searching neighbors
        /// * radius: cutoff radius distance
        pub fn neighbors(&self, n: usize, radius: f64) -> impl Iterator<Item = Neighbor> + '_ {
            // the index of host node `n` in point list.
            let (_, _, &pt) = self.points.get_full(&n).expect("invalid key");

            // FIXME: think twice
            // excluding self from the list
            let epsilon = 1e-6;
            self.search(pt, radius).filter_map(move |m| {
                if m.node == n && m.distance < epsilon {
                    return None;
                }
                Some(m)
            })
        }

        /// Return neighbors of a particle `pt` within distance cutoff `radius`.
        pub fn search(&self, pt: Point, radius: f64) -> impl Iterator<Item = Neighbor> + '_ {
            // inspired by: https://stackoverflow.com/a/54728634
            let mut iter_periodic = None;
            let mut iter_aperiodic = None;
            match self.lattice {
                Some(lattice) => {
                    let iter = self.search_neighbors_periodic(pt, radius, lattice);
                    iter_periodic = Some(iter);
                }
                None => {
                    let iter = self.search_neighbors_aperiodic(pt, radius);
                    iter_aperiodic = Some(iter);
                }
            }
            iter_periodic
                .into_iter()
                .flatten()
                .chain(iter_aperiodic.into_iter().flatten())
        }

        /// Return current number of points.
        pub fn npoints(&self) -> usize {
            self.points.len()
        }

        /// Set lattice for applying periodic boundary conditions
        pub fn set_lattice(&mut self, mat: [[f64; 3]; 3]) {
            let lat = Lattice::new(mat);
            self.lattice = Some(lat);
        }
    }
}
// api:1 ends here

// pub

// [[file:~/Workspace/Programming/gchemol-rs/neighbors/neighbors.note::*pub][pub:1]]
pub use crate::api::*;
pub use crate::base::*;
// pub:1 ends here