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

use super::bragg;
use super::parameters;


// JCP 88, 2547 (1988), eq. 20
#[inline]
fn f3(x: f64, hardness: usize) -> f64 {
    let mut f = x;

    for _ in 0..hardness {
        f *= 1.5 - 0.5 * f * f;
    }

    f
}

fn distance(p1: &(f64, f64, f64), p2: &(f64, f64, f64)) -> f64 {
    let dx = p1.0 - p2.0;
    let dy = p1.1 - p2.1;
    let dz = p1.2 - p2.2;

    (dx * dx + dy * dy + dz * dz).sqrt()
}

// JCP 88, 2547 (1988)
pub fn partitioning_weight(
    center_index: usize,
    center_coordinates_bohr: &[(f64, f64, f64)],
    proton_charges: &[i32],
    grid_coordinates_bohr: (f64, f64, f64),
    hardness: usize,
) -> f64 {
    let num_centers = proton_charges.len();

    let mut pa = vec![1.0; num_centers];

    for ia in 0..num_centers {
        let dist_a = distance(&grid_coordinates_bohr, &center_coordinates_bohr[ia]);

        let r_a = bragg::get_bragg_angstrom(proton_charges[ia]);

        for ib in 0..ia {
            let dist_b = distance(&grid_coordinates_bohr, &center_coordinates_bohr[ib]);

            let r_b = bragg::get_bragg_angstrom(proton_charges[ib]);

            let dist_ab = distance(&center_coordinates_bohr[ia], &center_coordinates_bohr[ib]);

            // JCP 88, 2547 (1988), eq. 11
            let mu_ab = (dist_a - dist_b) / dist_ab;

            let mut nu_ab = mu_ab;
            if (r_a - r_b).abs() > parameters::SMALL {
                let u_ab = (r_a + r_b) / (r_b - r_a);
                let a_ab = u_ab / (u_ab * u_ab - 1.0);

                nu_ab += a_ab.min(0.5).max(-0.5) * (1.0 - mu_ab * mu_ab);
            }

            let f = f3(nu_ab, hardness);

            if (1.0 - f).abs() > parameters::SMALL {
                pa[ia] *= 0.5 * (1.0 - f);
                pa[ib] *= 0.5 * (1.0 + f);
            } else {
                // avoid numerical issues
                pa[ia] = 0.0;
            }
        }
    }

    let w: f64 = pa.iter().sum();

    if w.abs() > parameters::SMALL {
        pa[center_index] / w
    } else {
        1.0
    }
}