From 365e79fa10c168de091e19c7268f8229e72f6628 Mon Sep 17 00:00:00 2001 From: Andrea Grisafi Date: Fri, 6 Oct 2023 16:16:30 +0200 Subject: [PATCH] Exploited subtraction theorem in LODE structure factor (#242) --- .../calculators/lode/spherical_expansion.rs | 29 +++++++++++++++---- 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/rascaline/src/calculators/lode/spherical_expansion.rs b/rascaline/src/calculators/lode/spherical_expansion.rs index afaed528e..a10666e51 100644 --- a/rascaline/src/calculators/lode/spherical_expansion.rs +++ b/rascaline/src/calculators/lode/spherical_expansion.rs @@ -126,14 +126,31 @@ fn compute_structure_factors(positions: &[Vector3D], species: &[i32], k_vectors: .map(|s| (s, Array2::from_elem((n_atoms, n_k_vectors), 0.0))) .collect::>(); + // Precompute species-dependent sum over neighbours j of sines and cosines + let mut sumjcos = all_species.iter().copied() + .map(|s| (s, Array1::from_elem(n_k_vectors, 0.0))) + .collect::>(); + let mut sumjsin = all_species.iter().copied() + .map(|s| (s, Array1::from_elem(n_k_vectors, 0.0))) + .collect::>(); + for j in 0..n_atoms { + let sumjcos = sumjcos.get_mut(&species[j]).unwrap(); + let sumjsin = sumjsin.get_mut(&species[j]).unwrap(); + for k in 0..n_k_vectors { + sumjcos[k] += cosines[[j, k]]; + sumjsin[k] += sines[[j, k]]; + } + } + + // Compute Sum_j cos(k*r_ij) and Sum_j sin(k*r_ij) using the subtraction theorem for i in 0..n_atoms { - for j in 0..n_atoms { + for (species, real_per_center) in &mut real_per_center { + let sumjcos = sumjcos.get_mut(&species).unwrap(); + let sumjsin = sumjsin.get_mut(&species).unwrap(); + let imag_per_center = imag_per_center.get_mut(&species).unwrap(); for k in 0..n_k_vectors { - let real = cosines[[i, k]] * cosines[[j, k]] + sines[[i, k]] * sines[[j, k]]; - let imag = sines[[i, k]] * cosines[[j, k]] - cosines[[i, k]] * sines[[j, k]]; - - let real_per_center = real_per_center.get_mut(&species[j]).unwrap(); - let imag_per_center = imag_per_center.get_mut(&species[j]).unwrap(); + let real = cosines[[i, k]] * sumjcos[k] + sines[[i, k]] * sumjsin[k]; + let imag = sines[[i, k]] * sumjcos[k] - cosines[[i, k]] * sumjsin[k]; real_per_center[[i, k]] += 2.0 * real; imag_per_center[[i, k]] += 2.0 * imag; }