diff --git a/python/rascaline/rascaline/systems/ase.py b/python/rascaline/rascaline/systems/ase.py index ef8ca1c70..4a4bb308d 100644 --- a/python/rascaline/rascaline/systems/ase.py +++ b/python/rascaline/rascaline/systems/ase.py @@ -91,10 +91,23 @@ def compute_neighbors(self, cutoff): nl_result = neighborlist.neighbor_list("ijdDS", self._atoms, cutoff) for i, j, d, D, S in zip(*nl_result): + # we want a half neighbor list, so drop all duplicated neighbors if j < i: - # we want a half neighbor list, so drop all duplicated - # neighbors continue + elif i == j: + if S[0] == 0 and S[1] == 0 and S[2] == 0: + # only create pairs with the same atom twice if the pair spans more + # than one unit cell + continue + elif S[0] + S[1] + S[2] < 0 or ( + (S[0] + S[1] + S[2] == 0) and (S[2] < 0 or (S[2] == 0 and S[1] < 0)) + ): + # When creating pairs between an atom and one of its periodic + # images, the code generate multiple redundant pairs (e.g. with + # shifts 0 1 1 and 0 -1 -1); and we want to only keep one of these. + # We keep the pair in the positive half plane of shifts. + continue + self._pairs.append((i, j, d, D, S)) self._pairs_by_center = [] diff --git a/python/rascaline/tests/systems/ase.py b/python/rascaline/tests/systems/ase.py index 0ca19795d..892060dd7 100644 --- a/python/rascaline/tests/systems/ase.py +++ b/python/rascaline/tests/systems/ase.py @@ -1,6 +1,7 @@ import numpy as np import pytest +from rascaline import SphericalExpansion from rascaline.systems import AseSystem @@ -145,3 +146,56 @@ def test_no_pbc_cell(): ) with pytest.warns(Warning, match=message): AseSystem(atoms) + + +def test_same_spherical_expansion(): + system = ase.Atoms( + "CaC6", + positions=[ + (0.0, 0.0, 0.0), + (1.88597, 1.92706, 0.0113749), + (2.66157, 3.55479, 7.7372), + (2.35488, 3.36661, 3.88), + (2.19266, 2.11524, 3.86858), + (2.52936, 4.62777, 3.87771), + (2.01817, 0.85408, 3.87087), + ], + cell=[[3.57941, 0, 0], [0.682558, 5.02733, 0], [0.285565, 0.454525, 7.74858]], + pbc=True, + ) + + calculator = SphericalExpansion( + # make sure to choose a cutoff larger then the cell to test for pairs crossing + # multiple periodic boundaries + cutoff=9, + max_radial=5, + max_angular=5, + atomic_gaussian_width=0.3, + radial_basis={"Gto": {}}, + center_atom_weight=1.0, + cutoff_function={"Step": {}}, + ) + + rascaline_nl = calculator.compute( + system, gradients=["positions", "cell"], use_native_system=True + ) + + ase_nl = calculator.compute( + system, gradients=["positions", "cell"], use_native_system=False + ) + + for key, block in rascaline_nl.items(): + ase_block = ase_nl.block(key) + + assert ase_block.samples == block.samples + # Since the pairs are in a different order, the values are slightly different + assert np.allclose(ase_block.values, block.values, atol=1e-16, rtol=1e-9) + + for parameter in ["cell", "positions"]: + gradient = block.gradient(parameter) + ase_gradient = ase_block.gradient(parameter) + + assert gradient.samples == ase_gradient.samples + assert np.allclose( + ase_gradient.values, gradient.values, atol=1e-16, rtol=1e-6 + ) diff --git a/rascaline/src/calculators/neighbor_list.rs b/rascaline/src/calculators/neighbor_list.rs index 0834713c0..8537a679d 100644 --- a/rascaline/src/calculators/neighbor_list.rs +++ b/rascaline/src/calculators/neighbor_list.rs @@ -39,7 +39,7 @@ pub struct NeighborList { pub full_neighbor_list: bool, /// Should individual atoms be considered their own neighbor? Setting this /// to `true` will add "self pairs", i.e. pairs between an atom and itself, - /// with the distance 0. The `pair_id` of such pairs is set to -1. + /// with the distance 0. pub self_pairs: bool, } @@ -423,7 +423,8 @@ impl FullNeighborList { let cell_c = pair.cell_shift_indices[2]; if species_first == species_second { - // same species for both atoms in the pair + // same species for both atoms in the pair, add the pair + // twice in both directions. if species[pair.first] == species_first.i32() && species[pair.second] == species_second.i32() { builder.add(&[ LabelValue::from(system_i), @@ -434,18 +435,14 @@ impl FullNeighborList { LabelValue::from(cell_c), ]); - if pair.first != pair.second { - // if the pair is between two different atoms, - // also add the reversed (second -> first) pair. - builder.add(&[ - LabelValue::from(system_i), - LabelValue::from(pair.second), - LabelValue::from(pair.first), - LabelValue::from(-cell_a), - LabelValue::from(-cell_b), - LabelValue::from(-cell_c), - ]); - } + builder.add(&[ + LabelValue::from(system_i), + LabelValue::from(pair.second), + LabelValue::from(pair.first), + LabelValue::from(-cell_a), + LabelValue::from(-cell_b), + LabelValue::from(-cell_c), + ]); } } else { // different species, find the right order for the pair @@ -501,6 +498,11 @@ impl FullNeighborList { let species = system.species()?; for pair in system.pairs()? { + if pair.first == pair.second { + // self pairs should not be part of the neighbors list + assert_ne!(pair.cell_shift_indices, [0, 0, 0]); + } + let first_block_i = descriptor.keys().position(&[ species[pair.first].into(), species[pair.second].into() ]); @@ -565,11 +567,6 @@ impl FullNeighborList { } } - if pair.first == pair.second { - // do not duplicate self pairs - continue; - } - // then the pair second -> first if let Some(second_block_i) = second_block_i { let mut block = descriptor.block_mut_by_id(second_block_i); @@ -764,6 +761,75 @@ mod tests { assert_relative_eq!(array, expected, max_relative=1e-6); } + #[test] + fn periodic_neighbor_list() { + let mut calculator = Calculator::from(Box::new(NeighborList{ + cutoff: 12.0, + full_neighbor_list: false, + self_pairs: false, + }) as Box); + + let mut systems = test_systems(&["CH"]); + + let descriptor = calculator.compute(&mut systems, Default::default()).unwrap(); + assert_eq!(*descriptor.keys(), Labels::new( + ["species_first_atom", "species_second_atom"], + &[[1, 1], [1, 6], [6, 6]] + )); + + // H-H block + let block = descriptor.block_by_id(0); + assert_eq!(block.samples(), Labels::new( + ["structure", "first_atom", "second_atom", "cell_shift_a", "cell_shift_b", "cell_shift_c"], + // the pairs only differ in cell shifts + &[[0, 1, 1, 0, 0, 1], [0, 1, 1, 0, 1, 0], [0, 1, 1, 1, 0, 0]] + )); + + let array = block.values().to_array(); + let expected = &ndarray::arr3(&[ + [[0.0], [0.0], [10.0]], + [[0.0], [10.0], [0.0]], + [[10.0], [0.0], [0.0]], + ]).into_dyn(); + assert_relative_eq!(array, expected, max_relative=1e-6); + + // now a full NL + let mut calculator = Calculator::from(Box::new(NeighborList{ + cutoff: 12.0, + full_neighbor_list: true, + self_pairs: false, + }) as Box); + + let descriptor = calculator.compute(&mut systems, Default::default()).unwrap(); + assert_eq!(*descriptor.keys(), Labels::new( + ["species_first_atom", "species_second_atom"], + &[[1, 1], [1, 6], [6, 1], [6, 6]] + )); + + // H-H block + let block = descriptor.block_by_id(0); + assert_eq!(block.samples(), Labels::new( + ["structure", "first_atom", "second_atom", "cell_shift_a", "cell_shift_b", "cell_shift_c"], + // twice as many pairs + &[ + [0, 1, 1, 0, 0, 1], [0, 1, 1, 0, 0, -1], + [0, 1, 1, 0, 1, 0], [0, 1, 1, 0, -1, 0], + [0, 1, 1, 1, 0, 0], [0, 1, 1, -1, 0, 0], + ] + )); + + let array = block.values().to_array(); + let expected = &ndarray::arr3(&[ + [[0.0], [0.0], [10.0]], + [[0.0], [0.0], [-10.0]], + [[0.0], [10.0], [0.0]], + [[0.0], [-10.0], [0.0]], + [[10.0], [0.0], [0.0]], + [[-10.0], [0.0], [0.0]], + ]).into_dyn(); + assert_relative_eq!(array, expected, max_relative=1e-6); + } + #[test] fn finite_differences_positions() { // half neighbor list diff --git a/rascaline/src/calculators/soap/spherical_expansion.rs b/rascaline/src/calculators/soap/spherical_expansion.rs index 63a8ddf19..aac10d3ca 100644 --- a/rascaline/src/calculators/soap/spherical_expansion.rs +++ b/rascaline/src/calculators/soap/spherical_expansion.rs @@ -241,12 +241,6 @@ impl SphericalExpansion { } } - if pair.first == pair.second { - // do not compute for the reversed pair if the pair is - // between an atom and its image - continue; - } - if let Some(mapped_center) = result.centers_mapping[pair.second] { // add the pair contribution to the atomic environnement // corresponding to the **second** atom in the pair @@ -778,7 +772,7 @@ mod tests { fn parameters() -> SphericalExpansionParameters { SphericalExpansionParameters { - cutoff: 3.5, + cutoff: 7.8, max_radial: 6, max_angular: 6, atomic_gaussian_width: 0.3, diff --git a/rascaline/src/calculators/soap/spherical_expansion_pair.rs b/rascaline/src/calculators/soap/spherical_expansion_pair.rs index 737b68797..898a31523 100644 --- a/rascaline/src/calculators/soap/spherical_expansion_pair.rs +++ b/rascaline/src/calculators/soap/spherical_expansion_pair.rs @@ -755,13 +755,7 @@ impl CalculatorBase for SphericalExpansionByPair { } } - // also check for the block with a reversed pair, except if - // we are handling a pair between an atom and it's own - // periodic image - if pair.first == pair.second { - continue; - } - + // also check for the block with a reversed pair contribution.inverse_pair(&self.m_1_pow_l); for spherical_harmonics_l in 0..=self.parameters.max_angular { @@ -817,7 +811,7 @@ mod tests { fn parameters() -> SphericalExpansionParameters { SphericalExpansionParameters { - cutoff: 3.5, + cutoff: 7.3, max_radial: 6, max_angular: 6, atomic_gaussian_width: 0.3, diff --git a/rascaline/src/systems/neighbors.rs b/rascaline/src/systems/neighbors.rs index fca5f24b8..e72141787 100644 --- a/rascaline/src/systems/neighbors.rs +++ b/rascaline/src/systems/neighbors.rs @@ -132,9 +132,9 @@ impl CellList { // number of cells to search in each direction to make sure all possible // pairs below the cutoff are accounted for. let mut n_search = [ - f64::trunc(cutoff * n_cells[0] / distances_between_faces[0]) as i32, - f64::trunc(cutoff * n_cells[1] / distances_between_faces[1]) as i32, - f64::trunc(cutoff * n_cells[2] / distances_between_faces[2]) as i32, + f64::ceil(cutoff * n_cells[0] / distances_between_faces[0]) as i32, + f64::ceil(cutoff * n_cells[1] / distances_between_faces[1]) as i32, + f64::ceil(cutoff * n_cells[2] / distances_between_faces[2]) as i32, ]; let n_cells = [ @@ -174,7 +174,7 @@ impl CellList { let n_cells = self.cells.shape(); let n_cells = [n_cells[0], n_cells[1], n_cells[2]]; - // find the cell in which this atom should go + // find the subcell in which this atom 'should go' let cell_index = [ f64::floor(fractional[0] * n_cells[0] as f64) as i32, f64::floor(fractional[1] * n_cells[1] as f64) as i32, @@ -250,10 +250,39 @@ impl CellList { let shift = CellShift(cell_shift) + atom_i.shift - atom_j.shift; let shift_is_zero = shift[0] == 0 && shift[1] == 0 && shift[2] == 0; - if atom_i.index == atom_j.index && shift_is_zero { - // only create pair with the same atom twice - // if the pair spans more than one unit cell - continue; + if atom_i.index == atom_j.index { + if shift_is_zero { + // only create pairs with the same atom twice + // if the pair spans more than one unit cell + continue; + } + + // When creating pairs between an atom and one of its + // periodic images, the code generate multiple redundant + // pairs (e.g. with shifts 0 1 1 and 0 -1 -1); and we + // want to only keep one of these. + if shift[0] + shift[1] + shift[2] < 0 { + // drop shifts on the negative half-space + continue; + } + + if (shift[0] + shift[1] + shift[2] == 0) + && (shift[2] < 0 || (shift[2] == 0 && shift[1] < 0)) { + + // drop shifts in the negative half plane or the + // negative shift[1] axis. See below for a graphical + // representation: we are keeping the shifts indicated + // with `O` and dropping the ones indicated with `X` + // + // O O O │ O O O + // O O O │ O O O + // O O O │ O O O + // ─X─X─X─┼─O─O─O─ + // X X X │ X X X + // X X X │ X X X + // X X X │ X X X + continue; + } } if self.unit_cell.is_infinite() && !shift_is_zero { @@ -378,8 +407,7 @@ impl NeighborsList { mod tests { use approx::assert_ulps_eq; - use crate::Matrix3; - + use crate::types::Matrix3; use super::*; #[test] @@ -426,21 +454,15 @@ mod tests { let neighbors = NeighborsList::new(&positions, cell, 3.0); let expected = [ - (Vector3D::new(0.0, -1.0, -1.0), [-1, 0, 0]), (Vector3D::new(1.0, 0.0, -1.0), [-1, 0, 1]), (Vector3D::new(1.0, -1.0, 0.0), [-1, 1, 0]), - (Vector3D::new(-1.0, 0.0, -1.0), [0, -1, 0]), (Vector3D::new(0.0, 1.0, -1.0), [0, -1, 1]), - (Vector3D::new(-1.0, -1.0, 0.0), [0, 0, -1]), (Vector3D::new(1.0, 1.0, 0.0), [0, 0, 1]), - (Vector3D::new(0.0, -1.0, 1.0), [0, 1, -1]), (Vector3D::new(1.0, 0.0, 1.0), [0, 1, 0]), - (Vector3D::new(-1.0, 1.0, 0.0), [1, -1, 0]), - (Vector3D::new(-1.0, 0.0, 1.0), [1, 0, -1]), (Vector3D::new(0.0, 1.0, 1.0), [1, 0, 0]), ]; - assert_eq!(neighbors.pairs.len(), 12); + assert_eq!(neighbors.pairs.len(), 6); for (pair, (vector, shifts)) in neighbors.pairs.iter().zip(&expected) { assert_eq!(pair.first, 0); assert_eq!(pair.second, 0); @@ -480,4 +502,61 @@ mod tests { assert_ulps_eq!(pair.distance, 2.0); } } + + #[test] + fn small_cell_large_cutoffs() { + let cell = UnitCell::cubic(0.5); + let positions = [Vector3D::new(0.0, 0.0, 0.0)]; + let neighbors = NeighborsList::new(&positions, cell, 0.6); + + let expected = [ + (Vector3D::new(0.0, 0.0, 0.5), [0, 0, 1]), + (Vector3D::new(0.0, 0.5, 0.0), [0, 1, 0]), + (Vector3D::new(0.5, 0.0, 0.0), [1, 0, 0]), + ]; + + assert_eq!(neighbors.pairs.len(), 3); + for (pair, (vector, shifts)) in neighbors.pairs.iter().zip(&expected) { + assert_eq!(pair.first, 0); + assert_eq!(pair.second, 0); + assert_ulps_eq!(pair.distance, 0.5); + assert_ulps_eq!(pair.vector, vector); + assert_eq!(&pair.cell_shift_indices, shifts); + } + } + + #[test] + fn non_cubic_cell() { + let cell = UnitCell::from(Matrix3::new([ + [4.26, -2.45951215, 0.0], + [2.13, 1.22975607, 0.0], + [0.0, 0.0, 50.0], + ])); + let positions = [ + Vector3D::new(1.42, 0.0, 0.0), + Vector3D::new(2.84, 0.0, 0.0), + Vector3D::new(3.55, -1.22975607, 0.0), + Vector3D::new(4.97, -1.22975607, 0.0), + ]; + let neighbors = NeighborsList::new(&positions, cell, 6.4); + + assert_eq!(neighbors.pairs.len(), 90); + + let previously_missing = [ + (0, 3, [-2, 0, 0]), + (0, 3, [-2, 1, 0]), + (0, 3, [-2, 2, 0]), + ]; + + for missing in previously_missing { + let mut found = false; + for pair in &neighbors.pairs { + if pair.first == missing.0 && pair.second == missing.1 + && pair.cell_shift_indices == missing.2 { + found = true; + } + } + assert!(found, "could not find pair {:?}", missing); + } + } }