Atoms that contribute to the density not correctly found

[pfebrer]

Describe the bug

This part of the code: https://github.com/zerothi/sisl/blob/115479f6e54591de16196986fe01526380497ca7/src/sisl/physics/densitymatrix.py#L913-L926

which finds atoms with an orbital that reaches the unit cell, does not work properly.

This is critical because it gives wrong densities, since atoms that should contribute are not taken into account.

Code to reproduce problem

The simplest reproducing code that I could come up with is the following. I build a cell with two atoms inside and no periodic boundary conditions. Clearly, it should find both atoms to contribute to the unit cell.

With a cell [10, 10, 20], if the second atom is at z=10 it is correctly found, while if it is at z=15 it is not found.

I have no idea what causes this problem, but it is pretty nasty :(

(I didn't change anything in this part when I implemented the new algorithm for the density)

import sisl
import numpy as np
from sisl import Lattice

# This is what is being used to find atoms that contribute to the density
def get_contributing_atoms(geom):
    self = geom
    # Instead of looping all atoms in the supercell we find the exact atoms
    # and their supercell indices.
    add_R = np.full(3, self.maxR())
    # Calculate the required additional vectors required to increase the fictitious
    # supercell by add_R in each direction.
    # For extremely skewed lattices this will be way too much, hence we make
    # them square.
    o = self.lattice.to.Cuboid(True)
    lattice = Lattice(o._v + np.diag(2 * add_R), origin=o.origin - add_R)

    # Retrieve all atoms within the grid supercell
    # (and the neighbours that connect into the cell)
    IA, XYZ, ISC = self.within_inf(lattice, periodic=self.pbc)
    XYZ -= self.lattice.origin.reshape(1, 3)

    return IA, XYZ, ISC

# With the second atom at z=10, it determines that it contributes. 
geom = sisl.Geometry([[1, 1, 2], [1, 1, 10]], lattice=[10, 10, 20], atoms=sisl.Atom("C", R=3))
geom.lattice.set_boundary_condition("")

print("THIS IS FINE")
print(get_contributing_atoms(geom))

# With the second atom at z=15, it determines that it doesn't contribute
geom = sisl.Geometry([[1, 1, 2], [1, 1, 15]], lattice=[10, 10, 20], atoms=sisl.Atom("C", R=3))
geom.lattice.set_boundary_condition("")

print("THIS IS NOT")
print(get_contributing_atoms(geom))

THIS IS FINE
(array([0, 1], dtype=int32), array([[ 1.,  1.,  2.],
       [ 1.,  1., 10.]]), array([[0, 0, 0],
       [0, 0, 0]], dtype=int32))
THIS IS NOT
(array([0], dtype=int32), array([[1., 1., 2.]]), array([[0, 0, 0]], dtype=int32))

[zerothi]

This might be related to #511.

[pfebrer]

No, I don't think within_inf is the problem here.

The computed lattice to pass to within_inf seems wrong. If I print it I get:

Lattice{nsc: [1 1 1],
 origin=[-7.0000, -7.0000, -12.0000],
 A=[16.0000, 0.0000, 0.0000],
 B=[0.0000, 16.0000, 0.0000],
 C=[0.0000, 0.0000, 26.0000],
 bc=[Periodic,
     Periodic,
     Periodic]
}

so it only includes up to z=14.

I don't know what toCuboid is suposed to do :sweat_smile:

[pfebrer]

Aaaah ok, I see what is happening. I think previously orthogonal was the first argument of toCuboid and now it is not.

So toCuboid(True) must be changed to toCuboid(orthogonal=True)

[zerothi]

Thanks for backtracking, and fixing! Have been hung up all day!