molcrafts-molpack 0.1.0

Faithful Rust port of Packmol molecular packing (split off from molrs)
Documentation 
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

//! 3D linked cell list implementation.
//! Exact port of `cell_indexing.f90`, `resetcells.f90` and `pbc.f90` from Packmol.

use molrs::types::F;
/// Wrap coordinate into periodic box, returning value in [pbc_min, pbc_min + pbc_length).
#[inline]
fn in_box(x: F, pbc_min: F, pbc_length: F) -> F {
    let mut v = x - pbc_min;
    v -= (v / pbc_length).floor() * pbc_length;
    v + pbc_min
}

/// Compute cell index (0-based) for a given position along one axis.
///
/// When `periodic` is `true` the coordinate is wrapped into
/// `[pbc_min, pbc_min + pbc_length)` before indexing; when `false`,
/// positions outside the box are clamped to the nearest edge cell
/// (so atoms drifting past the confinement still land in a valid cell).
/// Port of Fortran `setcell`.
#[inline]
pub fn axis_cell(
    x: F,
    pbc_min: F,
    pbc_length: F,
    cell_length: F,
    ncells: usize,
    periodic: bool,
) -> usize {
    let xt = if periodic {
        in_box(x, pbc_min, pbc_length)
    } else {
        x
    };
    let raw = ((xt - pbc_min) / cell_length).floor();
    if raw < 0.0 {
        0
    } else {
        (raw as usize).min(ncells - 1)
    }
}

/// Compute the 3D cell indices (0-based) for a position, honoring
/// per-axis periodicity.
#[inline(always)]
pub fn setcell(
    pos: &[F; 3],
    pbc_min: &[F; 3],
    pbc_length: &[F; 3],
    cell_length: &[F; 3],
    ncells: &[usize; 3],
    periodic: &[bool; 3],
) -> [usize; 3] {
    [
        axis_cell(
            pos[0],
            pbc_min[0],
            pbc_length[0],
            cell_length[0],
            ncells[0],
            periodic[0],
        ),
        axis_cell(
            pos[1],
            pbc_min[1],
            pbc_length[1],
            cell_length[1],
            ncells[1],
            periodic[1],
        ),
        axis_cell(
            pos[2],
            pbc_min[2],
            pbc_length[2],
            cell_length[2],
            ncells[2],
            periodic[2],
        ),
    ]
}

/// Convert a 3D cell index (0-based) to a flat index.
/// Port of `index_cell` (Fortran 1-based → Rust 0-based).
#[inline]
pub fn index_cell(cell: &[usize; 3], ncells: &[usize; 3]) -> usize {
    cell[0] * ncells[1] * ncells[2] + cell[1] * ncells[2] + cell[2]
}

/// Convert a flat index to a 3D cell index (0-based).
/// Port of `icell_to_cell` (converted to 0-based).
pub fn icell_to_cell(icell: usize, ncells: &[usize; 3]) -> [usize; 3] {
    let k = icell % ncells[2];
    let rem = icell / ncells[2];
    let j = rem % ncells[1];
    let i = rem / ncells[1];
    [i, j, k]
}

/// Wrap a cell index with periodic boundary conditions (0-based).
/// Port of Fortran `cell_ind` (which uses 1-based, here we use 0-based).
#[inline]
pub fn cell_ind(idx: isize, ncells: usize) -> usize {
    idx.rem_euclid(ncells as isize) as usize
}

/// Compute PBC-corrected difference vector with per-axis wrapping.
/// Port of `delta_vector` from `pbc.f90`.
#[inline]
pub fn delta_vector(xi: &[F; 3], xj: &[F; 3], pbc_length: &[F; 3], periodic: &[bool; 3]) -> [F; 3] {
    let mut d = [xi[0] - xj[0], xi[1] - xj[1], xi[2] - xj[2]];
    for k in 0..3 {
        if periodic[k] {
            d[k] -= (d[k] / pbc_length[k]).round() * pbc_length[k];
        }
    }
    d
}