(
alpha_min: HashMap<usize, f64>,
alpha_max: f64,
radial_precision: f64,
min_num_angular_points: usize,
max_num_angular_points: usize,
proton_charges: Vec<i32>,
center_inde
| 40 | |
| 41 | #[pyfunction] |
| 42 | pub fn atom_grid( |
| 43 | alpha_min: HashMap<usize, f64>, |
| 44 | alpha_max: f64, |
| 45 | radial_precision: f64, |
| 46 | min_num_angular_points: usize, |
| 47 | max_num_angular_points: usize, |
| 48 | proton_charges: Vec<i32>, |
| 49 | center_index: usize, |
| 50 | center_coordinates_bohr: Vec<(f64, f64, f64)>, |
| 51 | hardness: usize, |
| 52 | ) -> (Vec<(f64, f64, f64)>, Vec<f64>) { |
| 53 | let (rs, weights_radial) = radial::radial_grid_lmg( |
| 54 | alpha_min, |
| 55 | alpha_max, |
| 56 | radial_precision, |
| 57 | proton_charges[center_index], |
| 58 | ); |
| 59 | |
| 60 | // factors match DIRAC code |
| 61 | let rb = bragg::get_bragg_angstrom(proton_charges[center_index]) / (5.0 * 0.529177249); |
| 62 | |
| 63 | let mut coordinates = Vec::new(); |
| 64 | let mut weights = Vec::new(); |
| 65 | |
| 66 | let pi = std::f64::consts::PI; |
| 67 | |
| 68 | let cx = center_coordinates_bohr[center_index].0; |
| 69 | let cy = center_coordinates_bohr[center_index].1; |
| 70 | let cz = center_coordinates_bohr[center_index].2; |
| 71 | |
| 72 | for (&r, &weight_radial) in rs.iter().zip(weights_radial.iter()) { |
| 73 | // we read the angular grid at each radial step because of pruning |
| 74 | // this can be optimized |
| 75 | let mut num_angular = max_num_angular_points; |
| 76 | if r < rb { |
| 77 | num_angular = ((max_num_angular_points as f64) * r / rb) as usize; |
| 78 | num_angular = lebedev::get_closest_num_angular(num_angular); |
| 79 | if num_angular < min_num_angular_points { |
| 80 | num_angular = min_num_angular_points; |
| 81 | } |
| 82 | } |
| 83 | let (coordinates_angular, weights_angular) = lebedev::angular_grid(num_angular); |
| 84 | |
| 85 | let wt = 4.0 * pi * weight_radial; |
| 86 | for (&xyz, &weight_angular) in coordinates_angular.iter().zip(weights_angular.iter()) { |
| 87 | let x = cx + r * xyz.0; |
| 88 | let y = cy + r * xyz.1; |
| 89 | let z = cz + r * xyz.2; |
| 90 | |
| 91 | coordinates.push((x, y, z)); |
| 92 | weights.push(wt * weight_angular); |
| 93 | } |
| 94 | } |
| 95 | |
| 96 | if center_coordinates_bohr.len() > 1 { |
| 97 | let w_partitioning: Vec<f64> = coordinates |
| 98 | .par_iter() |
| 99 | .map(|c| { |
no test coverage detected