subtraction theorem exploited in LODE structure factor

[andreagrisafi]

Subtraction theorem exploited when computing the LODE structure factor in order to remove the double sum over atoms.

---

📚 Documentation preview 📚: https://rascaline--242.org.readthedocs.build/en/242/

[Luthaf]

For context, when running on the first frame of molecular_crystals (176 atoms), the old code produced this profile:

╔════╦══════════════════════════════════════════════╦════════════╦═══════════╦══════════╦═════════════╗
║ id ║ span name                                    ║ call count ║ called by ║ total    ║ mean        ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  2 ║ Full calculation                             ║          1 ║         — ║    6.19s ║       6.19s ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  4 ║ LodeSphericalExpansion::compute              ║          1 ║         2 ║    6.17s ║       6.17s ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  1 ║ LodeRadialIntegralSpline::with_accuracy      ║          2 ║      4, 2 ║  24.31ms ║     12.16ms ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  0 ║ LodeRadialIntegralGto::compute               ║       5122 ║         1 ║  21.18ms ║      4.14µs ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  3 ║ Calculator::prepare                          ║          1 ║         2 ║ 926.31µs ║    926.31µs ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  5 ║ compute_k_vectors                            ║          1 ║         4 ║ 161.05µs ║    161.05µs ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  7 ║ project_k_to_nlm                             ║          1 ║         4 ║  14.74ms ║     14.74ms ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  6 ║ SplinedRadialIntegral::compute               ║      25662 ║         7 ║   2.90ms ║ 113.00ns ⚠️ ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  8 ║ SphericalHarmonics::compute                  ║      25662 ║         7 ║ 875.10µs ║  34.00ns ⚠️ ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  9 ║ compute_structure_factors                    ║          1 ║         4 ║    3.60s ║       3.60s ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║ 10 ║ compute_density_fourrier                     ║          1 ║         4 ║ 164.45µs ║    164.45µs ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║ 11 ║ LODE main loop                               ║          1 ║         4 ║    2.53s ║       2.53s ║

and the new code produce this profile:

╔════╦══════════════════════════════════════════════╦════════════╦═══════════╦══════════╦═════════════╗
║ id ║ span name                                    ║ call count ║ called by ║ total    ║ mean        ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  2 ║ Full calculation                             ║          1 ║         — ║    2.75s ║       2.75s ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  4 ║ LodeSphericalExpansion::compute              ║          1 ║         2 ║    2.73s ║       2.73s ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  1 ║ LodeRadialIntegralSpline::with_accuracy      ║          2 ║      4, 2 ║  25.36ms ║     12.68ms ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  0 ║ LodeRadialIntegralGto::compute               ║       5122 ║         1 ║  22.30ms ║      4.35µs ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  3 ║ Calculator::prepare                          ║          1 ║         2 ║   1.04ms ║      1.04ms ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  5 ║ compute_k_vectors                            ║          1 ║         4 ║ 216.40µs ║    216.40µs ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  7 ║ project_k_to_nlm                             ║          1 ║         4 ║  15.19ms ║     15.19ms ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  6 ║ SplinedRadialIntegral::compute               ║      25662 ║         7 ║   2.88ms ║ 112.00ns ⚠️ ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  8 ║ SphericalHarmonics::compute                  ║      25662 ║         7 ║ 957.23µs ║  37.00ns ⚠️ ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║  9 ║ compute_structure_factors                    ║          1 ║         4 ║ 119.99ms ║    119.99ms ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║ 10 ║ compute_density_fourrier                     ║          1 ║         4 ║ 129.11µs ║    129.11µs ║
╠════╬══════════════════════════════════════════════╬════════════╬═══════════╬══════════╬═════════════╣
║ 11 ║ LODE main loop                               ║          1 ║         4 ║    2.58s ║       2.58s ║
╚════╩══════════════════════════════════════════════╩════════════╩═══════════╩══════════╩═════════════╝

compute_structure_factors is down to 120ms from 3.6s, and the overall calculation is a bit more than twice as fast 🎉