fixed #633 due to weird combination of boundary lookups

[zerothi]

The problem of close arose when the searched sphere was just crossing a border of the unit-cell (or supercell). In those cases could the translated to unit-cell indices miss a neighbour.

The fix was rather simple; once the atoms closests has been found, only loop over those in the primary unit-cell.

There are many other small changes here, but they
are more cosmetic and should provide more stability when using orbital distances very close to atomic coordinates.

• Closes Bug in construct? #633
• Tests added
• Ranned isort . at top--level
• Documentation for functionality, docs/
• Changes documented, CHANGELOG.md

@tfrederiksen might you have a 2nd look here, I tested on some geometries, all seems fine, but this would warrant a 2nd look

[zerothi]

I have found the issue. In some rare cases the algorithm would search radii just crossing the boundaries which in the end could result in some atoms being part of the first neighbour circle, but its true neighbours would not be part of the 2nd circle.

This was a very subtle issue, and would generally only occur for somewhat large structures.

[pfebrer]

Maybe moving to #393 would fix it as well ;)

[zerothi]

Damn, not fixed yet...
@pfebrer yeah, but...

[zerothi]

I don't know if that will be stable given I have to prepare for the workshop...

[zerothi]

I'll play a bit more tonight...

[pfebrer]

I mean replacing the internals, while keeping the same geom.close API.

[zerothi]

I mean replacing the internals, while keeping the same geom.close API.

I know. I don't think that is smart to do just before a workshop... It's a big change, and it requires handling the supercell indices in a strict manner. I need more time.

[codecov]

Codecov Report

Attention: 14 lines in your changes are missing coverage. Please review.

Comparison is base (d76ae42) 87.73% compared to head (5d5753d) 87.73%.
Report is 14 commits behind head on main.

❗ Current head 5d5753d differs from pull request most recent head 748a026. Consider uploading reports for the commit 748a026 to get more accurate results

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #634   +/-   ##
=======================================
  Coverage   87.73%   87.73%           
=======================================
  Files         362      362           
  Lines       48348    48308   -40     
=======================================
- Hits        42417    42382   -35     
+ Misses       5931     5926    -5     

Files	Coverage Δ
src/sisl/geom/_common.py	100.00% <100.00%> (ø)
src/sisl/lattice.py	93.09% <ø> (ø)
src/sisl/orbital.py	90.78% <ø> (ø)
src/sisl/physics/_brillouinzone_apply.py	79.69% <100.00%> (ø)
src/sisl/io/_multiple.py	90.55% <92.85%> (-0.86%)	⬇️
src/sisl/sparse_geometry.py	90.28% <88.88%> (-0.09%)	⬇️
src/sisl/_dispatcher.py	69.26% <42.85%> (-1.03%)	⬇️
src/sisl/geometry.py	86.35% <84.44%> (-0.04%)	⬇️

... and 5 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[tfrederiksen]

I am afraid this is not a complete fix. Try this:

import sisl
import numpy as np

def test(g, runs=100):
    x = []
    for i in range(runs):
        H = sisl.Hamiltonian(g)
        H.construct([(0.1, 3), (1, -2.7)])
        x.append(H.nnz)
        if x[-1] != x[0]:
            print(f"Error encountered after {i+1} runs: {x}")
            return
    print(f"Success after {runs} runs")

np.random.seed(1234)
g1 = sisl.geom.zgnr(10) * (20, 1, 1)
g2 = sisl.geom.sc(1, sisl.Atom('H', R=1))
g = g1.add(g2)
g.set_nsc([1, 1, 1])

test(g)
print(sisl.__version__)

which returns:

Error encountered after 3 runs: [4749, 4749, 4747]
0.14.3.dev30+g32b3be34

[zerothi]

Yeah, I discovered, see my 2nd comment... :(

[zerothi]

@tfrederiksen could you try now...

Please do note, that the problem arises because of the way the search is built. It does a batched search space made of 2 radii. One that is an integer number times the self.maxR() and one just one bigger. Thus when you use an R that is bigger than maxR, then you end up in situations where the edges are not great... :(
Now I am trying to catch this, but it will only do so with pure construct (as the ones you have). If people start doing their own def construct_func, then I can't get that information... And thus resulting in incorrect usage.

For now there is a warning raised when sisl thinks it is used incorrectly. This is why your above script fails. It is very hard for the batched search to correct this as it does not know your R / maxR() which is what sisl would need to figure it out correctly.

I have tested ~6 geometries with correct input, and they don't fail... :)
But... As always, there may be some corner case I haven't looked into... :(

[tfrederiksen]

It looks good to me now. Despite several attempts I couldn't design a case with failure :)

Would it make sense to have a test with a critical comparison of close vs construct for some random structure? Perhaps the following script can serve as inspiration (although everything but the last assert would have passed before the fix)?

import sisl
import numpy as np

def test(g):
    T = 1 + np.arange(na)
    # close call first
    H = sisl.Hamiltonian(g)
    for ia in g:
        idx = g.close(ia, R=R)
        for ja, t in zip(idx, T):
            H[ia, ja] = t
    x = [H.nnz + np.sum(H)]
    # now compare with construct calls
    for i in range(10):
        for method in ['rand', 'sphere', 'cube']:
            H = sisl.Hamiltonian(g)
            H.construct([R, T], method=method)
            x.append(H.nnz + np.sum(H))
            if x[-1] != x[0]:
                return False # inconsistency detected
    return x[0] # no inconsistency detected

np.random.seed(1234)

# build a random geometry with different orbitals
na = 7
R = 1e-4 + np.arange(na)
atoms = [sisl.Atom('HBCNOFPISK'[i], R=R[i]) for i in range(na)]
g = sisl.geom.fcc(2.5, atoms=atoms[0]) * (4, 3, 2)
g.xyz[:] += 10 * np.random.rand(g.na, 3)
for ia in g:
    g.atoms[ia] = atoms[ia % na]
g.reduce()

g.set_nsc([1, 1, 1])
assert test(g) == 1224.0
g.set_nsc([1, 1, 9])
assert test(g) == 4518.0
g.set_nsc([3, 3, 3])
assert test(g) == 18286.0

g2 = g.add(g.move(1e-2)) # put atoms almost on top of each other!
g2.set_nsc([1, 1, 1])
assert test(g2) == 4960.0
g2.set_nsc([1, 1, 3])
assert test(g2) == 13038.0

g3 = g.scale(1e3) # separate atoms by a lot!
g3.set_nsc([1, 1, 3])
assert test(g3) == 2 * g3.na

g4 = g.scale(1e-12) # collapse!
g4.set_nsc([1, 1, 3])
assert test(g4) == 6 * g4.na ** 2

Btw, I noticed some debug lines you may want to remove:

sisl/src/sisl/geometry.py

Lines 974 to 975 in a6344c4

	print(not_passed.nonzero()[0])
	print(np.sum(not_passed), len(self))

[zerothi]

@tfrederiksen thanks for the tests!! I'll see if some of that test can be added, currently there is a test that ensures symmetric construct, it is always hard to know when to cut the line here...